15

计算点到 3D 三角形的最小距离的一种明显方法是将点投影到三角形平面上,确定结果点的重心坐标,并使用它们来确定投影点是否位于三角形内。如果不是,则将其重心坐标钳制在 [0,1] 范围内,这将为您提供位于三角形内的最近点。

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

4

4 回答 4

12

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

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

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

本文将两者的性能与 2D 方法获胜进行了比较

于 2010-05-27T22:16:06.520 回答
7

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

在此处输入图像描述

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

}
于 2017-11-27T07:33:31.500 回答
6

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

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

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

我会尝试查找有关该主题的一些论文,但您可能想在 Google 上搜索“Ray Mesh Intersection”。当人们努力优化这些东西时,这将带来 80 年代和 90 年代的所有伟大工作。

于 2010-05-27T20:57:08.070 回答
3

我不认为重心坐标法本身可以加快很多,但是如果您要针对同一个三角形测试许多点,您可以预先计算一些东西。

这是我为计算这个投影而编写的一些健壮的代码,基于这个答案这篇文章(这又是基于Realtime Collision Detection书)。

请注意,您可以预先计算不直接或间接依赖的任何内容p(如果您多次重用同一个三角形,则可以为每个测试节省大约一半的计算工作)。

null如果点在三角形平面上的正交投影p不落在三角形内,则代码设计为返回。如果重心坐标超出范围,您可以扩展它以找到三角形边缘上的最近点,方法是计算投影点在每个三角形边缘向量上的投影,然后检查投影是否在任何一对三角形的顶点。(如果不是,角顶点将是最近的点。)但就我的目的而言,我只想要正交投影,这就是null如果正交投影不在三角形内,则此代码返回的原因。

代码是在 Java 中,使用 JOML 线性代数库。

/**
 * Find the closest orthogonal projection of a point p onto a triangle given by three vertices
 * a, b and c. Returns either the projection point, or null if the projection is not within
 * the triangle.
 */
public static Vector3d closestPoint(Vector3d p, Vector3d a, Vector3d b, Vector3d c) {
    // Find the normal to the plane: n = (b - a) x (c - a)
    Vector3d n = b.sub(a, new Vector3d()).cross(c.sub(a, new Vector3d()));

    // Normalize normal vector
    double nLen = n.length();
    if (nLen < 1.0e-30) {
        return null;  // Triangle is degenerate
    } else {
        n.mul(1.0f / nLen);
    }

    //    Project point p onto the plane spanned by a->b and a->c.
    //
    //    Given a plane
    //
    //        a : point on plane
    //        n : *unit* normal to plane
    //
    //    Then the *signed* distance from point p to the plane
    //    (in the direction of the normal) is
    //
    //        dist = p . n - a . n
    //
    double dist = p.dot(n) - a.dot(n);

    // Project p onto the plane by stepping the distance from p to the plane
    // in the direction opposite the normal: proj = p - dist * n
    Vector3d proj = p.add(n.mul(-dist, new Vector3d()), new Vector3d());

    // Find out if the projected point falls within the triangle -- see:
    // http://blackpawn.com/texts/pointinpoly/default.html

    // Compute edge vectors        
    double v0x = c.x - a.x;
    double v0y = c.y - a.y;
    double v0z = c.z - a.z;
    double v1x = b.x - a.x;
    double v1y = b.y - a.y;
    double v1z = b.z - a.z;
    double v2x = proj.x - a.x;
    double v2y = proj.y - a.y;
    double v2z = proj.z - a.z;

    // Compute dot products
    double dot00 = v0x * v0x + v0y * v0y + v0z * v0z;
    double dot01 = v0x * v1x + v0y * v1y + v0z * v1z;
    double dot02 = v0x * v2x + v0y * v2y + v0z * v2z;
    double dot11 = v1x * v1x + v1y * v1y + v1z * v1z;
    double dot12 = v1x * v2x + v1y * v2y + v1z * v2z;

    // Compute barycentric coordinates (u, v) of projection point
    double denom = (dot00 * dot11 - dot01 * dot01);
    if (Math.abs(denom) < 1.0e-30) {
        return null; // Triangle is degenerate
    }
    double invDenom = 1.0 / denom;
    double u = (dot11 * dot02 - dot01 * dot12) * invDenom;
    double v = (dot00 * dot12 - dot01 * dot02) * invDenom;

    // Check barycentric coordinates
    if ((u >= 0) && (v >= 0) && (u + v < 1)) {
        // Nearest orthogonal projection point is in triangle
        return proj;
    } else {
        // Nearest orthogonal projection point is outside triangle
        return null;
    }
}
于 2020-03-04T17:04:52.493 回答