fix small face issue for ptmeshdist

Summary:
Fix small face issue for point_mesh distance computation.

The issue lies in the computation of `IsInsideTriangle` which is unstable and non-symmetrical when faces with small areas are given as input. This diff fixes the issue by returning `False` for `IsInsideTriangle` when small faces are given as input.

Reviewed By: bottler

Differential Revision: D29163052

fbshipit-source-id: be297002f26b5e6eded9394fde00553a37406bee

pytorch3d/csrc/utils/geometry_utils.cuh

@@ -461,6 +461,27 @@ PointTriangleDistanceBackward(
 //                          vec3 utils                           //
 // ************************************************************* //
 
+// Computes the area of a triangle (v0, v1, v2).
+//
+// Args:
+//     v0, v1, v2: vec3 coordinates of the triangle vertices
+//
+// Returns
+//     area: float: The area of the triangle
+//
+__device__ inline float
+AreaOfTriangle(const float3& v0, const float3& v1, const float3& v2) {
+  float3 p0 = v1 - v0;
+  float3 p1 = v2 - v0;
+
+  // compute the hypotenus of the scross product (p0 x p1)
+  float dd = hypot(
+      p0.y * p1.z - p0.z * p1.y,
+      hypot(p0.z * p1.x - p0.x * p1.z, p0.x * p1.y - p0.y * p1.x));
+
+  return dd / 2.0;
+}
+
 // Computes the barycentric coordinates of a point p relative
 // to a triangle (v0, v1, v2), i.e. p = w0 * v0 + w1 * v1 + w2 * v2
 // s.t. w0 + w1 + w2 = 1.0
@@ -503,6 +524,7 @@ __device__ inline float3 BarycentricCoords3Forward(
 // Checks whether the point p is inside the triangle (v0, v1, v2).
 // A point is inside the triangle, if all barycentric coordinates
 // wrt the triangle are >= 0 & <= 1.
+// If the triangle is degenerate, aka line or point, then return False.
 //
 // NOTE that this function assumes that p lives on the space spanned
 // by (v0, v1, v2).
@@ -521,11 +543,16 @@ __device__ inline bool IsInsideTriangle(
     const float3& v0,
     const float3& v1,
     const float3& v2) {
-  float3 bary = BarycentricCoords3Forward(p, v0, v1, v2);
-  bool x_in = 0.0f <= bary.x && bary.x <= 1.0f;
-  bool y_in = 0.0f <= bary.y && bary.y <= 1.0f;
-  bool z_in = 0.0f <= bary.z && bary.z <= 1.0f;
-  bool inside = x_in && y_in && z_in;
+  bool inside;
+  if (AreaOfTriangle(v0, v1, v2) < 1e-5) {
+    inside = 0;
+  } else {
+    float3 bary = BarycentricCoords3Forward(p, v0, v1, v2);
+    bool x_in = 0.0f <= bary.x && bary.x <= 1.0f;
+    bool y_in = 0.0f <= bary.y && bary.y <= 1.0f;
+    bool z_in = 0.0f <= bary.z && bary.z <= 1.0f;
+    inside = x_in && y_in && z_in;
+  }
   return inside;
 }
 

pytorch3d/csrc/utils/geometry_utils.h

@@ -560,6 +560,26 @@ PointTriangleDistanceBackward(
   return std::make_tuple(grad_p, grad_v0, grad_v1, grad_v2);
 }
 
+// Computes the area of a triangle (v0, v1, v2).
+// Args:
+//     v0, v1, v2: vec3 coordinates of the triangle vertices
+//
+// Returns:
+//     area: float: the area of the triangle
+//
+template <typename T>
+T AreaOfTriangle(const vec3<T>& v0, const vec3<T>& v1, const vec3<T>& v2) {
+  vec3<T> p0 = v1 - v0;
+  vec3<T> p1 = v2 - v0;
+
+  // compute the hypotenus of the scross product (p0 x p1)
+  float dd = std::hypot(
+      p0.y * p1.z - p0.z * p1.y,
+      std::hypot(p0.z * p1.x - p0.x * p1.z, p0.x * p1.y - p0.y * p1.x));
+
+  return dd / 2.0;
+}
+
 // Computes the squared distance of a point p relative to a triangle (v0, v1,
 // v2). If the point's projection p0 on the plane spanned by (v0, v1, v2) is
 // inside the triangle with vertices (v0, v1, v2), then the returned value is
@@ -604,6 +624,7 @@ vec3<T> BarycentricCoords3Forward(
 // Checks whether the point p is inside the triangle (v0, v1, v2).
 // A point is inside the triangle, if all barycentric coordinates
 // wrt the triangle are >= 0 & <= 1.
+// If the triangle is degenerate, aka line or point, then return False.
 //
 // NOTE that this function assumes that p lives on the space spanned
 // by (v0, v1, v2).
@@ -623,11 +644,16 @@ static bool IsInsideTriangle(
     const vec3<T>& v0,
     const vec3<T>& v1,
     const vec3<T>& v2) {
-  vec3<T> bary = BarycentricCoords3Forward(p, v0, v1, v2);
-  bool x_in = 0.0f <= bary.x && bary.x <= 1.0f;
-  bool y_in = 0.0f <= bary.y && bary.y <= 1.0f;
-  bool z_in = 0.0f <= bary.z && bary.z <= 1.0f;
-  bool inside = x_in && y_in && z_in;
+  bool inside;
+  if (AreaOfTriangle(v0, v1, v2) < 1e-5) {
+    inside = 0;
+  } else {
+    vec3<T> bary = BarycentricCoords3Forward(p, v0, v1, v2);
+    bool x_in = 0.0f <= bary.x && bary.x <= 1.0f;
+    bool y_in = 0.0f <= bary.y && bary.y <= 1.0f;
+    bool z_in = 0.0f <= bary.z && bary.z <= 1.0f;
+    inside = x_in && y_in && z_in;
+  }
   return inside;
 }
 

tests/test_point_mesh_distance.py

@@ -96,7 +96,7 @@ def _point_to_bary(point: torch.Tensor, tri: torch.Tensor) -> torch.Tensor:
         d20 = v2.dot(v0)
         d21 = v2.dot(v1)
 
-        denom = d00 * d11 - d01 * d01
+        denom = d00 * d11 - d01 * d01 + TestPointMeshDistance.eps()
         s2 = (d11 * d20 - d01 * d21) / denom
         s3 = (d00 * d21 - d01 * d20) / denom
         s1 = 1.0 - s2 - s3
@@ -117,6 +117,13 @@ def _is_inside_triangle(point: torch.Tensor, tri: torch.Tensor) -> torch.Tensor:
         Returns:
             inside: BoolTensor of shape (1)
         """
+        v0 = tri[1] - tri[0]
+        v1 = tri[2] - tri[0]
+        area = torch.cross(v0, v1).norm() / 2.0
+
+        # check if triangle is a line or a point. In that case, return False
+        if area < 1e-5:
+            return False
         bary = TestPointMeshDistance._point_to_bary(point, tri)
         inside = ((bary >= 0.0) * (bary <= 1.0)).all()
         return inside
@@ -836,6 +843,28 @@ def test_point_mesh_face_distance(self):
             )
             self.assertClose(pcls.points_list()[i].grad, pcls_op.points_list()[i].grad)
 
+    def test_small_faces_case(self):
+        for device in [torch.device("cpu"), torch.device("cuda:0")]:
+            mesh_vertices = torch.tensor(
+                [
+                    [-0.0021, -0.3769, 0.7146],
+                    [-0.0161, -0.3771, 0.7146],
+                    [-0.0021, -0.3771, 0.7147],
+                ],
+                dtype=torch.float32,
+                device=device,
+            )
+            mesh1_faces = torch.tensor([[0, 2, 1]], device=device)
+            mesh2_faces = torch.tensor([[2, 0, 1]], device=device)
+            pcd_points = torch.tensor([[-0.3623, -0.5340, 0.7727]], device=device)
+            mesh1 = Meshes(verts=[mesh_vertices], faces=[mesh1_faces])
+            mesh2 = Meshes(verts=[mesh_vertices], faces=[mesh2_faces])
+            pcd = Pointclouds(points=[pcd_points])
+
+            loss1 = point_mesh_face_distance(mesh1, pcd)
+            loss2 = point_mesh_face_distance(mesh2, pcd)
+            self.assertClose(loss1, loss2)
+
     @staticmethod
     def point_mesh_edge(N: int, V: int, F: int, P: int, device: str):
         device = torch.device(device)