Intersections between Geodesics (shortest distance paths) on the surface of a sphere

Question

I've searched far and wide but have yet to find a suitable answer to this problem. Given two lines on a sphere, each defined by their start and end points, determine whether or not and where they intersect. I've found this site (http://mathforum.org/library/drmath/view/62205.html) which runs through a good algorithm for the intersections of two great circles, although I'm stuck on determining whether the given point lies along the finite section of the great circles.

I've found several sites which claim they've implemented this, Including some questions here and on stackexchange, but they always seem to reduce back to the intersections of two great circles.

The python class I'm writing is as follows and seems to almost work:

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: I forgot to specify that if two curves are determined to intersect, I then need to have the point of intersection.

Answer 1

Intersection using plane trig can be calculated using the below code in UBasic.

5   'interx.ub adapted from code at

6   'https://rosettacode.org
7   '/wiki/Find_the_intersection_of_two_linesSinclair_ZX81_BASIC

8  'In U Basic by yuji kida https://en.wikipedia.org/wiki/UBASIC

10   XA=48.7815144526:'669595.708

20   YA=-117.2847245001:'2495736.332

30   XB=48.7815093807:'669533.412

40   YB=-117.2901673467:'2494425.458

50   XC=48.7824947147:'669595.708

60   YC=-117.28751374:'2495736.332

70   XD=48.77996737:'669331.214

80   YD=-117.2922957:'2494260.804

90   print "THE TWO LINES ARE:"

100   print "YAB=";YA-XA*((YB-YA)/(XB-XA));"+X*";((YB-YA)/(XB-XA))

110   print "YCD=";YC-XC*((YD-YC)/(XD-XC));"+X*";((YD-YC)/(XD-XC))

120   X=((YC-XC*((YD-YC)/(XD-XC)))-(YA-XA*((YB-YA)/(XB-XA))))/(((YB-YA)/(XB-XA))-((YD-YC)/(XD-XC)))

130   print "Lat =  ";X

140   Y=YA-XA*((YB-YA)/(XB-XA))+X*((YB-YA)/(XB-XA))

150   print "Lon = ";Y

160   'print "YCD=";YC-XC*((YD-YC)/(XD-XC))+X*((YD-YC)/(XD-XC))

Answer 2

A simpler approach is to express the problem in terms of geometric primitive operations like the dot product, the cross product, and the triple product. The sign of the determinant of u, v, and w tells you which side of the plane spanned by v and w contains u. This enables us to detect when two points are on opposite sites of a plane. That's equivalent to testing whether a great circle segment crosses another great circle. Performing this test twice tells us whether two great circle segments cross each other.

The implementation requires no trigonometric functions, no division, no comparisons with pi, and no special behavior around the poles!

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))))

If you want to initialize unit vectors in terms of angles theta and phi, you can do that, but I recommend immediately converting to Cartesian (x, y, z) coordinates to perform all subsequent calculations.