Skip to content

marching_cubes

pleat.marching_cubes

Marching-cubes mesh generation from signed-distance fields.

Thin wrapper around :mod:sdf for converting SDFs into triangle meshes that can be exported to STL via :mod:meshio. Used to produce 3D-printable thickened versions of crease patterns. Optional dependency: threed.

CartesianGrid 

CartesianGrid(
    ndim: int,
    origin: ndarray = None,
    cell_size: ndarray = None,
)

Axis aligned cartesian grid in n-dimensional euclidean space

Source code in pleat/marching_cubes.py 
35
36
37
38
39
40
41
42
def __init__(self, ndim: int, origin: np.ndarray = None, cell_size: np.ndarray = None) -> None:
    self.ndim = ndim
    self.origin = np.zeros(self.ndim) if origin is None else origin
    self.cell_size = np.ones(self.ndim, dtype=np.float32) if cell_size is None else cell_size
    assert (
        self.origin.shape == self.cell_size.shape == (self.ndim,)
    ), f"Inconsistent shapes: {self.origin.shape=}, {self.cell_size.shape=}, {self.ndim=}"
    self.cell_center_to_corner_distance = np.linalg.norm(self.cell_size) / 2

cell_centers 

cell_centers(indices: ndarray) -> np.ndarray

Get the centers of grid cells at specified indices

Parameters:

Name Type Description Default
indices ndarray

Array of shape (n, dim), of indices along the different axes.

 required

Returns:

Type Description
ndarray

center positions as array of shape (n, dim).
Source code in pleat/marching_cubes.py 
44
45
46
47
48
49
50
51
52
53
54
def cell_centers(self, indices: np.ndarray) -> np.ndarray:
    """Get the centers of grid cells at specified indices

    Args:
        indices: Array of shape (n, dim), of indices along the different axes.

    Returns:
        center positions as array of shape (n, dim).
    """

    return (indices + 0.5) * self.cell_size[None]

get_1d_axis_grids 

get_1d_axis_grids() -> list['CartesianGrid']

Get a list of 1D grids, corresponding to the axes of this grid

Source code in pleat/marching_cubes.py 
65
66
67
68
69
def get_1d_axis_grids(self) -> list["CartesianGrid"]:
    """Get a list of 1D grids, corresponding to the axes of this grid"""
    return [
        CartesianGrid(ndim=1, origin=o, cell_size=s) for o, s in zip(self.origin[:, None], self.cell_size[:, None])
    ]

subdivide 

subdivide(
    subdivisions: int | tuple | ndarray = 2,
) -> "CartesianGrid"

Subdivide the grid, resulting in a finer grid

Parameters:

Name Type Description Default
subdivisions int | tuple | ndarray

Number of subdivisions per dimension. If integer, each dimension will be subdivided by the given number. Default is 2.

 2

Returns:

Type Description
'CartesianGrid'

The subdivided :class:CartesianGrid
Source code in pleat/marching_cubes.py 
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
def subdivide(self, subdivisions: int | tuple | np.ndarray = 2) -> "CartesianGrid":
    """Subdivide the grid, resulting in a finer grid

    Args:
        subdivisions: Number of subdivisions per dimension. If integer, each dimension will be subdivided by the
            given number. Default is 2.

    Returns:
        The subdivided :class:`CartesianGrid`
    """
    subdivisions = self._parse_subdivisions(subdivisions)

    return CartesianGrid(
        ndim=self.ndim,
        origin=self.origin,
        cell_size=self.cell_size / subdivisions,
    )

coarsen 

coarsen(
    subdivisions: int | tuple | ndarray = 2,
) -> "CartesianGrid"

The inverse to :meth:subdivide.

Source code in pleat/marching_cubes.py 
89
90
91
92
93
94
95
96
97
def coarsen(self, subdivisions: int | tuple | np.ndarray = 2) -> "CartesianGrid":
    """The inverse to :meth:`subdivide`."""
    subdivisions = self._parse_subdivisions(subdivisions)

    return CartesianGrid(
        ndim=self.ndim,
        origin=self.origin,
        cell_size=self.cell_size * subdivisions,
    )

coarsen_indices 

coarsen_indices(
    indices: ndarray, subdivisions: int | tuple | ndarray
) -> np.ndarray

Get list of unique indices of cells in coarser grid, that the grid cells at given indices are part of.

Source code in pleat/marching_cubes.py 
 99
