计算交点的数值不准确答案

【问题标题】：Numerical inaccuracy calculating intersection计算交点的数值不准确
【发布时间】：2014-11-22 22:09:10
【问题描述】：

我想计算射线和线段之间的交点。为此，我形成线性方程并寻找交点。现在我遇到一个例子的数值问题。我的代码的缩写：

public class Test {
    public static void main(String[] args) {
        double rayAX = 443.19661703858895d;
        double rayAY = 666.3485960845833d;

        double rayBX = 443.196744279195d;
        double rayBY = 103.21654864924565d;

        double segAX = 450.0d;
        double segAY = 114.42801992127828d;

        double segBX = 443.196744279195d;
        double segBY = 103.21654864924565d;


        double a1 = (rayBY - rayAY) / (rayBX - rayAX);
        double t1 = rayAY - rayAX * a1;

        double a2 = (segBY - segAY) / (segBX - segAX);
        double t2 = segAY - segAX * a2;

        double x = (t2 - t1) / (a1 - a2);
        double y = a1 * x + t1;

        System.out.println(x);
        System.out.println(y);
    }
}

显然返回应该是 (443.196744279195, 103.21654864924565) 因为这个点在射线和线段上是相同的。但实际回报是在我的情况下（443.19674427919506, 103.21654844284058）

在第二个数字中，小数点后六位已经有错误。我猜这个错误是因为 rayAX 和 rayBX 的值非常接近。我的问题是：在计算交点时，我能得到更精确的结果吗？

【问题讨论】：

当您的方程组接近 singular 时，总是可能出现数值问题。

标签： java geometry intersection numeric

【解决方案1】：

这是一种在数值上更稳定的获取交点的方法（请注意，它实际上是两条线的交点……看起来您的原始代码也没有检查交点是否在线段内）：

double rX = rayBX - rayAX;                                                                                                                                           
double rY = rayBY - rayAY;                                                                                                                                           

double sAX = segAX - rayAX;                                                                                                                                          
double sAY = segAY - rayAY;                                                                                                                                          
double areaA = sAX * rY - sAY * rX;                                                                                                                                  

double sBX = segBX - rayAX;                                                                                                                                          
double sBY = segBY - rayAY;                                                                                                                                          
double areaB = sBX * rY - sBY * rX;                                                                                                                                  

double t = areaA / (areaA - areaB);
// if t is not between 0 and 1, intersection is not in segment                                                                                                 
double x = (1 - t) * segAX + t * segBX;                                                                                                                              
double y = (1 - t) * segAY + t * segBY;

粗略解释：设A和B为射线的端点，设X和Y为线段的端点。让P 成为我们正在寻找的交点。那么，PX 与PY 的比率等于@987654330@ 的面积与ABY 的面积之比。您可以使用叉积计算面积，这就是上面的代码所做的。请注意此过程如何仅使用一个除法，这有助于最大限度地减少数值不稳定性。

【讨论】：

这肯定更准确！你有一个链接可以让我获得更多信息吗？
如果我想相交两个线段，我想如果我在计算过程中通过切换“seg”和“ray”来计算一个新的 areaA 和一个新的 areaB，我可以找到另一个 t？对于这个新的 t，我还可以检查它是否在 0 和 1 之间。
您是否有任何指向论文的链接或通常描述这种方法的地方？还是你自己想出了这个方法？
抱歉，我没有参考资料。这只是我根据一般几何知识自己想出的东西。不过，这看起来是个不错的资源：@987654321@.

【解决方案2】：

据我所知，最好的数值稳定性是通过全旋转的高斯或 Gauss-Jordan 方法实现的。

您需要求解R 和S 的线性2x2 系统。

(Brx - Arx).R - (Bsx - Asx).S = Asx - Arx
(Bxy - Ary).R - (Bsx - Asx).S = Asy - Ary

总旋转告诉您选择具有最大模块的 LHS 系数。有四种可能的选择，因此您必须实现该算法的四个版本。

例如，假设左上角的系数在系统中占主导地位

A.X + B.Y = C
D.X + E.Y = F

然后

  X + (B/A).Y = (C/A)
D.X + E    .Y = F

(E - D.(B/A)) Y = F - D.(C/A)

和

Y = (F - D.(C/A)) / (E - D.(B/A))
X = (C/A) - (B/A).Y

使用精确算术，这确实等同于 Cramer 规则，但从数值的角度来看可能更好。

其他情况对称处理。

【讨论】：