Fastest way to compute point to triangle distance

问题:

One obvious method for computing the minimum distance from a point to a 3D triangle is to project the point onto the plane of the triangle, determine the barycentric coordinates of the resulting point, and use them to determine whether the projected point lies within the triangle. If not, clamp its the barycentric coordinates to be in the range [0,1], and that gives you the closest point that lies inside the triangle.

Is there a way to speed this up or simplify it somehow?

回答1:

There are different approaches to finding the distance from a point P0 to a triangle P1,P2,P3.

1. The 3D method. Project the point onto the plane of the triangle and use barycentric coordinates or some other means of finding the closest point in the triangle. The distance is found in the usual way.

2. The 2D method. Apply a translation/rotation to the points so that P1 is on the origin, P2 is on the z-axis, P3 in the yz plane. Projection is of the point P0 is trivial (neglect the x coordinate). This results in a 2D problem. Using the edge equation it's possible to determine the closest vertex or edge of the triangle. Calculating distance is then easy-peasy.

This paper compares the performance of both with the 2D method winning.

回答2:

Assuming you're using one of the known fast algorithms, the only way to speed it up is when you are doing a lot of measurements on a lot of triangles. In that case, you can keep a lot of quantities precomputed in "edge" or "winding" structures. Instead of storing the 3 points, you store meshes comprised of edge structures. Projection then becomes very quick and barycentric tests can be coded so that they are branch-predictable.

The real key is to just keep everything in cache. Processors can do MUL and DIV in nearly 1 clock cycle so memory is usually the bottleneck.

Also, consider writing the algo in SSE3 or something similar (such as Mono's SIMD support). It's work, but you can usually do a couple triangles at a time if you think hard enough about it.

I'll try to find some papers on the topic, but you might want to Google for "Ray Mesh Intersection". That will bring up all the great work from the 80s and 90s when people worked hard on optimizing this stuff.

回答3:

I'll lead with my test case results.

The test case code and implementation is in 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  );

    }

The implementation code is fiddly as I have a number of framework classes. Hopefully you can treat this as pseudo code and pull out the algorithm. The raw vector types are from https://www.nuget.org/packages/System.DoubleNumerics/.

Note that some of the properties of Triangle could be cached to improve performance.

Note that to return the closest point does not require any square roots and does not require transforming the problem to 2D.

The algorithm first quickly tests if the test point is closest to an end point region. If that is inconclusive it then tests the edge external regions one by one. If those tests fail then the point is inside the triangle. Note that for randomly selected points far from the triangle it is most likely that the closest point will be a corner point of the triangle.

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 );

    }
}

And the edge structure

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;

}

And the plane structure

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;

}