Tangent lines to a rotated ellipse from a point

Question

I need to calculate the equations of tangent lines to a rotated ellipse from a specific point outside of the ellipse. I found Ellipse.rotate() method in Sympy which is almost exactly what I need, but only rotations with increments of pi/2 are supported and I want to be able to do this with other rotation angles. Any suggestion on how to do this in Python?

Answer 1

You can define your ellipse as a sequence of transformations that are applied to the unit circle:

1. Stretch by a and b.
2. Rotate by angle.
3. Translate by center.

Then this sequence of transformations can be applied in reverse order and using the inverse of each transformation to the reference point, in order to solve the problem w.r.t. the unit circle. Here the angles of the two tangent points on the unit circle are given by the constraint of perpendicularity: 1 - p0*cos(t) - q0*sin(t) = 0 where p0, q0 is the reference point.

See this answer for mathematical details.

from functools import reduce
from typing import Tuple

import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
import numpy as np


class Transformation:
    def transform(self, point: Tuple[float, float]) -> Tuple[float, float]:
        raise NotImplementedError

    def inverse_transform(self, point: Tuple[float, float]) -> Tuple[float, float]:
        raise NotImplementedError


class Translation(Transformation):
    def __init__(self, offset: Tuple[float, float]):
        self.dx, self.dy = offset

    def transform(self, point):
        return point[0] + self.dx, point[1] + self.dy

    def inverse_transform(self, point):
        return point[0] - self.dx, point[1] - self.dy


class Stretch(Transformation):
    def __init__(self, scale: Tuple[float, float]):
        self.scale = scale

    def transform(self, point):
        return point[0]*self.scale[0], point[1]*self.scale[1]

    def inverse_transform(self, point):
        return point[0]/self.scale[0], point[1]/self.scale[1]


class Rotation(Transformation):
    def __init__(self, angle: float):
        self.c, self.s = np.cos(angle), np.sin(angle)

    def transform(self, point):
        return self.c*point[0] - self.s*point[1], self.s*point[0] + self.c*point[1]

    def inverse_transform(self, point):
        return self.c*point[0] + self.s*point[1], -self.s*point[0] + self.c*point[1]


class Transformations(Transformation):
    def __init__(self, *transformations: Transformation):
        self.transformations = transformations

    def transform(self, point):
        return reduce(lambda p,t: t.transform(p), self.transformations, point)

    def inverse_transform(self, point):
        return reduce(lambda p,t: t.inverse_transform(p), self.transformations[::-1], point)


def find_tangent_lines(
    center: Tuple[float, float],
    semi_axes: Tuple[float, float],
    rotation: float,
    reference_point: Tuple[float, float],
):
    """Find the Ellipse's two tangents that go through a reference point.

    Args:
        center: The center of the ellipse.
        semi_axes: The semi-major and semi-minor axes of the ellipse.
        rotation: The counter-clockwise rotation of the ellipse in radians.
        reference_point: The coordinates of the reference point.

    Returns:
        (m1, h1): Slope and intercept of the first tangent.
        (m2, h2): Slope and intercept of the second tangent.
    """
    transformations = Transformations(
        Stretch(semi_axes),
        Rotation(rotation),
        Translation(center),
    )
    p = transformations.inverse_transform(reference_point)
    p_norm = p[0]**2 + p[1]**2

    if p_norm < 1:
        raise ValueError('Reference point lies inside the ellipse')

    theta = (
        2*np.arctan((p[1] + np.sqrt(p_norm - 1))/(p[0]+1)),
        2*np.arctan((p[1] - np.sqrt(p_norm - 1))/(p[0]+1)),
    )
    p_ellipse = [
        transformations.transform((np.cos(t), np.sin(t)))
        for t in theta
    ]
    slope = [
        (e[1] - reference_point[1]) / (e[0] - reference_point[0])
        for e in p_ellipse
    ]
    intercept = [
        e[1] - s*e[0]
        for e, s in zip(p_ellipse, slope)
    ]
    return (
        (slope[0], intercept[0]),
        (slope[1], intercept[1]),
    )


