用3D计算点到三角形距离的最快方法?

问题描述:

计算从点到三角形三角形的最小距离的一个显而易见的方法是将点投影到三角形的平面上,确定所得点的重心坐标,并使用它们来确定投影点是否位于三角形。如果不是,则将其重心坐标限制在[0,1]范围内,并给出位于三角形内的最近点。用3D计算点到三角形距离的最快方法?

有没有办法加快或简化它?

+2

夹紧重心坐标不给正交投影到最近的边缘,因此不正确的距离,如果一个点(没有顶点)是最接近的。 – 2015-01-12 10:47:53

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

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

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

这个paper比较二者的性能与二维方法获胜。

+0

我不同意以下结论:如果正确地构建查询,则需要4个平方根来计算与3D一起工作时与三角形的距离。它们都可以通过交叉点和点积来完成。 – bradgonesurfing 2017-11-27 07:10:11

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

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

此外,请考虑在SSE3或类似的东西(例如Mono's SIMD support)中编写算法。这是工作,但如果你对此有足够的努力,通常可以一次做几个三角形。

我会尽力找到关于该主题的一些论文,但您可能希望Google针对“Ray Mesh Intersection”。当人们努力优化这些东西时,这将带来所有来自80年代和90年代的伟大工作。

+0

“处理器可以在接近1个时钟周期内完成... DIV,因此内存通常是瓶颈。” - 真的?我听到的最后一个数字接近17个周期。 – 2016-03-09 11:26:38

我会带着我的测试用例结果。

enter image description here

测试用例代码和实现是在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; 

}