从一点到旋转椭圆的切线

Tangent lines to a rotated ellipse from a point

我需要从椭圆外的特定点计算旋转椭圆的切线方程。我在 Sympy 中找到了 Ellipse.rotate() 方法,这几乎正是我所需要的,但只支持增量为 pi/2 的旋转,我希望能够使用其他旋转角度来执行此操作。关于如何在 Python?

中执行此操作的任何建议

您可以使用以下代码段来计算切线的斜率和截距。感谢 this answer 的数学输入。由于浮点精度有限,如果参考点正好在椭圆上,可能会导致计算失败。

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 = (-a**2*s**2 - b**2*c**2 + (y0-q0)**2)
    B = 2*(c*s*(a**2-b**2) - (x0-p0)*(y0-q0))
    C = (-a**2*c**2 - b**2*s**2 + (x0-p0)**2)

    if B**2 - 4*A*C < 0:
        raise ValueError('Reference point lies inside the ellipse')

    t1, t2 = (
        (-B + np.sqrt(B**2 - 4*A*C))/(2*A),
        (-B - np.sqrt(B**2 - 4*A*C))/(2*A),
    )
    return (
        (1/t1, q0 - p0/t1),
        (1/t2, q0 - p0/t2),
    )


# 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()

您可以将椭圆定义为应用于单位圆的一系列变换:

  1. 拉伸 ab
  2. 旋转angle
  3. center 翻译。

然后可以按相反的顺序应用这个变换序列,并使用每个变换的逆向参考点,以解决问题w.r.t。单位圆。这里单位圆上的两个切点的角度由垂直约束给出:1 - p0*cos(t) - q0*sin(t) = 0 其中p0, q0是参考点。

有关数学详细信息,请参阅 this answer

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()

另一种选择是考虑所有通过参考点椭圆任意一点的线。将椭圆角映射到这些线的斜率的函数将为椭圆的两条切线提供两个极值。因此,求解 d/dt slopes(t) = 0 将得到 t(椭圆角)的两个值。

有关数学详细信息,请参阅 this answer

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()