From 88f5d790886b26efb9f370fb9e1ea2fa17079d19 Mon Sep 17 00:00:00 2001 From: Georgia Gkioxari Date: Fri, 18 Jun 2021 09:29:01 -0700 Subject: [PATCH] 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 | 37 +++++++++++++++++++++---- pytorch3d/csrc/utils/geometry_utils.h | 36 ++++++++++++++++++++---- tests/test_point_mesh_distance.py | 31 ++++++++++++++++++++- 3 files changed, 93 insertions(+), 11 deletions(-) diff --git a/pytorch3d/csrc/utils/geometry_utils.cuh b/pytorch3d/csrc/utils/geometry_utils.cuh index 16746396..2ba60d27 100644 --- a/pytorch3d/csrc/utils/geometry_utils.cuh +++ b/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; } diff --git a/pytorch3d/csrc/utils/geometry_utils.h b/pytorch3d/csrc/utils/geometry_utils.h index fd95a815..ef1f437d 100644 --- a/pytorch3d/csrc/utils/geometry_utils.h +++ b/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 +T AreaOfTriangle(const vec3& v0, const vec3& v1, const vec3& v2) { + vec3 p0 = v1 - v0; + vec3 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 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& v0, const vec3& v1, const vec3& v2) { - vec3 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 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; } diff --git a/tests/test_point_mesh_distance.py b/tests/test_point_mesh_distance.py index bef5d9af..6f35af0f 100644 --- a/tests/test_point_mesh_distance.py +++ b/tests/test_point_mesh_distance.py @@ -96,7 +96,7 @@ class TestPointMeshDistance(TestCaseMixin, unittest.TestCase): 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 @@ class TestPointMeshDistance(TestCaseMixin, unittest.TestCase): 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 @@ class TestPointMeshDistance(TestCaseMixin, unittest.TestCase): ) 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)