Search code examples
pythonmathgeometryellipse

Find intersection point between line and ellipse in Python

I am trying to find the point which intersects with the line and the ellipse. I can find the point manually as seen in the code but how can I find it automatically. Also when the red point is inside the ellipse the green point should also be projected on to the ellipse when extending the line. How can I solve this? Code examples are appreciated! Thanks!

enter image description here

import cv2
import numpy as np

# Parameters
ellipse_center = np.array([238, 239])
ellipse_axes = np.array([150, 100])  # Semi-major and semi-minor axes
ellipse_angle = np.radians(70)  # Angle of rotation (in radians)

line_point1 = ellipse_center  # Point on the line
line_point2 = np.array([341, 125])  # Another point on the line (moving point)

# Initialize OpenCV window
window_name = "Projecting Point on Rotated Ellipse"
cv2.namedWindow(window_name)
canvas = np.zeros((400, 400, 3), dtype=np.uint8)
    
while True:
    # Clear canvas
    canvas.fill(0)

    # Draw the rotated ellipse
    cv2.ellipse(canvas, tuple(ellipse_center), tuple(ellipse_axes), np.degrees(ellipse_angle), 0, 360, (255, 255, 255), 2)

    # Draw the moving point
    cv2.circle(canvas, (int(line_point2[0]), int(line_point2[1])), 5, (0, 0, 255), -1)

    # Draw the center point of the ellipse
    cv2.circle(canvas, tuple(ellipse_center), 5, (0, 255, 255), -1)

    # Draw the line between moving point and the center of the ellipse
    cv2.line(canvas, (int(line_point2[0]), int(line_point2[1])), tuple(ellipse_center), (0, 0, 255), 1)

    intersection_points = (310, 159) # <- How to find this?
    cv2.circle(canvas, intersection_points, 5, (0, 255, 0), -1)

    # Show canvas
    cv2.imshow(window_name, canvas)

    # Exit if 'q' is pressed
    if cv2.waitKey(1) & 0xFF == ord('q'):
        break

# Close OpenCV window
cv2.destroyAllWindows()

When using the ellipse equation and the line equation I get this:

import cv2
import numpy as np
from scipy.optimize import fsolve

# Parameters
ellipse_center = np.array([238, 239])
ellipse_axes = np.array([150, 100])  # Semi-major and semi-minor axes
ellipse_angle = np.radians(70)  # Angle of rotation (in radians)

line_point1 = ellipse_center  # Point on the line
line_point2 = np.array([341, 125])  # Another point on the line (moving point)

# Ellipse equation: (x - h)^2 / a^2 + (y - k)^2 / b^2 = 1, where (h, k) is the center of the ellipse
def ellipse_eq(params):
    x, y = params
    h, k = ellipse_center
    a, b = ellipse_axes
    return ((x - h) ** 2 / a ** 2) + ((y - k) ** 2 / b ** 2) - 1

# Line equation: y = mx + c
def line_eq(x):
    m = (line_point2[1] - line_point1[1]) / (line_point2[0] - line_point1[0])
    c = line_point1[1] - m * line_point1[0]
    return m * x + c

# Intersection of the ellipse and the line
def intersection_point():
    x_intercept = fsolve(lambda x: ellipse_eq([x, line_eq(x)]), line_point2[0])
    y_intercept = line_eq(x_intercept)
    return x_intercept[0], y_intercept[0]

# Initialize OpenCV window
window_name = "Projecting Point on Rotated Ellipse"
cv2.namedWindow(window_name)
canvas = np.zeros((400, 400, 3), dtype=np.uint8)

while True:
    # Clear canvas
    canvas.fill(0)

    # Draw the rotated ellipse
    cv2.ellipse(canvas, tuple(ellipse_center), tuple(ellipse_axes), np.degrees(ellipse_angle), 0, 360, (255, 255, 255), 2)

    # Draw the moving point
    cv2.circle(canvas, (int(line_point2[0]), int(line_point2[1])), 5, (0, 0, 255), -1)

    # Draw the center point of the ellipse
    cv2.circle(canvas, tuple(ellipse_center), 5, (0, 255, 255), -1)

    # Draw the line between moving point and the center of the ellipse
    cv2.line(canvas, (int(line_point2[0]), int(line_point2[1])), tuple(ellipse_center), (0, 0, 255), 1)

    # Find and draw the intersection point
    intersection_points = intersection_point()
    cv2.circle(canvas, (int(intersection_points[0]), int(intersection_points[1])), 5, (0, 255, 0), -1)

    # Show canvas
    cv2.imshow(window_name, canvas)

    # Exit if 'q' is pressed
    if cv2.waitKey(1) & 0xFF == ord('q'):
        break

# Close OpenCV window
cv2.destroyAllWindows()

enter image description here

Which does not solve my problem. What am I doing wrong here?

EDIT: Using MBo's approach with the code:

import cv2
import numpy as np

# Parameters
ellipse_center = np.array([238, 239])
ellipse_axes = np.array([150, 100])  # Semi-major and semi-minor axes
ellipse_angle = np.radians(70)  # Angle of rotation (in radians)

line_point1 = ellipse_center  # Point on the line (ellipse center)
line_point2 = np.array([341, 125])  # Another point on the line (moving point)

# Initialize OpenCV window
window_name = "Projecting Point on Rotated Ellipse"
cv2.namedWindow(window_name)
canvas = np.zeros((400, 400, 3), dtype=np.uint8)

