Repository: mikedh/trimesh
Pull request #310 (2019-02-03 15:27:31)

closest_point() fix for (obtuse) triangles
Hai again,

I have noticed an issue with the signed_distance() and on_surface() methods and were able to trace it back to an incorrect computation of the closest_point() method of the triangle class.
The selection for case_vertex and case_edge based on the number of positive barycentric coordinates is only correct for certain points but not for all outside of a triangle.
For cases which belong in case_vertex but end up in case_edge, the clipping at the end (to the edge's extend) will fix this and return the correct vertex position.
However, for obtuse triangles there are scenarios that belong in case_edge but end up in case_vertex and will return a wrong point.
I extended the triangle class with a minimally invasive fix on the barycentric-coordinate's 'positive' mask to also handle these cases correctly.
Furthermore, I have added a simple test for one point with ground-truth result and a circle of points with a naive point-to-edge distance comparison to check the results.

The fix feels kind of hacky, but as far as I have checked, it seems to produce correct results now.
Do you think its worth to refactor the method? Let me know, I think I could come up with something.

Cheers,
Dennis

Also - here is a nice wysiwyg-tool which I used to play around with barycentric coordinates and to figure out when and why the current method failed:
https://www.geogebra.org/m/ZuvmPjmy


Patch:
diff --git a/trimesh/triangles.py b/trimesh/triangles.py
index 30683308..d16bc58c 100644
--- a/trimesh/triangles.py
+++ b/trimesh/triangles.py
@@ -73,6 +73,25 @@ def normals(triangles=None, crosses=None):
     return unit, valid
 
 
+def angles(triangles):
+    """
+    Calculates the angles of input triangles
+
+    Parameters
+    ------------
+    triangles:   (n, 3, 3) float, vertex positions
+
+    Returns
+    ------------
+    angles: (n, 3) float, angles at vertex positions (radian)
+    """
+    ABCs = util.unitize((triangles - np.roll(triangles, -1, axis=1)).reshape(-1,3))
+    inner1d = lambda us, vs: np.einsum('...i,...i', us, vs)
+    return np.arccos([inner1d(ABCs[::3], -ABCs[2::3]),      # AB o AC
+                      inner1d(-ABCs[::3], ABCs[1::3]),      # BA o BC
+                      inner1d(ABCs[2::3], -ABCs[1::3])]).T  # CA o CB
+
+
 def all_coplanar(triangles):
     """
     Check to see if a list of triangles are all coplanar
@@ -485,6 +504,24 @@ def closest_point(triangles, points):
 
     # signs of barycentric coordinates
     positive = barycentric > -tol.zero
+
+    # compute angles at the triangles corners and check for obtuse tringles
+    corner_angles = angles(triangles)
+    corners_obtuse = corner_angles >= np.pi/2 - tol.merge
+
+    # the case selection below is not really correct for points outside of a triangle
+    # but for most triangles the clipping at the end will fix this
+    # however obtuse triangles require special care
+    # if the obtuse corner is the only one with a positive barycentric coordinate
+    positive_and_obtuse = np.sum(positive * corners_obtuse, axis=1) == 1
+
+    # these cases do not belong in case_vertex so we hack the 'positive' mask
+    # the case_edge treatment requires another vertex (one neighbor of the pos. obtuse corner)
+    # we identify this neighbor as the one with the larger barycentric coordinate (both < 0)
+    # use argsort: the pos. corner is index 2, the larger of the two neg. neighbors is index 1
+    bary_coord_idxs = np.argsort(barycentric[positive_and_obtuse], axis=1)
+    positive[positive_and_obtuse, bary_coord_idxs[:,1]] = True
+    
     positive_sum = positive.sum(axis=1)
     # cases for signs of barycentric coordinates:
     # 2 negative, 1 positive: closest point is positive vertex


Test patch:
diff --git a/tests/test_triangles.py b/tests/test_triangles.py
index 0255560f..c44eb868 100644
--- a/tests/test_triangles.py
+++ b/tests/test_triangles.py
@@ -34,6 +34,51 @@ class TrianglesTest(g.unittest.TestCase):
         self.assertTrue((comparison < 1e-8).all())
         g.log.info('finished closest check on %d triangles', len(closest))
 
+    def test_closest_obtuse(self):
+        # simple triangle in the xy-plane with an obtuse corner at vertex A
+        ABC = g.np.float32([[0,0,0], [2,0,0], [-2,1,0]])
+        D = g.np.float32([1, -1, 0])
+
+        # ground truth: closest point from D is the center of the AB edge: (1,0,0)
+        gt_closest = g.np.float32([1,0,0])
+        tm_closest = g.trimesh.triangles.closest_point([ABC], [D])[0]
+
+        self.assertTrue(g.np.linalg.norm(gt_closest - tm_closest) < g.tol.merge)
+
+
+        # create a circle of points around the triangle
+        # with a radius so that all points are outside of the triangle
+        radius = 3
+        nPtsOnCircle = 100
+        alphas = g.np.linspace(g.np.pi/nPtsOnCircle, g.np.pi*2-g.np.pi/nPtsOnCircle, nPtsOnCircle)
+        ptsOnCircle = g.np.transpose([g.np.cos(alphas), g.np.sin(alphas), g.np.zeros(nPtsOnCircle)]) * radius
+
+        norm = lambda v: g.np.sqrt(g.np.einsum('...i,...i', v, v))
+        distToLine = lambda o, v, p: norm((o-p)-g.np.dot(o-p,v)*v)
+
+        def distPointToEdge(U, V, P): # edge [U, V], point P
+            UtoV = V - U
+            UtoP = P - U
+            VtoP = P - V
+    
+            if g.np.dot(UtoV, UtoP) <= 0:
+                # P is 'behind' U
+                return norm(UtoP)
+            elif g.np.dot(-UtoV, VtoP) <= 0:
+                # P is 'behind' V
+                return norm(VtoP)
+            else:
+                # P is 'between' U and V
+                return distToLine(U, UtoV/norm(UtoV), P)
+
+        # get closest points from trimesh and compute distances to the circle points
+        tm_dists = norm(ptsOnCircle - g.trimesh.triangles.closest_point([ABC]*nPtsOnCircle, ptsOnCircle))
+        # compute naive point-to-edge distances for all points and take the min of the three edges
+        gt_dists = g.np.float32([[distPointToEdge(ABC[i], ABC[(i+1)%3], pt) for i in range(3)] for pt in ptsOnCircle]).min(axis=1)
+        diff_dists = tm_dists - gt_dists
+
+        self.assertTrue(g.np.dot(diff_dists, diff_dists) < g.tol.merge)
+
     def test_degenerate(self):
         tri = [[[0, 0, 0],
                 [1, 0, 0],