100
101
102
103
def coarsen_indices(self, indices: np.ndarray, subdivisions: int | tuple | np.ndarray) -> np.ndarray:
    """Get list of unique indices of cells in coarser grid, that the grid cells at given indices are part of."""
    subdivisions = self._parse_subdivisions(subdivisions)
    assert len(indices.shape) == 2 and indices.shape[1] == self.ndim, f"{indices.shape=}, {self.ndim=}"
    return np.unique(indices // subdivisions, axis=0)

coarsen_block_indices 

coarsen_block_indices(
    indices_per_dim: list[ndarray],
    subdivisions: int | tuple | ndarray,
) -> list[np.ndarray]

Apply :meth:coarsen_indices independently per axis to a block of indices.

Source code in pleat/marching_cubes.py 
105
106
107
108
109
110
111
112
113
114
115
def coarsen_block_indices(
    self, indices_per_dim: list[np.ndarray], subdivisions: int | tuple | np.ndarray
) -> list[np.ndarray]:
    """Apply :meth:`coarsen_indices` independently per axis to a block of indices."""
    subdivisions = self._parse_subdivisions(subdivisions)
    assert len(indices_per_dim) == self.ndim, f"{len(indices_per_dim)=} != {self.ndim=}"
    axis_grids = self.get_1d_axis_grids()
    return [
        grid.coarsen_indices(indices[:, None], subdivisions=s)[:, 0]
        for grid, indices, s in zip(axis_grids, indices_per_dim, subdivisions)
    ]

subdivide_indices 

subdivide_indices(
    indices: ndarray, subdivisions: int | tuple | ndarray
) -> np.ndarray

Get list of indices of cells in finer grid that encompass the grid cells at given indices.

Source code in pleat/marching_cubes.py 
117
118
119
120
121
122
123
124
125
126
127
128
def subdivide_indices(self, indices: np.ndarray, subdivisions: int | tuple | np.ndarray) -> np.ndarray:
    """Get list of indices of cells in finer grid that encompass the grid cells at given indices."""
    subdivisions = self._parse_subdivisions(subdivisions)

    # Indices of the cells in the finer grid that make up the (0, 0, ..., 0) cell in the coarser grid.
    # Shape (np.prod(subdivisions), self.ndim)
    block_indices = np.stack(np.meshgrid(*[np.arange(s) for s in subdivisions], indexing="ij"), axis=-1).reshape(
        -1, self.ndim
    )

    fine_indices = ((indices * subdivisions)[:, None] + block_indices[None]).reshape(-1, self.ndim)
    return fine_indices

oct_tree_marching_cubes 

oct_tree_marching_cubes(
    f: callable,
    step: ndarray,
    bounds: ndarray | None = None,
) -> tuple[np.ndarray, CartesianGrid]

Locate boundary cells of an SDF using an oct-tree refinement strategy.

Builds successively coarser grids until the bounding box fits in one cell, then refines top-down, keeping only cells whose distance is within the cell radius. Returns the surviving cell indices and the finest grid.

Parameters:

Name Type Description Default
f callable

Signed-distance function mapping (n, 3) positions to (n, 1) distances.

 required
step ndarray

Grid spacing per axis (or scalar applied to all).

 required
bounds ndarray | None

Axis-aligned bounding box ((x0, y0, z0), (x1, y1, z1)) to sample. If None, estimated via :func:sdf.mesh._estimate_bounds.

 None

Returns:

Type Description
ndarray

(indices, grid) where indices is the array of surviving cell
CartesianGrid

indices in the finest grid and grid is that :class:CartesianGrid.
Source code in pleat/marching_cubes.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
def oct_tree_marching_cubes(
    f: callable, step: np.ndarray, bounds: np.ndarray | None = None
) -> tuple[np.ndarray, CartesianGrid]:
    """Locate boundary cells of an SDF using an oct-tree refinement strategy.

    Builds successively coarser grids until the bounding box fits in one
    cell, then refines top-down, keeping only cells whose distance is within
    the cell radius. Returns the surviving cell indices and the finest grid.

    Args:
        f: Signed-distance function mapping ``(n, 3)`` positions to ``(n, 1)``
            distances.
        step: Grid spacing per axis (or scalar applied to all).
        bounds: Axis-aligned bounding box ``((x0, y0, z0), (x1, y1, z1))`` to
            sample. If ``None``, estimated via :func:`sdf.mesh._estimate_bounds`.

    Returns:
        ``(indices, grid)`` where ``indices`` is the array of surviving cell
        indices in the finest grid and ``grid`` is that :class:`CartesianGrid`.
    """
    ndim = 3

    if bounds is None:
        bounds = sdf.mesh._estimate_bounds(f)
    bounds = np.asarray(bounds)

    if isinstance(step, Number):
        step = np.array((step,) * ndim)

    assert np.all(bounds[0] < bounds[1]), f"Bounds must be sorted. Got {bounds=}."
    assert np.all(step > 0), f"Got negative step size: {step=}"

    grids = [CartesianGrid(ndim=ndim, origin=bounds[0], cell_size=step)]

    box_size = bounds[1] - bounds[0]

    subdivisions = 2
    level = 0
    block_indices = [np.arange(n) for n in box_size // step]

    def block_indices_n_cells(block_indices):
        return np.product([len(indices) for indices in block_indices])

    n_cells_naive = block_indices_n_cells(block_indices)
    logger.info("n_cells_naive=%d", n_cells_naive)

    # coarsen until the whole box fits into a single grid cell
    while block_indices_n_cells(block_indices) > 1:
        logger.debug("level %d, n_cells=%d", level, block_indices_n_cells(block_indices))
        new_block_indices = grids[-1].coarsen_block_indices(block_indices, subdivisions=subdivisions)
        assert block_indices_n_cells(new_block_indices) < block_indices_n_cells(block_indices)
        block_indices = new_block_indices
        grids.append(grids[-1].coarsen(subdivisions=subdivisions))
        level += 1

    indices = np.zeros((1, ndim), dtype=np.int32)
    for i, grid in enumerate(reversed(grids)):
        # filter indices
        pts = grid.cell_centers(indices)
        dists = f(pts)[:, 0]  # for some reason, the distances are returned in an array of shape (n, 1)
        mask = np.abs(dists) <= grid.cell_center_to_corner_distance
        logger.debug("level %d, mask.mean=%.3f, mask.sum=%d", len(grids) - 1 - i, mask.mean(), mask.sum())
        indices = indices[mask]

        if i != len(grids) - 1:
            # get indices in finder grid
            indices = grid.subdivide_indices(indices, subdivisions=subdivisions)

    logger.info("n_cells=%d", len(indices))
    return indices, grids[0]