Самый быстрый способ вычислить расстояние от точки до треугольника в 3D?

Question

Один очевидный метод для вычисления минимального расстояния от точки до трехмерного треугольника состоит в том, чтобы спроецировать точку на плоскость треугольника, определить барицентрические координаты полученной точки и использовать их, чтобы определить, находится ли проецируемая точка в пределах треугольник Если нет, зажмите его барицентрические координаты в диапазоне [0,1], и это даст вам ближайшую точку, которая находится внутри треугольника.

Есть ли способ ускорить или упростить это?

Answer 1

Существуют разные подходы к нахождению расстояния от точки P0 до треугольника P1, P2, P3.

1. 3D метод. Спроецируйте точку на плоскость треугольника и используйте барицентрические координаты или другие способы нахождения ближайшей точки в треугольнике. Расстояние определяется обычным способом.

2. 2D метод. Примените перемещение / вращение к точкам так, чтобы P1 находился в начале координат, P2 - на оси z, P3 - в плоскости yz. Проекция точки P0 тривиальна (не учитывайте координату x). Это приводит к двумерной проблеме. Используя уравнение ребра, можно определить ближайшую вершину или ребро треугольника. Тогда вычисление расстояния легко-peasy.

Эта статья сравнивает эффективность обоих с победой в 2D-методе.

Answer 2

Предполагая, что вы используете один из известных быстрых алгоритмов, единственный способ ускорить его - это когда вы делаете много измерений на множестве треугольников.В этом случае вы можете сохранить много величин, предварительно вычисленных в «ребристых» или «извилистых» структурах .Вместо того, чтобы хранить 3 точки, вы сохраняете сетки, состоящие из краевых структур.Затем проекция становится очень быстрой, и барицентрические тесты можно кодировать так, чтобы они были предсказуемыми ветвями .

Реальный ключ - просто хранить все в кэше.Процессоры могут выполнять MUL и DIV почти за 1 такт, поэтому память обычно является узким местом.

Кроме того, рассмотрите возможность записи алгоритма в SSE3 или что-то подобное (например, Поддержка SIMD в Mono).).Это работа, но вы обычно можете сделать пару треугольников за раз, если вы достаточно серьезно об этом думаете.

Я постараюсь найти несколько статей по этой теме, но вы можете обратиться в Google за "Ray Mesh"Точка пересечения».Это поднимет всю замечательную работу 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/.

Обратите внимание, что некоторые свойства треугольника могут быть кэшированы для повышения производительности.

Обратите внимание, что для возврата ближайшей точки не требуются квадратные корни и не требуется преобразование задачи в 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;

}