Нахождение внутренней общей касательной к паре коник

Он сводится к поиску алгоритма решения системы двух квадратных уравнений с двумя переменными путем интерпретации ее как пучка коник проективной геометрии, а затем к нахождению трех вырожденных коник пучка вместе с проективным преобразованием, упрощающим эти три вырожденных коники. до такой степени, что вы сможете очень легко считывать решения в новой упрощающей системе координат, а затем преобразовывать их обратно в исходную систему координат.

Набросал алгоритм на питоне, думаю на вашем примере вроде работает... но торопился и не проверил как следует, так что могут быть баги...

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, и это единственные четыре общие точки. Вы можете доказать это, заметив, что если x.T * C1 * x = x.T * C1 * x = 0, то 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. Это полиномиальное уравнение степени 3 относительно t, поэтому существует три, скажем, разных решения 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 на 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 (точки пересечения являются вершинами прямоугольника, поэтому, если вы найдете одну из них, другие просто зеркально симметричны друг другу).