射线/段交流测试的并行性和共线性使Python的浮点精度失败

Question

我正在尝试实现一个函数，以在Gareth Rees的出色指示下找到python中的射线/段相交：https://stackoverflow.com/a/14318254/7235455和https://stackoverflow.com/a/565282/7235455

这是我的功能：

from math import radians, sin, cos
import numpy as np
def find_intersection(point0, theta, point1, point2):
    # convert arguments to arrays:
    p = np.array(point0, dtype=np.float) # ray origin
    q = np.array(point1, dtype=np.float) # segment point 1
    q2 = np.array(point2, dtype=np.float) # segment point 2
    r = np.array((cos(theta),sin(theta))) # theta as vector (= ray as vector)
    s = q2 - q # vector from point1 to point2
    rxs = np.cross(r,s)
    qpxs = np.cross(q-p,s)
    qpxr = np.cross(q-p,r)
    t = qpxs/rxs
    u = qpxr/rxs
    if rxs == 0 and qpxr == 0:
        t0 = np.dot(q-p,r)/np.dot(r,r)
        t1 = np.dot(t0+s,r)/np.dot(r,r)
        return "collinear"
    elif rxs == 0 and qpxr != 0:
        return "parallel"
    elif rxs != 0 and 0 <= t and 0 <= u and u <= 1: # removed t <= 1 since ray is inifinte
        intersection = p+t*r
        return "intersection is {0}".format(intersection)
    else:
        return None

有交集时，该函数正常工作。但是它无法识别并行性或共线性，因为条件rxs == 0和qpxr == 0未达到（曾经？）。运行，例如：

p0 = (0.0,0.0)
theta = radians(45.0)
p1 = (1.0,1.0) 
p2 = (3.0,3.0)
c = find_intersection(p0,theta,p1,p2)

什么都没有。在IF-Block给出

之前，添加RXS和QPXR的打印语句

rxs =  2.22044604925e-16 qpxr =  -1.11022302463e-16

我的结论是，由于浮点问题，该功能无法捕获第一个if statement的条件。2.22044604925E-16和-1.11022302463E-16很小，但不幸的是。

我的结论是正确的还是我错过了什么？有什么想法避免了这个问题吗？非常感谢！

Answer 1

是的，您的结论是正确的，问题是"并行"谓词的数值稳定性。

您可以将结果与少量数字进行比较（例如，eps=1.0E-9）。它的大小可能取决于坐标范围（请注意，交叉产品给出了两倍的三角形区域，因此MaxVecLen**2正常化的eps看起来很合理）。

更复杂但更精确的选项 - 使用像这样的强大几何谓词。也许用于计算几何形状的Python/Numpy库包含一些用于此类操作的实现。

Answer 2

有一种简单且安全的方法来解决此问题。

写下射线的隐式方程（S(X, Y) = a X + b Y + c = 0）。当您在函数S中插入片段端点的坐标时，您将获得两个值，让S0和S1。如果它们的符号相反，则射线的支撑线与段之间有一个相交。

在这种情况下，沿段的交叉位置由参数的值给出，该值等于

- S0 / (S1 - S0).

此表达式享有始终可计算的属性（前提是符号更改）和范围[0, 1]。它允许安全地计算交叉点。

仅选择所需的半行（射线）上的那些相交，您只需计算射线的起源 S(Xo, Yo)的符号。

此过程无法检测到平行线或共线射线，但这没关系。无论如何，它会产生声音结果。