diff --git a/pytorch3d/csrc/marching_cubes/marching_cubes_cpu.cpp b/pytorch3d/csrc/marching_cubes/marching_cubes_cpu.cpp index fa128e71..95e9e453 100644 --- a/pytorch3d/csrc/marching_cubes/marching_cubes_cpu.cpp +++ b/pytorch3d/csrc/marching_cubes/marching_cubes_cpu.cpp @@ -34,6 +34,18 @@ std::tuple MarchingCubesCpu( const int H = vol.size(1); const int W = vol.size(2); + // The edge-id hash packs two grid-point ids as v1_id * (W*H*D) + v2_id, whose + // maximum value is (W*H*D)^2 - 1. Guard against signed int64 overflow so an + // oversized volume fails loudly instead of silently producing colliding edge + // ids and a corrupted mesh. floor(sqrt(INT64_MAX)) == 3037000499 (~1448^3). + const int64_t num_grid_points = (int64_t)W * H * D; + TORCH_CHECK( + num_grid_points <= 3037000499LL, + "Volume too large for CPU marching cubes: W*H*D (", + num_grid_points, + ") exceeds 3037000499 (~1448^3), which would overflow the int64 edge-id " + "hash and corrupt the mesh."); + // Create tensor accessors auto vol_a = vol.accessor(); // edge_id_to_v maps from an edge id to a vertex position @@ -66,19 +78,23 @@ std::tuple MarchingCubesCpu( tri.push_back(edge); ps.push_back(interp_points[e]); - // Check if the triangle face is degenerate. A triangle face - // is degenerate if any of the two verices share the same 3D position - if ((j + 1) % 3 == 0 && ps[0] != ps[1] && ps[1] != ps[2] && - ps[2] != ps[0]) { - for (int k = 0; k < 3; k++) { - int64_t v = tri.at(k); - edge_id_to_v[v] = ps.at(k); - if (!uniq_edge_id.count(v)) { - uniq_edge_id[v] = verts.size(); - verts.push_back(edge_id_to_v[v]); + if (ps.size() == 3) { + // Check if the triangle face is degenerate. A triangle face + // is degenerate if any of the two vertices share the same 3D + // position + if (ps[0] != ps[1] && ps[1] != ps[2] && ps[2] != ps[0]) { + for (int k = 0; k < 3; k++) { + int64_t v = tri.at(k); + edge_id_to_v[v] = ps.at(k); + if (!uniq_edge_id.count(v)) { + uniq_edge_id[v] = verts.size(); + verts.push_back(edge_id_to_v[v]); + } + faces.push_back(uniq_edge_id[v]); } - faces.push_back(uniq_edge_id[v]); } + // Clear unconditionally - a rejected degenerate triangle must not + // corrupt the buffer for the next triangle in this cube. tri.clear(); ps.clear(); } // endif diff --git a/pytorch3d/csrc/marching_cubes/marching_cubes_utils.h b/pytorch3d/csrc/marching_cubes/marching_cubes_utils.h index 486e0339..00627fd0 100644 --- a/pytorch3d/csrc/marching_cubes/marching_cubes_utils.h +++ b/pytorch3d/csrc/marching_cubes/marching_cubes_utils.h @@ -126,8 +126,7 @@ struct Cube { // edge with an integer to address floating point precision issue. // // Args: - // v1_id: global id of vertex 1 - // v2_id: global id of vertex 2 + // edge: edge (ID) whose two endpoint vertices are hashed // W: width of the 3d grid // H: height of the 3d grid // D: depth of the 3d grid @@ -135,11 +134,19 @@ struct Cube { // Returns: // hashing for a pair of vertex ids // - int64_t HashVpair(const int edge, int W, int H, int D) { + int64_t HashVpair(const int edge, int64_t W, int64_t H, int64_t D) const { const int v1 = _EDGE_TO_VERTICES[edge][0]; const int v2 = _EDGE_TO_VERTICES[edge][1]; - const int v1_id = p[v1].x + p[v1].y * W + p[v1].z * W * H; - const int v2_id = p[v2].x + p[v2].y * W + p[v2].z * W * H; - return (int64_t)v1_id * (W + W * H + W * H * D) + (int64_t)v2_id; + // p[v].x/y/z hold integral corner coordinates, so casting to int64_t is + // exact. We cast before multiplication to avoid float32 precision loss at + // grids > 256^3. + const int64_t v1_id = + (int64_t)p[v1].x + (int64_t)p[v1].y * W + (int64_t)p[v1].z * W * H; + const int64_t v2_id = + (int64_t)p[v2].x + (int64_t)p[v2].y * W + (int64_t)p[v2].z * W * H; + // v_id max is (W*H*D - 1). A stride of W*H*D perfectly separates v1 and v2. + // Note: This hash fits in int64_t safely up to grids of ~1448^3. + const int64_t stride = W * H * D; + return v1_id * stride + v2_id; } }; diff --git a/pytorch3d/ops/marching_cubes.py b/pytorch3d/ops/marching_cubes.py index b6f90d8c..9a040ce9 100644 --- a/pytorch3d/ops/marching_cubes.py +++ b/pytorch3d/ops/marching_cubes.py @@ -159,7 +159,6 @@ def marching_cubes_naive( """ batched_verts, batched_faces = [], [] D, H, W = vol_batch.shape[1:] - # each edge is represented with its two endpoints (represented with global id) for i in range(len(vol_batch)): vol = vol_batch[i] @@ -184,7 +183,7 @@ def marching_cubes_naive( # triangle vertex IDs and positions tri = [] ps = [] - for i, edge in enumerate(edge_indices): + for edge in edge_indices: interp_points[edge] = cube.vert_interp(thresh, edge, vol) # Bind interpolated vertex with a global edge_id, which @@ -196,20 +195,19 @@ def marching_cubes_naive( ) tri.append(edge_id) ps.append(interp_points[edge]) - # when the isolevel are the same as the edge endpoints, the interploated - # vertices can share the same values, and lead to degenerate triangles. - if ( - (i + 1) % 3 == 0 - and ps[0] != ps[1] - and ps[1] != ps[2] - and ps[2] != ps[0] - ): - for j, edge_id in enumerate(tri): - edge_id_to_v[edge_id] = ps[j] - if edge_id not in uniq_edge_id: - uniq_edge_id[edge_id] = len(verts) - verts.append(edge_id_to_v[edge_id]) - faces.append([uniq_edge_id[tri[j]] for j in range(3)]) + if len(ps) == 3: + # when the isolevel are the same as the edge + # endpoints, the interpolated vertices can share + # the same values, and lead to degenerate + # triangles. + if ps[0] != ps[1] and ps[1] != ps[2] and ps[2] != ps[0]: + for j, edge_id in enumerate(tri): + edge_id_to_v[edge_id] = ps[j] + if edge_id not in uniq_edge_id: + uniq_edge_id[edge_id] = len(verts) + verts.append(edge_id_to_v[edge_id]) + faces.append([uniq_edge_id[tri[j]] for j in range(3)]) + tri = [] ps = [] diff --git a/tests/test_marching_cubes.py b/tests/test_marching_cubes.py index b48ed038..b874efef 100644 --- a/tests/test_marching_cubes.py +++ b/tests/test_marching_cubes.py @@ -780,6 +780,79 @@ class TestMarchingCubes(TestCaseMixin, unittest.TestCase): self.assertClose(faces[0], expected_faces) self.assertTrue(verts[0].ge(-1).all() and verts[0].le(1).all()) + def test_degenerate_triangle_keeps_later_faces(self): + # A cube whose first triangle is degenerate must not suppress the + # remaining triangles of the same cube. Here the isolevel coincides + # with the value of the outside corners, so every interpolated point + # snaps onto a corner and the cube's 4 candidate triangles collapse + # onto the 4 outside corners: triangles 1 and 4 become degenerate + # while triangles 2 and 3 stay valid. + volume_data = torch.ones(1, 2, 2, 2) # (B, W, H, D) + volume_data[0, 1, 0, 0] = 0 + volume_data[0, 1, 0, 1] = 0 + volume_data[0, 0, 1, 1] = 0 + volume_data[0, 1, 1, 1] = 0 + volume_data = volume_data.permute(0, 3, 2, 1) # (B, D, H, W) + + # The four inside corners are separated from the four outside corners + # (1, 1, 0), (0, 0, 1), (0, 1, 0) and (0, 0, 0), which the surface + # passes exactly through, giving a quad made of two triangles. + expected_verts = torch.tensor( + [ + [1.0, 1.0, 0.0], + [0.0, 0.0, 1.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 0.0], + ] + ) + expected_faces = torch.tensor([[0, 1, 2], [0, 3, 1]]) + + verts, faces = marching_cubes_naive(volume_data, 1, return_local_coords=False) + self.assertClose(verts[0], expected_verts) + self.assertClose(faces[0], expected_faces) + + verts, faces = marching_cubes(volume_data, 1, return_local_coords=False) + self.assertClose(verts[0], expected_verts) + self.assertClose(faces[0], expected_faces) + + def test_large_grid_edge_ids(self): + # The C++ implementation identifies a vertex by hashing the pair of + # grid-point ids of the edge it lies on. Grid-point ids run up to + # W * H * D - 1, so once the grid exceeds 2 ** 24 points the ids must + # not be computed in float32 or distinct edges collide and their + # vertices get merged. Use a long, thin volume to cross that bound + # with a volume small enough to allocate (~67MB). + W, H, D = 4_200_000, 2, 2 + # Isolated points on the (y=1, z=1) row, whose ids are x + W + W * H + # and therefore the largest in the grid. Spacing 3 keeps the blobs + # from sharing a cube while cycling through every id residue. + first_x = 2**24 - W * H - W + xs = [first_x + 3 * i for i in range(16)] + self.assertLess(xs[-1] + 1, W) + + vol = torch.ones(1, D, H, W) + for x in xs: + vol[0, 1, 1, x] = 0.0 + verts, faces = marching_cubes(vol, 0.5, return_local_coords=False) + + # Each isolated point cuts the four grid edges leading away from it, + # so every blob contributes 4 distinct vertices and 2 faces. + expected_verts = set() + for x in xs: + expected_verts.update( + [ + (x - 0.5, 1.0, 1.0), + (x + 0.5, 1.0, 1.0), + (float(x), 0.5, 1.0), + (float(x), 1.0, 0.5), + ] + ) + self.assertEqual( + {tuple(v) for v in verts[0].tolist()}, + expected_verts, + ) + self.assertEqual(faces[0].shape[0], 2 * len(xs)) + def test_sphere(self): # (B, W, H, D) volume = torch.Tensor(