找到一对圆锥曲线的内公切线

Question

我一直在尝试获得两个旋转椭圆的公切线。我按照Edward Doolittle在following thread中给出的方法。这两个椭圆以Wiki 中给出的等式给出。

在矩阵形式中，椭圆可以显示为：

第一个椭圆以 (0,0) 为中心旋转 45 度，半长轴和半短轴长度为 2,1。第二个椭圆以 (15,0) 为中心，旋转 120 度，半长轴和半短轴长度为 3,1

两个椭圆的伴随矩阵的线性组合是两个椭圆组合的每对偶
我得到了这个值
.

然后我试图找到 t 的值，它会使圆锥（矩阵上方）退化。

我发现 t 的值为 (-0.05,0.29,2.46)。但是，当我将这些值放回上述矩阵时，我无法将矩阵简化为两个变量的形式。我总是在处理 3 个变量。例如，如果我输入 t = -0.05，那么我会得到以下结果：

有人可以帮我解决这个问题吗？

Answer 1

它归结为找到一种算法来求解一个由两个变量组成的二次方程组，通过将其解释为圆锥曲线的投影几何铅笔，然后找到铅笔的三个退化圆锥曲线以及简化这些曲线的投影变换三个退化的圆锥曲线，您可以在新的简化坐标系中轻松读取解，然后将它们转换回原始坐标系。

我在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 = 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×3数组k和特征向量的3×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 将坐标系从一般坐标系更改为特殊坐标系，其中两个退化圆锥曲线由平行线对给出，并将它们组合成一个矩形。第三个圆锥由矩形的诊断表示。

在这张新图片中，相交问题可以很容易地解决，只需从矩阵(U.T)*(C1 - k[j]*C2)*U 的条目中读取一个解（相交点是矩形的顶点，所以如果你找到其中一个，其他的只是彼此镜像对称）。