Пересечения между геодезическими (кратчайшими дорожками) на поверхности сферы

Question

Я искал далеко и широко, но пока не нашел подходящего ответа на эту проблему. Учитывая две линии на сфере, каждая из которых определяется их начальной и конечной точками, определяют, пересекаются или нет. Я нашел этот сайт (http://mathforum.org/library/drmath/view/62205.html), который проходит через хороший алгоритм для пересечений двух больших кругов, хотя я застрял на определении того, лежит ли данная точка вдоль конечного участка больших кругов.

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

Класс питон Я пишу это так и кажется, что почти работа:

class Geodesic(Boundary): 
    def _SecondaryInitialization(self): 
    self.theta_1 = self.point1.theta 
    self.theta_2 = self.point2.theta 
    self.phi_1 = self.point1.phi 
    self.phi_2 = self.point2.phi 

    sines = math.sin(self.phi_1) * math.sin(self.phi_2) 
    cosines = math.cos(self.phi_1) * math.cos(self.phi_2) 
    self.d = math.acos(sines - cosines * math.cos(self.theta_2 - self.theta_1)) 

    self.x_1 = math.cos(self.theta_1) * math.cos(self.phi_1) 
    self.x_2 = math.cos(self.theta_2) * math.cos(self.phi_2) 
    self.y_1 = math.sin(self.theta_1) * math.cos(self.phi_1) 
    self.y_2 = math.sin(self.theta_2) * math.cos(self.phi_2) 
    self.z_1 = math.sin(self.phi_1) 
    self.z_2 = math.sin(self.phi_2) 

    self.theta_wraps = (self.theta_2 - self.theta_1 > PI) 
    self.phi_wraps = ((self.phi_1 < self.GetParametrizedCoords(0.01).phi and 
     self.phi_2 < self.GetParametrizedCoords(0.99).phi) or (
     self.phi_1 > self.GetParametrizedCoords(0.01).phi) and 
     self.phi_2 > self.GetParametrizedCoords(0.99)) 

    def Intersects(self, boundary): 
    A = self.y_1 * self.z_2 - self.z_1 * self.y_2 
    B = self.z_1 * self.x_2 - self.x_1 * self.z_2 
    C = self.x_1 * self.y_2 - self.y_1 * self.x_2 
    D = boundary.y_1 * boundary.z_2 - boundary.z_1 * boundary.y_2 
    E = boundary.z_1 * boundary.x_2 - boundary.x_1 * boundary.z_2 
    F = boundary.x_1 * boundary.y_2 - boundary.y_1 * boundary.x_2 

    try: 
     z = 1/math.sqrt(((B * F - C * E) ** 2/(A * E - B * D) ** 2) 
      + ((A * F - C * D) ** 2/(B * D - A * E) ** 2) + 1) 
    except ZeroDivisionError: 
     return self._DealWithZeroZ(A, B, C, D, E, F, boundary) 

    x = ((B * F - C * E)/(A * E - B * D)) * z 
    y = ((A * F - C * D)/(B * D - A * E)) * z 

    theta = math.atan2(y, x) 
    phi = math.atan2(z, math.sqrt(x ** 2 + y ** 2)) 

    if self._Contains(theta, phi): 
     return point.SPoint(theta, phi) 

    theta = (theta + 2* PI) % (2 * PI) - PI 
    phi = -phi 

    if self._Contains(theta, phi): 
     return spoint.SPoint(theta, phi) 

    return None 

    def _Contains(self, theta, phi): 
    contains_theta = False 
    contains_phi = False 

    if self.theta_wraps: 
     contains_theta = theta > self.theta_2 or theta < self.theta_1 
    else: 
     contains_theta = theta > self.theta_1 and theta < self.theta_2 

    phi_wrap_param = self._PhiWrapParam() 
    if phi_wrap_param <= 1.0 and phi_wrap_param >= 0.0: 
     extreme_phi = self.GetParametrizedCoords(phi_wrap_param).phi 
     if extreme_phi < self.phi_1: 
     contains_phi = (phi < max(self.phi_1, self.phi_2) and 
      phi > extreme_phi) 
     else: 
     contains_phi = (phi > min(self.phi_1, self.phi_2) and 
      phi < extreme_phi) 
    else: 
     contains_phi = (phi > min(self.phi_1, self.phi_2) and 
      phi < max(self.phi_1, self.phi_2)) 

    return contains_phi and contains_theta 

    def _PhiWrapParam(self): 
    a = math.sin(self.d) 
    b = math.cos(self.d) 
    c = math.sin(self.phi_2)/math.sin(self.phi_1) 
    param = math.atan2(c - b, a)/self.d 

    return param 

    def _DealWithZeroZ(self, A, B, C, D, E, F, boundary): 
    if (A - D) is 0: 
     y = 0 
     x = 1 
    elif (E - B) is 0: 
     y = 1 
     x = 0 
    else: 
     y = 1/math.sqrt(((E - B)/(A - D)) ** 2 + 1) 
     x = ((E - B)/(A - D)) * y 

    theta = (math.atan2(y, x) + PI) % (2 * PI) - PI 
    return point.SPoint(theta, 0) 

def GetParametrizedCoords(self, param_value): 
    A = math.sin((1 - param_value) * self.d)/math.sin(self.d) 
    B = math.sin(param_value * self.d)/math.sin(self.d) 

    x = A * math.cos(self.phi_1) * math.cos(self.theta_1) + (
    B * math.cos(self.phi_2) * math.cos(self.theta_2)) 
    y = A * math.cos(self.phi_1) * math.sin(self.theta_1) + (
     B * math.cos(self.phi_2) * math.sin(self.theta_2)) 
    z = A * math.sin(self.phi_1) + B * math.sin(self.phi_2) 

    new_phi = math.atan2(z, math.sqrt(x**2 + y**2)) 
    new_theta = math.atan2(y, x) 

    return point.SPoint(new_theta, new_phi) 

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

Answer 1

Простой подход заключается в том, чтобы выразить проблему с точки зрения геометрических примитивных операций, таких как dot product, cross product и triple product. Знак определителя у, v и ж говорит вам, с какой стороны плоскости, натянутой на об и ш содержит ˙U. Это позволяет нам обнаружить, когда две точки находятся на противоположных участках плоскости. Это эквивалентно проверке того, пересекает ли сегмент большого круга другой большой круг. Выполнение этого теста дважды говорит нам, пересекаются ли два больших круга друг с другом.

Реализация не требует тригонометрических функций, без деления, сравнений с pi и никакого специального поведения вокруг полюсов!

class Vector: 
    def __init__(self, x, y, z): 
     self.x = x 
     self.y = y 
     self.z = z 

def dot(v1, v2): 
    return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z 

def cross(v1, v2): 
    return Vector(v1.y * v2.z - v1.z * v2.y, 
        v1.z * v2.x - v1.x * v2.z, 
        v1.x * v2.y - v1.y * v2.x) 

def det(v1, v2, v3): 
    return dot(v1, cross(v2, v3)) 

class Pair: 
    def __init__(self, v1, v2): 
     self.v1 = v1 
     self.v2 = v2 

# Returns True if the great circle segment determined by s 
# straddles the great circle determined by l 
def straddles(s, l): 
    return det(s.v1, l.v1, l.v2) * det(s.v2, l.v1, l.v2) < 0 

# Returns True if the great circle segments determined by a and b 
# cross each other 
def intersects(a, b): 
    return straddles(a, b) and straddles(b, a) 

# Test. Note that we don't need to normalize the vectors. 
print(intersects(Pair(Vector(1, 0, 1), Vector(-1, 0, 1)), 
       Pair(Vector(0, 1, 1), Vector(0, -1, 1)))) 

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