From 6883d2da47b7501278cdd48ccd2d834ad0140e3c Mon Sep 17 00:00:00 2001 From: Jonathan Swartz Date: Wed, 24 Jun 2026 05:05:48 +0000 Subject: [PATCH 01/10] Add `dual_contour` (DC+QEF) mesh extraction op Add dual_contour, a VoxelBlockManager-based dual-contouring mesher that turns a narrow-band SDF on an OnIndex grid into a triangle mesh, alongside the existing marching_cubes. It returns (vertices, faces, normals) as three JaggedTensors jagged over the grid batch, where normals are the normalized SDF gradient at each vertex. Algorithm: one vertex per surface cell placed by minimizing a quadratic error function over the 12 edge zero-crossings and their interpolated normals (centroid-recentred + Tikhonov-regularized); dual connectivity (one quad per sign-changing minimal grid edge, triangulated and grad and optional cluster-collapse decimation via `reduce` (uniform F^3 blocks) or `adaptivity` (flat-block collapse). CUDA only; consumes a >= ~3 half-width band (e.g. from retopologize_sdf) for a watertight result. Follows the dual-contouring approach of OpenVDB's VolumeToMesh, but pla QEF rather than an averaged mass-point and decimates by clustering rather than a seam-stitched octree merge. Sources are cited in the code: Ju Garland & Heckbert 1997, Kobbelt et al. 2001, Lorensen & Cline 1987. API: Grid.dual_contour / GridBatch.dual_contour, functional dual_contour_{single,batch}, the pybind binding, and the .pyi s Factor the shared NanoVDB VBM scaffolding (grid/buffer aliases, the build-once VBMHelper, and the per-block decode + 6-face stencil helpers) into detail/utils/cuda/VoxelBlockManagerHelper.h, and refactor onto it (replacing its VBM_FACES_BEGIN macro with the shared device helpers). Also tidy the SDF op docstrings: describe `band` as the narrow- and use the repo's shape/dtype docstring style. Tests (tests/unit/test_dc.py): mesh validity + closed genus-0 (Euler characteristic) topology; analytic ground truth (a planar SDF r exactly, sphere radius/normals); composition with reinitialize/retopologize_sdf; batch-vs-single parity; decimation bounds; and the empty-surfac Signed-off-by: Jonathan Swartz --- fvdb/_fvdb_cpp.pyi | 7 + fvdb/functional/__init__.py | 20 +- fvdb/functional/_meshing.py | 63 + fvdb/functional/_sdf.py | 28 +- fvdb/grid.py | 39 +- fvdb/grid_batch.py | 41 +- src/CMakeLists.txt | 1 + src/fvdb/detail/ops/DualContour.cu | 1285 +++++++++++++++++ src/fvdb/detail/ops/DualContour.h | 49 + src/fvdb/detail/ops/ReinitializeSdf.cu | 80 +- src/fvdb/detail/ops/ReinitializeSdf.h | 6 +- .../utils/cuda/VoxelBlockManagerHelper.h | 131 ++ src/python/GridBatchOps.cpp | 9 + tests/unit/test_dc.py | 195 +++ 14 files changed, 1856 insertions(+), 98 deletions(-) create mode 100644 src/fvdb/detail/ops/DualContour.cu create mode 100644 src/fvdb/detail/ops/DualContour.h create mode 100644 src/fvdb/detail/utils/cuda/VoxelBlockManagerHelper.h create mode 100644 tests/unit/test_dc.py diff --git a/fvdb/_fvdb_cpp.pyi b/fvdb/_fvdb_cpp.pyi index f522d955d..063201379 100644 --- a/fvdb/_fvdb_cpp.pyi +++ b/fvdb/_fvdb_cpp.pyi @@ -847,6 +847,13 @@ def marching_cubes( field: JaggedTensor, level: float, ) -> list[JaggedTensor]: ... +def dual_contour( + grid: GridBatchData, + field: JaggedTensor, + iso: float, + reduce: int, + adaptivity: float, +) -> tuple[JaggedTensor, JaggedTensor, JaggedTensor]: ... def integrate_tsdf( grid: GridBatchData, truncation_margin: float, diff --git a/fvdb/functional/__init__.py b/fvdb/functional/__init__.py index 299a77527..c451b8592 100644 --- a/fvdb/functional/__init__.py +++ b/fvdb/functional/__init__.py @@ -85,6 +85,8 @@ # Meshing / TSDF from ._meshing import ( + dual_contour_batch, + dual_contour_single, integrate_tsdf_batch, integrate_tsdf_single, integrate_tsdf_with_features_batch, @@ -93,14 +95,6 @@ marching_cubes_single, ) -# Signed distance fields -from ._sdf import ( - reinitialize_sdf_batch, - reinitialize_sdf_single, - retopologize_sdf_batch, - retopologize_sdf_single, -) - # Pooling / refinement from ._pooling import ( avg_pool_batch, @@ -145,6 +139,14 @@ voxels_along_rays_single, ) +# Signed distance fields +from ._sdf import ( + reinitialize_sdf_batch, + reinitialize_sdf_single, + retopologize_sdf_batch, + retopologize_sdf_single, +) + # Grid topology from ._topology import ( clip_batch, @@ -264,6 +266,8 @@ # Meshing "marching_cubes_batch", "marching_cubes_single", + "dual_contour_batch", + "dual_contour_single", "integrate_tsdf_batch", "integrate_tsdf_single", "integrate_tsdf_with_features_batch", diff --git a/fvdb/functional/_meshing.py b/fvdb/functional/_meshing.py index 286ddee03..f448a32a8 100644 --- a/fvdb/functional/_meshing.py +++ b/fvdb/functional/_meshing.py @@ -70,6 +70,69 @@ def marching_cubes_single( return result[0].jdata, result[1].jdata, result[2].jdata +def dual_contour_batch( + grid: GridBatch, + field: JaggedTensor, + iso: float = 0.0, + reduce: int = 1, + adaptivity: float = 0.0, +) -> tuple[JaggedTensor, JaggedTensor, JaggedTensor]: + """Extract a dual-contouring (DC + QEF) mesh from an SDF on a grid batch. + + Places one vertex per surface cell by minimizing a quadratic error function from the cell's edge + zero-crossings and SDF-gradient normals (sharp-feature preserving), with optional cluster-collapse + decimation. Requires a narrow-band SDF with a >= ~3-voxel band for a watertight result (e.g. from + :func:`retopologize_sdf_batch`). + + Args: + grid (GridBatch): The grid batch defining the sparse topology. + field (JaggedTensor): Per-voxel signed distance values. + iso (float): Isovalue at which to extract the surface. Default ``0.0``. + reduce (int): Uniform ``F x F x F`` cluster-collapse decimation factor (``1`` = full detail). + adaptivity (float): Curvature-adaptive decimation in ``[0, 1.5]`` (``0`` = uniform/off): + collapse flat blocks while keeping detail at features. + + Returns: + vertices (JaggedTensor): Mesh vertex positions, shape ``(B, -1, 3)``. + faces (JaggedTensor): Triangle face indices, shape ``(B, -1, 3)``. + normals (JaggedTensor): Per-vertex (normalized SDF gradient) normals, shape ``(B, -1, 3)``. + + .. seealso:: :func:`dual_contour_single`, :func:`marching_cubes_batch` + """ + grid_data = grid.data + result = _fvdb_cpp.dual_contour(grid_data, field._impl, iso, reduce, adaptivity) + return JaggedTensor(impl=result[0]), JaggedTensor(impl=result[1]), JaggedTensor(impl=result[2]) + + +def dual_contour_single( + grid: Grid, + field: torch.Tensor, + iso: float = 0.0, + reduce: int = 1, + adaptivity: float = 0.0, +) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor]: + """Extract a dual-contouring (DC + QEF) mesh from an SDF on a single grid. + + Args: + grid (Grid): The single grid defining the sparse topology. + field (torch.Tensor): Per-voxel signed distance values. + iso (float): Isovalue at which to extract the surface. Default ``0.0``. + reduce (int): Uniform ``F x F x F`` cluster-collapse decimation factor (``1`` = full detail). + adaptivity (float): Curvature-adaptive decimation in ``[0, 1.5]`` (``0`` = uniform/off). + + Returns: + vertices (torch.Tensor): Vertex positions, shape ``(V, 3)``. + faces (torch.Tensor): Triangle face indices, shape ``(T, 3)``. + normals (torch.Tensor): Per-vertex (normalized SDF gradient) normals, shape ``(V, 3)``. + + .. seealso:: :func:`dual_contour_batch`, :func:`marching_cubes_single` + """ + grid_data = grid.data + field_jt = JaggedTensor(field) + result = _fvdb_cpp.dual_contour(grid_data, field_jt._impl, iso, reduce, adaptivity) + return result[0].jdata, result[1].jdata, result[2].jdata + + def integrate_tsdf_batch( grid: GridBatch, truncation_distance: float, diff --git a/fvdb/functional/_sdf.py b/fvdb/functional/_sdf.py index d39dc472a..932b64067 100644 --- a/fvdb/functional/_sdf.py +++ b/fvdb/functional/_sdf.py @@ -119,7 +119,7 @@ def retopologize_sdf_batch( """Retopologize a signed field into a clean narrow-band SDF on a (possibly pruned) grid batch. If ``pad`` is ``True`` the grid is first dilated by ``band`` voxels (so the eikonal solve has room - to propagate a full-width band), then :func:`reinitialize_sdf_batch` is run, and finally, if + to fill the full ``±band*vx`` band), then :func:`reinitialize_sdf_batch` is run, and finally, if ``prune`` is ``True``, the grid is pruned to the voxels strictly inside the band (``|phi| < band*vx*0.999``). The prune reuses :meth:`GridBatch.pruned_grid`; the resulting field is selected in the grid's canonical voxel order so it stays aligned with the pruned grid. @@ -135,12 +135,12 @@ def retopologize_sdf_batch( :attr:`~fvdb.SmoothingMode.TAUBIN` (volume-preserving). Only used when ``smooth > 0``. redistance_iters (int): Number of redistancing sweeps. ``<= 0`` uses the default. pad (bool): If ``True`` (default) dilate the grid by ``band`` voxels before redistancing so - the output narrow band is a full ``band`` voxels wide even if the input grid had a - thinner active region. Newly added voxels are seeded as *exterior* (``+band*vx``), which - is correct when the dilation extends outward -- i.e. when the grid's interior (the - ``phi < 0`` region) is already represented (the usual case for occupancy/TSDF/mesh-derived - fields). For a thin shell that does not fill its interior, pass ``pad=False`` and supply a - grid that already has an adequate band. + the output band reaches the full ``band`` voxels on each side of the surface (``±band*vx``) + even if the input grid had a thinner active region. Newly added voxels are seeded as + *exterior* (``+band*vx``), which is correct when the dilation extends outward -- i.e. when + the grid's interior (the ``phi < 0`` region) is already represented (the usual case for + occupancy/TSDF/mesh-derived fields). For a thin shell that does not fill its interior, pass + ``pad=False`` and supply a grid that already has an adequate band. prune (bool): If ``True`` prune to the narrow band; if ``False`` return the (possibly padded) grid and the re-initialized field unchanged. @@ -179,7 +179,7 @@ def retopologize_sdf_single( """Retopologize a signed field into a clean narrow-band SDF on a (possibly pruned) single grid. If ``pad`` is ``True`` the grid is first dilated by ``band`` voxels (so the eikonal solve has room - to propagate a full-width band), then :func:`reinitialize_sdf_single` is run, and finally, if + to fill the full ``±band*vx`` band), then :func:`reinitialize_sdf_single` is run, and finally, if ``prune`` is ``True``, the grid is pruned to the voxels strictly inside the band (``|phi| < band*vx*0.999``). @@ -194,12 +194,12 @@ def retopologize_sdf_single( :attr:`~fvdb.SmoothingMode.TAUBIN` (volume-preserving). Only used when ``smooth > 0``. redistance_iters (int): Number of redistancing sweeps. ``<= 0`` uses the default. pad (bool): If ``True`` (default) dilate the grid by ``band`` voxels before redistancing so - the output narrow band is a full ``band`` voxels wide even if the input grid had a - thinner active region. Newly added voxels are seeded as *exterior* (``+band*vx``), which - is correct when the dilation extends outward -- i.e. when the grid's interior (the - ``phi < 0`` region) is already represented (the usual case for occupancy/TSDF/mesh-derived - fields). For a thin shell that does not fill its interior, pass ``pad=False`` and supply a - grid that already has an adequate band. + the output band reaches the full ``band`` voxels on each side of the surface (``±band*vx``) + even if the input grid had a thinner active region. Newly added voxels are seeded as + *exterior* (``+band*vx``), which is correct when the dilation extends outward -- i.e. when + the grid's interior (the ``phi < 0`` region) is already represented (the usual case for + occupancy/TSDF/mesh-derived fields). For a thin shell that does not fill its interior, pass + ``pad=False`` and supply a grid that already has an adequate band. prune (bool): If ``True`` prune to the narrow band; if ``False`` return the (possibly padded) grid and the re-initialized field unchanged. diff --git a/fvdb/grid.py b/fvdb/grid.py index a1e8d85dc..c98ac21c6 100644 --- a/fvdb/grid.py +++ b/fvdb/grid.py @@ -1483,6 +1483,34 @@ def marching_cubes( return functional.marching_cubes_single(self, field, level) + def dual_contour( + self, + field: torch.Tensor, + iso: float = 0.0, + reduce: int = 1, + adaptivity: float = 0.0, + ) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor]: + """Extract a dual-contouring (DC + QEF) mesh from a per-voxel SDF on this grid. + + Places one vertex per surface cell via quadratic-error-function minimization (sharp-feature + preserving), with optional cluster-collapse decimation. Needs a narrow-band SDF with a + >= ~3-voxel band for a watertight result (e.g. from :meth:`retopologize_sdf`). + + Args: + field (torch.Tensor): Per-voxel signed distance values, shape ``(num_voxels,)``. + iso (float): Isovalue at which to extract the surface. + reduce (int): Uniform ``F x F x F`` cluster-collapse decimation factor (``1`` = full detail). + adaptivity (float): Curvature-adaptive decimation in ``[0, 1.5]`` (``0`` = uniform/off). + + Returns: + vertices (torch.Tensor): Vertex positions, shape ``(V, 3)``. + faces (torch.Tensor): Triangle face indices, shape ``(T, 3)``. + normals (torch.Tensor): Per-vertex (normalized SDF gradient) normals, shape ``(V, 3)``. + """ + from . import functional + + return functional.dual_contour_single(self, field, iso, reduce, adaptivity) + def reinitialize_sdf( self, field: torch.Tensor, @@ -1528,7 +1556,7 @@ def retopologize_sdf( """Retopologize a signed field into a clean narrow-band SDF on a (possibly pruned) grid. If ``pad`` is ``True`` this grid is first dilated by ``band`` voxels so the redistance has - room to build a full-width band, then :meth:`reinitialize_sdf` is run, and finally, if + room to fill the full ``±band*vx`` band, then :meth:`reinitialize_sdf` is run, and finally, if ``prune`` is ``True``, the grid is pruned to the voxels strictly inside the band (``|phi| < band*vx*0.999``). @@ -1541,10 +1569,11 @@ def retopologize_sdf( :attr:`~fvdb.SmoothingMode.MEAN_CURVATURE` (default) or :attr:`~fvdb.SmoothingMode.TAUBIN` (volume-preserving). Only used when ``smooth > 0``. redistance_iters (int): Number of redistancing sweeps; ``<= 0`` uses the default. - pad (bool): If ``True`` (default) dilate by ``band`` first so the output band is a full - ``band`` voxels wide even if the input grid was thinner. New voxels are seeded as - exterior (``+band*vx``), which is correct when the interior (``phi < 0``) is already - represented; for a hollow thin shell, pass ``pad=False`` with a pre-banded grid. + pad (bool): If ``True`` (default) dilate by ``band`` first so the output band reaches the + full ``band`` voxels on each side of the surface (``±band*vx``) even if the input grid + was thinner. New voxels are seeded as exterior (``+band*vx``), which is correct when the + interior (``phi < 0``) is already represented; for a hollow thin shell, pass + ``pad=False`` with a pre-banded grid. prune (bool): If ``True`` prune to the narrow band, else return the (padded) grid. Returns: diff --git a/fvdb/grid_batch.py b/fvdb/grid_batch.py index c7a4495a4..f71b6666c 100644 --- a/fvdb/grid_batch.py +++ b/fvdb/grid_batch.py @@ -1029,6 +1029,36 @@ def marching_cubes( return functional.marching_cubes_batch(self, field, level) + def dual_contour( + self, + field: JaggedTensor, + iso: float = 0.0, + reduce: int = 1, + adaptivity: float = 0.0, + ) -> tuple[JaggedTensor, JaggedTensor, JaggedTensor]: + """Extract a dual-contouring (DC + QEF) mesh from a per-voxel SDF on this grid batch. + + Places one vertex per surface cell via quadratic-error-function minimization (sharp-feature + preserving), with optional cluster-collapse decimation. Needs a narrow-band SDF with a + >= ~3-voxel band for a watertight result (e.g. from :meth:`retopologize_sdf`). + + Args: + field (JaggedTensor): Per-voxel signed distance values. + iso (float): Isovalue at which to extract the surface. + reduce (int): Uniform ``F x F x F`` cluster-collapse decimation factor (``1`` = full detail). + adaptivity (float): Curvature-adaptive decimation in ``[0, 1.5]`` (``0`` = uniform/off). + + Returns: + vertices (JaggedTensor): Mesh vertex positions, shape ``(batch_size, -1, 3)``. + faces (JaggedTensor): Triangle face indices, shape ``(batch_size, -1, 3)``. + normals (JaggedTensor): Per-vertex (normalized SDF gradient) normals, shape ``(batch_size, -1, 3)``. + + .. seealso:: :meth:`Grid.dual_contour` + """ + from . import functional + + return functional.dual_contour_batch(self, field, iso, reduce, adaptivity) + def reinitialize_sdf( self, field: JaggedTensor, @@ -1076,7 +1106,7 @@ def retopologize_sdf( """Retopologize a signed field into a clean narrow-band SDF on a (possibly pruned) grid batch. If ``pad`` is ``True`` the grid is first dilated by ``band`` voxels so the redistance has - room to build a full-width band, then :meth:`reinitialize_sdf` is run, and finally, if + room to fill the full ``±band*vx`` band, then :meth:`reinitialize_sdf` is run, and finally, if ``prune`` is ``True``, the grid is pruned to the voxels strictly inside the band (``|phi| < band*vx*0.999``). @@ -1089,10 +1119,11 @@ def retopologize_sdf( :attr:`~fvdb.SmoothingMode.MEAN_CURVATURE` (default) or :attr:`~fvdb.SmoothingMode.TAUBIN` (volume-preserving). Only used when ``smooth > 0``. redistance_iters (int): Number of redistancing sweeps; ``<= 0`` uses the default. - pad (bool): If ``True`` (default) dilate by ``band`` first so the output band is a full - ``band`` voxels wide even if the input grid was thinner. New voxels are seeded as - exterior (``+band*vx``), which is correct when the interior (``phi < 0``) is already - represented; for a hollow thin shell, pass ``pad=False`` with a pre-banded grid. + pad (bool): If ``True`` (default) dilate by ``band`` first so the output band reaches the + full ``band`` voxels on each side of the surface (``±band*vx``) even if the input grid + was thinner. New voxels are seeded as exterior (``+band*vx``), which is correct when the + interior (``phi < 0``) is already represented; for a hollow thin shell, pass + ``pad=False`` with a pre-banded grid. prune (bool): If ``True`` prune to the narrow band, else return the (padded) grid batch. Returns: diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 3060e5bec..c65eef869 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -80,6 +80,7 @@ set(FVDB_CU_FILES fvdb/detail/ops/ClipGrid.cu fvdb/detail/ops/CloneGrid.cu fvdb/detail/ops/ConcatenateGrids.cu + fvdb/detail/ops/DualContour.cu fvdb/detail/ops/IndexGrid.cu fvdb/detail/ops/MakeContiguous.cu fvdb/detail/ops/SerializeGrid.cu diff --git a/src/fvdb/detail/ops/DualContour.cu b/src/fvdb/detail/ops/DualContour.cu new file mode 100644 index 000000000..f985c3004 --- /dev/null +++ b/src/fvdb/detail/ops/DualContour.cu @@ -0,0 +1,1285 @@ +// Copyright Contributors to the OpenVDB Project +// SPDX-License-Identifier: Apache-2.0 +// +// Dual-contouring (DC + QEF) mesher: meshes an OnIndex grid carrying a narrow-band SDF into a +// triangle mesh, with QEF vertex placement and optional cluster-collapse decimation. The 3x3x3 box +// stencil is gathered once per voxel via the NanoVDB VoxelBlockManager. +// +// This follows the dual-contouring approach of OpenVDB's tools/VolumeToMesh.h, but places vertices +// by minimising a quadratic error function (classic DC) rather than averaging the edge crossings +// (the "mass point" VolumeToMesh uses by default), and decimates by simple cluster-collapse rather +// than VolumeToMesh's seam-stitched octree region merge. +// +// References: +// [Ju et al. 2002] T. Ju, F. Losasso, S. Schaefer, J. Warren, "Dual Contouring of Hermite +// Data", ACM TOG 21(3) (SIGGRAPH 2002), 339-346. The core method: one +// vertex per cell minimising a QEF over the edge Hermite data (crossings + +// normals), the dual connectivity, the numerically-stable centroid-biased +// QEF, and the octree-based adaptive simplification our decimation is a +// simplified form of. +// [Garland & Heckbert 1997] M. Garland, P. Heckbert, "Surface Simplification Using Quadric Error +// Metrics", SIGGRAPH 1997, 209-216. The quadric error metric (A = sum n +// n^T) that the vertex placement minimises. +// [Kobbelt et al. 2001] L. Kobbelt, M. Botsch, U. Schwanecke, H.-P. Seidel, "Feature Sensitive +// Surface Extraction from Volume Data", SIGGRAPH 2001. Placing the vertex +// at the intersection of the edge-crossing tangent planes (the feature +// point). +// [Lorensen & Cline 1987] W. Lorensen, H. Cline, "Marching Cubes: A High Resolution 3D Surface +// Construction Algorithm", SIGGRAPH 1987, 163-169. The sign-based +// isosurface extraction that dual contouring is the dual of. +// +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +namespace fvdb { +namespace detail { +namespace ops { + +namespace { + +// ------------------------- box-stencil geometry (compile-time) ------------------------- +// A 3x3x3 box-stencil neighbour is addressed by a spoke index in [0,27): +// spoke = (di+1)*9 + (dj+1)*3 + (dk+1); the centre (the voxel itself) is spoke 13. The mesher only +// touches a subset of the 27 spokes (8 cube corners, 6 faces, 3 connectivity fans), so the gather +// emits just those `numColumns` columns. All of this is static, so it is computed once at compile +// time into `kGeometry` and passed by value into the kernels (no __constant__ / per-call setup). +constexpr int +spokeIndex(int di, int dj, int dk) { + return (di + 1) * 9 + (dj + 1) * 3 + (dk + 1); +} + +/// @brief Column of `spoke` within the sorted-unique `usedSpokes[0..numColumns)`, or -1 if absent. +/// Binary search (std::lower_bound is constexpr in C++20) since usedSpokes is sorted. +constexpr int +columnOfSpoke(const int *usedSpokes, int numColumns, int spoke) { + const int *end = usedSpokes + numColumns; + const int *found = std::lower_bound(usedSpokes, end, spoke); + return (found != end && *found == spoke) ? int(found - usedSpokes) : -1; +} + +struct BoxStencilGeometry { + int numColumns{0}; + int usedSpokes[27]{}; // the `numColumns` gathered spokes (sorted, unique) + int cornerColumn[8]{}; // column of each cube corner in the gathered (numVoxels,numColumns) + // table + int faceColumn[6]{}; // column of each face neighbour (-x,+x,-y,+y,-z,+z) + int centerColumn{0}; // column of spoke 13 (the voxel itself) + int fanForwardColumn[3]{}; // forward-edge spoke column for the x/y/z minimal edges + int fanColumn[3][4]{}; // the 4 surrounding-cell spoke columns per edge fan + int edgeCornerA[12]{}; // the 12 cube edges as corner-index pairs (A,B) + int edgeCornerB[12]{}; + int cornerSpoke[8]{}; // raw spoke of each corner (the fused gather reads the stencil directly) + int faceSpoke[6]{}; + float cornerOffset[8][3]{}; // the 8 corner offsets in [0,1]^3 +}; + +constexpr BoxStencilGeometry +makeBoxStencilGeometry() { + BoxStencilGeometry geometry{}; + const int cornerOffset[8][3] = { + {0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {1, 1, 0}, {1, 0, 1}, {0, 1, 1}, {1, 1, 1}}; + const int cubeEdge[12][2] = {{0, 1}, + {2, 4}, + {3, 5}, + {6, 7}, + {0, 2}, + {1, 4}, + {3, 6}, + {5, 7}, + {0, 3}, + {1, 5}, + {2, 6}, + {4, 7}}; + for (int corner = 0; corner < 8; ++corner) { + for (int axis = 0; axis < 3; ++axis) { + geometry.cornerOffset[corner][axis] = (float)cornerOffset[corner][axis]; + } + geometry.cornerSpoke[corner] = + spokeIndex(cornerOffset[corner][0], cornerOffset[corner][1], cornerOffset[corner][2]); + } + for (int edge = 0; edge < 12; ++edge) { + geometry.edgeCornerA[edge] = cubeEdge[edge][0]; + geometry.edgeCornerB[edge] = cubeEdge[edge][1]; + } + geometry.faceSpoke[0] = spokeIndex(-1, 0, 0); + geometry.faceSpoke[1] = spokeIndex(1, 0, 0); + geometry.faceSpoke[2] = spokeIndex(0, -1, 0); + geometry.faceSpoke[3] = spokeIndex(0, 1, 0); + geometry.faceSpoke[4] = spokeIndex(0, 0, -1); + geometry.faceSpoke[5] = spokeIndex(0, 0, 1); + const int fanForward[3] = {spokeIndex(1, 0, 0), spokeIndex(0, 1, 0), spokeIndex(0, 0, 1)}; + const int fanSpoke[3][4] = { + {13, spokeIndex(0, 0, -1), spokeIndex(0, -1, -1), spokeIndex(0, -1, 0)}, + {13, spokeIndex(0, 0, -1), spokeIndex(-1, 0, -1), spokeIndex(-1, 0, 0)}, + {13, spokeIndex(0, -1, 0), spokeIndex(-1, -1, 0), spokeIndex(-1, 0, 0)}}; + + // sorted-unique union of all touched spokes -> usedSpokes[0..numColumns) + int allSpokes[64]{}; + int spokeCount = 0; + for (int i = 0; i < 6; ++i) + allSpokes[spokeCount++] = geometry.faceSpoke[i]; + for (int i = 0; i < 8; ++i) + allSpokes[spokeCount++] = geometry.cornerSpoke[i]; + for (int dir = 0; dir < 3; ++dir) { + allSpokes[spokeCount++] = fanForward[dir]; + for (int i = 0; i < 4; ++i) + allSpokes[spokeCount++] = fanSpoke[dir][i]; + } + // sort then drop duplicates (std::sort / std::unique are constexpr in C++20) + std::sort(allSpokes, allSpokes + spokeCount); + const int numColumns = int(std::unique(allSpokes, allSpokes + spokeCount) - allSpokes); + for (int i = 0; i < numColumns; ++i) + geometry.usedSpokes[i] = allSpokes[i]; + geometry.numColumns = numColumns; + + for (int corner = 0; corner < 8; ++corner) { + geometry.cornerColumn[corner] = + columnOfSpoke(geometry.usedSpokes, numColumns, geometry.cornerSpoke[corner]); + } + for (int face = 0; face < 6; ++face) { + geometry.faceColumn[face] = + columnOfSpoke(geometry.usedSpokes, numColumns, geometry.faceSpoke[face]); + } + geometry.centerColumn = columnOfSpoke(geometry.usedSpokes, numColumns, 13); + for (int dir = 0; dir < 3; ++dir) { + geometry.fanForwardColumn[dir] = + columnOfSpoke(geometry.usedSpokes, numColumns, fanForward[dir]); + for (int i = 0; i < 4; ++i) { + geometry.fanColumn[dir][i] = + columnOfSpoke(geometry.usedSpokes, numColumns, fanSpoke[dir][i]); + } + } + return geometry; +} + +static constexpr BoxStencilGeometry kGeometry = makeBoxStencilGeometry(); + +/// @brief Solve a symmetric 3x3 SPD system A x = b (A given by its 6 unique entries) via Cramer's +/// rule. +__device__ inline void +solveSymmetric3x3(const double A[6], const double b[3], double x[3]) { + const double a0 = A[0], a1 = A[1], a2 = A[2], a3 = A[3], a4 = A[4], a5 = A[5]; + const double cof00 = a3 * a5 - a4 * a4, cof01 = a4 * a2 - a1 * a5, cof02 = a1 * a4 - a3 * a2; + double determinant = a0 * cof00 + a1 * cof01 + a2 * cof02; + if (fabs(determinant) < 1e-30) { + x[0] = x[1] = x[2] = 0.0; + return; + } + const double invDet = 1.0 / determinant; + const double cof11 = a0 * a5 - a2 * a2, cof12 = a1 * a2 - a0 * a4, cof22 = a0 * a3 - a1 * a1; + x[0] = invDet * (cof00 * b[0] + cof01 * b[1] + cof02 * b[2]); + x[1] = invDet * (cof01 * b[0] + cof11 * b[1] + cof12 * b[2]); + x[2] = invDet * (cof02 * b[0] + cof12 * b[1] + cof22 * b[2]); +} + +// ===================== fused VBM decode ==================================================== +/// @brief Decoding the VBM stencil is the expensive part, so each voxel does it ONCE here and +/// writes four value-indexed outputs that every later kernel then reads in O(1) (no re-walking the +/// tree): +/// - neighborTable [valueCount x numColumns] : the gathered box-stencil neighbour value-indices +/// - gradient : central-difference SDF gradient (the surface +/// normal) +/// - surfaceFlag : 1 iff the iso surface passes through this voxel's +/// cube +/// - voxelCoord : the voxel's index-space (i,j,k) +__global__ void +gatherFusedKernel(const OnIndexGridT *grid, + const uint32_t *firstLeafID, + const uint64_t *jumpMap, + uint64_t firstOffset, + const float *sdf, + BoxStencilGeometry geometry, + float iso, + int32_t *neighborTable, + float *gradient, + uint8_t *surfaceFlag, + int32_t *voxelCoord) { + __shared__ VbmBlockMaps maps; + if (!vbmDecodeBlock(grid, firstLeafID, jumpMap, firstOffset, maps)) + return; + const int threadId = threadIdx.x; + using VoxelBlockManagerT = nanovdb::tools::cuda::VoxelBlockManager; + uint64_t stencil[27]; + VoxelBlockManagerT::template computeBoxStencil( + grid, maps.leafIndex, maps.voxelOffset, stencil); + const uint64_t centerIndex = stencil[13]; // centre value index (>= 1 for active) + const int numColumns = geometry.numColumns; + + for (int column = 0; column < numColumns; ++column) { + neighborTable[centerIndex * numColumns + column] = + (int32_t)stencil[geometry.usedSpokes[column]]; + } + + // central-difference gradient (inactive face stencil==0 -> use the centre value) + const float centerSdf = sdf[centerIndex]; + auto faceValue = [&](int spoke) { + uint64_t neighbor = stencil[spoke]; + return neighbor > 0 ? sdf[neighbor] : centerSdf; + }; + gradient[centerIndex * 3 + 0] = + 0.5f * (faceValue(geometry.faceSpoke[1]) - faceValue(geometry.faceSpoke[0])); + gradient[centerIndex * 3 + 1] = + 0.5f * (faceValue(geometry.faceSpoke[3]) - faceValue(geometry.faceSpoke[2])); + gradient[centerIndex * 3 + 2] = + 0.5f * (faceValue(geometry.faceSpoke[5]) - faceValue(geometry.faceSpoke[4])); + + // surface-cell flag: the cube anchored at this voxel is a surface cell iff its 8 corners + // straddle iso (0 < #inside < 8) AND are all active -- a full active cube is needed to place a + // QEF vertex and to let the cell take part in quads. + int numInsideCorners = 0; + bool allCornersActive = true; +#pragma unroll + for (int corner = 0; corner < 8; ++corner) { + uint64_t cornerIndex = stencil[geometry.cornerSpoke[corner]]; + numInsideCorners += (sdf[cornerIndex] < iso) ? 1 : 0; + if (corner) + allCornersActive &= (cornerIndex > 0); + } + surfaceFlag[centerIndex] = + (allCornersActive && numInsideCorners > 0 && numInsideCorners < 8) ? 1u : 0u; + + const auto &leaf = grid->tree().getFirstNode<0>()[maps.leafIndex[threadId]]; + const nanovdb::Coord globalCoord = leaf.offsetToGlobalCoord(maps.voxelOffset[threadId]); + voxelCoord[centerIndex * 3 + 0] = globalCoord[0]; + voxelCoord[centerIndex * 3 + 1] = globalCoord[1]; + voxelCoord[centerIndex * 3 + 2] = globalCoord[2]; +} + +/// @brief (QEF) Place one dual-contouring vertex per surface cell. The vertex is the point x +/// minimising the quadratic error function E(x) = sum_i [ n_i . (x - p_i) ]^2 over the cell's +/// edge intersections, where p_i is the i-th edge zero-crossing and n_i the unit surface normal +/// there. Each term is the squared distance from x to the tangent plane at crossing i, so the +/// minimiser is the point that best fits all those planes -- it snaps onto sharp features +/// (edges/corners) that a plain crossing-average would round off. Setting dE/dx = 0 gives the 3x3 +/// normal equations +/// A x = b, A = sum_i n_i n_i^T, b = sum_i n_i (n_i . p_i), +/// solved below (re-centred on the crossing centroid and Tikhonov-regularised for stability). Also +/// stores the centre voxel's SDF gradient as the per-vertex reference normal (used to orient the +/// triangles, and emitted as the output normal once normalised). +/// Method: [Ju et al. 2002]; the QEF is the quadric error metric of [Garland & Heckbert 1997]; the +/// tangent-plane feature point is [Kobbelt et al. 2001]. (See the reference list at the top.) +__global__ void +placeQefVerticesKernel(const int32_t *surfaceCells, + int64_t numSurfaceCells, + const int32_t *neighborTable, + const float *sdf, + const float *gradient, + const int32_t *voxelCoord, + BoxStencilGeometry geometry, + float iso, + VoxelCoordTransform transform, + float *outVertices, + float *outRefNormal) { + const int64_t cellOrdinal = int64_t(blockIdx.x) * blockDim.x + threadIdx.x; + if (cellOrdinal >= numSurfaceCells) + return; + const int numColumns = geometry.numColumns; + const int32_t cellValueIndex = surfaceCells[cellOrdinal]; + // gather the 8 cube-corner SDF values and gradients (a surface cell's corners are all active) + float cornerSdf[8]; + float cornerGradient[8][3]; +#pragma unroll + for (int corner = 0; corner < 8; ++corner) { + int32_t cornerIndex = + neighborTable[int64_t(cellValueIndex) * numColumns + geometry.cornerColumn[corner]]; + cornerSdf[corner] = sdf[cornerIndex]; + cornerGradient[corner][0] = gradient[int64_t(cornerIndex) * 3 + 0]; + cornerGradient[corner][1] = gradient[int64_t(cornerIndex) * 3 + 1]; + cornerGradient[corner][2] = gradient[int64_t(cornerIndex) * 3 + 2]; + } + double normalMatrix[6] = {0, 0, 0, 0, 0, 0}; // sum of n n^T (6 unique entries) + double normalRhs[3] = {0, 0, 0}; // sum of n (n . crossing) + double crossingSum[3] = {0, 0, 0}; // sum of edge crossings + int numCrossings = 0; + // accumulate A and b over the cell's 12 edges: each edge that crosses the iso surface + // contributes one tangent-plane constraint (its crossing point p_i and the normal n_i there) +#pragma unroll + for (int edge = 0; edge < 12; ++edge) { + const int cornerA = geometry.edgeCornerA[edge], cornerB = geometry.edgeCornerB[edge]; + const float sdfA = cornerSdf[cornerA], sdfB = cornerSdf[cornerB]; + if ((sdfA < iso) == (sdfB < iso)) + continue; // endpoints on the same side of iso -> this edge does not cross the surface + // crossingT = fraction along A->B where the linearly-interpolated sdf reaches iso (the edge + // root); fall back to the midpoint if the two endpoints are (numerically) equal + float sdfDelta = sdfB - sdfA; + float crossingT = (fabsf(sdfDelta) > 1e-12f) ? (iso - sdfA) / sdfDelta : 0.5f; + crossingT = fminf(fmaxf(crossingT, 0.f), 1.f); + // crossing point p_i and normal n_i, both linearly interpolated between the two corners + // (the corner gradients supply the surface normal) + float crossing[3], edgeNormal[3]; +#pragma unroll + for (int axis = 0; axis < 3; ++axis) { + crossing[axis] = geometry.cornerOffset[cornerA][axis] + + crossingT * (geometry.cornerOffset[cornerB][axis] - + geometry.cornerOffset[cornerA][axis]); + edgeNormal[axis] = + cornerGradient[cornerA][axis] + + crossingT * (cornerGradient[cornerB][axis] - cornerGradient[cornerA][axis]); + } + // normalise n_i (the eps keeps a (near-)zero gradient finite) + float invLen = rsqrtf(edgeNormal[0] * edgeNormal[0] + edgeNormal[1] * edgeNormal[1] + + edgeNormal[2] * edgeNormal[2] + 1e-24f); + double nx = edgeNormal[0] * invLen, ny = edgeNormal[1] * invLen, + nz = edgeNormal[2] * invLen; + double normalDotCrossing = nx * crossing[0] + ny * crossing[1] + nz * crossing[2]; + // A += n n^T (6 unique entries), b-term normalRhs += n (n.p), and track the crossing + // centroid + normalMatrix[0] += nx * nx; + normalMatrix[1] += nx * ny; + normalMatrix[2] += nx * nz; + normalMatrix[3] += ny * ny; + normalMatrix[4] += ny * nz; + normalMatrix[5] += nz * nz; + normalRhs[0] += nx * normalDotCrossing; + normalRhs[1] += ny * normalDotCrossing; + normalRhs[2] += nz * normalDotCrossing; + crossingSum[0] += crossing[0]; + crossingSum[1] += crossing[1]; + crossingSum[2] += crossing[2]; + ++numCrossings; + } + double localVertex[3]; + double invCount = numCrossings > 0 ? 1.0 / numCrossings : 0.0; + double crossingCentroid[3] = { + crossingSum[0] * invCount, crossingSum[1] * invCount, crossingSum[2] * invCount}; + if (numCrossings == 0) { + // no bipolar edges (a real surface cell always has >= 1; defensive) -> just use the + // centroid + localVertex[0] = crossingCentroid[0]; + localVertex[1] = crossingCentroid[1]; + localVertex[2] = crossingCentroid[2]; + } else { + // Solve the normal equations re-centred on the crossing centroid c (much better + // conditioned): with y = x - c, (A + lambda I) y = b - A c = normalRhs - A.centroid, then + // x = c + y. The lambda = 0.05 Tikhonov term keeps the system solvable when A is + // rank-deficient (flat cells or too few independent normals leave the surface position + // underconstrained) and nudges those free directions gently back toward the centroid. + double matTimesCentroid[3] = { + normalMatrix[0] * crossingCentroid[0] + normalMatrix[1] * crossingCentroid[1] + + normalMatrix[2] * crossingCentroid[2], + normalMatrix[1] * crossingCentroid[0] + normalMatrix[3] * crossingCentroid[1] + + normalMatrix[4] * crossingCentroid[2], + normalMatrix[2] * crossingCentroid[0] + normalMatrix[4] * crossingCentroid[1] + + normalMatrix[5] * crossingCentroid[2]}; + double rhs[3] = {normalRhs[0] - matTimesCentroid[0], + normalRhs[1] - matTimesCentroid[1], + normalRhs[2] - matTimesCentroid[2]}; + double regularizedMatrix[6] = {normalMatrix[0] + 0.05, + normalMatrix[1], + normalMatrix[2], + normalMatrix[3] + 0.05, + normalMatrix[4], + normalMatrix[5] + 0.05}; + double solution[3]; + solveSymmetric3x3(regularizedMatrix, rhs, solution); + // x = centroid + y, clamped into the unit cell so the vertex never leaves its own cell + for (int axis = 0; axis < 3; ++axis) + localVertex[axis] = fmin(fmax(crossingCentroid[axis] + solution[axis], 0.0), 1.0); + } + // cell-local vertex -> world: localVertex is in [0,1]^3 within the cell, so add the cell's + // index-space origin and map index space -> world with the grid transform. + const float originX = (float)voxelCoord[int64_t(cellValueIndex) * 3 + 0]; + const float originY = (float)voxelCoord[int64_t(cellValueIndex) * 3 + 1]; + const float originZ = (float)voxelCoord[int64_t(cellValueIndex) * 3 + 2]; + auto worldPos = transform.applyInv(originX + (float)localVertex[0], + originY + (float)localVertex[1], + originZ + (float)localVertex[2]); + outVertices[cellOrdinal * 3 + 0] = worldPos[0]; + outVertices[cellOrdinal * 3 + 1] = worldPos[1]; + outVertices[cellOrdinal * 3 + 2] = worldPos[2]; + // reference normal = the centre voxel's SDF gradient; orients the triangles and (once + // normalised) becomes the emitted per-vertex normal. + outRefNormal[cellOrdinal * 3 + 0] = gradient[int64_t(cellValueIndex) * 3 + 0]; + outRefNormal[cellOrdinal * 3 + 1] = gradient[int64_t(cellValueIndex) * 3 + 1]; + outRefNormal[cellOrdinal * 3 + 2] = gradient[int64_t(cellValueIndex) * 3 + 2]; +} + +/// @brief (connectivity) Dual-contouring connectivity is the dual of marching cubes: wherever a +/// grid edge changes sign (crosses the surface), the 4 cells sharing that edge each hold a vertex, +/// and joining those 4 vertices makes one quad straddling the edge. To emit each quad exactly once, +/// a voxel owns only its 3 "minimal" forward edges (+x/+y/+z, selected by `dir`). Returns true and +/// fills `quadVertices[4]` with the 4 surrounding cells' vertices when this voxel's `dir` edge is +/// bipolar and all 4 of those cells carry a vertex; false otherwise. +/// Dual connectivity is [Ju et al. 2002]; the sign-based isosurface it is dual to is [Lorensen & +/// Cline 1987]. +__device__ inline bool +tryBuildQuad(int64_t voxelIndex, + int dir, + const int32_t *neighborTable, + const float *sdf, + const int32_t *cellVertex, + BoxStencilGeometry geometry, + float iso, + int32_t quadVertices[4]) { + if (voxelIndex == 0) + return false; + const int numColumns = geometry.numColumns; + int32_t forwardIndex = neighborTable[voxelIndex * numColumns + geometry.fanForwardColumn[dir]]; + if (forwardIndex <= 0) + return false; + if ((sdf[voxelIndex] < iso) == (sdf[forwardIndex] < iso)) + return false; // edge is not bipolar +#pragma unroll + for (int i = 0; i < 4; ++i) { + int column = geometry.fanColumn[dir][i]; + int32_t neighborCellIndex = (column == geometry.centerColumn) + ? int32_t(voxelIndex) + : neighborTable[voxelIndex * numColumns + column]; + int32_t vertex = cellVertex[neighborCellIndex]; // cellVertex[0] == -1 (background) + if (vertex < 0) + return false; + quadVertices[i] = vertex; + } + return true; +} + +/// @brief Diagonal-aware split of the cyclic quad; returns false if triangle `whichTriangle` is +/// degenerate (a repeated vertex -- happens only when cells contract onto one cluster). +__device__ inline bool +splitQuadTriangle(const int32_t quadVertices[4], int whichTriangle, int32_t outTriangle[3]) { + bool useDiagonal13 = + (quadVertices[0] == quadVertices[2]) && (quadVertices[1] != quadVertices[3]); + if (whichTriangle == 0) { + outTriangle[0] = quadVertices[0]; + outTriangle[1] = quadVertices[1]; + outTriangle[2] = useDiagonal13 ? quadVertices[3] : quadVertices[2]; + } else { + outTriangle[0] = useDiagonal13 ? quadVertices[1] : quadVertices[0]; + outTriangle[1] = quadVertices[2]; + outTriangle[2] = quadVertices[3]; + } + return outTriangle[0] != outTriangle[1] && outTriangle[1] != outTriangle[2] && + outTriangle[0] != outTriangle[2]; +} + +/// @brief (connectivity, pass 1) Count the triangles each surface cell emits; a cumsum of these +/// counts then gives each cell a contiguous, deterministic write range (no atomics, unlike a single +/// shared counter). Only surface cells can own a quad: every quad requires its 4 fan cells (one of +/// which is the owning voxel itself) to carry a vertex, so iterating surface cells is complete -- +/// and far cheaper than sweeping every active voxel. +__global__ void +countTrianglesKernel(const int32_t *surfaceCells, + int64_t numSurfaceCells, + const int32_t *neighborTable, + const float *sdf, + const int32_t *cellVertex, + BoxStencilGeometry geometry, + float iso, + int64_t *triangleCounts) { + const int64_t i = int64_t(blockIdx.x) * blockDim.x + threadIdx.x; + if (i >= numSurfaceCells) + return; + const int64_t voxelIndex = surfaceCells[i]; + int32_t quadVertices[4], triangle[3]; + int localCount = 0; +#pragma unroll + for (int dir = 0; dir < 3; ++dir) { + if (tryBuildQuad( + voxelIndex, dir, neighborTable, sdf, cellVertex, geometry, iso, quadVertices)) { + for (int whichTriangle = 0; whichTriangle < 2; ++whichTriangle) { + if (splitQuadTriangle(quadVertices, whichTriangle, triangle)) + ++localCount; + } + } + } + triangleCounts[i] = localCount; +} + +/// @brief (connectivity, pass 2) Write each surface cell's triangles into its own range, starting +/// at the exclusive prefix sum `triangleOffsets[i]`. +__global__ void +writeTrianglesKernel(const int32_t *surfaceCells, + int64_t numSurfaceCells, + const int32_t *neighborTable, + const float *sdf, + const int32_t *cellVertex, + BoxStencilGeometry geometry, + float iso, + const int64_t *triangleOffsets, + int32_t *triangles) { + const int64_t i = int64_t(blockIdx.x) * blockDim.x + threadIdx.x; + if (i >= numSurfaceCells) + return; + const int64_t voxelIndex = surfaceCells[i]; + int32_t quadVertices[4], triangle[3]; + int64_t slot = triangleOffsets[i]; +#pragma unroll + for (int dir = 0; dir < 3; ++dir) { + if (!tryBuildQuad( + voxelIndex, dir, neighborTable, sdf, cellVertex, geometry, iso, quadVertices)) + continue; + for (int whichTriangle = 0; whichTriangle < 2; ++whichTriangle) { + if (!splitQuadTriangle(quadVertices, whichTriangle, triangle)) + continue; + triangles[slot * 3 + 0] = triangle[0]; + triangles[slot * 3 + 1] = triangle[1]; + triangles[slot * 3 + 2] = triangle[2]; + ++slot; + } + } +} + +/// @brief (orient) Flip each triangle so its geometric normal agrees with the summed reference +/// normal. +__global__ void +orientTrianglesKernel(int64_t numTriangles, + const float *vertices, + const float *refNormal, + int32_t *triangles) { + const int64_t face = int64_t(blockIdx.x) * blockDim.x + threadIdx.x; + if (face >= numTriangles) + return; + using Vec3f = nanovdb::math::Vec3; + const int32_t i0 = triangles[face * 3 + 0], i1 = triangles[face * 3 + 1], + i2 = triangles[face * 3 + 2]; + auto vertex = [&](int32_t i) { + return Vec3f(vertices[i * 3], vertices[i * 3 + 1], vertices[i * 3 + 2]); + }; + auto normal = [&](int32_t i) { + return Vec3f(refNormal[i * 3], refNormal[i * 3 + 1], refNormal[i * 3 + 2]); + }; + const Vec3f p0 = vertex(i0); + const Vec3f faceNormal = (vertex(i1) - p0).cross(vertex(i2) - p0); + const Vec3f refSum = normal(i0) + normal(i1) + normal(i2); + if (faceNormal.dot(refSum) < 0.f) { + triangles[face * 3 + 1] = i2; + triangles[face * 3 + 2] = i1; + } +} + +// ===================== decimation (reduce / adaptivity) ==================================== +// Cluster-collapse decimation: instead of one vertex per surface cell, group cells into clusters +// and emit one vertex per cluster. `reduce = F` uses uniform F x F x F index-space blocks; +// `adaptivity` keeps feature cells at full detail and only collapses "flat" coarse blocks (those +// whose cells' edge normals are well aligned). Each cell's QEF is summed into its cluster +// (clusterQefAccumKernel) and a single merged vertex is solved per cluster (clusterSolveKernel); +// connectivity then merges for free, since quads referencing cells that collapsed to the same +// vertex become degenerate and are dropped. NOTE: the per-cluster accumulation uses double +// atomicAdd, so the summation order -- and hence the exact decimated vertex positions -- are not +// bit-reproducible run to run. +// This is a simplified, single-level form of the octree-based adaptive DC simplification of [Ju et +// al. 2002]; OpenVDB's VolumeToMesh instead does the full seam-stitched octree region merge. +__device__ inline void +loadCellCorners(int32_t cellValueIndex, + const int32_t *neighborTable, + const float *sdf, + const float *gradient, + BoxStencilGeometry geometry, + float cornerSdf[8], + float cornerGradient[8][3]) { + const int numColumns = geometry.numColumns; +#pragma unroll + for (int corner = 0; corner < 8; ++corner) { + int32_t cornerIndex = + neighborTable[int64_t(cellValueIndex) * numColumns + geometry.cornerColumn[corner]]; + cornerSdf[corner] = sdf[cornerIndex]; + cornerGradient[corner][0] = gradient[int64_t(cornerIndex) * 3 + 0]; + cornerGradient[corner][1] = gradient[int64_t(cornerIndex) * 3 + 1]; + cornerGradient[corner][2] = gradient[int64_t(cornerIndex) * 3 + 2]; + } +} + +/// @brief Per-cell flatness signal for adaptive decimation: sum the cell's unit edge-crossing +/// normals and normalize. If the crossings share a direction (a flat patch) the sum stays near unit +/// length; if they disagree (a feature) it partly cancels. Block flatness is then the mean of these +/// over a block. +__global__ void +cellNormalKernel(const int32_t *surfaceCells, + int64_t numSurfaceCells, + const int32_t *neighborTable, + const float *sdf, + const float *gradient, + BoxStencilGeometry geometry, + float iso, + float *cellNormals) { + int64_t i = int64_t(blockIdx.x) * blockDim.x + threadIdx.x; + if (i >= numSurfaceCells) + return; + float cornerSdf[8], cornerGradient[8][3]; + loadCellCorners( + surfaceCells[i], neighborTable, sdf, gradient, geometry, cornerSdf, cornerGradient); + float normalSum[3] = {0, 0, 0}; +#pragma unroll + for (int edge = 0; edge < 12; ++edge) { + int cornerA = geometry.edgeCornerA[edge], cornerB = geometry.edgeCornerB[edge]; + float sdfA = cornerSdf[cornerA], sdfB = cornerSdf[cornerB]; + if ((sdfA < iso) == (sdfB < iso)) + continue; + float sdfDelta = sdfB - sdfA; + float crossingT = (fabsf(sdfDelta) > 1e-12f) ? (iso - sdfA) / sdfDelta : 0.5f; + crossingT = fminf(fmaxf(crossingT, 0.f), 1.f); + float nx = cornerGradient[cornerA][0] + + crossingT * (cornerGradient[cornerB][0] - cornerGradient[cornerA][0]), + ny = cornerGradient[cornerA][1] + + crossingT * (cornerGradient[cornerB][1] - cornerGradient[cornerA][1]), + nz = cornerGradient[cornerA][2] + + crossingT * (cornerGradient[cornerB][2] - cornerGradient[cornerA][2]); + float invLen = rsqrtf(nx * nx + ny * ny + nz * nz + 1e-24f); + normalSum[0] += nx * invLen; + normalSum[1] += ny * invLen; + normalSum[2] += nz * invLen; + } + float invLen = rsqrtf(normalSum[0] * normalSum[0] + normalSum[1] * normalSum[1] + + normalSum[2] * normalSum[2] + 1e-24f); + cellNormals[i * 3 + 0] = normalSum[0] * invLen; + cellNormals[i * 3 + 1] = normalSum[1] * invLen; + cellNormals[i * 3 + 2] = normalSum[2] * invLen; +} + +__device__ inline long long +packKey(int x, int y, int z, int spanY, int spanZ) { + return ((long long)x * spanY + y) * spanZ + z; +} + +/// @brief Uniform-block cluster key (also the coarse-block key reused in the adaptive pass). +__global__ void +clusterKeyKernel(const int32_t *surfaceCells, + int64_t numSurfaceCells, + const int32_t *voxelCoord, + int originX, + int originY, + int originZ, + int spanY, + int spanZ, + int blockSize, + long long *clusterKeys) { + int64_t i = int64_t(blockIdx.x) * blockDim.x + threadIdx.x; + if (i >= numSurfaceCells) + return; + int32_t cellValueIndex = surfaceCells[i]; + int localX = voxelCoord[int64_t(cellValueIndex) * 3 + 0] - originX, + localY = voxelCoord[int64_t(cellValueIndex) * 3 + 1] - originY, + localZ = voxelCoord[int64_t(cellValueIndex) * 3 + 2] - originZ; + clusterKeys[i] = + packKey(localX / blockSize, localY / blockSize, localZ / blockSize, spanY, spanZ); +} + +/// @brief Adaptive cluster key: a cell in a flat block takes its coarse-block id (so the whole +/// block collapses to one vertex); a cell in a feature block takes its own per-voxel key (so it +/// stays full detail). The flat-block keys are offset by fineKeyCount so the fine and coarse key +/// spaces can never collide. +__global__ void +clusterKeyAdaptiveKernel(const int32_t *surfaceCells, + int64_t numSurfaceCells, + const int32_t *voxelCoord, + int originX, + int originY, + int originZ, + int spanY, + int spanZ, + int blockSize, + const int32_t *blockOfCell, + const uint8_t *blockIsFlat, + long long fineKeyCount, + long long *clusterKeys) { + int64_t i = int64_t(blockIdx.x) * blockDim.x + threadIdx.x; + if (i >= numSurfaceCells) + return; + int32_t cellValueIndex = surfaceCells[i]; + int localX = voxelCoord[int64_t(cellValueIndex) * 3 + 0] - originX, + localY = voxelCoord[int64_t(cellValueIndex) * 3 + 1] - originY, + localZ = voxelCoord[int64_t(cellValueIndex) * 3 + 2] - originZ; + if (blockIsFlat[blockOfCell[i]]) + clusterKeys[i] = + fineKeyCount + + packKey(localX / blockSize, localY / blockSize, localZ / blockSize, spanY, spanZ); + else + clusterKeys[i] = packKey(localX, localY, localZ, spanY, spanZ); +} + +/// @brief Accumulate the per-cluster QEF: each surface cell builds the same 12-edge normal-equation +/// terms as placeQefVerticesKernel (A = sum n n^T, b = sum n (n.p), and the crossing centroid), but +/// expressed in the cluster's shared origin frame and atomic-added into its cluster's running +/// totals. clusterSolveKernel then solves one vertex from each cluster's summed system. (double +/// atomicAdd -> order-dependent; see the decimation note above.) +__global__ void +clusterQefAccumKernel(const int32_t *surfaceCells, + const int32_t *clusterIds, + int64_t numSurfaceCells, + const int32_t *neighborTable, + const float *sdf, + const float *gradient, + const int32_t *voxelCoord, + BoxStencilGeometry geometry, + float iso, + int originX, + int originY, + int originZ, + double *normalMatrix, + double *normalRhs, + double *crossingSum, + double *crossingCount, + double *normalSum) { + int64_t i = int64_t(blockIdx.x) * blockDim.x + threadIdx.x; + if (i >= numSurfaceCells) + return; + int32_t cellValueIndex = surfaceCells[i]; + int clusterId = clusterIds[i]; + float cornerSdf[8], cornerGradient[8][3]; + loadCellCorners( + cellValueIndex, neighborTable, sdf, gradient, geometry, cornerSdf, cornerGradient); + double cellOriginX = voxelCoord[int64_t(cellValueIndex) * 3 + 0] - originX, + cellOriginY = voxelCoord[int64_t(cellValueIndex) * 3 + 1] - originY, + cellOriginZ = voxelCoord[int64_t(cellValueIndex) * 3 + 2] - originZ; + double localMatrix[6] = {0, 0, 0, 0, 0, 0}, localRhs[3] = {0, 0, 0}, + localCrossingSum[3] = {0, 0, 0}, localNormalSum[3] = {0, 0, 0}; +#pragma unroll + for (int edge = 0; edge < 12; ++edge) { + int cornerA = geometry.edgeCornerA[edge], cornerB = geometry.edgeCornerB[edge]; + float sdfA = cornerSdf[cornerA], sdfB = cornerSdf[cornerB]; + if ((sdfA < iso) == (sdfB < iso)) + continue; + float sdfDelta = sdfB - sdfA; + float crossingT = (fabsf(sdfDelta) > 1e-12f) ? (iso - sdfA) / sdfDelta : 0.5f; + crossingT = fminf(fmaxf(crossingT, 0.f), 1.f); + float crossingX = + geometry.cornerOffset[cornerA][0] + + crossingT * (geometry.cornerOffset[cornerB][0] - geometry.cornerOffset[cornerA][0]); + float crossingY = + geometry.cornerOffset[cornerA][1] + + crossingT * (geometry.cornerOffset[cornerB][1] - geometry.cornerOffset[cornerA][1]); + float crossingZ = + geometry.cornerOffset[cornerA][2] + + crossingT * (geometry.cornerOffset[cornerB][2] - geometry.cornerOffset[cornerA][2]); + float nx = cornerGradient[cornerA][0] + + crossingT * (cornerGradient[cornerB][0] - cornerGradient[cornerA][0]), + ny = cornerGradient[cornerA][1] + + crossingT * (cornerGradient[cornerB][1] - cornerGradient[cornerA][1]), + nz = cornerGradient[cornerA][2] + + crossingT * (cornerGradient[cornerB][2] - cornerGradient[cornerA][2]); + float invLen = rsqrtf(nx * nx + ny * ny + nz * nz + 1e-24f); + nx *= invLen; + ny *= invLen; + nz *= invLen; + double px = cellOriginX + crossingX, py = cellOriginY + crossingY, + pz = cellOriginZ + crossingZ, normalDotCrossing = nx * px + ny * py + nz * pz; + localMatrix[0] += nx * nx; + localMatrix[1] += nx * ny; + localMatrix[2] += nx * nz; + localMatrix[3] += ny * ny; + localMatrix[4] += ny * nz; + localMatrix[5] += nz * nz; + localRhs[0] += nx * normalDotCrossing; + localRhs[1] += ny * normalDotCrossing; + localRhs[2] += nz * normalDotCrossing; + localCrossingSum[0] += px; + localCrossingSum[1] += py; + localCrossingSum[2] += pz; + localNormalSum[0] += nx; + localNormalSum[1] += ny; + localNormalSum[2] += nz; + } + for (int j = 0; j < 6; ++j) + atomicAdd(&normalMatrix[clusterId * 6 + j], localMatrix[j]); + for (int j = 0; j < 3; ++j) { + atomicAdd(&normalRhs[clusterId * 3 + j], localRhs[j]); + atomicAdd(&crossingSum[clusterId * 3 + j], localCrossingSum[j]); + atomicAdd(&normalSum[clusterId * 3 + j], localNormalSum[j]); + } + int numCrossings = 0; +#pragma unroll + for (int edge = 0; edge < 12; ++edge) { + int cornerA = geometry.edgeCornerA[edge], cornerB = geometry.edgeCornerB[edge]; + if ((cornerSdf[cornerA] < iso) != (cornerSdf[cornerB] < iso)) + ++numCrossings; + } + atomicAdd(&crossingCount[clusterId], (double)numCrossings); +} + +/// @brief Solve one vertex per cluster from its accumulated QEF -- same centroid re-centring + 0.05 +/// Tikhonov regularisation as placeQefVerticesKernel -- clamp it within +/- blockSize of the +/// cluster centroid, then write the world-space vertex and its reference normal (the summed cluster +/// normal). +__global__ void +clusterSolveKernel(int64_t numClusters, + const double *normalMatrix, + const double *normalRhs, + const double *crossingSum, + const double *crossingCount, + const double *normalSum, + int originX, + int originY, + int originZ, + int blockSize, + VoxelCoordTransform transform, + float *vertices, + float *refNormal) { + int64_t clusterId = int64_t(blockIdx.x) * blockDim.x + threadIdx.x; + if (clusterId >= numClusters) + return; + double count = crossingCount[clusterId], invCount = count > 0 ? 1.0 / count : 0.0; + double crossingCentroid[3] = {crossingSum[clusterId * 3 + 0] * invCount, + crossingSum[clusterId * 3 + 1] * invCount, + crossingSum[clusterId * 3 + 2] * invCount}; + double regularizedMatrix[6] = {normalMatrix[clusterId * 6 + 0] + 0.05, + normalMatrix[clusterId * 6 + 1], + normalMatrix[clusterId * 6 + 2], + normalMatrix[clusterId * 6 + 3] + 0.05, + normalMatrix[clusterId * 6 + 4], + normalMatrix[clusterId * 6 + 5] + 0.05}; + double rhs[3] = {normalRhs[clusterId * 3 + 0] + 0.05 * crossingCentroid[0], + normalRhs[clusterId * 3 + 1] + 0.05 * crossingCentroid[1], + normalRhs[clusterId * 3 + 2] + 0.05 * crossingCentroid[2]}; + double solution[3]; + solveSymmetric3x3(regularizedMatrix, rhs, solution); + for (int axis = 0; axis < 3; ++axis) + solution[axis] = fmin(fmax(solution[axis], crossingCentroid[axis] - blockSize), + crossingCentroid[axis] + blockSize); + auto worldPos = transform.applyInv((float)(solution[0] + originX), + (float)(solution[1] + originY), + (float)(solution[2] + originZ)); + vertices[clusterId * 3 + 0] = worldPos[0]; + vertices[clusterId * 3 + 1] = worldPos[1]; + vertices[clusterId * 3 + 2] = worldPos[2]; + refNormal[clusterId * 3 + 0] = (float)normalSum[clusterId * 3 + 0]; + refNormal[clusterId * 3 + 1] = (float)normalSum[clusterId * 3 + 1]; + refNormal[clusterId * 3 + 2] = (float)normalSum[clusterId * 3 + 2]; +} + +// ------------------------- host helpers ------------------------- +/// @brief Dense [0,numUnique) relabel of `keys` in original order: expressed via torch::unique_dim +/// exactly as MarchingCubes does. Returns the per-element int32 labels (sorted-key order) and +/// writes the unique count to `numUnique`. +static torch::Tensor +denseRelabel(const torch::Tensor &keys, int64_t &numUnique) { + auto result = torch::unique_dim(keys, 0, /*sorted=*/true, /*return_inverse=*/true); + numUnique = std::get<0>(result).size(0); + return std::get<1>(result).to(torch::kInt32); +} + +/// @brief Per-dtype TensorOptions for this op's scratch/output tensors (all on the same device). +/// Bundled so the per-grid mesher takes a single argument instead of one option per dtype. +struct ScratchOptions { + torch::TensorOptions f32, i32, i64, u8, f64; + explicit ScratchOptions(torch::Device device) + : f32(torch::TensorOptions().dtype(torch::kFloat32).device(device)), + i32(torch::TensorOptions().dtype(torch::kInt32).device(device)), + i64(torch::TensorOptions().dtype(torch::kInt64).device(device)), + u8(torch::TensorOptions().dtype(torch::kUInt8).device(device)), + f64(torch::TensorOptions().dtype(torch::kFloat64).device(device)) {} +}; + +/// @brief Mesh ONE grid into a triangle mesh. `sdf` is value-indexed (length valueCount, slot 0 = +/// the "outside" sentinel used for inactive neighbours). Pipeline: +/// 1. gatherFusedKernel -- one VBM decode per voxel -> neighbour table + gradient + surface flag +/// + index coord +/// 2. compact the surface cells (torch::nonzero) -- the cells the iso surface passes through +/// 3. place one vertex per surface cell: QEF (placeQefVerticesKernel), or when decimating group +/// the +/// cells into clusters and solve one vertex per cluster +/// 4. connectivity: one quad per sign-changing minimal grid edge -> triangles (count -> cumsum -> +/// write), each oriented outward by the SDF gradient +/// 5. prune unreferenced vertices and normalise the per-vertex reference normals +/// Returns {vertices (V,3) f32, normals (V,3) f32, triangles (T,3) int32} with triangle indices +/// local to this grid. +static std::tuple +meshOneGrid(OnIndexGridT *grid, + const VBMHelper &vbm, + const float *sdf, + int64_t valueCount, + VoxelCoordTransform transform, + float iso, + int reduce, + double adaptivity, + cudaStream_t stream, + const ScratchOptions &opts) { + const bool decimate = (reduce != 1) || (adaptivity > 0.0); + const int numColumns = kGeometry.numColumns; + auto emptyVertices = [&] { return torch::empty({0, 3}, opts.f32); }; + auto emptyTriangles = [&] { return torch::empty({0, 3}, opts.i32); }; + + if (valueCount <= 1 || vbm.blockCount == 0) + return {emptyVertices(), emptyVertices(), emptyTriangles()}; + + torch::Tensor neighborTableBuf = torch::empty({valueCount * numColumns}, opts.i32); + torch::Tensor gradientBuf = torch::empty({valueCount * 3}, opts.f32); + torch::Tensor surfaceFlagBuf = torch::empty({valueCount}, opts.u8); + torch::Tensor voxelCoordBuf = torch::empty({valueCount * 3}, opts.i32); + int32_t *neighborTable = neighborTableBuf.data_ptr(); + float *gradient = gradientBuf.data_ptr(); + uint8_t *surfaceFlag = surfaceFlagBuf.data_ptr(); + int32_t *voxelCoord = voxelCoordBuf.data_ptr(); + C10_CUDA_CHECK(cudaMemsetAsync(gradient, 0, 3 * sizeof(float), stream)); // background slot 0 + C10_CUDA_CHECK(cudaMemsetAsync(surfaceFlag, 0, sizeof(uint8_t), stream)); + + gatherFusedKernel<<>>(grid, + vbm.firstLeafID(), + vbm.jumpMap(), + vbm.firstOffset, + sdf, + kGeometry, + iso, + neighborTable, + gradient, + surfaceFlag, + voxelCoord); + C10_CUDA_KERNEL_LAUNCH_CHECK(); + + // compact surface cells: torch::nonzero returns the set value indices in ascending order, + // giving a stable vertex order (== the source's copy_if over a counting iterator). + torch::Tensor surfaceCellBuf = torch::nonzero(surfaceFlagBuf).flatten().to(torch::kInt32); + const int32_t *surfaceCells = surfaceCellBuf.data_ptr(); + const int64_t numSurfaceCells = surfaceCellBuf.size(0); + if (numSurfaceCells == 0) + return {emptyVertices(), emptyVertices(), emptyTriangles()}; + // int64 view of the surface-cell indices, reused by the index_select / index_put_ scatters + // below + torch::Tensor surfaceCellsLong = surfaceCellBuf.to(torch::kInt64); + + torch::Tensor cellVertexBuf = torch::full({valueCount}, -1, opts.i32); + int32_t *cellVertex = cellVertexBuf.data_ptr(); + + int64_t numVertices = 0; + torch::Tensor vertexBuf, refNormalBuf; + if (!decimate) { + numVertices = numSurfaceCells; + // cellVertex[surfaceCells[i]] = i (the compaction ordinal); background slots stay -1. + cellVertexBuf.index_put_({surfaceCellsLong}, torch::arange(numSurfaceCells, opts.i32)); + vertexBuf = torch::empty({numSurfaceCells * 3}, opts.f32); + refNormalBuf = torch::empty({numSurfaceCells * 3}, opts.f32); + placeQefVerticesKernel<<>>(surfaceCells, + numSurfaceCells, + neighborTable, + sdf, + gradient, + voxelCoord, + kGeometry, + iso, + transform, + vertexBuf.data_ptr(), + refNormalBuf.data_ptr()); + C10_CUDA_KERNEL_LAUNCH_CHECK(); + } else { + const int blockSize = (adaptivity > 0.0 && reduce <= 1) ? 8 : std::max(1, reduce); + // surface-cell coord bounds (for the mixed-radix cluster keys): a gather + amin/amax + // reduction in place of a hand-rolled atomic-min/max kernel. + torch::Tensor surfaceCoord = + voxelCoordBuf.view({valueCount, 3}).index_select(0, surfaceCellsLong); + torch::Tensor coordBounds = torch::stack({surfaceCoord.amin(0), surfaceCoord.amax(0)}, 0) + .cpu(); // (2,3): row 0 min, row 1 max + auto bounds = coordBounds.accessor(); + const int originX = bounds[0][0], originY = bounds[0][1], originZ = bounds[0][2]; + const int spanY = bounds[1][1] - originY + 1, spanZ = bounds[1][2] - originZ + 1; + const long long fineKeyCount = (long long)(bounds[1][0] - originX + 1) * spanY * spanZ; + + torch::Tensor clusterKeyBuf = torch::empty({numSurfaceCells}, opts.i64); + auto *clusterKeys = reinterpret_cast(clusterKeyBuf.data_ptr()); + torch::Tensor clusterIdBuf; + int64_t numClusters = 0; + if (adaptivity <= 0.0) { // uniform blocks + clusterKeyKernel<<>>(surfaceCells, + numSurfaceCells, + voxelCoord, + originX, + originY, + originZ, + spanY, + spanZ, + blockSize, + clusterKeys); + C10_CUDA_KERNEL_LAUNCH_CHECK(); + clusterIdBuf = denseRelabel(clusterKeyBuf, numClusters); + } else { // flat blocks collapse, feature cells stay fine + torch::Tensor cellNormalBuf = torch::empty({numSurfaceCells * 3}, opts.f32); + cellNormalKernel<<>>(surfaceCells, + numSurfaceCells, + neighborTable, + sdf, + gradient, + kGeometry, + iso, + cellNormalBuf.data_ptr()); + C10_CUDA_KERNEL_LAUNCH_CHECK(); + clusterKeyKernel<<>>(surfaceCells, + numSurfaceCells, + voxelCoord, + originX, + originY, + originZ, + spanY, + spanZ, + blockSize, + clusterKeys); + C10_CUDA_KERNEL_LAUNCH_CHECK(); + int64_t numBlocks = 0; + torch::Tensor blockOfCellBuf = denseRelabel(clusterKeyBuf, numBlocks); + // a block is "flat" when the mean magnitude of its cells' unit normals (i.e. how + // aligned they are) clears cos(threshold); flat blocks then collapse to one vertex. The + // scatter-add and per-block reduction are index_add_ + bincount + a vectorized + // comparison. + torch::Tensor blockNormalSum = + torch::zeros({numBlocks, 3}, opts.f32) + .index_add_(0, blockOfCellBuf, cellNormalBuf.view({-1, 3})); + torch::Tensor blockCellCount = + torch::bincount(blockOfCellBuf, /*weights=*/{}, numBlocks) + .clamp_min(1) + .to(torch::kFloat32); + const float flatCosThreshold = + (float)std::cos(adaptivity * 60.0 * 3.14159265358979323846 / 180.0); + torch::Tensor blockFlatBuf = + (blockNormalSum.norm(2, /*dim=*/1) / blockCellCount >= flatCosThreshold) + .to(torch::kUInt8); + clusterKeyAdaptiveKernel<<>>(surfaceCells, + numSurfaceCells, + voxelCoord, + originX, + originY, + originZ, + spanY, + spanZ, + blockSize, + blockOfCellBuf.data_ptr(), + blockFlatBuf.data_ptr(), + fineKeyCount, + clusterKeys); + C10_CUDA_KERNEL_LAUNCH_CHECK(); + clusterIdBuf = denseRelabel(clusterKeyBuf, numClusters); + } + int32_t *clusterIds = clusterIdBuf.data_ptr(); + numVertices = numClusters; + // cellVertex[surfaceCells[i]] = clusterIds[i] (background slots stay -1). + cellVertexBuf.index_put_({surfaceCellsLong}, clusterIdBuf); + if (numClusters == 0) + return {emptyVertices(), emptyVertices(), emptyTriangles()}; + torch::Tensor normalMatrixBuf = torch::zeros({numClusters * 6}, opts.f64), + normalRhsBuf = torch::zeros({numClusters * 3}, opts.f64), + crossingSumBuf = torch::zeros({numClusters * 3}, opts.f64), + crossingCountBuf = torch::zeros({numClusters}, opts.f64), + normalSumBuf = torch::zeros({numClusters * 3}, opts.f64); + clusterQefAccumKernel<<>>(surfaceCells, + clusterIds, + numSurfaceCells, + neighborTable, + sdf, + gradient, + voxelCoord, + kGeometry, + iso, + originX, + originY, + originZ, + normalMatrixBuf.data_ptr(), + normalRhsBuf.data_ptr(), + crossingSumBuf.data_ptr(), + crossingCountBuf.data_ptr(), + normalSumBuf.data_ptr()); + C10_CUDA_KERNEL_LAUNCH_CHECK(); + vertexBuf = torch::empty({numClusters * 3}, opts.f32); + refNormalBuf = torch::empty({numClusters * 3}, opts.f32); + clusterSolveKernel<<>>(numClusters, + normalMatrixBuf.data_ptr(), + normalRhsBuf.data_ptr(), + crossingSumBuf.data_ptr(), + crossingCountBuf.data_ptr(), + normalSumBuf.data_ptr(), + originX, + originY, + originZ, + blockSize, + transform, + vertexBuf.data_ptr(), + refNormalBuf.data_ptr()); + C10_CUDA_KERNEL_LAUNCH_CHECK(); + } + + // connectivity: per-surface-cell triangle counts -> cumsum into contiguous write ranges -> + // write. This is the MarchingCubes pattern; the deterministic offsets replace a shared atomic + // append counter, and iterating surface cells (not every active voxel) keeps the launch tight. + torch::Tensor triangleCountBuf = torch::empty({numSurfaceCells}, opts.i64); + countTrianglesKernel<<>>(surfaceCells, + numSurfaceCells, + neighborTable, + sdf, + cellVertex, + kGeometry, + iso, + triangleCountBuf.data_ptr()); + C10_CUDA_KERNEL_LAUNCH_CHECK(); + torch::Tensor triangleOffsetBuf = torch::cumsum(triangleCountBuf, 0); // inclusive prefix sum + int64_t numTriangles = triangleOffsetBuf[-1].item(); + triangleOffsetBuf = torch::roll(triangleOffsetBuf, {1}); // shift -> exclusive prefix sum + triangleOffsetBuf[0] = 0; + torch::Tensor triangleBuf = torch::empty({numTriangles * 3}, opts.i32); + int32_t *triangles = triangleBuf.data_ptr(); + if (numTriangles > 0) { + writeTrianglesKernel<<>>(surfaceCells, + numSurfaceCells, + neighborTable, + sdf, + cellVertex, + kGeometry, + iso, + triangleOffsetBuf.data_ptr(), + triangles); + C10_CUDA_KERNEL_LAUNCH_CHECK(); + orientTrianglesKernel<<>>( + numTriangles, vertexBuf.data_ptr(), refNormalBuf.data_ptr(), triangles); + C10_CUDA_KERNEL_LAUNCH_CHECK(); + } + + // prune unreferenced (free-point) vertices and reindex the triangles. torch::unique_dim over + // the flattened triangle indices returns, in ascending order, the referenced vertex ids and + // (via return_inverse) the triangles already reindexed into [0, numKeptVertices) -- the same + // merge idiom MarchingCubes uses. Compacting the vertices/normals is then a single index_select + // gather. + torch::Tensor outVertices, outNormals; + int64_t numOutVertices = 0; + if (numTriangles > 0 && numVertices > 0) { + auto unique = torch::unique_dim(triangleBuf, 0, /*sorted=*/true, /*return_inverse=*/true); + torch::Tensor keptVertexIds = + std::get<0>(unique).to(torch::kInt64); // referenced ids, ascending + numOutVertices = keptVertexIds.size(0); + outVertices = vertexBuf.view({numVertices, 3}).index_select(0, keptVertexIds); + outNormals = refNormalBuf.view({numVertices, 3}).index_select(0, keptVertexIds); + triangleBuf = std::get<1>(unique).to(torch::kInt32); + } else { + outVertices = emptyVertices(); + outNormals = emptyVertices(); + numTriangles = 0; + } + + // normalize the per-vertex reference normal -> unit SDF-gradient normal + if (numOutVertices > 0) + outNormals = outNormals / outNormals.norm(2, /*dim=*/1, /*keepdim=*/true).clamp_min(1e-12); + return {outVertices.view({numOutVertices, 3}), + outNormals.view({numOutVertices, 3}), + triangleBuf.view({numTriangles, 3})}; +} + +} // namespace + +std::tuple +dualContour(const GridBatchData &batchHdl, + const JaggedTensor &field, + double iso, + int reduce, + double adaptivity) { + TORCH_CHECK_VALUE( + field.ldim() == 1, + "Expected field to have 1 list dimension (a single list of per-voxel values)"); + TORCH_CHECK_TYPE(field.is_floating_point(), "field must have a floating point type"); + TORCH_CHECK_VALUE(field.numel() == batchHdl.totalVoxels(), + "field value count does not match the number of voxels in the grid"); + TORCH_CHECK_VALUE(field.num_outer_lists() == batchHdl.batchSize(), + "field batch size does not match the grid batch size"); + batchHdl.checkDevice(field); + TORCH_CHECK(field.device().is_cuda(), "dual_contour currently requires a CUDA device"); + adaptivity = std::clamp(adaptivity, 0.0, 1.5); // beyond 1.5 the flatness bar saturates + reduce = std::max(1, reduce); + + cudaStream_t stream = at::cuda::getCurrentCUDAStream(batchHdl.device().index()).stream(); + const auto device = field.device(); + ScratchOptions opts(device); + auto jidxOpts = torch::TensorOptions().dtype(fvdb::JIdxScalarType).device(device); + + torch::Tensor fieldData = field.jdata().contiguous(); + if (fieldData.dim() != 1) + fieldData = fieldData.view({-1}); + if (fieldData.scalar_type() != torch::kFloat32) + fieldData = + fieldData.to(torch::kFloat32); // QEF accumulators are double; the SDF stays float32 + const float *fieldPtr = fieldData.data_ptr(); + const float isoFloat = (float)iso; + + std::vector vertexList, normalList, faceList, vertexBatchList, faceBatchList; + for (int64_t b = 0; b < batchHdl.batchSize(); ++b) { + const int64_t numVoxels = batchHdl.numVoxelsAt(b); + torch::Tensor vertices, normals, triangles; + if (numVoxels == 0) { + vertices = torch::empty({0, 3}, opts.f32); + normals = torch::empty({0, 3}, opts.f32); + triangles = torch::empty({0, 3}, opts.i32); + } else { + OnIndexGridT *grid = batchHdl.mGridHdl->deviceGrid((uint32_t)b); + const int64_t voxelOffset = batchHdl.cumVoxelsAt(b); + VoxelCoordTransform transform = batchHdl.primalTransformAt(b); + VBMHelper vbm(grid, stream); + const int64_t valueCount = (int64_t)vbm.valueCount; // numVoxels + 1 + + // gather the field into a value-indexed SDF buffer: slot 0 = "outside" sentinel (>= + // iso) so inactive corners classify as outside, and slots [1..numVoxels] = this grid's + // field. The memcpy fills everything except slot 0, so only slot 0 needs the sentinel + // write. + torch::Tensor sdfBuf = torch::empty({valueCount}, opts.f32); + float *sdf = sdfBuf.data_ptr(); + sdfBuf.narrow(0, 0, 1).fill_(isoFloat + 1.0f); + C10_CUDA_CHECK(cudaMemcpyAsync(sdf + 1, + fieldPtr + voxelOffset, + numVoxels * sizeof(float), + cudaMemcpyDeviceToDevice, + stream)); + std::tie(vertices, normals, triangles) = meshOneGrid( + grid, vbm, sdf, valueCount, transform, isoFloat, reduce, adaptivity, stream, opts); + C10_CUDA_CHECK(cudaStreamSynchronize(stream)); // per-grid scratch freed at scope exit + } + const int64_t numGridVertices = vertices.size(0), numGridTriangles = triangles.size(0); + vertexList.push_back(vertices); + normalList.push_back(normals); + faceList.push_back(triangles.to(torch::kInt64)); + vertexBatchList.push_back(torch::full({numGridVertices}, b, jidxOpts)); + faceBatchList.push_back(torch::full({numGridTriangles}, b, jidxOpts)); + } + + // concatenate per-grid meshes (already grouped by batch) and build the jagged outputs + torch::Tensor allVertices = torch::cat(vertexList, 0); + torch::Tensor allNormals = torch::cat(normalList, 0); + torch::Tensor allFaces = torch::cat(faceList, 0); + torch::Tensor vertexBatchIdx = torch::cat(vertexBatchList, 0); + torch::Tensor faceBatchIdx = torch::cat(faceBatchList, 0); + + JaggedTensor retVertices = JaggedTensor::from_data_indices_and_list_ids( + allVertices, vertexBatchIdx, batchHdl.jlidx(), batchHdl.batchSize()); + JaggedTensor retFaces = JaggedTensor::from_data_indices_and_list_ids( + allFaces, faceBatchIdx, batchHdl.jlidx(), batchHdl.batchSize()); + JaggedTensor retNormals = JaggedTensor::from_data_indices_and_list_ids( + allNormals, vertexBatchIdx, batchHdl.jlidx(), batchHdl.batchSize()); + return {retVertices, retFaces, retNormals}; +} + +} // namespace ops +} // namespace detail +} // namespace fvdb diff --git a/src/fvdb/detail/ops/DualContour.h b/src/fvdb/detail/ops/DualContour.h new file mode 100644 index 000000000..af4e31eea --- /dev/null +++ b/src/fvdb/detail/ops/DualContour.h @@ -0,0 +1,49 @@ +// Copyright Contributors to the OpenVDB Project +// SPDX-License-Identifier: Apache-2.0 +// +#ifndef FVDB_DETAIL_OPS_DUALCONTOUR_H +#define FVDB_DETAIL_OPS_DUALCONTOUR_H + +#include +#include + +#include + +#include + +namespace fvdb { +namespace detail { +namespace ops { + +/// @brief Dual-contour (DC + QEF) mesh an OnIndex grid carrying a narrow-band signed distance +/// field. +/// +/// For each cube cell whose 8 corners are all active and whose corner signs straddle @p iso, places +/// one vertex by minimising a quadratic error function built from the 12 edge zero-crossings and +/// their interpolated SDF-gradient normals (regularised toward the crossing centroid, clamped to +/// the cell). One quad per bipolar minimal grid edge joins the 4 surrounding cells; quads are +/// triangulated and oriented outward by the SDF gradient. Optional cluster-collapse decimation +/// (@p reduce / @p adaptivity) contracts the connectivity onto fewer vertices. +/// Requires a >= ~3-voxel band for a watertight result. +/// +/// @param batchHdl Grid batch defining the sparse topology. +/// @param field Per-voxel signed field: a floating-point JaggedTensor of shape [B, -1]. +/// @param iso Isovalue of the surface to extract. +/// @param reduce Uniform F x F x F cluster-collapse decimation factor (1 = full resolution). +/// @param adaptivity Curvature-adaptive decimation in [0, 1.5] (0 = uniform/off): collapse flat +/// blocks while keeping full detail at features. +/// @return A (vertices, faces, normals) tuple, each jagged over the grid batch: vertices and +/// normals are float32 JaggedTensors of shape [B, -1, 3] and faces is an int64 JaggedTensor +/// of shape [B, -1, 3]. faces holds grid-local triangle vertex indices; normals is the +/// normalized SDF gradient at each vertex. +std::tuple dualContour(const GridBatchData &batchHdl, + const JaggedTensor &field, + double iso, + int reduce, + double adaptivity); + +} // namespace ops +} // namespace detail +} // namespace fvdb + +#endif // FVDB_DETAIL_OPS_DUALCONTOUR_H diff --git a/src/fvdb/detail/ops/ReinitializeSdf.cu b/src/fvdb/detail/ops/ReinitializeSdf.cu index 2780efecd..b124fd9b9 100644 --- a/src/fvdb/detail/ops/ReinitializeSdf.cu +++ b/src/fvdb/detail/ops/ReinitializeSdf.cu @@ -3,6 +3,7 @@ // #include #include +#include #include #include @@ -23,43 +24,6 @@ namespace ops { namespace { -using OnIndexGridT = nanovdb::NanoGrid; -using VbmBuffer = nanovdb::cuda::DeviceBuffer; - -// log2 of the VoxelBlockManager block width: each VBM block spans 2^9 = 512 active voxels. -static constexpr int kLog2BlockWidth = 9; - -// ------------------------- VBM fused 6-face stencil preamble ------------------------- -// The VBM decode gives the centre coord / value-index for free; we then read just the 6 FACE -// neighbours through a cached ReadAccessor. Yields `centerIndex` (the centre voxel's value index) -// and `faceIndex[6]` (the 6 face-neighbour value indices, -x,+x,-y,+y,-z,+z; 0 = -// inactive/background). Kernels using it take the grid/firstLeafID/jumpMap/firstOffset parameters -// by these exact names. -#define VBM_FACES_BEGIN() \ - constexpr int blockWidth = 1 << kLog2BlockWidth, jumpMapWordCount = blockWidth / 64; \ - using VoxelBlockManagerT = nanovdb::tools::cuda::VoxelBlockManager; \ - __shared__ uint32_t sharedLeafIndex[blockWidth]; \ - __shared__ uint16_t sharedVoxelOffset[blockWidth]; \ - VoxelBlockManagerT::template decodeInverseMaps( \ - grid, \ - firstLeafID[blockIdx.x], \ - &jumpMap[uint64_t(blockIdx.x) * jumpMapWordCount], \ - firstOffset + uint64_t(blockIdx.x) * blockWidth, \ - sharedLeafIndex, \ - sharedVoxelOffset); \ - if (sharedLeafIndex[threadIdx.x] == VoxelBlockManagerT::UnusedLeafIndex) \ - return; \ - const auto &leaf = grid->tree().template getFirstNode<0>()[sharedLeafIndex[threadIdx.x]]; \ - const nanovdb::Coord centerCoord = leaf.offsetToGlobalCoord(sharedVoxelOffset[threadIdx.x]); \ - const uint64_t centerIndex = leaf.getValue(sharedVoxelOffset[threadIdx.x]); \ - auto accessor = grid->getAccessor(); \ - const uint64_t faceIndex[6] = {accessor.getValue(centerCoord.offsetBy(-1, 0, 0)), \ - accessor.getValue(centerCoord.offsetBy(1, 0, 0)), \ - accessor.getValue(centerCoord.offsetBy(0, -1, 0)), \ - accessor.getValue(centerCoord.offsetBy(0, 1, 0)), \ - accessor.getValue(centerCoord.offsetBy(0, 0, -1)), \ - accessor.getValue(centerCoord.offsetBy(0, 0, 1))}; - // ===================== fused stencil kernels ==================================================== // frozen Peng smoothed sign from a field's value + central-difference gradient. template @@ -71,7 +35,11 @@ signFusedKernel(const OnIndexGridT *grid, const ScalarT *field, ScalarT voxelSize, ScalarT *sign) { - VBM_FACES_BEGIN(); + VbmFaceStencil faces; + if (!vbmDecodeFaceStencil(grid, firstLeafID, jumpMap, firstOffset, faces)) + return; + const uint64_t centerIndex = faces.centerIndex; + const uint64_t *faceIndex = faces.faceIndex; ScalarT xm = field[faceIndex[0]], xp = field[faceIndex[1]], ym = field[faceIndex[2]], yp = field[faceIndex[3]], zm = field[faceIndex[4]], zp = field[faceIndex[5]]; ScalarT gradX = (xp - xm) / (2 * voxelSize), gradY = (yp - ym) / (2 * voxelSize), @@ -109,7 +77,11 @@ godunovFusedKernel(const OnIndexGridT *grid, const ScalarT *sign, ScalarT voxelSize, ScalarT *rhs) { - VBM_FACES_BEGIN(); + VbmFaceStencil faces; + if (!vbmDecodeFaceStencil(grid, firstLeafID, jumpMap, firstOffset, faces)) + return; + const uint64_t centerIndex = faces.centerIndex; + const uint64_t *faceIndex = faces.faceIndex; ScalarT center = field[centerIndex], sgn = sign[centerIndex]; ScalarT xm = field[faceIndex[0]], xp = field[faceIndex[1]], ym = field[faceIndex[2]], yp = field[faceIndex[3]], zm = field[faceIndex[4]], zp = field[faceIndex[5]]; @@ -131,8 +103,12 @@ smoothFusedKernel(const OnIndexGridT *grid, const ScalarT *in, ScalarT weight, ScalarT *out) { - VBM_FACES_BEGIN(); - ScalarT center = in[centerIndex]; + VbmFaceStencil faces; + if (!vbmDecodeFaceStencil(grid, firstLeafID, jumpMap, firstOffset, faces)) + return; + const uint64_t centerIndex = faces.centerIndex; + const uint64_t *faceIndex = faces.faceIndex; + ScalarT center = in[centerIndex]; ScalarT faceMean = (in[faceIndex[0]] + in[faceIndex[1]] + in[faceIndex[2]] + in[faceIndex[3]] + in[faceIndex[4]] + in[faceIndex[5]]) * (ScalarT(1) / ScalarT(6)); @@ -189,28 +165,6 @@ heunKernel(ScalarT *out, out[i] = nanovdb::math::Min(nanovdb::math::Max(value, -bandWidth), bandWidth); } -// small VBM helper: build once, expose the block count + the firstLeafID/jumpMap device pointers. -struct VBMHelper { - nanovdb::tools::VoxelBlockManagerHandle handle; - uint32_t blockCount{0}; - uint64_t firstOffset{0}, valueCount{1}; - VBMHelper(OnIndexGridT *grid, cudaStream_t stream) { - handle = nanovdb::tools::cuda::buildVoxelBlockManager( - grid, 0, 0, 0, stream); - blockCount = (uint32_t)handle.blockCount(); - firstOffset = handle.firstOffset(); - valueCount = handle.lastOffset() + 1; - } - const uint32_t * - firstLeafID() const { - return handle.deviceFirstLeafID(); - } - const uint64_t * - jumpMap() const { - return handle.deviceJumpMap(); - } -}; - // Redistance (|grad phi| = 1) + optional de-staircase one grid's value-indexed buffer, in place. // `phi`/scratch are length `valueCount` with slot 0 holding the +bandWidth background; the // stencil/combiner kernels never write slot 0 (so inactive-neighbour reads always see the boundary diff --git a/src/fvdb/detail/ops/ReinitializeSdf.h b/src/fvdb/detail/ops/ReinitializeSdf.h index 4d363f747..721d68f3b 100644 --- a/src/fvdb/detail/ops/ReinitializeSdf.h +++ b/src/fvdb/detail/ops/ReinitializeSdf.h @@ -30,8 +30,7 @@ enum class SmoothingMode : int32_t { /// same per-voxel ordering as the input. /// /// @param batchHdl Grid batch defining the sparse topology. -/// @param field Per-voxel signed field (a single list of scalar values, numel == -/// totalVoxels). +/// @param field Per-voxel signed field: a floating-point JaggedTensor of shape [B, -1]. /// @param band Narrow-band half-width in voxels. The field is clamped to [-band*vx, band*vx] /// each sweep and the "outside" Dirichlet value for inactive neighbours is /// +band*vx. @@ -40,7 +39,8 @@ enum class SmoothingMode : int32_t { /// @param order TVD-RK order: 1 (forward Euler), 2 (Heun), or 3 (Shu-Osher). /// @param smooth Number of smoothing passes (0 disables smoothing). /// @param smoothing Which Laplacian flow each smoothing pass applies (mean-curvature or Taubin). -/// @return A new per-voxel SDF JaggedTensor on the same grid/ordering as @p field. +/// @return A per-voxel SDF: a JaggedTensor of shape [B, -1] with the same dtype, grid, and ordering +/// as @p field. JaggedTensor reinitializeSdf(const GridBatchData &batchHdl, const JaggedTensor &field, int band, diff --git a/src/fvdb/detail/utils/cuda/VoxelBlockManagerHelper.h b/src/fvdb/detail/utils/cuda/VoxelBlockManagerHelper.h new file mode 100644 index 000000000..dc2427adf --- /dev/null +++ b/src/fvdb/detail/utils/cuda/VoxelBlockManagerHelper.h @@ -0,0 +1,131 @@ +// Copyright Contributors to the OpenVDB Project +// SPDX-License-Identifier: Apache-2.0 +// +#ifndef FVDB_DETAIL_UTILS_CUDA_VOXELBLOCKMANAGERHELPER_H +#define FVDB_DETAIL_UTILS_CUDA_VOXELBLOCKMANAGERHELPER_H + +#include +#include +#include +#include + +#include + +#include + +namespace fvdb { +namespace detail { +namespace ops { + +// Shared NanoVDB VoxelBlockManager (VBM) scaffolding used by the VBM-based ops (ReinitializeSdf, +// DualContour, ...). Each op still writes its own per-voxel neighbour read -- they differ in which +// neighbours they touch (6 faces vs the full 3x3x3 box stencil) -- but the grid type, buffer type, +// block width, the build-once handle wrapper, and the per-block decode below are identical across +// all of them. + +using OnIndexGridT = nanovdb::NanoGrid; +using VbmBuffer = nanovdb::cuda::DeviceBuffer; + +// log2 of the VoxelBlockManager block width: each VBM block spans 2^9 = 512 active voxels. +static constexpr int kLog2BlockWidth = 9; + +// small VBM helper: build once, expose the block count + the firstLeafID/jumpMap device pointers. +struct VBMHelper { + nanovdb::tools::VoxelBlockManagerHandle handle; + uint32_t blockCount{0}; + uint64_t firstOffset{0}, valueCount{1}; + VBMHelper(OnIndexGridT *grid, cudaStream_t stream) { + handle = nanovdb::tools::cuda::buildVoxelBlockManager( + grid, 0, 0, 0, stream); + blockCount = (uint32_t)handle.blockCount(); + firstOffset = handle.firstOffset(); + valueCount = handle.lastOffset() + 1; + } + const uint32_t * + firstLeafID() const { + return handle.deviceFirstLeafID(); + } + const uint64_t * + jumpMap() const { + return handle.deviceJumpMap(); + } +}; + +// ------------------------- per-block decode helpers (device) ------------------------- +// Per-VBM-block shared decode scratch: the inverse leaf/offset maps for one 2^Log2BlockWidth block. +// Declare one as __shared__ in the kernel and hand it to vbmDecodeBlock. +template struct VbmBlockMaps { + uint32_t leafIndex[1 << Log2BlockWidth]; + uint16_t voxelOffset[1 << Log2BlockWidth]; +}; + +// Cooperatively decode block blockIdx.x into `maps` (ALL threads in the block must call). Returns +// true if this thread's slot maps to an active voxel; the caller returns on false. The kernel must +// have parameters named grid / firstLeafID / jumpMap / firstOffset (a VBMHelper supplies the latter +// three). After this, maps.leafIndex[threadIdx.x] / maps.voxelOffset[threadIdx.x] locate the voxel. +template +__device__ inline bool +vbmDecodeBlock(const OnIndexGridT *grid, + const uint32_t *firstLeafID, + const uint64_t *jumpMap, + uint64_t firstOffset, + VbmBlockMaps &maps) { + using VoxelBlockManagerT = nanovdb::tools::cuda::VoxelBlockManager; + constexpr uint64_t blockWidth = uint64_t(1) << Log2BlockWidth, + jumpMapWordCount = blockWidth / 64; + VoxelBlockManagerT::template decodeInverseMaps( + grid, + firstLeafID[blockIdx.x], + &jumpMap[blockIdx.x * jumpMapWordCount], + firstOffset + blockIdx.x * blockWidth, + maps.leafIndex, + maps.voxelOffset); + return maps.leafIndex[threadIdx.x] != VoxelBlockManagerT::UnusedLeafIndex; +} + +// The centre voxel's value-index plus its 6 face-neighbour value-indices (-x,+x,-y,+y,-z,+z; +// 0 = inactive/background). +struct VbmFaceStencil { + uint64_t centerIndex; + uint64_t faceIndex[6]; +}; + +// Read this thread's centre + 6-face value indices through a cached accessor. Call only for active +// threads (i.e. after vbmDecodeBlock returned true for this thread). +template +__device__ inline VbmFaceStencil +vbmReadFaceStencil(const OnIndexGridT *grid, const VbmBlockMaps &maps) { + const auto &leaf = grid->tree().template getFirstNode<0>()[maps.leafIndex[threadIdx.x]]; + const nanovdb::Coord centerCoord = leaf.offsetToGlobalCoord(maps.voxelOffset[threadIdx.x]); + auto accessor = grid->getAccessor(); + return {leaf.getValue(maps.voxelOffset[threadIdx.x]), + {accessor.getValue(centerCoord.offsetBy(-1, 0, 0)), + accessor.getValue(centerCoord.offsetBy(1, 0, 0)), + accessor.getValue(centerCoord.offsetBy(0, -1, 0)), + accessor.getValue(centerCoord.offsetBy(0, 1, 0)), + accessor.getValue(centerCoord.offsetBy(0, 0, -1)), + accessor.getValue(centerCoord.offsetBy(0, 0, 1))}}; +} + +// Full face-stencil kernel preamble: decode this block, and for active threads read the 6-face +// stencil into `out`. Returns false (caller should return) for unused slots. Owns its __shared__ +// scratch, so the kernel needs no shared declaration of its own. +template +__device__ inline bool +vbmDecodeFaceStencil(const OnIndexGridT *grid, + const uint32_t *firstLeafID, + const uint64_t *jumpMap, + uint64_t firstOffset, + VbmFaceStencil &out) { + __shared__ VbmBlockMaps maps; + if (!vbmDecodeBlock(grid, firstLeafID, jumpMap, firstOffset, maps)) + return false; + out = vbmReadFaceStencil(grid, maps); + return true; +} + +} // namespace ops +} // namespace detail +} // namespace fvdb + +#endif // FVDB_DETAIL_UTILS_CUDA_VOXELBLOCKMANAGERHELPER_H diff --git a/src/python/GridBatchOps.cpp b/src/python/GridBatchOps.cpp index 5fbcf821a..4fd5a832d 100644 --- a/src/python/GridBatchOps.cpp +++ b/src/python/GridBatchOps.cpp @@ -50,6 +50,7 @@ #include // Meshing / TSDF +#include #include #include #include @@ -459,6 +460,14 @@ bind_grid_batch_ops(py::module &m) { m.def( "marching_cubes", &ops::marchingCubes, py::arg("grid"), py::arg("field"), py::arg("level")); + m.def("dual_contour", + &ops::dualContour, + py::arg("grid"), + py::arg("field"), + py::arg("iso"), + py::arg("reduce"), + py::arg("adaptivity")); + m.def("reinitialize_sdf", &ops::reinitializeSdf, py::arg("grid"), diff --git a/tests/unit/test_dc.py b/tests/unit/test_dc.py new file mode 100644 index 000000000..456c76424 --- /dev/null +++ b/tests/unit/test_dc.py @@ -0,0 +1,195 @@ +# Copyright Contributors to the OpenVDB Project +# SPDX-License-Identifier: Apache-2.0 +# +import unittest + +import torch + +import fvdb + + +def _dense_cube_grid(vx: float, half: int, device: torch.device) -> "fvdb.Grid": + """A dense (2*half+1)^3 cube grid built by placing one point at each voxel centre.""" + rng = torch.arange(-half, half + 1, device=device, dtype=torch.float32) + ii, jj, kk = torch.meshgrid(rng, rng, rng, indexing="ij") + ijk = torch.stack([ii, jj, kk], dim=-1).reshape(-1, 3) + return fvdb.Grid.from_points(ijk * vx, voxel_size=vx) + + +def _edge_stats(faces: torch.Tensor, num_verts: int): + """(num_edges, boundary_edges, nonmanifold_edges) from a triangle list, on-device.""" + f = faces.long() + e = torch.cat([f[:, [0, 1]], f[:, [1, 2]], f[:, [2, 0]]], dim=0) + key = torch.minimum(e[:, 0], e[:, 1]) * num_verts + torch.maximum(e[:, 0], e[:, 1]) + counts = key.unique(return_counts=True)[1] + return int(counts.shape[0]), int((counts == 1).sum()), int((counts > 2).sum()) + + +class DualContourTests(unittest.TestCase): + def setUp(self): + """Build the shared fixture: a dense cube grid sampling an analytic sphere SDF, then a clean + narrow-band SDF on a pruned band grid (``mesh_grid`` / ``phi``) -- the kind of input dual + contouring is meant to consume.""" + if not torch.cuda.is_available(): + self.skipTest("dual_contour requires a CUDA device") + torch.manual_seed(0) + self.device = torch.device("cuda:0") + self.vx = 0.05 + self.R = 0.3 + self.band = 3 + self.bw = self.band * self.vx + self.grid = _dense_cube_grid(self.vx, half=12, device=self.device) + self.analytic = (self.grid.ijk.float() * self.vx).norm(dim=1) - self.R + # a clean narrow-band SDF on a pruned band grid (what DC is meant to consume) + self.mesh_grid, self.phi = self.grid.retopologize_sdf(self.analytic.clamp(-self.bw, self.bw), band=self.band) + + def test_dc_composes_with_reinitialize_and_retopologize(self): + """Dual contouring composes with the SDF ops that produce its input: meshing the output of + ``reinitialize_sdf`` (topology preserved) and of ``retopologize_sdf(pad=True)`` (re-banded / + pruned grid) both produce non-empty meshes.""" + # reinitialize_sdf keeps topology; the dense cube has a thick band -> DC meshes it + phi_re = self.grid.reinitialize_sdf(self.analytic.clamp(-self.bw, self.bw), band=self.band) + v, f, _ = self.grid.dual_contour(phi_re, iso=0.0) + self.assertGreater(v.shape[0], 0) + self.assertGreater(f.shape[0], 0) + # retopologize with padding then DC + g2, phi2 = self.grid.retopologize_sdf(self.analytic.clamp(-self.bw, self.bw), band=self.band, pad=True) + v2, f2, _ = g2.dual_contour(phi2, iso=0.0) + self.assertGreater(v2.shape[0], 0) + self.assertGreater(f2.shape[0], 0) + + def test_dc_extent_matches_marching_cubes(self): + """Dual contouring traces the same isosurface as marching cubes: the two meshes' bounding + boxes agree to within a voxel or two. A cross-implementation sanity check -- the precise + surface location is pinned analytically by the planar and sphere-radius tests.""" + v_dc, _, _ = self.mesh_grid.dual_contour(self.phi, iso=0.0) + v_mc, _, _ = self.mesh_grid.marching_cubes(self.phi, level=0.0) + # different algorithms, but both trace the same isosurface -> same bounding box (within a voxel) + self.assertTrue(torch.allclose(v_dc.amin(0), v_mc.amin(0), atol=2 * self.vx)) + self.assertTrue(torch.allclose(v_dc.amax(0), v_mc.amax(0), atol=2 * self.vx)) + + def test_dc_batch_matches_single(self): + """A batched dual_contour matches running each grid on its own: for every grid in the batch + the per-grid vertex/face counts and vertex positions match the single-grid result, and face + indices stay grid-local (max index < that grid's vertex count). Exercises the non-decimated + path, whose vertex order is deterministic (stable surface-cell compaction).""" + vx = self.vx + gb = fvdb.GridBatch.from_points( + fvdb.JaggedTensor([self._cube_points(vx, 12), self._cube_points(vx, 10)]), + voxel_sizes=vx, + ) + analytic = (gb.ijk.jdata.float() * vx).norm(dim=1) - self.R + gb2, phi_b = gb.retopologize_sdf(gb.jagged_like(analytic.clamp(-self.bw, self.bw)), band=self.band) + vb, fb, nb = gb2.dual_contour(phi_b, iso=0.0) + for i, half in enumerate([12, 10]): + single = fvdb.Grid.from_points(self._cube_points(vx, half), voxel_size=vx) + a = (single.ijk.float() * vx).norm(dim=1) - self.R + g2, phi_s = single.retopologize_sdf(a.clamp(-self.bw, self.bw), band=self.band) + vs, fs, _ = g2.dual_contour(phi_s, iso=0.0) + # vertex order is deterministic (stable surface-cell compaction); positions match + self.assertEqual(vb[i].jdata.shape[0], vs.shape[0]) + self.assertEqual(fb[i].jdata.shape[0], fs.shape[0]) + self.assertTrue(torch.allclose(vb[i].jdata, vs, atol=1e-4)) + self.assertLess(int(fb[i].jdata.max()), vb[i].jdata.shape[0]) # faces grid-local + + def test_dc_decimation_reduces(self): + """Decimation reduces the mesh: uniform ``reduce=4`` yields strictly fewer (still valid, + in-range) vertices than full detail, and curvature-adaptive ``adaptivity=0.5`` yields a + non-empty mesh with no more vertices than full detail. The decimated paths use double + atomics, so these are count/validity bounds rather than exact-equality checks.""" + v0, _, _ = self.mesh_grid.dual_contour(self.phi, iso=0.0) + # uniform reduce collapses ~F^2 fewer vertices + vr, fr, _ = self.mesh_grid.dual_contour(self.phi, iso=0.0, reduce=4) + self.assertGreater(vr.shape[0], 0) + self.assertLess(vr.shape[0], v0.shape[0]) + self.assertLess(int(fr.max()), vr.shape[0]) + # curvature-adaptive: non-empty, no more verts than full detail + va, fa, _ = self.mesh_grid.dual_contour(self.phi, iso=0.0, adaptivity=0.5) + self.assertGreater(va.shape[0], 0) + self.assertLessEqual(va.shape[0], v0.shape[0]) + self.assertLess(int(fa.max()), va.shape[0]) + + def test_dc_no_surface_is_empty(self): + """A field with no sign crossing (all-positive) has no surface cells, so dual_contour returns + empty but correctly shaped ``(0, 3)`` vertices, faces, and normals rather than erroring.""" + # an all-positive field has no sign crossing -> no surface cells -> empty mesh + field = torch.full((self.grid.num_voxels,), self.bw, device=self.device) + v, f, n = self.grid.dual_contour(field, iso=0.0) + self.assertEqual(v.shape[0], 0) + self.assertEqual(f.shape[0], 0) + self.assertEqual(tuple(v.shape), (0, 3)) + self.assertEqual(tuple(f.shape), (0, 3)) + self.assertEqual(tuple(n.shape), (0, 3)) + + def test_dc_planar_sdf_is_exact(self): + """A planar (linear) SDF is the one input with an exact analytic output: QEF reproduces a + plane exactly, so each interior vertex must lie on the plane ``n·x = d + iso`` and carry the + plane normal ``n``, to float precision. The field is the unclamped exact plane SDF (so + gradients are exact in the interior); vertices near the open mesh boundary are excluded since + there the central-difference gradient falls back to the centre value (the surface SDF has no + neighbour outside the finite grid). Covers axis-aligned and tilted planes plus a nonzero iso + (which shifts the plane).""" + grid = _dense_cube_grid(self.vx, half=12, device=self.device) + world = grid.ijk.float() * self.vx + + def unit(*components: float) -> torch.Tensor: + t = torch.tensor(components, device=self.device, dtype=torch.float32) + return t / t.norm() + + # (plane normal n, offset d, iso level) -> surface is the plane n·x = d + iso + for n, d, iso in [ + (unit(1.0, 0.0, 0.0), 0.0, 0.0), + (unit(0.0, 1.0, 0.0), 0.1, 0.0), + (unit(1.0, 1.0, 1.0), 0.0, 0.0), + (unit(2.0, -1.0, 0.0), -0.05, 0.05), + ]: + phi = world @ n - d # exact SDF of the plane (|grad| = 1); unclamped -> no band artifacts + v, _, nrm = grid.dual_contour(phi, iso=iso) + interior = v.abs().amax(dim=1) < 0.45 # exclude the open boundary (block spans +/-0.6) + self.assertGreater(int(interior.sum()), 0) + v, nrm = v[interior], nrm[interior] + # each interior vertex lies exactly on the plane n·x = d + iso + self.assertLess(float((v @ n - (d + iso)).abs().max()), 1e-4) + # each interior normal equals the constant plane normal + self.assertTrue(torch.allclose(nrm, n.expand_as(nrm), atol=1e-4)) + + def test_dc_sphere_radius_and_normals_accurate(self): + """For a sphere SDF the surface is analytically ``‖x‖ = R``: DC must place every vertex within + a sub-voxel band of that radius, and (the SDF gradient being exactly radial) every normal must + point along ``x/‖x‖``.""" + v, _, n = self.mesh_grid.dual_contour(self.phi, iso=0.0) + radius = v.norm(dim=1) + self.assertLess(float((radius - self.R).abs().max()), 1.0 * self.vx) + self.assertLess(float((radius - self.R).abs().mean()), 0.25 * self.vx) + radial = v / radius.clamp_min(1e-9).unsqueeze(1) + cos = (n * radial).sum(dim=1) # +1 when the normal is the outward radial direction + self.assertGreater(float(cos.min()), 0.8) + self.assertGreater(float(cos.mean()), 0.99) + + def test_dc_sphere_is_valid_closed_genus0(self): + """A closed sphere mesh is a well-formed genus-0 manifold: non-empty, face indices in range, + every triangle non-degenerate, normals aligned with the vertices, no boundary or non-manifold + edges, and Euler characteristic ``V - E + F == 2`` -- an exact topological invariant DC should + satisfy for a well-resolved closed SDF. (Subsumes the generic mesh-validity checks.)""" + v, f, n = self.mesh_grid.dual_contour(self.phi, iso=0.0) + self.assertGreater(v.shape[0], 0) + self.assertGreater(f.shape[0], 0) + self.assertEqual(n.shape, v.shape) + # face indices in range, no degenerate triangles + self.assertGreaterEqual(int(f.min()), 0) + self.assertLess(int(f.max()), v.shape[0]) + self.assertTrue(bool(((f[:, 0] != f[:, 1]) & (f[:, 1] != f[:, 2]) & (f[:, 0] != f[:, 2])).all())) + # closed genus-0 manifold: no boundary / non-manifold edges, Euler characteristic == 2 + num_edges, boundary, nonmanifold = _edge_stats(f, v.shape[0]) + self.assertEqual(boundary, 0) + self.assertEqual(nonmanifold, 0) + self.assertEqual(v.shape[0] - num_edges + f.shape[0], 2) + + def _cube_points(self, vx: float, half: int) -> torch.Tensor: + rng = torch.arange(-half, half + 1, device=self.device, dtype=torch.float32) + ii, jj, kk = torch.meshgrid(rng, rng, rng, indexing="ij") + return torch.stack([ii, jj, kk], dim=-1).reshape(-1, 3) * vx + + +if __name__ == "__main__": + unittest.main() From a48eb8f1cc363095a7d23c698cd36b0b0568caf0 Mon Sep 17 00:00:00 2001 From: Jonathan Swartz Date: Wed, 24 Jun 2026 05:14:40 +0000 Subject: [PATCH 02/10] docs update and use pi from nano math Signed-off-by: Jonathan Swartz --- src/fvdb/detail/ops/DualContour.cu | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/fvdb/detail/ops/DualContour.cu b/src/fvdb/detail/ops/DualContour.cu index f985c3004..e85701b35 100644 --- a/src/fvdb/detail/ops/DualContour.cu +++ b/src/fvdb/detail/ops/DualContour.cu @@ -573,6 +573,9 @@ orientTrianglesKernel(int64_t numTriangles, // bit-reproducible run to run. // This is a simplified, single-level form of the octree-based adaptive DC simplification of [Ju et // al. 2002]; OpenVDB's VolumeToMesh instead does the full seam-stitched octree region merge. + +/// @brief Gather a cell's 8 cube-corner SDF values and gradients from the neighbour table into +/// `cornerSdf` / `cornerGradient` (shared by the decimation kernels' per-cell QEF work). __device__ inline void loadCellCorners(int32_t cellValueIndex, const int32_t *neighborTable, @@ -1041,7 +1044,7 @@ meshOneGrid(OnIndexGridT *grid, .clamp_min(1) .to(torch::kFloat32); const float flatCosThreshold = - (float)std::cos(adaptivity * 60.0 * 3.14159265358979323846 / 180.0); + (float)std::cos(adaptivity * 60.0 * nanovdb::math::pi() / 180.0); torch::Tensor blockFlatBuf = (blockNormalSum.norm(2, /*dim=*/1) / blockCellCount >= flatCosThreshold) .to(torch::kUInt8); From a33dabb781c345f307065459dd9c3f1e1f982e43 Mon Sep 17 00:00:00 2001 From: Jonathan Swartz Date: Thu, 25 Jun 2026 00:48:32 +0000 Subject: [PATCH 03/10] Address review: guard int32 value indexing, document isotropic-voxel assumption Add a per-grid TORCH_CHECK in dualContour so a grid with more than 2^31 voxels fails clearly instead of overflowing the int32 value-index buffers (neighbour table / surface cells / coords). Document that vertex placement and the emitted normals are computed in index space, which is exact for isotropic voxel sizes -- the assumption of the reinitialize_sdf/retopologize_sdf pipeline this op consumes. (Scaling the gradient by the voxel size, as suggested in review, is a no-op for isotropic voxels and would break the index-space QEF for anisotropic ones; full anisotropic support would need a world-space QEF.) Co-Authored-By: Claude Opus 4.8 Signed-off-by: Jonathan Swartz --- src/fvdb/detail/ops/DualContour.cu | 13 ++++++++++++- src/fvdb/detail/ops/DualContour.h | 5 +++++ 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/src/fvdb/detail/ops/DualContour.cu b/src/fvdb/detail/ops/DualContour.cu index e85701b35..383026b9c 100644 --- a/src/fvdb/detail/ops/DualContour.cu +++ b/src/fvdb/detail/ops/DualContour.cu @@ -43,6 +43,7 @@ #include #include +#include #include namespace fvdb { @@ -224,7 +225,9 @@ gatherFusedKernel(const OnIndexGridT *grid, (int32_t)stencil[geometry.usedSpokes[column]]; } - // central-difference gradient (inactive face stencil==0 -> use the centre value) + // central-difference gradient (inactive face stencil==0 -> use the centre value). Computed in + // index space; normalised, this is the correct world normal for isotropic voxels (the uniform + // scale cancels) -- the assumption of the SDF ops this consumes. See the note in DualContour.h. const float centerSdf = sdf[centerIndex]; auto faceValue = [&](int spoke) { uint64_t neighbor = stencil[spoke]; @@ -1242,6 +1245,14 @@ dualContour(const GridBatchData &batchHdl, VoxelCoordTransform transform = batchHdl.primalTransformAt(b); VBMHelper vbm(grid, stream); const int64_t valueCount = (int64_t)vbm.valueCount; // numVoxels + 1 + // Value indices are stored as int32 in the neighbour table / surface-cell / coord + // buffers; fail loudly rather than silently overflow on an enormous single grid. + TORCH_CHECK(valueCount <= std::numeric_limits::max(), + "dual_contour: grid ", + b, + " has too many voxels (", + valueCount, + ") for int32 value indexing"); // gather the field into a value-indexed SDF buffer: slot 0 = "outside" sentinel (>= // iso) so inactive corners classify as outside, and slots [1..numVoxels] = this grid's diff --git a/src/fvdb/detail/ops/DualContour.h b/src/fvdb/detail/ops/DualContour.h index af4e31eea..f0473b218 100644 --- a/src/fvdb/detail/ops/DualContour.h +++ b/src/fvdb/detail/ops/DualContour.h @@ -26,6 +26,11 @@ namespace ops { /// (@p reduce / @p adaptivity) contracts the connectivity onto fewer vertices. /// Requires a >= ~3-voxel band for a watertight result. /// +/// @note Vertex placement and the emitted normals are computed in index space (the QEF and the +/// central-difference gradient). This is exact for isotropic voxel sizes -- the assumption of +/// the narrow-band SDF ops (reinitialize_sdf / retopologize_sdf) this op consumes, which use a +/// scalar voxel size; for strongly anisotropic voxels the emitted normals are approximate. +/// /// @param batchHdl Grid batch defining the sparse topology. /// @param field Per-voxel signed field: a floating-point JaggedTensor of shape [B, -1]. /// @param iso Isovalue of the surface to extract. From fc31dbfea520fd4042f36d4902ad8e35e7c0c4c6 Mon Sep 17 00:00:00 2001 From: Jonathan Swartz Date: Thu, 25 Jun 2026 00:51:43 +0000 Subject: [PATCH 04/10] format fix Signed-off-by: Jonathan Swartz --- src/fvdb/detail/ops/DualContour.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/fvdb/detail/ops/DualContour.h b/src/fvdb/detail/ops/DualContour.h index f0473b218..82782712e 100644 --- a/src/fvdb/detail/ops/DualContour.h +++ b/src/fvdb/detail/ops/DualContour.h @@ -28,8 +28,8 @@ namespace ops { /// /// @note Vertex placement and the emitted normals are computed in index space (the QEF and the /// central-difference gradient). This is exact for isotropic voxel sizes -- the assumption of -/// the narrow-band SDF ops (reinitialize_sdf / retopologize_sdf) this op consumes, which use a -/// scalar voxel size; for strongly anisotropic voxels the emitted normals are approximate. +/// the narrow-band SDF ops (reinitialize_sdf / retopologize_sdf) this op consumes, which use +/// a scalar voxel size; for strongly anisotropic voxels the emitted normals are approximate. /// /// @param batchHdl Grid batch defining the sparse topology. /// @param field Per-voxel signed field: a floating-point JaggedTensor of shape [B, -1]. From 7bf3003d0a43439065a93f715375a16ca4740901 Mon Sep 17 00:00:00 2001 From: Jonathan Swartz Date: Thu, 25 Jun 2026 01:07:38 +0000 Subject: [PATCH 05/10] Address review: clarify VBM decode-helper usage docs and fix a test docstring Clarify VoxelBlockManagerHelper.h: document the two decode paths together (vbmDecodeBlock with a caller-declared __shared__ VbmBlockMaps vs vbmDecodeFaceStencil which owns its own shared scratch), and warn that whichever is used must be called at most once per kernel since each allocates a block-sized __shared__ array. Fix the test_dc_sphere_is_valid_closed_genus0 docstring: it claimed to check normal/vertex alignment, but the body only checks n.shape == v.shape. Reword to 'one normal per vertex' and point to test_dc_sphere_radius_and_normals_accurate, which asserts normal direction. Co-Authored-By: Claude Opus 4.8 Signed-off-by: Jonathan Swartz --- src/fvdb/detail/utils/cuda/VoxelBlockManagerHelper.h | 12 ++++++++++-- tests/unit/test_dc.py | 7 ++++--- 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/src/fvdb/detail/utils/cuda/VoxelBlockManagerHelper.h b/src/fvdb/detail/utils/cuda/VoxelBlockManagerHelper.h index dc2427adf..0bf66ba78 100644 --- a/src/fvdb/detail/utils/cuda/VoxelBlockManagerHelper.h +++ b/src/fvdb/detail/utils/cuda/VoxelBlockManagerHelper.h @@ -52,8 +52,14 @@ struct VBMHelper { }; // ------------------------- per-block decode helpers (device) ------------------------- +// Two ways to decode a block (pick one, and invoke it at most once per kernel since each path +// allocates a block-sized __shared__ array): +// * vbmDecodeBlock(.., maps) -- you declare `__shared__ VbmBlockMaps maps;` in the kernel and own +// it afterwards (use this when the kernel needs the raw maps, e.g. for a custom stencil read). +// * vbmDecodeFaceStencil(.., out) -- the 6-face preamble; it owns its own __shared__ VbmBlockMaps +// internally, so the kernel declares no shared scratch of its own. +// // Per-VBM-block shared decode scratch: the inverse leaf/offset maps for one 2^Log2BlockWidth block. -// Declare one as __shared__ in the kernel and hand it to vbmDecodeBlock. template struct VbmBlockMaps { uint32_t leafIndex[1 << Log2BlockWidth]; uint16_t voxelOffset[1 << Log2BlockWidth]; @@ -109,7 +115,9 @@ vbmReadFaceStencil(const OnIndexGridT *grid, const VbmBlockMaps // Full face-stencil kernel preamble: decode this block, and for active threads read the 6-face // stencil into `out`. Returns false (caller should return) for unused slots. Owns its __shared__ -// scratch, so the kernel needs no shared declaration of its own. +// scratch (a VbmBlockMaps), so the kernel needs no shared declaration of its own -- but for that +// reason it must be called at most once per kernel. For any other neighbour pattern, use +// vbmDecodeBlock with your own `__shared__ VbmBlockMaps` instead. template __device__ inline bool vbmDecodeFaceStencil(const OnIndexGridT *grid, diff --git a/tests/unit/test_dc.py b/tests/unit/test_dc.py index 456c76424..9b347c0f5 100644 --- a/tests/unit/test_dc.py +++ b/tests/unit/test_dc.py @@ -168,9 +168,10 @@ def test_dc_sphere_radius_and_normals_accurate(self): def test_dc_sphere_is_valid_closed_genus0(self): """A closed sphere mesh is a well-formed genus-0 manifold: non-empty, face indices in range, - every triangle non-degenerate, normals aligned with the vertices, no boundary or non-manifold - edges, and Euler characteristic ``V - E + F == 2`` -- an exact topological invariant DC should - satisfy for a well-resolved closed SDF. (Subsumes the generic mesh-validity checks.)""" + every triangle non-degenerate, one normal per vertex (n.shape == v.shape), no boundary or + non-manifold edges, and Euler characteristic ``V - E + F == 2`` -- an exact topological + invariant DC should satisfy for a well-resolved closed SDF. (Normal *direction* is checked in + test_dc_sphere_radius_and_normals_accurate; this subsumes the generic mesh-validity checks.)""" v, f, n = self.mesh_grid.dual_contour(self.phi, iso=0.0) self.assertGreater(v.shape[0], 0) self.assertGreater(f.shape[0], 0) From 1df595475c867fb4365d6b4d273947de31ea232f Mon Sep 17 00:00:00 2001 From: Jonathan Swartz Date: Thu, 25 Jun 2026 01:10:02 +0000 Subject: [PATCH 06/10] Format fix Signed-off-by: Jonathan Swartz --- src/fvdb/detail/utils/cuda/VoxelBlockManagerHelper.h | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/fvdb/detail/utils/cuda/VoxelBlockManagerHelper.h b/src/fvdb/detail/utils/cuda/VoxelBlockManagerHelper.h index 0bf66ba78..27ec5eb25 100644 --- a/src/fvdb/detail/utils/cuda/VoxelBlockManagerHelper.h +++ b/src/fvdb/detail/utils/cuda/VoxelBlockManagerHelper.h @@ -54,10 +54,11 @@ struct VBMHelper { // ------------------------- per-block decode helpers (device) ------------------------- // Two ways to decode a block (pick one, and invoke it at most once per kernel since each path // allocates a block-sized __shared__ array): -// * vbmDecodeBlock(.., maps) -- you declare `__shared__ VbmBlockMaps maps;` in the kernel and own -// it afterwards (use this when the kernel needs the raw maps, e.g. for a custom stencil read). +// * vbmDecodeBlock(.., maps) -- you declare `__shared__ VbmBlockMaps maps;` in the kernel and +// own it afterwards (use this when the kernel needs the raw maps, e.g. for a custom stencil +// read). // * vbmDecodeFaceStencil(.., out) -- the 6-face preamble; it owns its own __shared__ VbmBlockMaps -// internally, so the kernel declares no shared scratch of its own. +// internally, so the kernel declares no shared scratch of its own. // // Per-VBM-block shared decode scratch: the inverse leaf/offset maps for one 2^Log2BlockWidth block. template struct VbmBlockMaps { From b1bd0a88bf3451d4891010955a743e402635ddd4 Mon Sep 17 00:00:00 2001 From: Jonathan Swartz Date: Thu, 25 Jun 2026 01:23:17 +0000 Subject: [PATCH 07/10] Address review: clarify reduce=1 is full detail only when adaptivity==0 The reduce parameter docs said reduce=1 means full resolution, but the adaptive path (adaptivity>0) collapses flat blocks using a default coarse block size (8) even at reduce=1. Clarify across DualContour.h and the Python docstrings (grid.py, grid_batch.py, functional/_meshing.py) that reduce is the decimation block size, that reduce=1 is full detail only when adaptivity==0, and that adaptivity>0 simplifies the mesh even at reduce=1. Co-Authored-By: Claude Opus 4.8 Signed-off-by: Jonathan Swartz --- fvdb/functional/_meshing.py | 9 ++++----- fvdb/grid.py | 4 ++-- fvdb/grid_batch.py | 4 ++-- src/fvdb/detail/ops/DualContour.h | 9 ++++++--- 4 files changed, 14 insertions(+), 12 deletions(-) diff --git a/fvdb/functional/_meshing.py b/fvdb/functional/_meshing.py index f448a32a8..4f485d4bf 100644 --- a/fvdb/functional/_meshing.py +++ b/fvdb/functional/_meshing.py @@ -88,9 +88,8 @@ def dual_contour_batch( grid (GridBatch): The grid batch defining the sparse topology. field (JaggedTensor): Per-voxel signed distance values. iso (float): Isovalue at which to extract the surface. Default ``0.0``. - reduce (int): Uniform ``F x F x F`` cluster-collapse decimation factor (``1`` = full detail). - adaptivity (float): Curvature-adaptive decimation in ``[0, 1.5]`` (``0`` = uniform/off): - collapse flat blocks while keeping detail at features. + reduce (int): Uniform ``F x F x F`` cluster-collapse decimation block size. ``reduce=1`` is full detail only when ``adaptivity == 0``; when ``adaptivity > 0`` it instead sets the coarse-block size for the adaptive collapse (default ``8`` when ``reduce <= 1``). + adaptivity (float): Curvature-adaptive decimation in ``[0, 1.5]`` (``0`` = off). When ``> 0``, flat blocks are collapsed while features keep full detail, so the mesh is simplified even at ``reduce=1``. Returns: vertices (JaggedTensor): Mesh vertex positions, shape ``(B, -1, 3)``. @@ -117,8 +116,8 @@ def dual_contour_single( grid (Grid): The single grid defining the sparse topology. field (torch.Tensor): Per-voxel signed distance values. iso (float): Isovalue at which to extract the surface. Default ``0.0``. - reduce (int): Uniform ``F x F x F`` cluster-collapse decimation factor (``1`` = full detail). - adaptivity (float): Curvature-adaptive decimation in ``[0, 1.5]`` (``0`` = uniform/off). + reduce (int): Uniform ``F x F x F`` cluster-collapse decimation block size. ``reduce=1`` is full detail only when ``adaptivity == 0``; when ``adaptivity > 0`` it instead sets the coarse-block size for the adaptive collapse (default ``8`` when ``reduce <= 1``). + adaptivity (float): Curvature-adaptive decimation in ``[0, 1.5]`` (``0`` = off). When ``> 0``, flat blocks are collapsed while features keep full detail, so the mesh is simplified even at ``reduce=1``. Returns: vertices (torch.Tensor): Vertex positions, shape ``(V, 3)``. diff --git a/fvdb/grid.py b/fvdb/grid.py index c98ac21c6..3a81037c4 100644 --- a/fvdb/grid.py +++ b/fvdb/grid.py @@ -1499,8 +1499,8 @@ def dual_contour( Args: field (torch.Tensor): Per-voxel signed distance values, shape ``(num_voxels,)``. iso (float): Isovalue at which to extract the surface. - reduce (int): Uniform ``F x F x F`` cluster-collapse decimation factor (``1`` = full detail). - adaptivity (float): Curvature-adaptive decimation in ``[0, 1.5]`` (``0`` = uniform/off). + reduce (int): Uniform ``F x F x F`` cluster-collapse decimation block size. ``reduce=1`` is full detail only when ``adaptivity == 0``; when ``adaptivity > 0`` it instead sets the coarse-block size for the adaptive collapse (default ``8`` when ``reduce <= 1``). + adaptivity (float): Curvature-adaptive decimation in ``[0, 1.5]`` (``0`` = off). When ``> 0``, flat blocks are collapsed while features keep full detail, so the mesh is simplified even at ``reduce=1``. Returns: vertices (torch.Tensor): Vertex positions, shape ``(V, 3)``. diff --git a/fvdb/grid_batch.py b/fvdb/grid_batch.py index f71b6666c..bcaa07a37 100644 --- a/fvdb/grid_batch.py +++ b/fvdb/grid_batch.py @@ -1045,8 +1045,8 @@ def dual_contour( Args: field (JaggedTensor): Per-voxel signed distance values. iso (float): Isovalue at which to extract the surface. - reduce (int): Uniform ``F x F x F`` cluster-collapse decimation factor (``1`` = full detail). - adaptivity (float): Curvature-adaptive decimation in ``[0, 1.5]`` (``0`` = uniform/off). + reduce (int): Uniform ``F x F x F`` cluster-collapse decimation block size. ``reduce=1`` is full detail only when ``adaptivity == 0``; when ``adaptivity > 0`` it instead sets the coarse-block size for the adaptive collapse (default ``8`` when ``reduce <= 1``). + adaptivity (float): Curvature-adaptive decimation in ``[0, 1.5]`` (``0`` = off). When ``> 0``, flat blocks are collapsed while features keep full detail, so the mesh is simplified even at ``reduce=1``. Returns: vertices (JaggedTensor): Mesh vertex positions, shape ``(batch_size, -1, 3)``. diff --git a/src/fvdb/detail/ops/DualContour.h b/src/fvdb/detail/ops/DualContour.h index 82782712e..da891a2e1 100644 --- a/src/fvdb/detail/ops/DualContour.h +++ b/src/fvdb/detail/ops/DualContour.h @@ -34,9 +34,12 @@ namespace ops { /// @param batchHdl Grid batch defining the sparse topology. /// @param field Per-voxel signed field: a floating-point JaggedTensor of shape [B, -1]. /// @param iso Isovalue of the surface to extract. -/// @param reduce Uniform F x F x F cluster-collapse decimation factor (1 = full resolution). -/// @param adaptivity Curvature-adaptive decimation in [0, 1.5] (0 = uniform/off): collapse flat -/// blocks while keeping full detail at features. +/// @param reduce Uniform F x F x F cluster-collapse decimation block size. reduce=1 is full +/// resolution only when adaptivity==0; when adaptivity>0 it instead sets the +/// coarse-block size for the adaptive collapse (defaulting to 8 when reduce<=1). +/// @param adaptivity Curvature-adaptive decimation in [0, 1.5] (0 = off). When >0, flat blocks are +/// collapsed while features keep full detail -- so the mesh is simplified even at +/// reduce=1. /// @return A (vertices, faces, normals) tuple, each jagged over the grid batch: vertices and /// normals are float32 JaggedTensors of shape [B, -1, 3] and faces is an int64 JaggedTensor /// of shape [B, -1, 3]. faces holds grid-local triangle vertex indices; normals is the From 2a2a2e9d8e6b2e4d4a2c112cecbf693422752791 Mon Sep 17 00:00:00 2001 From: Jonathan Swartz Date: Thu, 25 Jun 2026 01:25:22 +0000 Subject: [PATCH 08/10] format Signed-off-by: Jonathan Swartz --- src/fvdb/detail/ops/DualContour.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/fvdb/detail/ops/DualContour.h b/src/fvdb/detail/ops/DualContour.h index da891a2e1..a0adc7ef0 100644 --- a/src/fvdb/detail/ops/DualContour.h +++ b/src/fvdb/detail/ops/DualContour.h @@ -38,8 +38,8 @@ namespace ops { /// resolution only when adaptivity==0; when adaptivity>0 it instead sets the /// coarse-block size for the adaptive collapse (defaulting to 8 when reduce<=1). /// @param adaptivity Curvature-adaptive decimation in [0, 1.5] (0 = off). When >0, flat blocks are -/// collapsed while features keep full detail -- so the mesh is simplified even at -/// reduce=1. +/// collapsed while features keep full detail -- so the mesh is simplified even +/// at reduce=1. /// @return A (vertices, faces, normals) tuple, each jagged over the grid batch: vertices and /// normals are float32 JaggedTensors of shape [B, -1, 3] and faces is an int64 JaggedTensor /// of shape [B, -1, 3]. faces holds grid-local triangle vertex indices; normals is the From 4c42fcddaf60b57fbdfece6e26d80be61d3415fa Mon Sep 17 00:00:00 2001 From: Jonathan Swartz Date: Thu, 25 Jun 2026 03:47:50 +0000 Subject: [PATCH 09/10] DRY up the shared edge crossing code Improve some organization and comments Signed-off-by: Jonathan Swartz --- src/fvdb/detail/ops/DualContour.cu | 201 ++++++++++++++--------------- 1 file changed, 98 insertions(+), 103 deletions(-) diff --git a/src/fvdb/detail/ops/DualContour.cu b/src/fvdb/detail/ops/DualContour.cu index 383026b9c..180fa4e92 100644 --- a/src/fvdb/detail/ops/DualContour.cu +++ b/src/fvdb/detail/ops/DualContour.cu @@ -38,6 +38,7 @@ #include #include +#include #include #include @@ -53,25 +54,15 @@ namespace ops { namespace { // ------------------------- box-stencil geometry (compile-time) ------------------------- -// A 3x3x3 box-stencil neighbour is addressed by a spoke index in [0,27): -// spoke = (di+1)*9 + (dj+1)*3 + (dk+1); the centre (the voxel itself) is spoke 13. The mesher only -// touches a subset of the 27 spokes (8 cube corners, 6 faces, 3 connectivity fans), so the gather -// emits just those `numColumns` columns. All of this is static, so it is computed once at compile -// time into `kGeometry` and passed by value into the kernels (no __constant__ / per-call setup). -constexpr int -spokeIndex(int di, int dj, int dk) { - return (di + 1) * 9 + (dj + 1) * 3 + (dk + 1); -} - -/// @brief Column of `spoke` within the sorted-unique `usedSpokes[0..numColumns)`, or -1 if absent. -/// Binary search (std::lower_bound is constexpr in C++20) since usedSpokes is sorted. -constexpr int -columnOfSpoke(const int *usedSpokes, int numColumns, int spoke) { - const int *end = usedSpokes + numColumns; - const int *found = std::lower_bound(usedSpokes, end, spoke); - return (found != end && *found == spoke) ? int(found - usedSpokes) : -1; -} +/// @brief Compile-time description of the 3x3x3 box stencil used by the mesher. +/// +/// A box-stencil neighbour is addressed by a spoke index in [0,27): +/// `spoke = (di+1)*9 + (dj+1)*3 + (dk+1)`; the centre (the voxel itself) is spoke 13. The mesher +/// only touches a subset of the 27 spokes (8 cube corners, 6 faces, 3 connectivity fans), so the +/// gather emits just those `numColumns` columns. All of this is static, so it is computed once at +/// compile time into `kGeometry` and passed by value into the kernels (no `__constant__` / +/// per-call setup). struct BoxStencilGeometry { int numColumns{0}; int usedSpokes[27]{}; // the `numColumns` gathered spokes (sorted, unique) @@ -88,6 +79,20 @@ struct BoxStencilGeometry { float cornerOffset[8][3]{}; // the 8 corner offsets in [0,1]^3 }; +constexpr int +spokeIndex(int di, int dj, int dk) { + return (di + 1) * 9 + (dj + 1) * 3 + (dk + 1); +} + +/// @brief Column of `spoke` within the sorted-unique `usedSpokes[0..numColumns)`, or -1 if absent. +/// Binary search (std::lower_bound is constexpr in C++20) since usedSpokes is sorted. +constexpr int +columnOfSpoke(const int *usedSpokes, int numColumns, int spoke) { + const int *end = usedSpokes + numColumns; + const int *found = std::lower_bound(usedSpokes, end, spoke); + return (found != end && *found == spoke) ? int(found - usedSpokes) : -1; +} + constexpr BoxStencilGeometry makeBoxStencilGeometry() { BoxStencilGeometry geometry{}; @@ -262,6 +267,54 @@ gatherFusedKernel(const OnIndexGridT *grid, voxelCoord[centerIndex * 3 + 2] = globalCoord[2]; } +/// @brief One cube edge's iso-surface crossing: the linearly interpolated crossing point (in the +/// cell-local [0,1]^3 frame) and the unit surface normal there. +struct EdgeCrossing { + float point[3]; // crossing position in the cell's [0,1]^3 frame + float normal[3]; // unit surface normal at the crossing (from the interpolated corner gradients) +}; + +/// @brief Compute cube edge `edge`'s iso crossing from the cell's 8 corner SDFs/gradients. Returns +/// false (leaving `out` untouched) when both endpoints lie on the same side of `iso`, i.e. the edge +/// does not cross the surface. This is the per-edge primitive the three QEF kernels share +/// (placeQefVerticesKernel, cellNormalKernel, clusterQefAccumKernel), so the 12-edge +/// crossing/normal interpolation lives in exactly one place. +__device__ inline bool +computeEdgeCrossing(int edge, + const float cornerSdf[8], + const float cornerGradient[8][3], + BoxStencilGeometry geometry, + float iso, + EdgeCrossing &out) { + const int cornerA = geometry.edgeCornerA[edge], cornerB = geometry.edgeCornerB[edge]; + const float sdfA = cornerSdf[cornerA], sdfB = cornerSdf[cornerB]; + if ((sdfA < iso) == (sdfB < iso)) + return false; // endpoints on the same side of iso -> this edge does not cross the surface + // crossingT = fraction along A->B where the linearly-interpolated sdf reaches iso (the edge + // root); fall back to the midpoint if the two endpoints are (numerically) equal + const float sdfDelta = sdfB - sdfA; + float crossingT = (fabsf(sdfDelta) > 1e-12f) ? (iso - sdfA) / sdfDelta : 0.5f; + crossingT = fminf(fmaxf(crossingT, 0.f), 1.f); + // crossing point and (un-normalised) normal, both linearly interpolated between the two corners + float edgeNormal[3]; +#pragma unroll + for (int axis = 0; axis < 3; ++axis) { + out.point[axis] = geometry.cornerOffset[cornerA][axis] + + crossingT * (geometry.cornerOffset[cornerB][axis] - + geometry.cornerOffset[cornerA][axis]); + edgeNormal[axis] = + cornerGradient[cornerA][axis] + + crossingT * (cornerGradient[cornerB][axis] - cornerGradient[cornerA][axis]); + } + // normalise n_i (the eps keeps a (near-)zero gradient finite) + const float invLen = rsqrtf(edgeNormal[0] * edgeNormal[0] + edgeNormal[1] * edgeNormal[1] + + edgeNormal[2] * edgeNormal[2] + 1e-24f); + out.normal[0] = edgeNormal[0] * invLen; + out.normal[1] = edgeNormal[1] * invLen; + out.normal[2] = edgeNormal[2] * invLen; + return true; +} + /// @brief (QEF) Place one dual-contouring vertex per surface cell. The vertex is the point x /// minimising the quadratic error function E(x) = sum_i [ n_i . (x - p_i) ]^2 over the cell's /// edge intersections, where p_i is the i-th edge zero-crossing and n_i the unit surface normal @@ -312,33 +365,11 @@ placeQefVerticesKernel(const int32_t *surfaceCells, // contributes one tangent-plane constraint (its crossing point p_i and the normal n_i there) #pragma unroll for (int edge = 0; edge < 12; ++edge) { - const int cornerA = geometry.edgeCornerA[edge], cornerB = geometry.edgeCornerB[edge]; - const float sdfA = cornerSdf[cornerA], sdfB = cornerSdf[cornerB]; - if ((sdfA < iso) == (sdfB < iso)) - continue; // endpoints on the same side of iso -> this edge does not cross the surface - // crossingT = fraction along A->B where the linearly-interpolated sdf reaches iso (the edge - // root); fall back to the midpoint if the two endpoints are (numerically) equal - float sdfDelta = sdfB - sdfA; - float crossingT = (fabsf(sdfDelta) > 1e-12f) ? (iso - sdfA) / sdfDelta : 0.5f; - crossingT = fminf(fmaxf(crossingT, 0.f), 1.f); - // crossing point p_i and normal n_i, both linearly interpolated between the two corners - // (the corner gradients supply the surface normal) - float crossing[3], edgeNormal[3]; -#pragma unroll - for (int axis = 0; axis < 3; ++axis) { - crossing[axis] = geometry.cornerOffset[cornerA][axis] + - crossingT * (geometry.cornerOffset[cornerB][axis] - - geometry.cornerOffset[cornerA][axis]); - edgeNormal[axis] = - cornerGradient[cornerA][axis] + - crossingT * (cornerGradient[cornerB][axis] - cornerGradient[cornerA][axis]); - } - // normalise n_i (the eps keeps a (near-)zero gradient finite) - float invLen = rsqrtf(edgeNormal[0] * edgeNormal[0] + edgeNormal[1] * edgeNormal[1] + - edgeNormal[2] * edgeNormal[2] + 1e-24f); - double nx = edgeNormal[0] * invLen, ny = edgeNormal[1] * invLen, - nz = edgeNormal[2] * invLen; - double normalDotCrossing = nx * crossing[0] + ny * crossing[1] + nz * crossing[2]; + EdgeCrossing ec; + if (!computeEdgeCrossing(edge, cornerSdf, cornerGradient, geometry, iso, ec)) + continue; + const double nx = ec.normal[0], ny = ec.normal[1], nz = ec.normal[2]; + const double normalDotCrossing = nx * ec.point[0] + ny * ec.point[1] + nz * ec.point[2]; // A += n n^T (6 unique entries), b-term normalRhs += n (n.p), and track the crossing // centroid normalMatrix[0] += nx * nx; @@ -350,9 +381,9 @@ placeQefVerticesKernel(const int32_t *surfaceCells, normalRhs[0] += nx * normalDotCrossing; normalRhs[1] += ny * normalDotCrossing; normalRhs[2] += nz * normalDotCrossing; - crossingSum[0] += crossing[0]; - crossingSum[1] += crossing[1]; - crossingSum[2] += crossing[2]; + crossingSum[0] += ec.point[0]; + crossingSum[1] += ec.point[1]; + crossingSum[2] += ec.point[2]; ++numCrossings; } double localVertex[3]; @@ -621,23 +652,12 @@ cellNormalKernel(const int32_t *surfaceCells, float normalSum[3] = {0, 0, 0}; #pragma unroll for (int edge = 0; edge < 12; ++edge) { - int cornerA = geometry.edgeCornerA[edge], cornerB = geometry.edgeCornerB[edge]; - float sdfA = cornerSdf[cornerA], sdfB = cornerSdf[cornerB]; - if ((sdfA < iso) == (sdfB < iso)) + EdgeCrossing ec; + if (!computeEdgeCrossing(edge, cornerSdf, cornerGradient, geometry, iso, ec)) continue; - float sdfDelta = sdfB - sdfA; - float crossingT = (fabsf(sdfDelta) > 1e-12f) ? (iso - sdfA) / sdfDelta : 0.5f; - crossingT = fminf(fmaxf(crossingT, 0.f), 1.f); - float nx = cornerGradient[cornerA][0] + - crossingT * (cornerGradient[cornerB][0] - cornerGradient[cornerA][0]), - ny = cornerGradient[cornerA][1] + - crossingT * (cornerGradient[cornerB][1] - cornerGradient[cornerA][1]), - nz = cornerGradient[cornerA][2] + - crossingT * (cornerGradient[cornerB][2] - cornerGradient[cornerA][2]); - float invLen = rsqrtf(nx * nx + ny * ny + nz * nz + 1e-24f); - normalSum[0] += nx * invLen; - normalSum[1] += ny * invLen; - normalSum[2] += nz * invLen; + normalSum[0] += ec.normal[0]; + normalSum[1] += ec.normal[1]; + normalSum[2] += ec.normal[2]; } float invLen = rsqrtf(normalSum[0] * normalSum[0] + normalSum[1] * normalSum[1] + normalSum[2] * normalSum[2] + 1e-24f); @@ -743,36 +763,17 @@ clusterQefAccumKernel(const int32_t *surfaceCells, cellOriginZ = voxelCoord[int64_t(cellValueIndex) * 3 + 2] - originZ; double localMatrix[6] = {0, 0, 0, 0, 0, 0}, localRhs[3] = {0, 0, 0}, localCrossingSum[3] = {0, 0, 0}, localNormalSum[3] = {0, 0, 0}; + int numCrossings = 0; #pragma unroll for (int edge = 0; edge < 12; ++edge) { - int cornerA = geometry.edgeCornerA[edge], cornerB = geometry.edgeCornerB[edge]; - float sdfA = cornerSdf[cornerA], sdfB = cornerSdf[cornerB]; - if ((sdfA < iso) == (sdfB < iso)) + EdgeCrossing ec; + if (!computeEdgeCrossing(edge, cornerSdf, cornerGradient, geometry, iso, ec)) continue; - float sdfDelta = sdfB - sdfA; - float crossingT = (fabsf(sdfDelta) > 1e-12f) ? (iso - sdfA) / sdfDelta : 0.5f; - crossingT = fminf(fmaxf(crossingT, 0.f), 1.f); - float crossingX = - geometry.cornerOffset[cornerA][0] + - crossingT * (geometry.cornerOffset[cornerB][0] - geometry.cornerOffset[cornerA][0]); - float crossingY = - geometry.cornerOffset[cornerA][1] + - crossingT * (geometry.cornerOffset[cornerB][1] - geometry.cornerOffset[cornerA][1]); - float crossingZ = - geometry.cornerOffset[cornerA][2] + - crossingT * (geometry.cornerOffset[cornerB][2] - geometry.cornerOffset[cornerA][2]); - float nx = cornerGradient[cornerA][0] + - crossingT * (cornerGradient[cornerB][0] - cornerGradient[cornerA][0]), - ny = cornerGradient[cornerA][1] + - crossingT * (cornerGradient[cornerB][1] - cornerGradient[cornerA][1]), - nz = cornerGradient[cornerA][2] + - crossingT * (cornerGradient[cornerB][2] - cornerGradient[cornerA][2]); - float invLen = rsqrtf(nx * nx + ny * ny + nz * nz + 1e-24f); - nx *= invLen; - ny *= invLen; - nz *= invLen; - double px = cellOriginX + crossingX, py = cellOriginY + crossingY, - pz = cellOriginZ + crossingZ, normalDotCrossing = nx * px + ny * py + nz * pz; + // shift the cell-local crossing into the cluster's shared origin frame before accumulating + const float nx = ec.normal[0], ny = ec.normal[1], nz = ec.normal[2]; + const double px = cellOriginX + ec.point[0], py = cellOriginY + ec.point[1], + pz = cellOriginZ + ec.point[2], + normalDotCrossing = nx * px + ny * py + nz * pz; localMatrix[0] += nx * nx; localMatrix[1] += nx * ny; localMatrix[2] += nx * nz; @@ -788,22 +789,16 @@ clusterQefAccumKernel(const int32_t *surfaceCells, localNormalSum[0] += nx; localNormalSum[1] += ny; localNormalSum[2] += nz; + ++numCrossings; } for (int j = 0; j < 6; ++j) - atomicAdd(&normalMatrix[clusterId * 6 + j], localMatrix[j]); + gpuAtomicAddNoReturn(&normalMatrix[clusterId * 6 + j], localMatrix[j]); for (int j = 0; j < 3; ++j) { - atomicAdd(&normalRhs[clusterId * 3 + j], localRhs[j]); - atomicAdd(&crossingSum[clusterId * 3 + j], localCrossingSum[j]); - atomicAdd(&normalSum[clusterId * 3 + j], localNormalSum[j]); - } - int numCrossings = 0; -#pragma unroll - for (int edge = 0; edge < 12; ++edge) { - int cornerA = geometry.edgeCornerA[edge], cornerB = geometry.edgeCornerB[edge]; - if ((cornerSdf[cornerA] < iso) != (cornerSdf[cornerB] < iso)) - ++numCrossings; + gpuAtomicAddNoReturn(&normalRhs[clusterId * 3 + j], localRhs[j]); + gpuAtomicAddNoReturn(&crossingSum[clusterId * 3 + j], localCrossingSum[j]); + gpuAtomicAddNoReturn(&normalSum[clusterId * 3 + j], localNormalSum[j]); } - atomicAdd(&crossingCount[clusterId], (double)numCrossings); + gpuAtomicAddNoReturn(&crossingCount[clusterId], (double)numCrossings); } /// @brief Solve one vertex per cluster from its accumulated QEF -- same centroid re-centring + 0.05 From 2478de1b2f72ac0207d4ca20036058de7a0de1d6 Mon Sep 17 00:00:00 2001 From: Jonathan Swartz Date: Thu, 25 Jun 2026 04:23:23 +0000 Subject: [PATCH 10/10] Address review: clarify corner-0 active check and the adaptive 8^3 block size Add a comment in gatherFusedKernel explaining that the allCornersActive loop skips corner 0 because it is the cell's anchor voxel (cornerSpoke[0]==spoke 13==centerIndex), which is always active in this kernel -- so the guard skips a tautological check, not a real corner. Also add a comment in meshOneGrid explaining the adaptive path's 8^3 coarse-block size (it only sets how coarse flat regions become; feature blocks stay full resolution). Co-Authored-By: Claude Opus 4.8 Signed-off-by: Jonathan Swartz --- src/fvdb/detail/ops/DualContour.cu | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/fvdb/detail/ops/DualContour.cu b/src/fvdb/detail/ops/DualContour.cu index 180fa4e92..a5ce2a843 100644 --- a/src/fvdb/detail/ops/DualContour.cu +++ b/src/fvdb/detail/ops/DualContour.cu @@ -254,6 +254,8 @@ gatherFusedKernel(const OnIndexGridT *grid, for (int corner = 0; corner < 8; ++corner) { uint64_t cornerIndex = stencil[geometry.cornerSpoke[corner]]; numInsideCorners += (sdf[cornerIndex] < iso) ? 1 : 0; + // corner 0 is offset (0,0,0) -- the anchor voxel itself (== centerIndex), which is always + // active here (this thread only runs for an active voxel), so skip its redundant check. if (corner) allCornersActive &= (cornerIndex > 0); } @@ -968,6 +970,10 @@ meshOneGrid(OnIndexGridT *grid, refNormalBuf.data_ptr()); C10_CUDA_KERNEL_LAUNCH_CHECK(); } else { + // Cluster block size in voxels. Uniform decimation uses `reduce` directly. The curvature- + // adaptive path (adaptivity > 0, reduce left at its default) groups cells into 8^3 coarse + // blocks: a flat block collapses to one vertex while feature blocks stay full resolution, + // so 8 sets only how coarse the *flat* regions become. const int blockSize = (adaptivity > 0.0 && reduce <= 1) ? 8 : std::max(1, reduce); // surface-cell coord bounds (for the mixed-radix cluster keys): a gather + amin/amax // reduction in place of a hand-rolled atomic-min/max kernel.