|
| 1 | +"""scikit-image reference helpers for marching-cubes development checks. |
| 2 | +
|
| 3 | +This module is intentionally kept outside the package and test dependency set. |
| 4 | +It imports scikit-image lazily, mirrors the public ``pad`` extension, and is |
| 5 | +used by the parity and benchmark scripts in this directory. |
| 6 | +""" |
| 7 | + |
| 8 | +from __future__ import annotations |
| 9 | + |
| 10 | +from collections.abc import Sequence |
| 11 | + |
| 12 | +import numpy as np |
| 13 | + |
| 14 | + |
| 15 | +def _skimage_marching_cubes(): |
| 16 | + try: |
| 17 | + from skimage.measure import marching_cubes |
| 18 | + except ImportError as error: # pragma: no cover - development only |
| 19 | + raise ImportError( |
| 20 | + "scikit-image is required for the marching-cubes reference " |
| 21 | + "(`pip install scikit-image`)" |
| 22 | + ) from error |
| 23 | + return marching_cubes |
| 24 | + |
| 25 | + |
| 26 | +def reference_marching_cubes( |
| 27 | + volume: np.ndarray, |
| 28 | + level: float | None = None, |
| 29 | + *, |
| 30 | + spacing: Sequence[float] = (1.0, 1.0, 1.0), |
| 31 | + gradient_direction: str = "descent", |
| 32 | + step_size: int = 1, |
| 33 | + allow_degenerate: bool = True, |
| 34 | + method: str = "lewiner", |
| 35 | + mask: np.ndarray | None = None, |
| 36 | + pad: bool = False, |
| 37 | +): |
| 38 | + """Call scikit-image with the same public contract as ``bic.mesh``.""" |
| 39 | + image = np.ascontiguousarray(volume, dtype=np.float32) |
| 40 | + if level is None: |
| 41 | + level = 0.5 * (float(image.min()) + float(image.max())) |
| 42 | + mask_array = None if mask is None else np.ascontiguousarray(np.asarray(mask) != 0) |
| 43 | + spacing_array = np.asarray(spacing, dtype=np.float64) |
| 44 | + if pad: |
| 45 | + image = np.pad(image, 1, mode="constant", constant_values=0) |
| 46 | + if mask_array is not None: |
| 47 | + mask_array = np.pad(mask_array, 1, mode="constant", constant_values=True) |
| 48 | + result = _skimage_marching_cubes()( |
| 49 | + image, |
| 50 | + level, |
| 51 | + spacing=spacing_array, |
| 52 | + gradient_direction=gradient_direction, |
| 53 | + step_size=step_size, |
| 54 | + allow_degenerate=allow_degenerate, |
| 55 | + method=method, |
| 56 | + mask=mask_array, |
| 57 | + ) |
| 58 | + if pad: |
| 59 | + vertices = result[0] - spacing_array.astype(result[0].dtype, copy=False) |
| 60 | + result = (vertices, *result[1:]) |
| 61 | + return result |
| 62 | + |
| 63 | + |
| 64 | +def _sorted_rows(array: np.ndarray) -> tuple[np.ndarray, np.ndarray]: |
| 65 | + if len(array) == 0: |
| 66 | + return array, np.empty(0, dtype=np.int64) |
| 67 | + order = np.lexsort(tuple(array[:, axis] for axis in range(array.shape[1] - 1, -1, -1))) |
| 68 | + return array[order], order |
| 69 | + |
| 70 | + |
| 71 | +def _canonical_faces(vertices: np.ndarray, faces: np.ndarray) -> tuple[np.ndarray, np.ndarray]: |
| 72 | + unique_vertices, inverse = np.unique(vertices, axis=0, return_inverse=True) |
| 73 | + canonical = np.sort(inverse[faces], axis=1) |
| 74 | + canonical, _ = _sorted_rows(canonical) |
| 75 | + return unique_vertices, canonical |
| 76 | + |
| 77 | + |
| 78 | +def _surface_area(vertices: np.ndarray, faces: np.ndarray) -> float: |
| 79 | + total = 0.0 |
| 80 | + for begin in range(0, len(faces), 100_000): |
| 81 | + triangles = vertices[faces[begin : begin + 100_000]].astype(np.float64, copy=False) |
| 82 | + cross = np.cross(triangles[:, 1] - triangles[:, 0], triangles[:, 2] - triangles[:, 0]) |
| 83 | + total += float(np.linalg.norm(cross, axis=1).sum()) * 0.5 |
| 84 | + return total |
| 85 | + |
| 86 | + |
| 87 | +def assert_mesh_matches(actual, reference, *, normal_atol: float = 1e-5) -> None: |
| 88 | + """Assert geometry/topology parity independent of output ordering. |
| 89 | +
|
| 90 | + Normals and local-range values are compared after sorting complete vertex |
| 91 | + records by coordinates. Face winding is tested separately by the package |
| 92 | + tests because canonical triangle comparison intentionally ignores it. |
| 93 | + """ |
| 94 | + actual_vertices, actual_faces, actual_normals, actual_values = actual |
| 95 | + reference_vertices, reference_faces, reference_normals, reference_values = reference |
| 96 | + |
| 97 | + assert actual_vertices.shape == reference_vertices.shape |
| 98 | + assert actual_faces.shape == reference_faces.shape |
| 99 | + assert actual_normals.shape == reference_normals.shape |
| 100 | + assert actual_values.shape == reference_values.shape |
| 101 | + |
| 102 | + actual_unique, actual_triangles = _canonical_faces(actual_vertices, actual_faces) |
| 103 | + reference_unique, reference_triangles = _canonical_faces(reference_vertices, reference_faces) |
| 104 | + np.testing.assert_allclose(actual_unique, reference_unique, rtol=0.0, atol=1e-6) |
| 105 | + np.testing.assert_array_equal(actual_triangles, reference_triangles) |
| 106 | + np.testing.assert_allclose( |
| 107 | + _surface_area(actual_vertices, actual_faces), |
| 108 | + _surface_area(reference_vertices, reference_faces), |
| 109 | + rtol=1e-6, |
| 110 | + atol=1e-8, |
| 111 | + ) |
| 112 | + |
| 113 | + actual_records = np.column_stack( |
| 114 | + (actual_vertices, actual_normals, actual_values) |
| 115 | + ).astype(np.float64, copy=False) |
| 116 | + reference_records = np.column_stack( |
| 117 | + (reference_vertices, reference_normals, reference_values) |
| 118 | + ).astype(np.float64, copy=False) |
| 119 | + actual_records, _ = _sorted_rows(actual_records) |
| 120 | + reference_records, _ = _sorted_rows(reference_records) |
| 121 | + np.testing.assert_allclose( |
| 122 | + actual_records[:, :3], reference_records[:, :3], rtol=0.0, atol=1e-6 |
| 123 | + ) |
| 124 | + np.testing.assert_allclose( |
| 125 | + actual_records[:, 3:6], reference_records[:, 3:6], rtol=0.0, atol=normal_atol |
| 126 | + ) |
| 127 | + np.testing.assert_allclose( |
| 128 | + actual_records[:, 6], reference_records[:, 6], rtol=0.0, atol=1e-6 |
| 129 | + ) |
0 commit comments