Fix CPU marching cubes precision and topology at high resolutions (#1934) (#2043)

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
This commit is contained in:
erald ceni
2026-08-17 06:49:38 -07:00
committed by meta-codesync[bot]
parent 3143b3baf8
commit a5093b158b
4 changed files with 127 additions and 33 deletions

View File

@@ -34,6 +34,18 @@ std::tuple<at::Tensor, at::Tensor, at::Tensor> 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<float, 3>();
// edge_id_to_v maps from an edge id to a vertex position
@@ -66,19 +78,23 @@ std::tuple<at::Tensor, at::Tensor, at::Tensor> 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

View File

@@ -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;
}
};

View File

@@ -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 = []

View File

@@ -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(