在3D中计算点到三角距离的最快方法？

Question

计算从点到3D三角形的最小距离的一种显而易见的方法是将点投影到三角形的平面上,确定结果点的重心坐标,并使用它们来确定投影点是否位于三角形内.如果没有,将其重心坐标夹在[0,1]范围内,这样就可以得到三角形内最近的点.

有没有办法加快速度或以某种方式简化它？

Answer 1

找到从点P0到三角形P1,P2,P3的距离有不同的方法.

1. 3D方法.将点投影到三角形平面上,并使用重心坐标或其他方法找到三角形中的最近点.距离以通常的方式找到.

2. 2D方法.对点应用平移/旋转,使P1在原点上,P2在z轴上,P3在yz平面上.投影是P0是微不足道的(忽略x坐标).这导致2D问题.使用边缘方程,可以确定三角形的最近顶点或边缘.然后计算距离很容易.

这纸都与2D法获胜的性能进行比较.

Answer 2

假设您正在使用其中一种已知的快速算法,加速它的唯一方法是在很多三角形上进行大量测量.在这种情况下,您可以在"边缘"或"缠绕"结构中预先计算大量的数量.您可以存储由边缘结构组成的网格,而不是存储3个点.然后投影变得非常快,并且可以对重心测试进行编码,以使它们具有分支可预测性.

真正的关键是将所有内容保存在缓存中.处理器可以在近1个时钟周期内完成MUL和DIV,因此内存通常是瓶颈.

另外,考虑在SSE3中编写算法或类似的东西(例如Mono的SIMD支持).这是有用的,但如果你认真思考,你通常可以一次做几个三角形.

我会尝试找一些关于这个主题的论文,但你可能想要谷歌"Ray Mesh Intersection".当人们努力优化这些东西时,这将带来80年代和90年代的所有伟大工作.

Answer 3

我将带领我的测试用例结果。

测试用例代码和实现在C＃中

    public void ClosestPointToShouldWork()
    {
        var r = new Random(0);
        double next() => r.NextDouble() * 5 - 1;
        var t = new Triangle(new Vector3(0,0,0), new Vector3(3.5,2,0), new Vector3(3,0.0,0));

        DrawTriangle(t);

        var hash = new Vector3( 0, 0, 0 );
        for (int i = 0; i < 800; i++)
        {
            var pt = new Vector3( next(), next(), 0 );
            var pc = t.ClosestPointTo( pt );
            hash += pc;

            DrawLine(pc,pt);
        }

        // Test the hash
        // If it doesn't match then eyeball the visualization
        // and see what has gone wrong

        hash.ShouldBeApproximately( new Vector3(1496.28118561104,618.196568578824,0),1e-5  );

    }

由于我有许多框架类，因此实现代码很奇怪。希望您可以将其视为伪代码并提取算法。原始向量类型来自https://www.nuget.org/packages/System.DoubleNumerics/。

请注意，可以缓存Triangle的某些属性以提高性能。

请注意，返回最接近的点不需要任何平方根，也不需要将问题转换为2D。

该算法首先快速测试测试点是否最接近端点区域。如果不确定，则一个接一个地测试边缘外部区域。如果这些测试失败，则该点位于三角形内。请注意，对于远离三角形的随机选择的点，最有可能的点是最接近的点将是三角形的角点。

public class Triangle
{
    public Vector3 A => EdgeAb.A;
    public Vector3 B => EdgeBc.A;
    public Vector3 C => EdgeCa.A;

    public readonly Edge3 EdgeAb;
    public readonly Edge3 EdgeBc;
    public readonly Edge3 EdgeCa;

    public Triangle(Vector3 a, Vector3 b, Vector3 c)
    {
        EdgeAb = new Edge3( a, b );
        EdgeBc = new Edge3( b, c );
        EdgeCa = new Edge3( c, a );
        TriNorm = Vector3.Cross(a - b, a - c);
    }

    public Vector3[] Verticies => new[] {A, B, C};

    public readonly Vector3 TriNorm;

    private static readonly RangeDouble ZeroToOne = new RangeDouble(0,1);

    public Plane TriPlane => new Plane(A, TriNorm);

    // The below three could be pre-calculated to
    // trade off space vs time

    public Plane PlaneAb => new Plane(EdgeAb.A, Vector3.Cross(TriNorm, EdgeAb.Delta  ));
    public Plane PlaneBc => new Plane(EdgeBc.A, Vector3.Cross(TriNorm, EdgeBc.Delta  ));
    public Plane PlaneCa => new Plane(EdgeCa.A, Vector3.Cross(TriNorm, EdgeCa.Delta  ));

    public static readonly  RangeDouble Zero1 = new RangeDouble(0,1);

    public Vector3 ClosestPointTo(Vector3 p)
    {
        // Find the projection of the point onto the edge

        var uab = EdgeAb.Project( p );
        var uca = EdgeCa.Project( p );

        if (uca > 1 && uab < 0)
            return A;

        var ubc = EdgeBc.Project( p );

        if (uab > 1 && ubc < 0)
            return B;

        if (ubc > 1 && uca < 0)
            return C;

        if (ZeroToOne.Contains( uab ) && !PlaneAb.IsAbove( p ))
            return EdgeAb.PointAt( uab );

        if (ZeroToOne.Contains( ubc ) && !PlaneBc.IsAbove( p ))
            return EdgeBc.PointAt( ubc );

        if (ZeroToOne.Contains( uca ) && !PlaneCa.IsAbove( p ))
            return EdgeCa.PointAt( uca );

        // The closest point is in the triangle so 
        // project to the plane to find it
        return TriPlane.Project( p );

    }
}

和边缘结构

public struct Edge3
{

    public readonly Vector3 A;
    public readonly Vector3 B;
    public readonly Vector3 Delta;

    public Edge3(Vector3 a, Vector3 b)
    {
        A = a;
        B = b;
        Delta = b -a;
    }

    public Vector3 PointAt(double t) => A + t * Delta;
    public double LengthSquared => Delta.LengthSquared();

    public double Project(Vector3 p) => (p - A).Dot( Delta ) / LengthSquared;

}

和平面结构

public struct Plane
{
    public Vector3 Point;
    public Vector3 Direction;

    public Plane(Vector3 point, Vector3 direction )
    {
            Point = point;
            Direction = direction;
    }

    public bool IsAbove(Vector3 q) => Direction.Dot(q - Point) > 0;

}