我已经进行了广泛的搜索，但还没有找到这个问题的合适答案。给定球体上的两条线，每条线由它们的起点和终点定义，确定它们是否相交以及在哪里相交。我发现这个站点 ( 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)
编辑:我忘了指定如果确定两条曲线相交，那么我需要有交点。
最佳答案
一种更简单的方法是用几何基本运算(如 dot product 、 cross product 和 triple product )来表达问题。 u、v 和 w 的行列式的符号告诉您 v 和 w 跨越的平面的哪一侧包含 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))))
如果您想根据角度 theta 和 phi 初始化单位向量，您可以这样做，但我建议立即转换为笛卡尔 (x, y, z) 坐标以执行所有后续计算。
关于python - 球体表面测地线(最短距离路径)之间的交点，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/26668041/