我一直试图让两个旋转椭圆的公切线。我在以下线程中遵循Edward Doolittle给出的方法。这两个椭圆以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 x 3 对称矩阵表示，我将其表示为 C1 和 C2。由 C1 和 C2 生成的称为圆锥曲线的线性组合是 t 参数化的圆锥曲线族C1 - t*C2，其中 t 只是一个数字。关键是每个圆锥都C1 - tC2通过 C1 和 C2 的交点，这是它们唯一的四个共同点。您可以通过观察 if x.T * C1 * x = x.T * C1 * x = 0then
来证明这一点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 次多项式方程，因此k[0], k[1], k[2],由于 C1 和 C2 的好坏，存在三个不同的解。投影而言，每个退化二次曲线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 3U = [ 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将坐标系从一般坐标系更改为特殊坐标系，其中两个退化圆锥曲线由平行线对给出，并将它们组合成一个矩形。第三个圆锥由矩形的诊断表示。
在这张新图片中，交集问题可以很容易地解决，只需从矩阵的条目中读取一个解(U.T)*(C1 - k[j]*C2)*U（交点是矩形的顶点，所以如果你找到其中一个，其他的只是镜像对称彼此的）。