求一对二次曲线的内公切线
Finding inner common tangent to a pair of conics
我一直在尝试求出两个旋转椭圆的公切线。我正在按照 Edward Doolittle in the . The two ellipses are given as the equation given in Wiki.
给出的方法进行操作
矩阵形式的椭圆可以表示为：
第一个椭圆以 (0,0) 为中心旋转 45 度，半长轴和半短轴长度为 2,1。第二个椭圆以 (15,0) 为中心，旋转 120 度，半长轴和半短轴长度为 3,1
两个椭圆的伴随矩阵的线性组合是两个椭圆组合的对偶
我得到这个值
.
然后我试图找到 t 的值，这将使圆锥曲线（矩阵上方）退化。
我发现 t 的值为 (-0.05,0.29,2.46)。但是，当我将这些值放回上述矩阵时，我无法将矩阵简化为两个变量形式。我总是在处理 3 个变量。例如，如果我输入 t = -0.05，那么我会得到以下结果：
有人可以帮我解决这个问题吗？
它归结为找到一种算法来求解两个变量的两个二次方程组，通过将其解释为二次曲线的射影几何铅笔，然后找到铅笔的三个退化二次曲线以及射影变换将这三个退化二次曲线简化到可以在新的简化坐标系中非常容易地读出解，然后将它们转换回原始坐标系。
我在python中勾勒了一个算法，我认为它似乎适用于你的例子......但是我很匆忙并且没有正确检查它，所以可能存在错误......
import numpy as np
import math
# technical module, functions, details
def homogenize(x):
return np.array([x[0], x[1], 1])
def cos_sin(angle_deg):
return math.cos(angle_deg*math.pi/180), math.sin(angle_deg*math.pi/180)
def rotation(cs_sn):
return np.array([[cs_sn[0], -cs_sn[1]],
[cs_sn[1], cs_sn[0]]])
# defining the isometry (the rotation plus translation) transformation
# between the coordinate system aligned with the conic and a general (world)
# coordinate system
def isom_inverse(angle, translation):
'''
isometry from conic-aligned coordinate system (conic attached)
to global coordinate system (world system)
'''
cos_, sin_ = cos_sin(angle)
return np.array([[cos_, -sin_, translation[0]],
[sin_, cos_, translation[1]],
[ 0, 0, 1]])
def isom(angle, translation):
'''
isometry from global coordinate system (world system)
to conic-aligned coordinate system (conic attached)
'''
cos_, sin_ = cos_sin(-angle)
tr = - rotation((cos_, sin_)).dot(translation)
return np.array([[ cos_, -sin_, tr[0]],
[ sin_, cos_, tr[1]],
[ 0, 0, 1 ]])
# calculating the coinc defined by a pair of axes' lengts,
# axes rotation angle and center of the conic
def Conic(major, minor, angle, center):
D = np.array([[minor**2, 0, 0],
[ 0, major**2, 0],
[ 0, 0, -(minor*major)**2]])
U = isom(angle, center)
return (U.T).dot(D.dot(U))
# calculating the coinc dual to the conic defined by a pair of axes' lengths,
# axes rotation angle and center of the conic
def dual_Conic(major, minor, angle, center):
D_1 = np.array([[major**2, 0, 0],
[ 0, minor**2, 0],
[ 0, 0, -1]])
U_1 = isom_inverse(angle, center)
return (U_1).dot(D_1.dot(U_1.T))
# transforming the matrix of a conic into a vector of six coefficients
# of a quadratic equation with two variables
def conic_to_equation(C):
'''
c[0]*x**2 + c[1]*x*y + c[2]*y**2 + c[3]*x + c[4]*y + c[5] = 0
'''
return np.array([C[0,0], 2*C[0,1], C[1,1], 2*C[0,2], 2*C[1,2], C[2,2]])
# transforming the vector of six coefficients
# of a quadratic equation with two variables into a matrix of
# the corresponding conic
def equation_to_conic(eq):
'''
eq[0]*x**2 + eq[1]*x*y + eq[2]*y**2 + eq[3]*x + eq[4]*y + eq[5] = 0
'''
return np.array([[2*eq[0], eq[1], eq[3]],
[ eq[1], 2*eq[2], eq[4]],
[ eq[3], eq[4], 2*eq[5]]]) / 2
# given a point (x,y) define the vector (x^2, xy, y^2, x, y, 1)
def argument(x):
return np.array([x[0]**2, x[0]*x[1], x[1]**2, x[0], x[1], 1])
# given x = (x[0],x[1]) calculate the value of the quadratic equation with
# six coefficients coeff
def quadratic_equation(x, coeff):
'''
coeff[0]*x**2 + coeff[1]*x*y + coeff[2]*y**2 + coeff[3]*x + coeff[4]*y + coeff[5] = 0
'''
return coeff.dot( argument(x) )
# given a pair of conics, as a pair of symmetric matrices,
# calculate the vector k = (k[0], k[1], k[2]) of values for each of which
# the conic c1 - k[i]*c2 from the pencil of conics c1 - t*c2
# is a degenerate conic (the anti-symmetric product of a pair of linear forms)
# and also find the matrix U
# of the projective transformation that simplifies the geometry of
# the pair of conics, the geometry of the pencil c1 - t*c2 in general,
# as well as the geometry of the three degenerate conics in particular
def transform(c1, c2):
'''
c1 and c2 are 3 by 3 symmetric matrices of the two conics
'''
c21 = np.linalg.inv(c2).dot(c1)
k, U = np.linalg.eig(c21)
return k, U
# the same as before, but for a pair of equations instead of matrices of conics
def eq_transform(eq1, eq2):
'''
eq1 and eq2 = np.array([eq[0], eq[1], eq[2], eq[3], eq[4], eq[5]])
'''
C1 = equation_to_conic(eq1)
C2 = equation_to_conic(eq2)
return transform(C1, C2)
# realizing the matrix U as a projective transformation
def proj(U, x):
if len(x) == 2:
x = homogenize(x)
y = U.dot(x)
y = y / y[2]
return y[0:2]
# find the common points, i.e. points of intersection of a pair of conics
# represented by a pair of symmetric matrices
def find_common_points(c1, c2):
k, U = transform(c1, c2)
L1 = (U.T).dot((c1 - k[0]*c2).dot(U))
L2 = (U.T).dot((c1 - k[1]*c2).dot(U))
sol = np.empty((4,3), dtype=float)
for i in range(2):
for j in range(2):
sol[i+2*j,0:2] = np.array([math.sqrt(abs(L2[2,2] / L2[0,0]))*(-1)**i, math.sqrt(abs(L1[2,2] / L1[1,1]))*(-1)**j])
sol[i+2*j,0:2] = proj(U, sol[i+2*j,0:2])
sol[:,2] = np.ones(4)
return sol
# find the solutions, i.e. the points x=(x[0],x[1]) saisfying the pair
# of quadratic equations
# represented by a pair of vectors eq1 and eq2 of 6 coefficients
def solve_eq(eq1, eq2):
conic1 = equation_to_conic(eq1)
conic2 = equation_to_conic(eq2)
return find_common_points(conic1, conic2)
'''
Esample of finding the common tangents of a pair of conics:
conic 1: major axis = 2, minor axis = 1, angle = 45, center = (0,0)
conic 2: major axis = 3, minor axis = 1, angle = 120, center = (15,0)
'''
a = 2
b = 1
cntr = np.array([0,0])
w = 45
Q1 = Conic(a, b, w, cntr)
dQ1 = dual_Conic(a, b, w, cntr)
a = 3
b = 1
cntr = np.array([15,0])
w = 120
Q2 = Conic(a, b, w, cntr)
dQ2 = dual_Conic(a, b, w, cntr)
R = find_common_points(dQ1, dQ2)
print('')
print(R)
print('')
print('checking that the output forms common tangent lines: ')
print('')
print('conic 1: ')
print(np.diagonal(R.dot(dQ1.dot(R.T))) )
print('')
print('conic 2: ')
print(np.diagonal(R.dot(dQ2.dot(R.T))) )
#conic_to_equation(dQ1)
一些解释：假设你想找到两个二次曲线C1和C2的交点。为简单起见，我们假设它们是相交于四个不同点的真实椭圆（以避免复数）
在寻找一对二次曲线的公切线的情况下，
只需将两个二次曲线转换为两个对应的对偶，然后找到交点
双胞胎的。这些交点就是原二次曲线的切线的方程系数。
这个问题可能有几种不同的几何解释，但让我们用圆锥曲线来解释。两个二次曲线 C1 和 C2 由具有非零行列式的 3×3 对称矩阵表示，我已将其表示为 C1 和 C2。线性组合，称为由 C1 和 C2 生成的圆锥曲线，是 t 参数化的圆锥曲线族 C1 - t*C2 ，其中 t 只是一个数字。关键是每个圆锥曲线 C1 - tC2 都通过 C1 和 C2 的交点，这是它们仅有的四个共同点。您可以通过观察 if x.T * C1 * x = x.T * C1 * x = 0 then 来证明这一点
x.T * (C1 - t*C2) * x = x.T * C1 * x - t * x.T * C2 * x = 0 - t*0 = 0。此外，如果取 C1 和 C1 - t*C2 的交点，则 C2 = C1 - t*C2 + s*C2 可以在 s = t.
时应用相同的参数
在这个二次曲线族中，有三个退化二次曲线，它们在几何上是三对直线。它们恰好在 t 满足 det( C1 - t*C2 ) = 0 时发生。这是关于 t 的 3 次多项式方程，因此由于 C1 和 C2 的良好性，存在三个不同的解 k[0], k[1], k[2],。从投影上讲，每个退化二次曲线 C1 - k[j]*C2 是一对直线，它们有一个公共交点 u[:,j] = [ u[0,j] : u[1,j] : u[2,j] ]。而且，rank(C1 - k[j]*C2) = 2，所以ker(C1 - k[j]*C2) = 1。这个点u[:,j]被表征为方程的解
(C1 - k[j]*C2) * u[:,j] = 0。
由于C2是可逆的（非退化的），将方程两边乘以inverse(C2)得到等价的方程( (inverse(C2) * C1) - k[j]*Identity ) * u[:,j] = 0是一个特征值方程，以k[j]为特征值u[:,j] 作为特征向量。函数 transform() 的输出是特征值的 1 x 3 数组 k 和特征向量的 3 x 3 矩阵 U = [ u[:,0], u[:,1], u[:,2] ]。
共轭 C1 - k[j]*C2 by U, i.e. (U.T)*(C1 - k[j]*C2)*U 在几何上等同于执行投影变换，将 u[:,0] 和 u[:,1] 发送到无穷远，将 u[:,2] 发送到原点。因此，U 将坐标系从一般坐标系更改为特殊坐标系，其中两个退化二次曲线由成对的平行线给出，并将它们相交在一个矩形中。第三个圆锥曲线由矩形的 diagnols 表示。
在这张新图中，相交问题很容易解决，只需从矩阵的条目中读取一个解决方案即可(U.T)*(C1 - k[j]*C2)*U（交点是矩形的顶点，所以如果你找到其中一个，其他只是彼此镜像对称）。
我一直在尝试求出两个旋转椭圆的公切线。我正在按照 Edward Doolittle in the
矩阵形式的椭圆可以表示为：
第一个椭圆以 (0,0) 为中心旋转 45 度，半长轴和半短轴长度为 2,1。第二个椭圆以 (15,0) 为中心，旋转 120 度，半长轴和半短轴长度为 3,1
两个椭圆的伴随矩阵的线性组合是两个椭圆组合的对偶
我得到这个值
然后我试图找到 t 的值，这将使圆锥曲线（矩阵上方）退化。
我发现 t 的值为 (-0.05,0.29,2.46)。但是，当我将这些值放回上述矩阵时，我无法将矩阵简化为两个变量形式。我总是在处理 3 个变量。例如，如果我输入 t = -0.05，那么我会得到以下结果：
有人可以帮我解决这个问题吗？
它归结为找到一种算法来求解两个变量的两个二次方程组，通过将其解释为二次曲线的射影几何铅笔，然后找到铅笔的三个退化二次曲线以及射影变换将这三个退化二次曲线简化到可以在新的简化坐标系中非常容易地读出解，然后将它们转换回原始坐标系。
我在python中勾勒了一个算法，我认为它似乎适用于你的例子......但是我很匆忙并且没有正确检查它，所以可能存在错误......
import numpy as np
import math
# technical module, functions, details
def homogenize(x):
return np.array([x[0], x[1], 1])
def cos_sin(angle_deg):
return math.cos(angle_deg*math.pi/180), math.sin(angle_deg*math.pi/180)
def rotation(cs_sn):
return np.array([[cs_sn[0], -cs_sn[1]],
[cs_sn[1], cs_sn[0]]])
# defining the isometry (the rotation plus translation) transformation
# between the coordinate system aligned with the conic and a general (world)
# coordinate system
def isom_inverse(angle, translation):
'''
isometry from conic-aligned coordinate system (conic attached)
to global coordinate system (world system)
'''
cos_, sin_ = cos_sin(angle)
return np.array([[cos_, -sin_, translation[0]],
[sin_, cos_, translation[1]],
[ 0, 0, 1]])
def isom(angle, translation):
'''
isometry from global coordinate system (world system)
to conic-aligned coordinate system (conic attached)
'''
cos_, sin_ = cos_sin(-angle)
tr = - rotation((cos_, sin_)).dot(translation)
return np.array([[ cos_, -sin_, tr[0]],
[ sin_, cos_, tr[1]],
[ 0, 0, 1 ]])
# calculating the coinc defined by a pair of axes' lengts,
# axes rotation angle and center of the conic
def Conic(major, minor, angle, center):
D = np.array([[minor**2, 0, 0],
[ 0, major**2, 0],
[ 0, 0, -(minor*major)**2]])
U = isom(angle, center)
return (U.T).dot(D.dot(U))
# calculating the coinc dual to the conic defined by a pair of axes' lengths,
# axes rotation angle and center of the conic
def dual_Conic(major, minor, angle, center):
D_1 = np.array([[major**2, 0, 0],
[ 0, minor**2, 0],
[ 0, 0, -1]])
U_1 = isom_inverse(angle, center)
return (U_1).dot(D_1.dot(U_1.T))
# transforming the matrix of a conic into a vector of six coefficients
# of a quadratic equation with two variables
def conic_to_equation(C):
'''
c[0]*x**2 + c[1]*x*y + c[2]*y**2 + c[3]*x + c[4]*y + c[5] = 0
'''
return np.array([C[0,0], 2*C[0,1], C[1,1], 2*C[0,2], 2*C[1,2], C[2,2]])
# transforming the vector of six coefficients
# of a quadratic equation with two variables into a matrix of
# the corresponding conic
def equation_to_conic(eq):
'''
eq[0]*x**2 + eq[1]*x*y + eq[2]*y**2 + eq[3]*x + eq[4]*y + eq[5] = 0
'''
return np.array([[2*eq[0], eq[1], eq[3]],
[ eq[1], 2*eq[2], eq[4]],
[ eq[3], eq[4], 2*eq[5]]]) / 2
# given a point (x,y) define the vector (x^2, xy, y^2, x, y, 1)
def argument(x):
return np.array([x[0]**2, x[0]*x[1], x[1]**2, x[0], x[1], 1])
# given x = (x[0],x[1]) calculate the value of the quadratic equation with
# six coefficients coeff
def quadratic_equation(x, coeff):
'''
coeff[0]*x**2 + coeff[1]*x*y + coeff[2]*y**2 + coeff[3]*x + coeff[4]*y + coeff[5] = 0
'''
return coeff.dot( argument(x) )
# given a pair of conics, as a pair of symmetric matrices,
# calculate the vector k = (k[0], k[1], k[2]) of values for each of which
# the conic c1 - k[i]*c2 from the pencil of conics c1 - t*c2
# is a degenerate conic (the anti-symmetric product of a pair of linear forms)
# and also find the matrix U
# of the projective transformation that simplifies the geometry of
# the pair of conics, the geometry of the pencil c1 - t*c2 in general,
# as well as the geometry of the three degenerate conics in particular
def transform(c1, c2):
'''
c1 and c2 are 3 by 3 symmetric matrices of the two conics
'''
c21 = np.linalg.inv(c2).dot(c1)
k, U = np.linalg.eig(c21)
return k, U
# the same as before, but for a pair of equations instead of matrices of conics
def eq_transform(eq1, eq2):
'''
eq1 and eq2 = np.array([eq[0], eq[1], eq[2], eq[3], eq[4], eq[5]])
'''
C1 = equation_to_conic(eq1)
C2 = equation_to_conic(eq2)
return transform(C1, C2)
# realizing the matrix U as a projective transformation
def proj(U, x):
if len(x) == 2:
x = homogenize(x)
y = U.dot(x)
y = y / y[2]
return y[0:2]
# find the common points, i.e. points of intersection of a pair of conics
# represented by a pair of symmetric matrices
def find_common_points(c1, c2):
k, U = transform(c1, c2)
L1 = (U.T).dot((c1 - k[0]*c2).dot(U))
L2 = (U.T).dot((c1 - k[1]*c2).dot(U))
sol = np.empty((4,3), dtype=float)
for i in range(2):
for j in range(2):
sol[i+2*j,0:2] = np.array([math.sqrt(abs(L2[2,2] / L2[0,0]))*(-1)**i, math.sqrt(abs(L1[2,2] / L1[1,1]))*(-1)**j])
sol[i+2*j,0:2] = proj(U, sol[i+2*j,0:2])
sol[:,2] = np.ones(4)
return sol
# find the solutions, i.e. the points x=(x[0],x[1]) saisfying the pair
# of quadratic equations
# represented by a pair of vectors eq1 and eq2 of 6 coefficients
def solve_eq(eq1, eq2):
conic1 = equation_to_conic(eq1)
conic2 = equation_to_conic(eq2)
return find_common_points(conic1, conic2)
'''
Esample of finding the common tangents of a pair of conics:
conic 1: major axis = 2, minor axis = 1, angle = 45, center = (0,0)
conic 2: major axis = 3, minor axis = 1, angle = 120, center = (15,0)
'''
a = 2
b = 1
cntr = np.array([0,0])
w = 45
Q1 = Conic(a, b, w, cntr)
dQ1 = dual_Conic(a, b, w, cntr)
a = 3
b = 1
cntr = np.array([15,0])
w = 120
Q2 = Conic(a, b, w, cntr)
dQ2 = dual_Conic(a, b, w, cntr)
R = find_common_points(dQ1, dQ2)
print('')
print(R)
print('')
print('checking that the output forms common tangent lines: ')
print('')
print('conic 1: ')
print(np.diagonal(R.dot(dQ1.dot(R.T))) )
print('')
print('conic 2: ')
print(np.diagonal(R.dot(dQ2.dot(R.T))) )
#conic_to_equation(dQ1)
一些解释：假设你想找到两个二次曲线C1和C2的交点。为简单起见，我们假设它们是相交于四个不同点的真实椭圆（以避免复数）
在寻找一对二次曲线的公切线的情况下， 只需将两个二次曲线转换为两个对应的对偶，然后找到交点 双胞胎的。这些交点就是原二次曲线的切线的方程系数。
这个问题可能有几种不同的几何解释，但让我们用圆锥曲线来解释。两个二次曲线 C1 和 C2 由具有非零行列式的 3×3 对称矩阵表示，我已将其表示为 C1 和 C2。线性组合，称为由 C1 和 C2 生成的圆锥曲线，是 t 参数化的圆锥曲线族 C1 - t*C2 ，其中 t 只是一个数字。关键是每个圆锥曲线 C1 - tC2 都通过 C1 和 C2 的交点，这是它们仅有的四个共同点。您可以通过观察 if x.T * C1 * x = x.T * C1 * x = 0 then 来证明这一点
x.T * (C1 - t*C2) * x = x.T * C1 * x - t * x.T * C2 * x = 0 - t*0 = 0。此外，如果取 C1 和 C1 - t*C2 的交点，则 C2 = C1 - t*C2 + s*C2 可以在 s = t.
在这个二次曲线族中，有三个退化二次曲线，它们在几何上是三对直线。它们恰好在 t 满足 det( C1 - t*C2 ) = 0 时发生。这是关于 t 的 3 次多项式方程，因此由于 C1 和 C2 的良好性，存在三个不同的解 k[0], k[1], k[2],。从投影上讲，每个退化二次曲线 C1 - k[j]*C2 是一对直线，它们有一个公共交点 u[:,j] = [ u[0,j] : u[1,j] : u[2,j] ]。而且，rank(C1 - k[j]*C2) = 2，所以ker(C1 - k[j]*C2) = 1。这个点u[:,j]被表征为方程的解
(C1 - k[j]*C2) * u[:,j] = 0。
由于C2是可逆的（非退化的），将方程两边乘以inverse(C2)得到等价的方程( (inverse(C2) * C1) - k[j]*Identity ) * u[:,j] = 0是一个特征值方程，以k[j]为特征值u[:,j] 作为特征向量。函数 transform() 的输出是特征值的 1 x 3 数组 k 和特征向量的 3 x 3 矩阵 U = [ u[:,0], u[:,1], u[:,2] ]。
共轭 C1 - k[j]*C2 by U, i.e. (U.T)*(C1 - k[j]*C2)*U 在几何上等同于执行投影变换，将 u[:,0] 和 u[:,1] 发送到无穷远，将 u[:,2] 发送到原点。因此，U 将坐标系从一般坐标系更改为特殊坐标系，其中两个退化二次曲线由成对的平行线给出，并将它们相交在一个矩形中。第三个圆锥曲线由矩形的 diagnols 表示。
在这张新图中，相交问题很容易解决，只需从矩阵的条目中读取一个解决方案即可(U.T)*(C1 - k[j]*C2)*U（交点是矩形的顶点，所以如果你找到其中一个，其他只是彼此镜像对称）。