From a5093b158b3dcce80ac4f69bc34afb5104fbbeab Mon Sep 17 00:00:00 2001 From: erald ceni Date: Mon, 17 Aug 2026 06:49:38 -0700 Subject: [PATCH] Fix CPU marching cubes precision and topology at high resolutions (#1934) (#2043) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Summary: Fixes https://github.com/facebookresearch/pytorch3d/issues/1934. Two independent bugs in the CPU backend. ### 1. Degenerate-triangle filter discards valid faces `marching_cubes_cpu.cpp`, `marching_cubes.py` `tri.clear()` and `ps.clear()` sit inside the degeneracy check, so the buffers only reset when a triangle is *accepted*. Once a cube's first triangle is degenerate, `ps[0..2]` stay frozen on it, and every subsequent triangle in that cube fails the same stale check and is dropped. Fixed by gating on `ps.size() == 3` and clearing unconditionally. The old code could only ever drop faces, never emit incorrect ones, so this is strictly additive. ### 2. Edge hash computed in float32 `marching_cubes_utils.h` `p[v].x/y/z` hold integral coordinates but are stored as `float`, so `x + y*W + z*W*H` evaluates entirely in float32 before truncating to `int`. float32 is exact only to 2²⁴ − 1 = 16,777,215 — and 256³ maxes out at exactly that value. At 512³ the maximum id is 134,217,727, where float32 spacing is 8, so distinct vertices collide on one id and `uniq_edge_id` merges them. Fixed by widening `W/H/D` to `int64_t` and casting each coordinate before multiplying. Also tightens the stride from `(W + W*H + W*H*D)` to `W*H*D`. Raises the CPU ceiling from 256³ to 1448³. Scope is the CPU path only. `marching_cubes_naive` was never affected (Python ints are arbitrary-precision), and neither was CUDA: `hashVpair` there computes ids in `uint` rather than `float`, so it has no 2²⁴ cliff, and `MarchingCubes` already rejects volumes above 1024³ before the CUDA kernel runs. The new `TORCH_CHECK` bound and the "~1448³" note in the new comments describe `MarchingCubesCpu` only. ### Verification Ellipsoid SDF (0.1, 1, 1), `isolevel=0.0`, identical input tensors on both devices. | Resolution | CUDA V | CUDA F | CPU V | CPU F | Degenerate dropped | | -- | -- | -- | -- | -- | -- | | 32³ | 1,664 | 3,324 | 1,664 | 3,324 | 0 | | 64³ | 7,312 | 14,620 | 7,312 | 14,620 | 0 | | 128³ | 30,168 | 60,332 | 30,168 | 60,332 | 0 | | 256³ | 122,448 | 244,892 | 122,448 | 243,996 | 896 | | 512³ | 491,944 | 983,884 | 491,944 | 976,140 | 7,744 | Machine: Arch Linux, RTX 4080, Ryzen 7 7800X3D Vertex counts now match CUDA exactly at every resolution; 512³ previously produced 176,121. The remaining face gap is entirely degenerate geometry — the CUDA mesh at 512³ contains exactly 7,744 zero-area triangles. ### Tests `test_degenerate_triangle_keeps_later_faces` — a 2×2×2 volume at `isolevel=1` chosen so the cube's four candidate triangles collapse onto its four outside corners: triangles 1 and 4 become degenerate, 2 and 3 stay valid. Pre-fix, the first degeneracy suppresses the rest and the mesh comes back empty; post-fix it is the expected quad. Asserts both `marching_cubes_naive` and the C++ extension. `test_large_grid_edge_ids` — a 2×2×4,200,000 volume (~67MB, ~0.1s) holding 16 isolated interior points on the highest-id grid row, positioned so grid-point ids straddle 2²⁴. Each point cuts exactly four grid edges, so the 64-vertex expectation is derived geometrically rather than copied from output. Pre-fix, a hash collision merges two edges and one vertex is lost. All 26 pre-existing tests in `test_marching_cubes.py` pass **unchanged**. That includes `test_cube_no_duplicate_verts` (`isolevel=1`) and `test_sphere` (`isolevel=64`), which both exercise the degenerate path but whose output is identical before and after the fix — so no existing expectation was edited and `sphere_level64.pickle` does not need regenerating. The imported diff contained no test file. The two tests above were written during import and differ from the tests described in the upstream PR description. *Analysis and write-up done collaboratively with AI, figures from testing are done on my own machine and have been checked.* Pull Request resolved: https://github.com/facebookresearch/pytorch3d/pull/2043 Test Plan: ``` buck2 test fbcode//vision/fair/pytorch3d:tests -- --regex 'test_marching_cubes' ``` `Pass 28. Fail 0.` — 26 pre-existing tests plus the 2 new regression tests. Reverting all three source hunks to their pre-fix state and re-running the same command: both new tests fail (`test_degenerate_triangle_keeps_later_faces` returns an empty mesh instead of 4 verts / 2 faces; `test_large_grid_edge_ids` returns 31 verts instead of 32) and all 26 pre-existing tests still pass. The new tests are therefore pinned to exactly this change, and the change breaks nothing that was already covered. Reviewed By: MichaelRamamonjisoa Differential Revision: D115424433 Pulled By: bottler fbshipit-source-id: 547a260010b94253f52f3a3c223d4c4fa78a7ee6 --- .../marching_cubes/marching_cubes_cpu.cpp | 38 +++++++--- .../marching_cubes/marching_cubes_utils.h | 19 +++-- pytorch3d/ops/marching_cubes.py | 30 ++++---- tests/test_marching_cubes.py | 73 +++++++++++++++++++ 4 files changed, 127 insertions(+), 33 deletions(-) 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(