球面上的测地线(最短距离路径)之间的交点 [英] Intersections between Geodesics (shortest distance paths) on the surface of a sphere

问题描述

我已经进行了广泛的搜索，但是还没有找到解决该问题的合适方法.给定球体上的两条线(分别由它们的起点和终点定义)，确定它们是否相交以及在何处相交.我找到了这个网站( http://mathforum.org/library/drmath/view/62205.html )，它通过一个很好的算法来计算两个大圆的交点，尽管我坚持确定给定点是否位于大圆的有限部分上.

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.

我发现有几个站点声称已实现此目的，包括此处和关于stackexchange的一些问题，但它们似乎总是减少到两个大圆圈的交集.

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.

我正在编写的python类如下，并且似乎几乎可以正常工作:

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)  

我忘记指定如果确定两条曲线要相交，那么我需要有交点.

I forgot to specify that if two curves are determined to intersect, I then need to have the point of intersection.

推荐答案

一种更简单的方法是按照几何基本操作(例如点产品，跨产品和三重产品. u ， v 和 w 的行列式的符号告诉您 v 跨越平面的哪一侧并且 w 包含 u .这使我们能够检测到两个点何时位于平面的相对位置.这等效于测试一个大圆弧段是否与另一个大圆弧相交.进行两次测试可以告诉我们两个大圆弧段是否彼此交叉.

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.

该实现不需要三角函数，除法，与pi的比较以及极点附近的特殊行为！

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

如果要根据角度theta和phi初始化单位矢量，可以这样做，但是我建议立即转换为笛卡尔(x，y，z)坐标以执行所有后续计算.

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

这篇关于球面上的测地线(最短距离路径)之间的交点的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！