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

2024-04-19

我进行了广泛的搜索,但尚未找到该问题的合适答案。给定球体上的两条线,每条线由起点和终点定义,确定它们是否相交以及相交的位置。我找到了这个网站(http://mathforum.org/library/drmath/view/62205.html 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)  

编辑:我忘记指定如果两条曲线确定相交,那么我需要有交点。


一种更简单的方法是用几何基元运算来表达问题,例如点积 https://en.wikipedia.org/wiki/Dot_product, the 叉积 https://en.wikipedia.org/wiki/Cross_product,以及三重积 https://en.wikipedia.org/wiki/Triple_product。行列式的符号u, v, and w告诉你飞机的哪一侧跨越v and 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) 坐标以执行所有后续计算。

本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

球体表面上测地线(最短距离路径)之间的交点 的相关文章

  • 我正在用 python 编写一个电报机器人

    我想通过Python编写一个电报机器人 但它不起作用 import telebot bot telebot TeleBot my token bot message handler content types text def sendin
  • Python Numpy TypeError:输入类型不支持 ufunc 'isfinite'

    这是我的代码 def topK dataMat sensitivity meanVals np mean dataMat axis 0 meanRemoved dataMat meanVals covMat np cov meanRemov
  • 将 python scikit learn 模型导出到 pmml

    我想将 python scikit learn 模型导出到 PMML 中 什么 python 包最适合 我读到Augustus https github com opendatagroup augustus 但我找不到任何使用 scikit
  • 如何在 for 循环中跳过一些迭代

    在 python 中 我通常简单地循环遍历范围 for i in range 100 do something 但现在我想跳过循环中的几个步骤 更具体地说 我想要类似的东西continue 10 这样它就会跳过整个循环并将计数器增加 10
  • 使用 Python 连接从 FTP 检索文件

    我构建了这个简单的工具来暴力破解并连接到 ftp 服务器 import socket import ftplib from ftplib import FTP port 21 ip 192 168 1 108 file1 passwords
  • python 正则表达式中括号的奇怪行为

    我正在编写一个 python 正则表达式 它可以在文本文档中查找引用的字符串 从黑匣子中记录的航空公司飞行员的引用 我首先尝试编写具有以下规则的正则表达式 返回引号之间的内容 如果以 single 打开 则仅在以 single 关闭时返回
  • 为什么通过selenium切换到alert不稳定?

    为什么通过selenium切换到alert不稳定 例如 1 运行代码 一切顺利 一切都很顺利 但如果这段代码在几分钟内运行 那么可能会出现错误 例如 没有可以单击的元素 等等 2 在一个站点上有一个警报窗口 alert driver swi
  • 使用Python mysql.connector远程连接MySQL

    以下代码 在同一 LAN 内与 mysql 服务器不同的机器上运行 使用 Python3 和 mysql connector 本地连接到 MySQL 数据库 import mysql connector cnx mysql connecto
  • 如何从 python 脚本更改 python 文件中的变量值

    我目前有一个 python 文件 其中包含一堆带有值的全局变量 我想从一个单独的 python 脚本永久更改这些值 我尝试过 setattr 等 但似乎不起作用 有没有办法做到这一点 简短的回答是 不 不值得这么麻烦 听起来您正在尝试创建一
  • Matplotlib 动画未显示

    当我在家里的电脑上尝试这个时 它可以工作 但在工作的电脑上却不行 这是代码 import numpy as np import matplotlib pyplot as plt import matplotlib animation as
  • 多级QTreeView

    我很难理解如何使用 QTreeView 和 QStandardItemModel 设置多级 QTree 这是我所拥有的 from PySide QtGui import import sys class MainFrame QWidget
  • 如何将 NaN 数组插入 numpy 二维数组

    我试图在二维数组中的特定位置插入任意数量的 NaN 值行 我正在将来自微控制器的一些数据记录在 csv 文件中并使用 python 进行解析 数据存储在 3 列 2D 数组中 如下所示 122 0 1 0 47 0 123 0 1 0 47
  • scipy 的 curve_fit 函数的尺寸问题

    我对 python 中的曲线拟合以及一般的 python 都很陌生 目前 我正在尝试使用 scipy 中的 curve fit 模块来拟合 4 个光谱峰 简而言之 我的文本文件中有两列数据 所以我的第一步是将数据导入到两个数组中 一个包含
  • Python lmfit:拟合 2D 模型

    我正在尝试将二维高斯拟合到一些灰度图像数据 该数据由一个二维数组给出 lmfit 库实现了一个易于使用的模型类 它应该能够做到这一点 不幸的是文档 http lmfit github io lmfit py model html http
  • 如何使用增量值向 Pyspark 中的 DataFrame 添加列?

    我有一个名为 df 的 DataFrame 如下所示 Atr1 Atr2 Atr3 A A A B A A C A A 我想向其中添加一个具有增量值的新列并获取以下更新的 DataFrame Atr1 Atr2 Atr3
  • 有没有更快的方法将数字转换为名称?

    以下代码定义了映射到数字的名称序列 它的设计目的是获取一个号码并检索一个特定的名称 该类通过确保名称存在于其缓存中来进行操作 然后通过索引到其缓存中来返回名称 问题在这 如何在不存储缓存的情况下根据数字计算出名称 该名称可以被认为是一个以
  • numpy.polyfit 没有关键字“cov”

    我试图使用 polyfit 来找到一组数据的最佳拟合直线 但我还需要知道参数的不确定性 所以我也想要协方差矩阵 在线文档建议我写 polyfit x y 2 cov True 但这给出了错误 类型错误 polyfit 得到了意外的关键字参数
  • 了解 Tensorflow 中的 while 循环

    我正在使用用于 Tensorflow 的 Python API https www tensorflow org api docs python 我正在努力实施罗森布罗克函数 https www sfu ca ssurjano rosen
  • NumPy 中 exp(-x^2) 的快速傅立叶变换

    I have to calculate numerically the 2nd derivative of a Gaussian function I ve read every question on this topic here bu
  • 交响二阶颂歌

    我有一个简单的二阶 ODE 的齐次解 当我尝试使用 Sympy 求解初始值时 它返回相同的解 它应该替代 y 0 和 y 0 并产生一个没有常数的解 但事实并非如此 这是建立方程的代码 它是一个弹簧平衡方程 k 弹簧常数 m 质量 我在其他

随机推荐