数值精度计算交点

Question

数值精度计算交点

3

我希望能计算一条光线和一个线段之间的交点。为此，需要构建线性方程并寻找交点。现在我在处理一个案例时遇到了数值问题。以下是我的代码缩写：

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的值非常接近所致。我的问题是：在计算交点时，我能否得到更精确的结果？

- Geminus

当你的方程组接近于奇异矩阵时，数值问题往往会出现。 - Oliver Charlesworth

2个回答

0

据我所知，高斯或高斯-约旦方法配合完全主元选取可以实现最佳的数值稳定性。

您需要解决这个线性2x2系统，求出R和S的值。

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

全主元高斯消元法告诉你选择模数最大的左侧系数。有四种可能的选择，因此您将需要实现四个版本的算法。

例如，假设系统中左上角系数占优势。

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

使用精确算术，这确实等同于克莱姆法则，但从数值角度来看可能更好。

其他情况对称处理。

- user1196549

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- arghbleargh · Accepted Answer

以下是一种更加数值稳定的获取交点的方法（注意，实际上是两条线段的交点...似乎您原始的代码没有检查交点是否在线段内）：

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的比值等于ABX面积与ABY面积的比值。你可以使用叉积来计算面积，这就是上面代码所做的。请注意，此过程仅使用了一个除法，有助于最小化数值不稳定性。