从一点到旋转椭圆的切线
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()
您可以将椭圆定义为应用于单位圆的一系列变换:
- 拉伸
a
和 b
。
- 旋转
angle
。
- 由
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()
我需要从椭圆外的特定点计算旋转椭圆的切线方程。我在 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()
您可以将椭圆定义为应用于单位圆的一系列变换:
- 拉伸
a
和b
。 - 旋转
angle
。 - 由
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()