求解两个三维线段交点的算法

问题描述 投票:19回答:6

找到两个2D线段的交点很容易; the formula is straight forward。但我担心找到两条3D线段的交点不是。

算法是什么,在C#中最好找到两个3D线段的交点?

我发现了一个C++ implementation here。但我不相信解决方案,因为它优先考虑某个平面(看看perp在实现部分下实现的方式,它假设优先考虑z plane。任何通用算法都不能假设任何平面方向或偏好)。

有更好的解决方案吗?

c# math
6个回答
14
投票

大多数3D线不相交。一种可靠的方法是找到两条3D线之间的最短线。如果最短的线的长度为零(或距离小于您指定的任何公差),则您知道两条原始线相交。

enter image description here

找到the shortest line between two 3D lines, written by Paul Bourke的方法总结/解释如下:

在下文中,一条线将由位于其上的两个点定义,由点P1和P2定义的线“a”上的点具有等式

Pa = P1 + mua (P2 - P1)

类似地,由点P4和P4定义的第二行“b”上的点将被写为

Pb = P3 + mub (P4 - P3)

mua和mub的值从负无穷大到正无穷大。 P1 P2和P3 P4之间的线段具有0到1之间的对应μ。

有两种方法可以找到“a”和“b”行之间的最短线段。

方法一:

第一个是记下连接两条线的线段的长度,然后找到最小值。也就是说,最小化以下内容

|| Pb - Pa ||^2

用线的方程代替

|| P1 - P3 + mua (P2 - P1) - mub (P4 - P3) ||^2

然后可以在(x,y,z)分量中扩展上述内容。

有条件要满足,mua和mub的衍生物必须为零。 ...上述函数只有一个最小值,没有其他最小值或最大值。然后可以为mua和mub求解这两个方程,通过将mu的值代入线的原始方程中得到的实际交点。

方法二:

另一种给出完全相同方程的方法是认识到两条线之间的最短线段将垂直于两条线。这允许我们为点积写出两个方程式

(Pa - Pb) dot (P2 - P1) = 0
(Pa - Pb) dot (P4 - P3) = 0

根据线的等式扩展这些

( P1 - P3 + mua (P2 - P1) - mub (P4 - P3) ) dot (P2 - P1) = 0
( P1 - P3 + mua (P2 - P1) - mub (P4 - P3) ) dot (P4 - P3) = 0

根据坐标(x,y,z)扩展这些......结果如下

d1321 + mua d2121 - mub d4321 = 0
d1343 + mua d4321 - mub d4343 = 0

哪里

dmnop = (xm - xn)(xo - xp) + (ym - yn)(yo - yp) + (zm - zn)(zo - zp)

注意dmnop = dopmn

最后,解决mua给出了

mua = ( d1343 d4321 - d1321 d4343 ) / ( d2121 d4343 - d4321 d4321 )

并且代替代替mub

mub = ( d1343 + mua d4321 ) / d4343

这个方法被发现on Paul Bourke's website,这是一个很好的几何资源。该网站已重新组织,因此向下滚动以查找该主题。


4
投票
// This code in C++ works for me in 2d and 3d

// assume Coord has members x(), y() and z() and supports arithmetic operations
// that is Coord u + Coord v = u.x() + v.x(), u.y() + v.y(), u.z() + v.z()

inline Point
dot(const Coord& u, const Coord& v) 
{
return u.x() * v.x() + u.y() * v.y() + u.z() * v.z();   
}

inline Point
norm2( const Coord& v )
{
return v.x() * v.x() + v.y() * v.y() + v.z() * v.z();
}

inline Point
norm( const Coord& v ) 
{
return sqrt(norm2(v));
}

inline
Coord
cross( const Coord& b, const Coord& c) // cross product
{
return Coord(b.y() * c.z() - c.y() * b.z(), b.z() * c.x() - c.z() * b.x(), b.x() *  c.y() - c.x() * b.y());
}

bool 
intersection(const Line& a, const Line& b, Coord& ip)
// http://mathworld.wolfram.com/Line-LineIntersection.html
// in 3d; will also work in 2d if z components are 0
{
Coord da = a.second - a.first; 
Coord db = b.second - b.first;
    Coord dc = b.first - a.first;

if (dot(dc, cross(da,db)) != 0.0) // lines are not coplanar
    return false;

Point s = dot(cross(dc,db),cross(da,db)) / norm2(cross(da,db));
if (s >= 0.0 && s <= 1.0)
{
    ip = a.first + da * Coord(s,s,s);
    return true;
}

return false;
}

3
投票

我找到了一个解决方案:it's here

我的想法是利用矢量代数,使用dotcross来简化问题,直到这个阶段:

a (V1 X V2) = (P2 - P1) X V2

并计算a

请注意,此实现不需要任何平面或轴作为参考。


2
投票

我试过@Bill answer,它实际上每次都不起作用,我可以解释一下。基于link in his code.Let,例如这两个线段AB和CD。

A =(2,1,5),B =(1,2,5)和C =(2,1,3)和D =(2,1,2)

当你试图得到交叉点它可能会告诉你它是A点(不正确)或没有交叉点(正确)。根据您放置这些细分的顺序。

x = A +(B-A)s x = C +(D-C)t

比尔解决了s但从未解决过。并且由于您希望该交叉点位于两个线段上,因此s和t必须来自区间<0,1>。在我的例子中实际发生的是,只有从该间隔开始且t为-2。位于由C和D定义的线上,但不在线段CD上。

var s = Vector3.Dot(Vector3.Cross(dc, db), Vector3.Cross(da, db)) / Norm2(Vector3.Cross(da, db));

var t = Vector3.Dot(Vector3.Cross(dc, da), Vector3.Cross(da, db)) / Norm2(Vector3.Cross(da, db));

其中da是B-A,db是D-C,dc是C-A,我只保留了Bill提供的名字。

然后正如我所说,你必须检查s和t是否来自<0,1>并且你可以计算结果。基于上面的公式。

if ((s >= 0 && s <= 1) && (k >= 0 && k <= 1))
{
   Vector3 res = new Vector3(this.A.x + da.x * s, this.A.y + da.y * s, this.A.z + da.z * s);
}

比尔答案的另一个问题是当两条线共线并且有一个以上的交叉点时。会有零除。你想避免这种情况。


1
投票

但我担心找到两条3D线段的交点不是。

我觉得是这样的。您可以使用与2d(或任何其他维度)完全相同的方式找到交点。唯一的区别是,得到的线性方程组更可能没有解(意味着线不相交)。

您可以手动解决一般方程式,只需使用您的解决方案,或使用例如编程方法解决它。 Gaussian elemination


1
投票

你提到的原始资料仅适用于2d案例。 perp的实现很好。 x和y的使用只是变量而不是特定平面的偏好的指示。

同一个网站有这个3d案例:“http://geomalgorithms.com/a07-_distance.html

看起来Eberly撰写了回应:“https://www.geometrictools.com/Documentation/DistanceLine3Line3.pdf

把这些东西放在这里因为谷歌指向地理算法和这篇文章。

© www.soinside.com 2019 - 2024. All rights reserved.