球面上测地线(最短距离路径)的交点

4
我已经搜索了很多地方,但仍然没有找到一个合适的答案来解决这个问题。给定两条线在球面上,每条线由其起始点和结束点定义,确定它们是否相交以及它们相交的位置。我发现了这个网站(http://mathforum.org/library/drmath/view/62205.html),它介绍了一个很好的算法来计算两个大圆的交点,但我无法确定给定点是否位于大圆弧的有限部分上。
我找到了一些声称已经实现了这个功能的网站,包括一些在此处和stackexchange上的问题,但它们总是回到计算两个大圆的交点上。
我正在编写的python类如下,似乎几乎可以工作:
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)  

编辑:我忘记指定了,如果确定两条曲线相交,则需要得到相交点。

备注:本段内容涉及it技术。请仔细查看并确保翻译准确无误,并在翻译过程中保留所有html标签。
2个回答

10
一种更简单的方法是使用几何基本操作来表达问题,例如点积叉积三重积uvw的行列式符号告诉您由vw张成的平面的哪一侧包含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))))

如果你想用角度θ和φ来初始化单位向量,是可以做到的,但我建议立即转换为笛卡尔坐标系(x,y,z)进行所有后续计算。

2
我对这个解决方案持怀疑态度。这难道不是在欧几里得坐标系中计算交叉吗?但我们的坐标表示球面上的点。例如,线段(0,0)到(50,50)在欧几里得坐标系中不会与(26,25)到(80,-30)相交,但如果将这些点解释为(纬度,经度),它们在球面坐标系中会相交。 - Jim Newton
2
@JoeAmenta 对我来说它返回True:def polar(lat, long): return Vector(math.cos(long)*math.cos(lat), math.sin(long)*math.cos(lat), math.sin(lat)) deg = math.pi / 180 intersects(Pair(polar(0,0), polar(50*deg,50*deg)), Pair(polar(26*deg,25*deg), polar(80*deg,-30*deg))) True - Chris Culter
2
@TomMcLean 很高兴分享链接。希望它有所帮助。 https://gitlab.lrde.epita.fr/jnewton/earthmap/-/blob/041c4d262cec19d4bd7796842f7e0b66e6e0c372/globe/src/main/scala/globe/Earth.scala#L73 - Jim Newton
1
代码有点问题。它的功能是: - undefined
1
它的功能: 这两对点跨越两个大圆,通常在两个交点处相交,确定了4个半圆。代码检查这四个点是否位于四个不同的半圆上。 这是判断测地线是否相交的必要条件,但不是充分条件。还需要检查测地线是否击中了半圆的相同交点。 - undefined
显示剩余9条评论

-3

使用平面三角函数可以在 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))

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接