Я хочу рассчитать точки пересечения эллипса и линии с помощью Python. Однако как мне это сделать, используя только математическое квадратное уравнение, если заданы как конечные точки линии, так и центральная точка и радиусы эллипса?
Ни одно из решений в stackoverflow не описывает решение, которое я ищу. Они используют либо numpy, scipy и другие библиотеки, либо используют тригонометрию — ни одна из них не использует чисто математический процесс.
Сравнивая решение в этом посте, например:
У вас есть формула эллипса:
(x - center_x)^2 / a^2 + (y - center_y)^2 / b^2 = 1
И формула строки:
у = м*х + с
Решение в связанном сообщении описывает квадратичную формулу, полученную из него при перестановке:
А=а2м2+b2
B=2a2m(c-k)-2b2h
C=b2h2+a2(c-k)2-a2b2
Я хочу решить это уравнение, используя только математические вычисления, находя решение для x, а затем решая для y, используя формулу линии.
Однако при использовании приведенной выше формулы полученные точки пересечения ошибочны и не лежат на эллипсе или прямой.
Это код, который я пробовал:
from math import *
def GetIntersectionPoints(line_point_a, line_point_b, center_x, center_y, rad_x, rad_y):
# use line formula to calculate m and c:
x1, y1, x2, y2 = *line_point_a, *line_point_b
m = (y1 - y2) / (x1 - x2) if not x1 - x2 == 0 else 0
c = y2 - (m * x2)
# Using quadratic equation provided in linked solution:
A = (rad_y ** 2) + ((rad_x ** 2) * (m ** 2))
B = (2 * (rad_x ** 2) * m * (c - center_y)) - 2 * (rad_y ** 2) * center_x
C = (rad_y ** 2) * (center_x ** 2) + (rad_x ** 2) * (c - center_y) * 2 - (rad_x ** 2) * (rad_y ** 2)
# Solving quadratic equation:
xi1 = (-B - sqrt((B ** 2) - (4 * A * C))) / (2 * A)
xi2 = (-B + sqrt((B ** 2) - (4 * A * C))) / (2 * A)
# Solving for y:
yi1 = m * xi1 + c
yi2 = m * xi2 + c
return [[xi1, yi1], [xi2, yi2]]
Как мне правильно решить уравнение и формулу в Python, не используя scipy, numpy и т. д. (или тригонометрию и другие решения, описанные в stackoverflow), чисто математически, чтобы получить правильные точки пересечения, лежащие на эллипсе?
Ваш код для C
неправильный. Проверьте это внимательно. У вас есть термин (c - center_y) * 2
, и он должен быть (c - center_y) ** 2
(т. е. умножение на 2 должно быть «в степени 2»). Один из простых способов это выяснить — отметить, что если бы все термины имели физические единицы, то этот бит был бы несовместимым по размерности.
См. пояснения в коде, как вычисляются точки пересечения. Код не использует никаких модулей и выполняет все вычисления на чистом Python:
def get_intersection_points(line_point_a, line_point_b, center_x, center_y, rad_x, rad_y):
# Line endpoints
x1, y1 = line_point_a
x2, y2 = line_point_b
# Calculate line slope (m) and intercept (c)
if x1 != x2:
m = (y2 - y1) / (x2 - x1)
else:
m = float('inf') # Vertical line case
c = y1 - m * x1 if m != float('inf') else x1 # If line is vertical, c is x1
# Ellipse parameters
h = center_x
k = center_y
a = rad_x
b = rad_y
if m != float('inf'):
# Coefficients of the quadratic equation
A = b**2 + a**2 * m**2
B = 2 * a**2 * m * (c - k) - 2 * b**2 * h
C = b**2 * h**2 + a**2 * (c - k)**2 - a**2 * b**2
# Discriminant
D = B**2 - 4 * A * C
if D < 0:
return [] # No intersection points
# Calculate the x-coordinates of the intersection points
sqrt_D = D**0.5 # math.sqrt(D)
xi1 = (-B - sqrt_D) / (2 * A)
xi2 = (-B + sqrt_D) / (2 * A)
# Calculate the y-coordinates of the intersection points
yi1 = m * xi1 + c
yi2 = m * xi2 + c
return [[xi1, yi1], [xi2, yi2]]
else:
# Vertical line case (x = c)
x = c
A = 1 / a**2
B = -2 * k / b**2
C = (h**2 / a**2) + (k**2 / b**2) - 1
# Quadratic equation for y
A_y = 1 / b**2
B_y = -2 * k / b**2
C_y = (x - h)**2 / a**2 + k**2 / b**2 - 1
# Discriminant
D_y = B_y**2 - 4 * A_y * C_y
if D_y < 0:
return [] # No intersection points
# Calculate the y-coordinates of the intersection points
sqrt_D_y = D_y**0.5 # math.sqrt(D_y)
yi1 = (-B_y - sqrt_D_y) / (2 * A_y)
yi2 = (-B_y + sqrt_D_y) / (2 * A_y)
return [ [x , yi1], [ x , yi2] ]
# Example usage:
line_point_a = (-7,-3)
line_point_b = (7, 3)
center_x = 0
center_y = 0
rad_x = 7
rad_y = 3
print(get_intersection_points(line_point_a, line_point_b, center_x, center_y, rad_x, rad_y))
from turtle import *
import math
def draw_ellipse(center_x, center_y, rad_x, rad_y):
# Move the turtle to the starting position without drawing
up()
goto(center_x + rad_x, center_y)
down()
# Number of steps to make the ellipse smooth
steps = 360
# Draw the ellipse
for step in range(steps + 1):
# Calculate the angle in radians
angle = math.radians(step)
# Parametric equations for the ellipse
x = center_x + rad_x * math.cos(angle)
y = center_y + rad_y * math.sin(angle)
# Move the turtle to the calculated position
goto(x, y)
from turtle import *
Screen().setworldcoordinates(-10, -10, 10, 10)
Screen().bgcolor("lightgrey")
color("blue")
shape("circle")
pensize(5)
draw_ellipse(center_x, center_y, rad_x, rad_y)
up(); goto( center_x, center_y - rad_x ); color("cyan")
down(); circle(rad_x, steps=36)
up(); goto( center_x, center_y - rad_y );
down(); circle(rad_y, steps=36); color("green")
up(); goto( *line_point_a)
down(); goto(*line_point_b)
up()
color("red")
for point in get_intersection_points(line_point_a, line_point_b, center_x, center_y, rad_x, rad_y) :
goto(*point)
stamp()
shape("turtle")
color("black")
goto(-8,-8)
done()
Чтобы проверить корректность работы кода, Черепаха рисует в левом нижнем углу изображения эллипс, две окружности, линию и найденные точки пересечения, отдыхающие после проделанной работы:
Страница, на которую вы указываете выше, буквально показывает вам нужный код. Он использует различные математические функции из
numpy
, но они также доступны вmath
.