while True:
    # Clear canvas
    canvas.fill(0)

    # Draw the rotated ellipse
    cv2.ellipse(canvas, tuple(ellipse_center), tuple(ellipse_axes), np.degrees(ellipse_angle), 0, 360, (255, 255, 255), 2)

    # Draw the moving point
    cv2.circle(canvas, (int(line_point2[0]), int(line_point2[1])), 5, (0, 0, 255), -1)

    # Draw the center point of the ellipse
    cv2.circle(canvas, tuple(ellipse_center), 5, (0, 255, 255), -1)

    # Draw the line between moving point and the center of the ellipse
    cv2.line(canvas, (int(line_point2[0]), int(line_point2[1])), tuple(ellipse_center), (0, 0, 255), 1)

    # Calculate the intersection point
    t = np.arctan2(line_point2[1] - ellipse_center[1], line_point2[0] - ellipse_center[0]) - np.arctan2(np.sin(ellipse_angle) * ellipse_axes[0], np.cos(ellipse_angle) * ellipse_axes[1])
    x_intersection = ellipse_center[0] + ellipse_axes[0] * np.cos(t) * np.cos(ellipse_angle) - ellipse_axes[1] * np.sin(t) * np.sin(ellipse_angle)
    y_intersection = ellipse_center[1] + ellipse_axes[0] * np.cos(t) * np.sin(ellipse_angle) + ellipse_axes[1] * np.sin(t) * np.cos(ellipse_angle)
    intersection_point = (int(x_intersection), int(y_intersection))

    # Draw the intersection point
    cv2.circle(canvas, intersection_point, 5, (0, 255, 0), -1)

    # Show canvas
    cv2.imshow(window_name, canvas)

    # Exit if 'q' is pressed
    if cv2.waitKey(1) & 0xFF == ord('q'):
        break

# Close OpenCV window
cv2.destroyAllWindows()

Gets me:

enter image description here

Solution

  • First rotate the line back by ellipse_angle so that it corresponds to the unrotated ellipse.

    Substitute (for y) the equation of the line (y=mx+c) into the equation of the ellipse (x-h)2/a2+(y-k)2/b2=1. This will give a quadratic in x which can be solved by the usual quadratic-equation formula.

    The quadratic is Ax2+Bx+C=0, where, on rearrangement,

    A=a2m2+b2

    B=2a2m(c-k)-2b2h

    C=b2h2+a2(c-k)2-a2b2

    Solve for x and get y from the equation of the line. Be careful to get the correct quadratic root.

    Then rotate the intersection point back again by ellipse_angle.

    As with any quadratic equation you can run into occasions with a single root (line just touches the ellipse) and no roots (line misses the ellipse).

    Code below - I have indicated the modified section of your code.

    import cv2
import numpy as np

# Parameters
ellipse_center = np.array([238, 239])
ellipse_axes = np.array([150, 100])  # Semi-major and semi-minor axes
ellipse_angle = np.radians(70)  # Angle of rotation (in radians)

line_point1 = ellipse_center  # Point on the line (ellipse center)
line_point2 = np.array([341, 125])  # Another point on the line (moving point)

# Initialize OpenCV window
window_name = "Projecting Point on Rotated Ellipse"
cv2.namedWindow(window_name)
canvas = np.zeros((400, 400, 3), dtype=np.uint8)

while True:
    # Clear canvas
    canvas.fill(0)

    # Draw the rotated ellipse
    cv2.ellipse(canvas, tuple(ellipse_center), tuple(ellipse_axes), np.degrees(ellipse_angle), 0, 360, (255, 255, 255), 2)

    # Draw the moving point
    cv2.circle(canvas, (int(line_point2[0]), int(line_point2[1])), 5, (0, 0, 255), -1)

    # Draw the center point of the ellipse
    cv2.circle(canvas, tuple(ellipse_center), 5, (0, 255, 255), -1)

    # Draw the line between moving point and the center of the ellipse
    cv2.line(canvas, (int(line_point2[0]), int(line_point2[1])), tuple(ellipse_center), (0, 0, 255), 1)


    ############## MODIFIED SECTION ############
    # Rotate the line back to an unrotated ellipse
    theta = np.arctan2( line_point2[1] - line_point1[1], line_point2[0] - line_point1[0] ) - ellipse_angle
    m = np.tan( theta )
    c = line_point1[1] - m * line_point1[0]
    a, b, h, k = ellipse_axes[0], ellipse_axes[1], ellipse_center[0], ellipse_center[1]
    A = ( a * m ) ** 2 + b ** 2
    B = 2 * ( a ** 2 * m * ( c - k ) - b ** 2 * h )
    C = ( b * h ) ** 2 + ( a * ( c - k ) ) ** 2 - ( a * b ) ** 2
    discriminant = B ** 2 - 4 * A * C
    # Intersection of back-rotated line with back-rotated ellipse (need correct root)
    theta = np.mod( theta, 2 * np.pi )
    x1 = ( -B + np.sqrt( discriminant ) ) / ( 2 * A )                 # quadratic solution (larger root)
    if theta > 0.5 * np.pi and theta < 1.5 * np.pi: x1 = -B/A - x1    # need other root
    y1 = m * x1 + c
    # Rotate the line forward again
    x1, y1 = h + ( x1 - h ) * np.cos( ellipse_angle ) - ( y1 - k ) * np.sin( ellipse_angle ),  k + ( x1 - h ) * np.sin( ellipse_angle ) + ( y1 - k ) * np.cos( ellipse_angle )
    intersection_point = (int(x1), int(y1))
    ############## END OF MODIFIED SECTION ############

    # Draw the intersection point
    cv2.circle(canvas, intersection_point, 5, (0, 255, 0), -1)

    # Show canvas
    cv2.imshow(window_name, canvas)

    # Exit if 'q' is pressed
    if cv2.waitKey(1) & 0xFF == ord('q'):
        break

# Close OpenCV window
cv2.destroyAllWindows()

    Output: enter image description here