From 5b7b00c4d98bb5e44234732a14c37318f31bddf4 Mon Sep 17 00:00:00 2001 From: uermel Date: Mon, 26 Oct 2026 15:11:33 -0700 Subject: [PATCH] fix: seg2mesh writes closed x/y/z meshes, and mesh2seg fills them back exactly MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit seg2mesh (volume_to_mesh) ran marching cubes on the (z, y, x) array and stored the vertices in that order, so its meshes were mirrored (x and z swapped, triangles facing inward) relative to the tomogram and to every other copick mesh: seg2slab, validbox, picks2mesh and mesh2seg all use (x, y, z) in Å with voxel index i at i * spacing. It also did not pad the volume, so a segmentation touching the volume's edge gave an open mesh. Vertices are now reversed to (x, y, z), which also turns the triangles outward, and the volume is padded by a voxel of background. Fixes #132. seg2mesh --individual-meshes was documented but never reached the converter (one mesh was stored under a literal "{instance_id}" session). It now stores one mesh per label (label_meshes, each label's mask cropped to its bounding box), with {instance_id} replaced by the label. mesh2seg (watertight) cast its rays through the voxel centres, exactly where marching-cubes meshes have their vertices and edges, so crossings were counted twice or missed: a correct seg2mesh mesh of 10521 run 35962's valid-sample segmentation came back with IoU 0.65. Ray origins are now offset by ~0.015 voxel across the ray. Three more fixes in the same code: - a crossing was assigned to its nearest voxel (np.round), which counts the voxel before the surface in where a ray enters and leaves out the last one inside where it exits; a voxel is now inside when its centre is (the first centre after each crossing); - two crossings in the same voxel were applied once (fancy-index assignment with repeated indices); crossings now toggle by parity; - rays started one voxel before the volume, so a mesh reaching further lost its entry crossing; they now start before the mesh. On 10521 run 35962 (boundary-clean valid-sample), seg2mesh then mesh2seg now gives back the segmentation exactly (IoU 1.0, was 0.11). The published pipeline mesh still round-trips, now with voxels decided by their centres: +242,057 / -739 voxels (0.5%) next to its surface. Co-Authored-By: Claude Opus 5.5 --- src/copick_utils/cli/seg2mesh.py | 1 + .../converters/mesh_from_segmentation.py | 92 +++++++++- .../converters/segmentation_from_mesh.py | 95 ++++++---- tests/test_seg2mesh_mesh2seg.py | 163 ++++++++++++++++++ 4 files changed, 312 insertions(+), 39 deletions(-) create mode 100644 tests/test_seg2mesh_mesh2seg.py diff --git a/src/copick_utils/cli/seg2mesh.py b/src/copick_utils/cli/seg2mesh.py index 3c2aa7c..889c1b7 100644 --- a/src/copick_utils/cli/seg2mesh.py +++ b/src/copick_utils/cli/seg2mesh.py @@ -129,6 +129,7 @@ def seg2mesh( workers=workers, level=level, step_size=step_size, + individual_meshes=individual_meshes, ) successful = sum(1 for result in results.values() if result and result.get("processed", 0) > 0) diff --git a/src/copick_utils/converters/mesh_from_segmentation.py b/src/copick_utils/converters/mesh_from_segmentation.py index d2980b8..d8d9901 100644 --- a/src/copick_utils/converters/mesh_from_segmentation.py +++ b/src/copick_utils/converters/mesh_from_segmentation.py @@ -21,6 +21,10 @@ def volume_to_mesh(volume: np.ndarray, voxel_spacing: float, level: float = 0.5, """ Convert a binary volume to a mesh using marching cubes. + Vertices are (x, y, z) in physical units, with voxel index i at i * voxel_spacing, like copick's other meshes + and picks. The volume is padded with a voxel of background first, so the surface is closed where the + segmentation touches the edge of the volume. + Args: volume: Binary volume array with shape (z, y, x) voxel_spacing: Spacing between voxels in physical units @@ -28,24 +32,69 @@ def volume_to_mesh(volume: np.ndarray, voxel_spacing: float, level: float = 0.5, step_size: Step size for marching cubes (higher = coarser mesh) Returns: - Trimesh object representing the mesh + Trimesh object representing the mesh (outward-facing triangles) """ from skimage import measure # Generate mesh using marching cubes verts, faces, normals, values = measure.marching_cubes( - volume.astype(float), + np.pad(volume.astype(float), 1), level=level, step_size=step_size, spacing=(voxel_spacing, voxel_spacing, voxel_spacing), ) - # Create trimesh object - mesh = tm.Trimesh(vertices=verts, faces=faces, vertex_normals=normals) + # marching_cubes works in the array's (z, y, x) order: reverse to (x, y, z) and undo the padding. Reversing the + # axes mirrors the mesh, which also turns its triangles to face outward, so the faces stay as they are. + mesh = tm.Trimesh(vertices=verts[:, ::-1] - voxel_spacing, faces=faces, vertex_normals=normals[:, ::-1]) return mesh +def label_meshes( + volume: np.ndarray, + voxel_spacing: float, + level: float = 0.5, + step_size: int = 1, +) -> Dict[int, tm.Trimesh]: + """ + Convert each non-zero label of a volume to its own mesh (see :func:`volume_to_mesh`). + + Each label is meshed from its own mask, cropped to the label's bounding box. + + Args: + volume: Label volume with shape (z, y, x) + voxel_spacing: Spacing between voxels in physical units + level: Isosurface level for marching cubes (of each label's binary mask) + step_size: Step size for marching cubes + + Returns: + Meshes by label value, in increasing label order + """ + from scipy import ndimage + + if np.issubdtype(volume.dtype, np.integer) and volume.size and volume.min() >= 0: + labels = volume + values = None + else: + # Map the values to consecutive indices, 0 staying background. + values, labels = np.unique(volume, return_inverse=True) + labels = labels.reshape(volume.shape) + if values[0] != 0: + values = np.concatenate([[0], values]) + labels = labels + 1 + + meshes = {} + for index, box in enumerate(ndimage.find_objects(labels), start=1): + if box is None: + continue + mesh = volume_to_mesh(labels[box] == index, voxel_spacing, level, step_size) + # The crop starts at (z, y, x) = box starts: move the mesh there. + mesh.apply_translation(np.array([s.start for s in box][::-1]) * voxel_spacing) + meshes[int(index if values is None else values[index])] = mesh + return meshes + + def mesh_from_segmentation( segmentation: "CopickSegmentation", run: "CopickRun", @@ -55,6 +104,7 @@ def mesh_from_segmentation( level: float = 0.5, step_size: int = 1, voxel_spacing: Optional[float] = None, + individual_meshes: bool = False, **kwargs, ) -> Optional[Tuple["CopickMesh", Dict[str, int]]]: """ @@ -64,14 +114,16 @@ def mesh_from_segmentation( segmentation: CopickSegmentation object to convert run: CopickRun object object_name: Name for the output mesh object - session_id: Session ID for the output mesh + session_id: Session ID for the output mesh; with ``individual_meshes``, a template whose ``{instance_id}`` + is replaced by each label (without it, ``-{label}`` is appended) user_id: User ID for the output mesh level: Isosurface level for marching cubes step_size: Step size for marching cubes + individual_meshes: Store one mesh per non-zero label instead of one mesh of the whole segmentation Returns: - Tuple of (CopickMesh object, stats dict) or None if creation failed. - Stats dict contains 'vertices_created' and 'faces_created'. + Tuple of (CopickMesh object, stats dict) or None if creation failed. With ``individual_meshes``, the first + label's mesh and the totals. Stats dict contains 'vertices_created' and 'faces_created'. """ try: # Load the volume from the segmentation @@ -84,6 +136,32 @@ def mesh_from_segmentation( # Get voxel spacing from segmentation voxel_spacing = segmentation.voxel_size + if individual_meshes: + meshes = label_meshes(volume, voxel_spacing, level, step_size) + if not meshes: + logger.error("No labels to mesh") + return None + first = None + stats = {"vertices_created": 0, "faces_created": 0} + for label, mesh in meshes.items(): + if "{instance_id}" in session_id: + label_session_id = session_id.replace("{instance_id}", str(label)) + else: + label_session_id = f"{session_id}-{label}" + copick_mesh, mesh_stats = store_mesh_with_stats( + run, + mesh, + object_name, + label_session_id, + user_id, + "mesh", + ) + if first is None: + first = copick_mesh + for key in stats: + stats[key] += mesh_stats[key] + return first, stats + # Convert volume to mesh mesh = volume_to_mesh(volume, voxel_spacing, level, step_size) diff --git a/src/copick_utils/converters/segmentation_from_mesh.py b/src/copick_utils/converters/segmentation_from_mesh.py index e2f9d8a..e44edc8 100644 --- a/src/copick_utils/converters/segmentation_from_mesh.py +++ b/src/copick_utils/converters/segmentation_from_mesh.py @@ -41,62 +41,93 @@ def ensure_mesh(trimesh_object) -> Optional[tm.Trimesh]: raise ValueError("Input must be a Trimesh or Scene object") +#: Ray origins sit this far (in voxels) off the voxel centres, across the ray. Meshes made by marching cubes have +#: their vertices and edges on the voxel grid, where rays through the voxel centres would hit edges and vertices +#: exactly and count the crossing twice or not at all. Rounding back to voxel indices removes the offset; a voxel +#: centre within about this distance of a sloped surface may fall on either side of it. +_RAY_OFFSET = (0.01 * np.sqrt(2), 0.01 * np.sqrt(3)) + +#: A crossing this close (in voxels) to a voxel centre counts as before it, so a face through the centres (e.g. a +#: box at the volume's edge) includes them where the ray enters and excludes them where it leaves. +_CENTRE_TOLERANCE = 1e-6 + + +def _crossing_voxels(loc: np.ndarray, voxel_spacing: float, axis: int) -> np.ndarray: + """Voxel indices of ray crossings: the ray's own voxel across it, and along it the first voxel centre after the + crossing, so a voxel is inside when its centre is. Crossings before the volume move to its first voxel.""" + pos = loc / voxel_spacing + idx = np.round(pos).astype("int") + idx[:, axis] = np.maximum(np.ceil(pos[:, axis] - _CENTRE_TOLERANCE).astype("int"), 0) + return idx + + +def _fill_from_crossings( + int_loc: np.ndarray, + voxel_dims: Tuple[int, int, int], + axis: int, +) -> np.ndarray: + """Inside/outside of every voxel along rays on ``axis``, from the voxels of their crossings: a voxel is inside + after an odd number of crossings. Crossings in the same voxel cancel in pairs. Returns a (z, y, x) volume.""" + toggles = np.zeros((voxel_dims[2], voxel_dims[1], voxel_dims[0]), dtype=np.uint8) + zyx = int_loc[int_loc[:, axis] < voxel_dims[axis]][:, ::-1] + np.add.at(toggles, (zyx[:, 0], zyx[:, 1], zyx[:, 2]), 1) + toggles &= 1 + return np.bitwise_xor.accumulate(toggles, axis=2 - axis).astype(bool) + + +def _ray_start(mesh: tm.Trimesh, voxel_spacing: float, axis: int) -> float: + """Where rays along an axis start (in voxels): before the volume and before the mesh.""" + return min(-1.0, np.floor(mesh.bounds[0, axis] / voxel_spacing) - 1) + + def _onesmask_z(mesh: tm.Trimesh, voxel_dims: Tuple[int, int, int], voxel_spacing: float) -> np.ndarray: """Create mask by ray casting in Z direction.""" intersector = RayMeshIntersector(mesh) - # Create a grid of rays in XY plane, shooting in Z direction + # Create a grid of rays in XY plane (off the voxel centres, see _RAY_OFFSET), shooting in Z direction grid_x, grid_y = np.mgrid[0 : voxel_dims[0], 0 : voxel_dims[1]] - ray_grid = np.vstack([grid_x.ravel(), grid_y.ravel(), -np.ones((grid_x.size,))]).T * voxel_spacing + ray_grid = ( + np.vstack( + [ + grid_x.ravel() + _RAY_OFFSET[0], + grid_y.ravel() + _RAY_OFFSET[1], + np.full(grid_x.size, _ray_start(mesh, voxel_spacing, 2)), + ], + ).T + * voxel_spacing + ) ray_dir = np.zeros((ray_grid.shape[0], 3)) ray_dir[:, 2] = 1 loc, _, _ = intersector.intersects_location(ray_grid, ray_dir) - # Convert to voxel coordinates and sort by z - int_loc = np.round(loc / voxel_spacing).astype("int") - sort_idx = int_loc[:, 2].argsort() - int_loc = int_loc[sort_idx, :] - # Build volume by tracking crossings - img = np.zeros((voxel_dims[1], voxel_dims[0]), dtype="bool") - vol = np.zeros((voxel_dims[2], voxel_dims[1], voxel_dims[0]), dtype="bool") - - for z in range(voxel_dims[2]): - idx = int_loc[:, 2] == z - img[int_loc[idx, 1], int_loc[idx, 0]] = np.logical_not(img[int_loc[idx, 1], int_loc[idx, 0]]) - vol[z, :, :] = img - - return vol + return _fill_from_crossings(_crossing_voxels(loc, voxel_spacing, 2), voxel_dims, 2) def _onesmask_x(mesh: tm.Trimesh, voxel_dims: Tuple[int, int, int], voxel_spacing: float) -> np.ndarray: """Create mask by ray casting in X direction.""" intersector = RayMeshIntersector(mesh) - # Create a grid of rays in YZ plane, shooting in X direction + # Create a grid of rays in YZ plane (off the voxel centres, see _RAY_OFFSET), shooting in X direction grid_y, grid_z = np.mgrid[0 : voxel_dims[1], 0 : voxel_dims[2]] - ray_grid = np.vstack([-np.ones((grid_y.size,)), grid_y.ravel(), grid_z.ravel()]).T * voxel_spacing + ray_grid = ( + np.vstack( + [ + np.full(grid_y.size, _ray_start(mesh, voxel_spacing, 0)), + grid_y.ravel() + _RAY_OFFSET[0], + grid_z.ravel() + _RAY_OFFSET[1], + ], + ).T + * voxel_spacing + ) ray_dir = np.zeros((ray_grid.shape[0], 3)) ray_dir[:, 0] = 1 loc, _, _ = intersector.intersects_location(ray_grid, ray_dir) - # Convert to voxel coordinates and sort by x - int_loc = np.round(loc / voxel_spacing).astype("int") - sort_idx = int_loc[:, 0].argsort() - int_loc = int_loc[sort_idx, :] - # Build volume by tracking crossings - img = np.zeros((voxel_dims[2], voxel_dims[1]), dtype="bool") - vol = np.zeros((voxel_dims[2], voxel_dims[1], voxel_dims[0]), dtype="bool") - - for x in range(voxel_dims[0]): - idx = int_loc[:, 0] == x - img[int_loc[idx, 2], int_loc[idx, 1]] = np.logical_not(img[int_loc[idx, 2], int_loc[idx, 1]]) - vol[:, :, x] = img - - return vol + return _fill_from_crossings(_crossing_voxels(loc, voxel_spacing, 0), voxel_dims, 0) def mesh_to_volume(mesh: tm.Trimesh, voxel_dims: Tuple[int, int, int], voxel_spacing: float) -> np.ndarray: diff --git a/tests/test_seg2mesh_mesh2seg.py b/tests/test_seg2mesh_mesh2seg.py new file mode 100644 index 0000000..a73c58a --- /dev/null +++ b/tests/test_seg2mesh_mesh2seg.py @@ -0,0 +1,163 @@ +"""seg2mesh and mesh2seg: meshes are (x, y, z) in Å (voxel index i at i * spacing), closed, and convert back to the +voxels they came from.""" + +import copick +import numpy as np +import pytest +import trimesh as tm +from click.testing import CliRunner +from copick_utils.converters.mesh_from_segmentation import label_meshes, mesh_from_segmentation, volume_to_mesh +from copick_utils.converters.segmentation_from_mesh import mesh_to_volume + +VS = 10.0 + + +def box(shape_zyx, lo_xyz, hi_xyz, value=1, volume=None): + """A (z, y, x) volume with ``value`` in the voxels from ``lo_xyz`` to ``hi_xyz`` (inclusive).""" + vol = np.zeros(shape_zyx, np.uint16) if volume is None else volume + vol[lo_xyz[2] : hi_xyz[2] + 1, lo_xyz[1] : hi_xyz[1] + 1, lo_xyz[0] : hi_xyz[0] + 1] = value + return vol + + +def ball(shape_zyx, centre_xyz, radius): + z, y, x = np.indices(shape_zyx) + cx, cy, cz = centre_xyz + return ((x - cx) ** 2 + (y - cy) ** 2 + (z - cz) ** 2 < radius**2).astype(np.uint8) + + +def ring(shape_zyx): + """A thick ring around the z axis (a hole through it, so inside/outside alternates along rays).""" + z, y, x = np.indices(shape_zyx) + r = np.hypot(x - shape_zyx[2] / 2, y - shape_zyx[1] / 2) + return ((r > 3) & (r < 7) & (z > 1) & (z < shape_zyx[0] - 2)).astype(np.uint8) + + +# --- seg2mesh ---------------------------------------------------------------------------------------------------------- + + +def test_mesh_is_xyz_in_angstrom(): + # Voxels x 3-10, y 2-5, z 1-2 of a (z, y, x) = (6, 8, 14) volume: the surface lies half a voxel outside them. + mesh = volume_to_mesh(box((6, 8, 14), (3, 2, 1), (10, 5, 2)), VS) + np.testing.assert_allclose(mesh.bounds, [[2.5 * VS, 1.5 * VS, 0.5 * VS], [10.5 * VS, 5.5 * VS, 2.5 * VS]]) + + +def test_mesh_is_closed_and_faces_outward(): + mesh = volume_to_mesh(ball((14, 14, 14), (7, 7, 7), 5), VS) + assert mesh.is_watertight and mesh.is_winding_consistent + assert mesh.volume > 0 + outward = np.einsum("ij,ij->i", mesh.vertex_normals, mesh.vertices - 7 * VS) + assert (outward > 0).all() + + +def test_mesh_is_closed_where_the_segmentation_touches_the_volume_edge(): + mesh = volume_to_mesh(box((6, 6, 6), (0, 0, 0), (2, 5, 5)), VS) + assert mesh.is_watertight + np.testing.assert_allclose(mesh.bounds, [[-0.5 * VS] * 3, [2.5 * VS, 5.5 * VS, 5.5 * VS]]) + + +def test_label_meshes_one_per_label_in_place(): + vol = box((6, 8, 20), (1, 1, 1), (3, 3, 3), value=2) + box((6, 8, 20), (12, 2, 1), (15, 5, 4), value=7, volume=vol) + meshes = label_meshes(vol, VS) + assert list(meshes) == [2, 7] + np.testing.assert_allclose(meshes[7].bounds, [[11.5 * VS, 1.5 * VS, 0.5 * VS], [15.5 * VS, 5.5 * VS, 4.5 * VS]]) + assert all(m.is_watertight and m.volume > 0 for m in meshes.values()) + + +# --- mesh2seg ---------------------------------------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "volume", + [ + pytest.param(box((6, 8, 14), (3, 2, 1), (10, 5, 2)), id="box"), + pytest.param(ball((14, 16, 18), (9, 8, 7), 5), id="ball"), + pytest.param(ring((8, 16, 16)), id="ring"), + pytest.param(box((6, 6, 6), (0, 0, 0), (2, 5, 5)), id="at the volume edge"), + pytest.param(box((8, 10, 12), (2, 2, 4), (9, 7, 4)), id="one voxel thin"), + ], +) +def test_mesh_to_volume_gives_back_the_segmentation(volume): + # Marching-cubes meshes have their vertices and edges on the voxel grid, where the rays were cast. + mesh = volume_to_mesh(volume, VS) + np.testing.assert_array_equal(mesh_to_volume(mesh, volume.shape[::-1], VS), volume > 0) + + +def test_voxels_are_inside_when_their_centres_are(): + # x 22-61, y 31-79, z 47-93 Å: the voxel centres inside are x 30-60, y 40-70, z 50-90. + mesh = tm.creation.box(bounds=[[22, 31, 47], [61, 79, 93]]) + expected = np.zeros((12, 10, 10), bool) + expected[5:10, 4:8, 3:7] = True + np.testing.assert_array_equal(mesh_to_volume(mesh, (10, 10, 12), VS), expected) + + +def test_a_box_on_the_volume_edges_fills_the_volume(): + # Faces through the first voxel centres (0) and half a voxel past the last ones, like a valid-area box. + mesh = tm.creation.box(bounds=[[0, 0, 0], [100, 100, 120]]) + assert mesh_to_volume(mesh, (10, 10, 12), VS).all() + + +def test_a_mesh_reaching_beyond_the_volume(): + mesh = tm.creation.box(bounds=[[-55, -55, -55], [55, 45, 65]]) + expected = np.zeros((12, 10, 10), bool) + expected[:7, :5, :6] = True + np.testing.assert_array_equal(mesh_to_volume(mesh, (10, 10, 12), VS), expected) + + +# --- end to end ------------------------------------------------------------------------------------------------------ + + +def _invoke(command, args): + result = CliRunner().invoke(command, args, catch_exceptions=False) + assert result.exit_code == 0, result.output + return result + + +def test_seg2mesh_then_mesh2seg_round_trip(run, config_path): + from copick_utils.cli.mesh2seg import mesh2seg + from copick_utils.cli.seg2mesh import seg2mesh + + shape = (12, 16, 20) + volume = ball(shape, (12, 7, 5), 4) | box(shape, (0, 10, 6), (19, 15, 8)).astype(np.uint8) + run.new_voxel_spacing(VS, exist_ok=True) + run.get_voxel_spacing(VS).new_tomogram("wbp").from_numpy(np.zeros(shape, np.float32)) + run.new_segmentation(VS, "membrane", "1", is_multilabel=False, user_id="test").from_numpy(volume) + + _invoke(seg2mesh, ["-c", config_path, "-i", "membrane:test/1@10.0", "-o", "membrane:mesh/1", "-w", "1"]) + (mesh,) = run.get_meshes(object_name="membrane", user_id="mesh", session_id="1") + zyx = np.nonzero(volume) + expected = [[(zyx[a].min() - 0.5) * VS for a in (2, 1, 0)], [(zyx[a].max() + 0.5) * VS for a in (2, 1, 0)]] + np.testing.assert_allclose(mesh.mesh.bounds, expected) + + _invoke(mesh2seg, ["-c", config_path, "-i", "membrane:mesh/1", "-o", "membrane:back/1@10.0", "-w", "1"]) + # The CLI wrote through its own copick root: open the project again to see its output. + run = copick.from_file(config_path).get_run(run.name) + (back,) = run.get_segmentations(name="membrane", user_id="back", session_id="1", voxel_size=VS) + np.testing.assert_array_equal(back.numpy() > 0, volume > 0) + + +def test_seg2mesh_individual_meshes(run, config_path): + from copick_utils.cli.seg2mesh import seg2mesh + + vol = box((6, 8, 20), (1, 1, 1), (3, 3, 3), value=2) + box((6, 8, 20), (12, 2, 1), (15, 5, 4), value=7, volume=vol) + run.new_voxel_spacing(VS, exist_ok=True) + run.new_segmentation(VS, "membrane", "1", is_multilabel=False, user_id="test").from_numpy(vol) + + _invoke( + seg2mesh, + ["-c", config_path, "-i", "membrane:test/1@10.0", "-o", "membrane:mesh/inst-{instance_id}", "-im", "-w", "1"], + ) + meshes = {m.session_id: m for m in run.get_meshes(object_name="membrane", user_id="mesh")} + assert sorted(meshes) == ["inst-2", "inst-7"] + np.testing.assert_allclose(meshes["inst-2"].mesh.bounds, [[0.5 * VS] * 3, [3.5 * VS] * 3]) + + +def test_mesh_from_segmentation_individual_meshes_without_template(run): + vol = box((6, 8, 20), (1, 1, 1), (3, 3, 3), value=2) + run.new_voxel_spacing(VS, exist_ok=True) + seg = run.new_segmentation(VS, "membrane", "1", is_multilabel=False, user_id="test") + seg.from_numpy(vol) + mesh, stats = mesh_from_segmentation(seg, run, "membrane", "split", "mesh", individual_meshes=True) + assert mesh.session_id == "split-2" + assert stats["faces_created"] > 0