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

Question

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

Есть ли способ ускорить это или упростить его каким-то образом?

Answer 1

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

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

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

Этот paper сравнивает производительность как с методом 2D-метода.

Answer 2

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

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

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

Я попытаюсь найти некоторые статьи по этой теме, но вы можете захотеть в 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/.

Обратите внимание, что некоторые свойства 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; 

}