# Example:
CENTER = 1, 2
SEMI_AXES = 3, 1
ROTATION = np.pi/3
REFERENCE_POINT = -2, 3

(m1, h1), (m2, h2) = find_tangent_lines(
    center=CENTER,
    semi_axes=SEMI_AXES,
    rotation=ROTATION,
    reference_point=REFERENCE_POINT,
)

fig, ax = plt.subplots()
ax.add_patch(Ellipse(CENTER, 2*SEMI_AXES[0], 2*SEMI_AXES[1], ROTATION/np.pi*180, color='tab:blue', alpha=0.4))
ax.scatter([REFERENCE_POINT[0]], [REFERENCE_POINT[1]], color='tab:blue', s=50)
ax.axline((0, h1), slope=m1, color='tab:orange', lw=1)
ax.axline((0, h2), slope=m2, color='tab:orange', lw=1)
plt.show()

Answer 2

Another option is to consider all lines that go through the reference point and any point of the ellipse. The function that maps ellipse angles to the slopes of these lines will have two extrema for the two tangent lines of the ellipse. Hence, solving d/dt slopes(t) = 0 will yield the two values for t (the ellipse angles).

See this answer for mathematical details.

from typing import Tuple

import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
import numpy as np


def find_tangent_lines(
    center: Tuple[float, float],
    semi_axes: Tuple[float, float],
    rotation: float,
    reference_point: Tuple[float, float],
):
    """Find the Ellipse's two tangents that go through a reference point.

    Args:
        center: The center of the ellipse.
        semi_axes: The semi-major and semi-minor axes of the ellipse.
        rotation: The counter-clockwise rotation of the ellipse in radians.
        reference_point: The coordinates of the reference point.

    Returns:
        (m1, h1): Slope and intercept of the first tangent.
        (m2, h2): Slope and intercept of the second tangent.
    """
    x0, y0 = center
    a, b = semi_axes
    s, c = np.sin(rotation), np.cos(rotation)
    p0, q0 = reference_point

    A = y0 - q0
    D = x0 - p0
    B = a*s
    E = a*c
    C = b*c
    F = -b*s

    denominator = np.sqrt((C*D - A*F)**2 + (A*E - B*D)**2)
    if not (-1 <= (B*F - C*E) / denominator <= 1):
        raise ValueError('Reference point lies inside the ellipse')

    beta = np.arctan2(
        (C*D - A*F) / denominator,
        (A*E - B*D) / denominator,
    )
    theta = [
        -beta + np.arcsin((B*F - C*E) / denominator),
        -beta - np.arcsin((B*F - C*E) / denominator) + np.pi,
    ]
    p_ellipse = [
        (
            x0 + E*np.cos(t) + F*np.sin(t),
            y0 + B*np.cos(t) + C*np.sin(t),
        )
        for t in theta
    ]
    slope = [
        (e[1] - reference_point[1]) / (e[0] - reference_point[0])
        for e in p_ellipse
    ]
    intercept = [
        e[1] - s*e[0]
        for e, s in zip(p_ellipse, slope)
    ]
    return (
        (slope[0], intercept[0]),
        (slope[1], intercept[1]),
    )


# Example:
CENTER = 1, 2
SEMI_AXES = 3, 1
ROTATION = np.pi/3
REFERENCE_POINT = -2, 3

(m1, h1), (m2, h2) = find_tangent_lines(
    center=CENTER,
    semi_axes=SEMI_AXES,
    rotation=ROTATION,
    reference_point=REFERENCE_POINT,
)

fig, ax = plt.subplots()
ax.add_patch(Ellipse(CENTER, 2*SEMI_AXES[0], 2*SEMI_AXES[1], ROTATION/np.pi*180, color='tab:blue', alpha=0.4))
ax.scatter([REFERENCE_POINT[0]], [REFERENCE_POINT[1]], color='tab:blue', s=50)
ax.axline((0, h1), slope=m1, color='tab:orange', lw=1)
ax.axline((0, h2), slope=m2, color='tab:orange', lw=1)
plt.show()