diff --git a/docs/api/workflows.rst b/docs/api/workflows.rst index 3ad644f..e397cbc 100644 --- a/docs/api/workflows.rst +++ b/docs/api/workflows.rst @@ -57,7 +57,7 @@ Convert Image to USD workflow = WorkflowConvertImageToUSD( input_filenames=["cardiac_4d.nrrd"], output_directory="./results", - project_name="patient_001", + usd_project_name="patient_001", segmentation_method=SegmentChestTotalSegmentatorWithContrast(), registration_method=RegisterImagesICON(), ) @@ -77,6 +77,7 @@ Image to VTK import itk from physiomotion4d import ( + ContourTools, SegmentChestTotalSegmentatorWithContrast, WorkflowConvertImageToVTK, ) @@ -85,12 +86,12 @@ Image to VTK workflow = WorkflowConvertImageToVTK( segmentation_method=SegmentChestTotalSegmentatorWithContrast() ) - result = workflow.run_workflow( + result = workflow.process( input_image=image, anatomy_groups=["heart", "major_vessels"], ) - WorkflowConvertImageToVTK.save_combined_surface( + ContourTools.save_combined_surface( result["surfaces"], "./output", prefix="patient01", diff --git a/docs/troubleshooting.rst b/docs/troubleshooting.rst index 48c004d..b75c798 100644 --- a/docs/troubleshooting.rst +++ b/docs/troubleshooting.rst @@ -75,7 +75,10 @@ Poor Segmentation Quality .. code-block:: python - from physiomotion4d import SegmentChestTotalSegmentatorWithContrast + from physiomotion4d import ( + SegmentChestTotalSegmentatorWithContrast, + WorkflowConvertImageToUSD, + ) workflow = WorkflowConvertImageToUSD( ..., diff --git a/pyproject.toml b/pyproject.toml index ec10fe6..3b73d9e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -78,6 +78,7 @@ dependencies = [ "pyvista[all]>=0.47.0", "usd-core>=23.11", "trimesh>=4.0.0", + "netgen-mesher>=6.2.2606", # Utilities "ipykernel>=6.0.0", @@ -228,6 +229,8 @@ module = [ "itk.*", "matplotlib", "matplotlib.*", + "netgen", + "netgen.*", "nibabel", "nibabel.*", "nrrd", diff --git a/src/physiomotion4d/cli/convert_image_to_vtk.py b/src/physiomotion4d/cli/convert_image_to_vtk.py index df64822..bf21ee0 100644 --- a/src/physiomotion4d/cli/convert_image_to_vtk.py +++ b/src/physiomotion4d/cli/convert_image_to_vtk.py @@ -2,7 +2,8 @@ """Command-line interface for the image-to-VTK segmentation workflow. Segments a 3D image using a chosen backend and writes per-anatomy-group VTP -surfaces and VTU voxel meshes annotated with anatomy labels and colors. +surfaces and VTU tetrahedral volume meshes annotated with anatomy labels and +colors. """ import argparse @@ -37,12 +38,12 @@ def main() -> int: Output files — combined mode (default) --------------------------------------- {prefix}_surfaces.vtp all surfaces merged into one file - {prefix}_meshes.vtu all voxel meshes merged into one file + {prefix}_meshes.vtu all tetrahedral volume meshes merged into one file Output files — split mode (--split-files) ------------------------------------------ {prefix}_{group}.vtp one surface per anatomy group - {prefix}_{group}.vtu one voxel mesh per anatomy group + {prefix}_{group}.vtu one tetrahedral volume mesh per anatomy group Examples -------- @@ -108,6 +109,26 @@ def main() -> int: "Choices: " + " ".join(ANATOMY_GROUPS) ), ) + parser.add_argument( + "--surface-target-reduction", + type=float, + default=0.0, + help=( + "Fraction in [0, 1) of surface triangles to remove via " + "decimate_pro (default: 0.0, no decimation)." + ), + ) + parser.add_argument( + "--mesh-target-reduction", + type=float, + default=0.0, + help=( + "Fraction in [0, 1) of triangles to remove from the surface " + "(via decimate_pro) before it is meshed into a tetrahedral " + "volume mesh by netgen; a coarser input surface yields a " + "coarser volume mesh (default: 0.0, no decimation)." + ), + ) # ── Output ──────────────────────────────────────────────────────────── parser.add_argument( @@ -156,16 +177,18 @@ def main() -> int: print("=" * 70) try: - from .. import WorkflowConvertImageToVTK + from .. import ContourTools, WorkflowConvertImageToVTK workflow = WorkflowConvertImageToVTK( segmentation_method=build_segmentation_method( args.segmentation_method, contrast=args.contrast ), ) - result = workflow.run_workflow( + result = workflow.process( input_image=input_image, anatomy_groups=args.anatomy_groups, + surface_target_reduction=args.surface_target_reduction, + mesh_target_reduction=args.mesh_target_reduction, ) except (ValueError, RuntimeError, OSError) as exc: print(f"Error during workflow: {exc}") @@ -189,13 +212,13 @@ def main() -> int: if args.split_files: # One file per anatomy group if surfaces: - saved_surfaces = WorkflowConvertImageToVTK.save_surfaces( + saved_surfaces = ContourTools.save_surfaces( surfaces, args.output_dir, prefix=prefix ) for group, path in saved_surfaces.items(): print(f" Surface [{group:15s}] -> {path}") if meshes: - saved_meshes = WorkflowConvertImageToVTK.save_meshes( + saved_meshes = ContourTools.save_meshes( meshes, args.output_dir, prefix=prefix ) for group, path in saved_meshes.items(): @@ -203,12 +226,12 @@ def main() -> int: else: # Combined single-file output if surfaces: - surface_file = WorkflowConvertImageToVTK.save_combined_surface( + surface_file = ContourTools.save_combined_surface( surfaces, args.output_dir, prefix=prefix ) print(f" Combined surface -> {surface_file}") if meshes: - mesh_file = WorkflowConvertImageToVTK.save_combined_mesh( + mesh_file = ContourTools.save_combined_mesh( meshes, args.output_dir, prefix=prefix ) print(f" Combined mesh -> {mesh_file}") diff --git a/src/physiomotion4d/contour_tools.py b/src/physiomotion4d/contour_tools.py index f1ee433..cd8219d 100644 --- a/src/physiomotion4d/contour_tools.py +++ b/src/physiomotion4d/contour_tools.py @@ -5,7 +5,8 @@ from __future__ import annotations import logging -from typing import cast +import os +from typing import Optional, cast import itk import numpy as np @@ -70,6 +71,84 @@ def extract_contours( return contours + def extract_mesh( + self, surface: pv.PolyData, mesh_target_reduction: float = 0.0 + ) -> Optional[pv.UnstructuredGrid]: + """Generate a tetrahedral volume mesh (VTU) from a closed surface via netgen. + + Optionally decimates the surface with + :meth:`pyvista.PolyDataFilters.decimate_pro` before meshing — netgen + has no post-hoc decimation of its own, so a coarser input surface is + the only way to get a coarser tetrahedral mesh out. + + Builds netgen's surface mesh directly from *surface*'s indexed + points/triangles rather than round-tripping through an STL file. + STL has no shared-vertex topology, so every triangle corner is + written independently; float32 rounding during that round-trip can + make two corners that were exactly the same point diverge by ~1e-6, + which is well within netgen's own "identical point" merge tolerance + and makes its STL reader spin forever trying to reconcile them + (observed hang on Taubin-smoothed surfaces, e.g. from + :meth:`extract_contours`). Adding points/triangles directly preserves + the exact shared-vertex indices already present in *surface*, so + there is nothing for netgen to reconcile. + + Args: + surface: Closed, triangulated surface for one anatomy group (as + returned by :meth:`extract_contours`). + mesh_target_reduction: Fraction in ``[0, 1)`` of surface triangles + to remove via ``decimate_pro(mesh_target_reduction, + preserve_topology=True)`` before meshing. ``0.0`` (default) + skips decimation and meshes the surface as given. + + Returns: + :class:`pyvista.UnstructuredGrid` of tetrahedral cells, or + ``None`` if the surface is empty or netgen produced no volume + mesh from it. + """ + if surface.n_points == 0: + return None + + meshing_surface = surface.triangulate().clean() + if mesh_target_reduction > 0.0: + meshing_surface = meshing_surface.decimate_pro( + mesh_target_reduction, preserve_topology=True + ) + + import netgen.meshing as ngm # noqa: PLC0415 + + triangles = meshing_surface.faces.reshape(-1, 4)[:, 1:4] + ngmesh = ngm.Mesh() + ngmesh.dim = 3 + point_ids = [ + ngmesh.Add(ngm.MeshPoint(ngm.Pnt(*point))) + for point in meshing_surface.points + ] + face_descriptor = ngmesh.Add(ngm.FaceDescriptor(surfnr=1, domin=1, domout=0)) + for triangle in triangles: + ngmesh.Add(ngm.Element2D(face_descriptor, [point_ids[i] for i in triangle])) + ngmesh.GenerateVolumeMesh() + + elements = ngmesh.Elements3D() + if len(elements) == 0: + self.log_warning("netgen produced no volume mesh for this surface") + return None + + points = ngmesh.Coordinates() + tets = np.array( + [[vertex.nr - 1 for vertex in element.vertices] for element in elements], + dtype=np.int64, + ) + # netgen's tet vertex order is opposite VTK_TETRA's right-hand + # convention (confirmed by negative cell volumes); swapping the last + # two indices restores positive-volume orientation. + tets = tets[:, [0, 1, 3, 2]] + cells = np.hstack( + [np.full((tets.shape[0], 1), 4, dtype=np.int64), tets] + ).flatten() + cell_types = np.full(tets.shape[0], pv.CellType.TETRA, dtype=np.uint8) + return pv.UnstructuredGrid(cells, cell_types, points) + def transform_contours( self, contours: pv.PolyData, @@ -462,3 +541,130 @@ def create_deformation_field( ) return deformation_field_img + + # ─────────────────────────── I/O helpers ─────────────────────────────── + + @staticmethod + def save_surfaces( + surfaces: dict[str, pv.PolyData], + output_dir: str, + prefix: str = "", + ) -> dict[str, str]: + """Save each named surface to its own VTP file. + + Args: + surfaces: Mapping of name → surface (e.g. the ``'surfaces'`` + value from :meth:`WorkflowConvertImageToVTK.process`). + output_dir: Directory to write files into (created if absent). + prefix: Optional filename prefix. Each file is named + ``{prefix}_{name}.vtp`` (or ``{name}.vtp`` when *prefix* is empty). + + Returns: + Mapping of name → absolute path of the saved file. + """ + os.makedirs(output_dir, exist_ok=True) + saved: dict[str, str] = {} + for name, surface in surfaces.items(): + stem = f"{prefix}_{name}" if prefix else name + path = os.path.join(output_dir, f"{stem}.vtp") + surface.save(path) + saved[name] = path + return saved + + @staticmethod + def save_meshes( + meshes: dict[str, pv.UnstructuredGrid], + output_dir: str, + prefix: str = "", + ) -> dict[str, str]: + """Save each named volume mesh to its own VTU file. + + Args: + meshes: Mapping of name → mesh (e.g. the ``'meshes'`` value from + :meth:`WorkflowConvertImageToVTK.process`). + output_dir: Directory to write files into (created if absent). + prefix: Optional filename prefix. Each file is named + ``{prefix}_{name}.vtu`` (or ``{name}.vtu`` when *prefix* is empty). + + Returns: + Mapping of name → absolute path of the saved file. + """ + os.makedirs(output_dir, exist_ok=True) + saved: dict[str, str] = {} + for name, mesh in meshes.items(): + stem = f"{prefix}_{name}" if prefix else name + path = os.path.join(output_dir, f"{stem}.vtu") + mesh.save(path) + saved[name] = path + return saved + + @staticmethod + def save_combined_surface( + surfaces: dict[str, pv.PolyData], + output_dir: str, + prefix: str = "", + ) -> str: + """Merge all named surfaces into a single VTP file. + + The merged mesh retains per-cell ``Color`` (RGBA uint8) from each + surface's annotation, enabling colour-by-anatomy rendering in + Paraview, PyVista, etc. Per-object ``field_data`` is not preserved + in the merged file. + + Args: + surfaces: Mapping of name → surface. + output_dir: Directory to write the file into (created if absent). + prefix: Optional filename prefix. Output is ``{prefix}_surfaces.vtp`` + (or ``surfaces.vtp`` when *prefix* is empty). + + Returns: + Absolute path to the saved VTP file. + + Raises: + ValueError: If *surfaces* is empty. + """ + if not surfaces: + raise ValueError("No surfaces to save.") + os.makedirs(output_dir, exist_ok=True) + stem = f"{prefix}_surfaces" if prefix else "surfaces" + output_file = os.path.join(output_dir, f"{stem}.vtp") + merged = cast( + pv.PolyData, pv.merge(list(surfaces.values()), merge_points=False) + ) + merged.save(output_file) + return output_file + + @staticmethod + def save_combined_mesh( + meshes: dict[str, pv.UnstructuredGrid], + output_dir: str, + prefix: str = "", + ) -> str: + """Merge all named volume meshes into a single VTU file. + + The merged mesh retains per-cell ``Color`` (RGBA uint8) from each + mesh's annotation. Per-object ``field_data`` is not preserved in the + merged file. + + Args: + meshes: Mapping of name → volume mesh. + output_dir: Directory to write the file into (created if absent). + prefix: Optional filename prefix. Output is ``{prefix}_meshes.vtu`` + (or ``meshes.vtu`` when *prefix* is empty). + + Returns: + Absolute path to the saved VTU file. + + Raises: + ValueError: If *meshes* is empty. + """ + if not meshes: + raise ValueError("No meshes to save.") + os.makedirs(output_dir, exist_ok=True) + stem = f"{prefix}_meshes" if prefix else "meshes" + output_file = os.path.join(output_dir, f"{stem}.vtu") + merged = cast( + pv.UnstructuredGrid, pv.merge(list(meshes.values()), merge_points=False) + ) + merged.save(output_file) + return output_file diff --git a/src/physiomotion4d/segment_chest_total_segmentator_with_contrast.py b/src/physiomotion4d/segment_chest_total_segmentator_with_contrast.py index 9d1190a..96d5710 100644 --- a/src/physiomotion4d/segment_chest_total_segmentator_with_contrast.py +++ b/src/physiomotion4d/segment_chest_total_segmentator_with_contrast.py @@ -37,7 +37,7 @@ class SegmentChestTotalSegmentatorWithContrast(SegmentChestTotalSegmentator): >>> contrast_mask = result['contrast'] """ - def __init__(self, log_level: int | str = logging.INFO): + def __init__(self, log_level: int | str = logging.INFO) -> None: """Initialize the contrast-enhanced TotalSegmentator-based segmentation. Args: @@ -52,32 +52,32 @@ def _add_extra_taxonomy_groups(self) -> None: self.taxonomy.add_organ("contrast", 135, "contrast") def postprocess_after_labelmap( - self, input_image: itk.image, labelmap_image: itk.image - ) -> itk.image: + self, input_image: itk.Image, labelmap_image: itk.Image + ) -> itk.Image: """Run contrast-enhanced blood detection on the labelmap. Overrides :meth:`SegmentAnatomyBase.postprocess_after_labelmap`. Args: - input_image (itk.image): The original, unpreprocessed input image - labelmap_image (itk.image): The postprocessed segmentation labelmap + input_image (itk.Image): The original, unpreprocessed input image + labelmap_image (itk.Image): The postprocessed segmentation labelmap Returns: - itk.image: The labelmap, with contrast-enhanced regions labeled + itk.Image: The labelmap, with contrast-enhanced regions labeled """ return self.segment_contrast_agent(input_image, labelmap_image) def segment_connected_component( self, - preprocessed_image: itk.image, - labelmap_image: itk.image, + preprocessed_image: itk.Image, + labelmap_image: itk.Image, lower_threshold: int, upper_threshold: int, labelmap_ids: Optional[list[int]] = None, mask_id: int = 0, use_mid_slice: bool = True, hole_fill: int = 2, - ) -> itk.image: + ) -> itk.Image: """ Segment connected components based on intensity thresholding. @@ -87,8 +87,8 @@ def segment_connected_component( tissue types. Args: - preprocessed_image (itk.image): The preprocessed input image - labelmap_image (itk.image): Existing labelmap to constrain search + preprocessed_image (itk.Image): The preprocessed input image + labelmap_image (itk.Image): Existing labelmap to constrain search lower_threshold (int): Lower intensity threshold upper_threshold (int): Upper intensity threshold labelmap_ids (Optional[list[int]]): List of label IDs to search within. @@ -99,7 +99,7 @@ def segment_connected_component( hole_fill (int): Number of pixels to dilate/erode for hole filling Returns: - itk.image: Updated labelmap with new component labeled as mask_id + itk.Image: Updated labelmap with new component labeled as mask_id Example: >>> # Segment contrast-enhanced blood @@ -178,8 +178,8 @@ def segment_connected_component( return results_image def segment_contrast_agent( - self, preprocessed_image: itk.image, labelmap_image: itk.image - ) -> itk.image: + self, preprocessed_image: itk.Image, labelmap_image: itk.Image + ) -> itk.Image: """ Include contrast-enhanced blood in the labelmap. @@ -188,11 +188,11 @@ def segment_contrast_agent( focused on the middle slice where the heart is typically located. Args: - preprocessed_image (itk.image): The preprocessed CT image - labelmap_image (itk.image): Existing segmentation labelmap + preprocessed_image (itk.Image): The preprocessed CT image + labelmap_image (itk.Image): Existing segmentation labelmap Returns: - itk.image: Updated labelmap with contrast-enhanced regions labeled + itk.Image: Updated labelmap with contrast-enhanced regions labeled Note: Assumes the mid-z slice of the data contains the heart. diff --git a/src/physiomotion4d/workflow_convert_image_to_vtk.py b/src/physiomotion4d/workflow_convert_image_to_vtk.py index 4e240b5..0fcc739 100644 --- a/src/physiomotion4d/workflow_convert_image_to_vtk.py +++ b/src/physiomotion4d/workflow_convert_image_to_vtk.py @@ -1,9 +1,10 @@ """Workflow for segmenting a CT image and converting anatomy groups to VTK surfaces and meshes. The workflow segments a 3D CT image using a chosen backend, then extracts one VTP -(surface) and one VTU (voxel mesh) per non-empty anatomy group. Each output object -carries anatomy metadata and solid color from :class:`USDAnatomyTools` as field and -cell data so that downstream tools (PyVista, Paraview, USD pipeline) can use them +(surface) and one VTU (tetrahedral volume mesh, generated by netgen from that +surface) per non-empty anatomy group. Each output object carries anatomy +metadata and solid color from :class:`USDAnatomyTools` as field and cell data +so that downstream tools (PyVista, Paraview, USD pipeline) can use them directly. Typical usage:: @@ -17,20 +18,21 @@ ct = itk.imread('chest_ct.nii.gz') segmenter = SegmentChestTotalSegmentatorWithContrast() workflow = WorkflowConvertImageToVTK(segmentation_method=segmenter) - result = workflow.run_workflow(ct) + result = workflow.process( + ct, surface_target_reduction=0.5, mesh_target_reduction=0.7 + ) # Combined single-file output (default) - WorkflowConvertImageToVTK.save_combined_surface(result['surfaces'], './out', prefix='patient') - WorkflowConvertImageToVTK.save_combined_mesh(result['meshes'], './out', prefix='patient') + ContourTools.save_combined_surface(result['surfaces'], './out', prefix='patient') + ContourTools.save_combined_mesh(result['meshes'], './out', prefix='patient') # Per-group split output - WorkflowConvertImageToVTK.save_surfaces(result['surfaces'], './out', prefix='patient') - WorkflowConvertImageToVTK.save_meshes(result['meshes'], './out', prefix='patient') + ContourTools.save_surfaces(result['surfaces'], './out', prefix='patient') + ContourTools.save_meshes(result['meshes'], './out', prefix='patient') """ import logging -import os -from typing import Any, Optional, cast +from typing import Any, Optional import itk import numpy as np @@ -65,7 +67,7 @@ class WorkflowConvertImageToVTK(PhysioMotion4DBase): **VTK object annotation** Each :class:`pyvista.PolyData` surface and :class:`pyvista.UnstructuredGrid` mesh - returned by :meth:`run_workflow` carries: + returned by :meth:`process` carries: - ``field_data['AnatomyGroup']`` — anatomy group name, e.g. ``'heart'``. - ``field_data['SegmentationLabelNames']`` — individual structure names within the @@ -76,10 +78,12 @@ class WorkflowConvertImageToVTK(PhysioMotion4DBase): **I/O contract** - :meth:`run_workflow` performs *no* file I/O. Use the static helpers - :meth:`save_surfaces`, :meth:`save_meshes`, :meth:`save_combined_surface`, and - :meth:`save_combined_mesh` — or the CLI ``physiomotion4d-convert-image-to-vtk`` — to - write results to disk. + :meth:`process` performs *no* file I/O. Use + :class:`ContourTools`'s static helpers + :meth:`ContourTools.save_surfaces`, :meth:`ContourTools.save_meshes`, + :meth:`ContourTools.save_combined_surface`, and + :meth:`ContourTools.save_combined_mesh` — or the CLI + ``physiomotion4d-convert-image-to-vtk`` — to write results to disk. """ def __init__( @@ -183,38 +187,34 @@ def _extract_surface(self, mask_image: Any) -> Optional[pv.PolyData]: return None return self._contour_tools.extract_contours(mask_image) - def _extract_mesh(self, mask_image: Any) -> Optional[pv.UnstructuredGrid]: - """Extract a voxel-based volumetric mesh (VTU) from a binary mask image. + def _extract_mesh( + self, surface: pv.PolyData, mesh_target_reduction: float = 0.0 + ) -> Optional[pv.UnstructuredGrid]: + """Generate a tetrahedral volume mesh (VTU) from a closed surface. + + Delegates to :meth:`ContourTools.extract_mesh`. - Wraps the ITK image as a VTK ImageData and thresholds at 0.5 to obtain - hexahedral voxel cells for non-zero voxels. + Args: + surface: Closed, triangulated surface for one anatomy group (as + returned by :meth:`_extract_surface`). + mesh_target_reduction: Fraction in ``[0, 1)`` of surface triangles + to remove before meshing. ``0.0`` (default) skips decimation. Returns: - :class:`pyvista.UnstructuredGrid` of labeled voxels, or ``None`` if empty. + :class:`pyvista.UnstructuredGrid` of tetrahedral cells, or + ``None`` if the surface is empty or no volume mesh could be + generated from it. """ - arr = itk.GetArrayFromImage(mask_image) - if int(arr.sum()) == 0: - return None - - vtk_image = pv.wrap(itk.vtk_image_from_image(mask_image)) - if not isinstance(vtk_image, pv.ImageData): - self.log_warning( - "Expected pv.ImageData from vtk_image_from_image, got %s — skipping mesh", - type(vtk_image).__name__, - ) - return None - - thresholded = vtk_image.threshold(0.5) - if isinstance(thresholded, pv.UnstructuredGrid): - return thresholded - return cast(pv.UnstructuredGrid, thresholded.cast_to_unstructured_grid()) + return self._contour_tools.extract_mesh(surface, mesh_target_reduction) # ─────────────────────────── Main workflow ───────────────────────────── - def run_workflow( + def process( self, input_image: Any, anatomy_groups: Optional[list[str]] = None, + surface_target_reduction: float = 0.0, + mesh_target_reduction: float = 0.0, ) -> dict[str, Any]: """Segment the CT image and extract per-anatomy-group VTK objects. @@ -224,12 +224,22 @@ def run_workflow( processes all non-empty groups. Valid names are given by :attr:`ANATOMY_GROUPS`, derived from the active segmenter's taxonomy. + surface_target_reduction: Fraction in ``[0, 1)`` of surface + triangles to remove via ``decimate_pro(surface_target_reduction, + preserve_topology=True)``. ``0.0`` (default) skips decimation. + mesh_target_reduction: Fraction in ``[0, 1)`` of the surface's + triangles to remove (via the same ``decimate_pro`` call) + before it is handed to netgen for tetrahedral meshing. A + coarser input surface yields a coarser volume mesh — netgen + itself has no post-hoc decimation. ``0.0`` (default) skips + decimation and meshes the surface as extracted. Returns: ``dict`` with the following keys: - ``'surfaces'`` — ``dict[str, pv.PolyData]``: smoothed surface per group. - - ``'meshes'`` — ``dict[str, pv.UnstructuredGrid]``: voxel mesh per group. + - ``'meshes'`` — ``dict[str, pv.UnstructuredGrid]``: tetrahedral volume + mesh per group, generated by netgen from that group's surface. - ``'labelmap'`` — ``itk.Image``: detailed per-structure segmentation labelmap from the segmenter. - ``'segmentation_masks'`` — ``dict[str, itk.Image]``: per-group binary @@ -283,13 +293,20 @@ def run_workflow( color = self._anatomy_color_map.get(group, (0.7, 0.7, 0.7)) self.log_info(" Extracting surface for: %s", group) - surface = self._extract_surface(mask_image) - if surface is not None: - self._annotate(surface, group, label_names, label_ids, color) - surfaces[group] = surface + base_surface = self._extract_surface(mask_image) + if base_surface is None: + continue + + export_surface = base_surface + if surface_target_reduction > 0.0: + export_surface = export_surface.decimate_pro( + surface_target_reduction, preserve_topology=True + ) + self._annotate(export_surface, group, label_names, label_ids, color) + surfaces[group] = export_surface - self.log_info(" Extracting voxel mesh for: %s", group) - mesh = self._extract_mesh(mask_image) + self.log_info(" Extracting volume mesh for: %s", group) + mesh = self._extract_mesh(base_surface, mesh_target_reduction) if mesh is not None: self._annotate(mesh, group, label_names, label_ids, color) meshes[group] = mesh @@ -304,127 +321,3 @@ def run_workflow( "labelmap": seg_result["labelmap"], "segmentation_masks": seg_masks, } - - # ─────────────────────────── I/O helpers ─────────────────────────────── - - @staticmethod - def save_surfaces( - surfaces: dict[str, pv.PolyData], - output_dir: str, - prefix: str = "", - ) -> dict[str, str]: - """Save each group surface to its own VTP file. - - Args: - surfaces: Mapping of anatomy group name → surface (from - :meth:`run_workflow`). - output_dir: Directory to write files into (created if absent). - prefix: Optional filename prefix. Each file is named - ``{prefix}_{group}.vtp`` (or ``{group}.vtp`` when *prefix* is empty). - - Returns: - Mapping of anatomy group name → absolute path of the saved file. - """ - os.makedirs(output_dir, exist_ok=True) - saved: dict[str, str] = {} - for name, surface in surfaces.items(): - stem = f"{prefix}_{name}" if prefix else name - path = os.path.join(output_dir, f"{stem}.vtp") - surface.save(path) - saved[name] = path - return saved - - @staticmethod - def save_meshes( - meshes: dict[str, pv.UnstructuredGrid], - output_dir: str, - prefix: str = "", - ) -> dict[str, str]: - """Save each group voxel mesh to its own VTU file. - - Args: - meshes: Mapping of anatomy group name → mesh (from :meth:`run_workflow`). - output_dir: Directory to write files into (created if absent). - prefix: Optional filename prefix. Each file is named - ``{prefix}_{group}.vtu`` (or ``{group}.vtu`` when *prefix* is empty). - - Returns: - Mapping of anatomy group name → absolute path of the saved file. - """ - os.makedirs(output_dir, exist_ok=True) - saved: dict[str, str] = {} - for name, mesh in meshes.items(): - stem = f"{prefix}_{name}" if prefix else name - path = os.path.join(output_dir, f"{stem}.vtu") - mesh.save(path) - saved[name] = path - return saved - - @staticmethod - def save_combined_surface( - surfaces: dict[str, pv.PolyData], - output_dir: str, - prefix: str = "", - ) -> str: - """Merge all group surfaces into a single VTP file. - - The merged mesh retains per-cell ``Color`` (RGBA uint8) from each group's - annotation, enabling colour-by-anatomy rendering in Paraview, PyVista, etc. - Per-object ``field_data`` is not preserved in the merged file. - - Args: - surfaces: Mapping of anatomy group name → surface. - output_dir: Directory to write the file into (created if absent). - prefix: Optional filename prefix. Output is ``{prefix}_surfaces.vtp`` - (or ``surfaces.vtp`` when *prefix* is empty). - - Returns: - Absolute path to the saved VTP file. - - Raises: - ValueError: If *surfaces* is empty. - """ - if not surfaces: - raise ValueError("No surfaces to save.") - os.makedirs(output_dir, exist_ok=True) - stem = f"{prefix}_surfaces" if prefix else "surfaces" - output_file = os.path.join(output_dir, f"{stem}.vtp") - merged = cast( - pv.PolyData, pv.merge(list(surfaces.values()), merge_points=False) - ) - merged.save(output_file) - return output_file - - @staticmethod - def save_combined_mesh( - meshes: dict[str, pv.UnstructuredGrid], - output_dir: str, - prefix: str = "", - ) -> str: - """Merge all group meshes into a single VTU file. - - The merged mesh retains per-cell ``Color`` (RGBA uint8) from each group's - annotation. Per-object ``field_data`` is not preserved in the merged file. - - Args: - meshes: Mapping of anatomy group name → voxel mesh. - output_dir: Directory to write the file into (created if absent). - prefix: Optional filename prefix. Output is ``{prefix}_meshes.vtu`` - (or ``meshes.vtu`` when *prefix* is empty). - - Returns: - Absolute path to the saved VTU file. - - Raises: - ValueError: If *meshes* is empty. - """ - if not meshes: - raise ValueError("No meshes to save.") - os.makedirs(output_dir, exist_ok=True) - stem = f"{prefix}_meshes" if prefix else "meshes" - output_file = os.path.join(output_dir, f"{stem}.vtu") - merged = cast( - pv.UnstructuredGrid, pv.merge(list(meshes.values()), merge_points=False) - ) - merged.save(output_file) - return output_file diff --git a/src/physiomotion4d/workflow_fit_statistical_model_to_patient.py b/src/physiomotion4d/workflow_fit_statistical_model_to_patient.py index 07de1b8..bbedfc5 100644 --- a/src/physiomotion4d/workflow_fit_statistical_model_to_patient.py +++ b/src/physiomotion4d/workflow_fit_statistical_model_to_patient.py @@ -190,7 +190,7 @@ def __init__( segmentation_method=segmentation_method, log_level=log_level, ) - patient_models_data = convert_image_to_vtk.run_workflow( + patient_models_data = convert_image_to_vtk.process( input_image=patient_image, anatomy_groups=["heart"], ) diff --git a/tests/test_workflow_fit_statistical_model_to_patient.py b/tests/test_workflow_fit_statistical_model_to_patient.py index 5cf28ce..9842f19 100644 --- a/tests/test_workflow_fit_statistical_model_to_patient.py +++ b/tests/test_workflow_fit_statistical_model_to_patient.py @@ -100,8 +100,8 @@ class _FakeConvertImageToVTK: def __init__(self, **kwargs: Any) -> None: captured["init_kwargs"] = kwargs - def run_workflow(self, **kwargs: Any) -> dict[str, Any]: - captured["run_kwargs"] = kwargs + def process(self, **kwargs: Any) -> dict[str, Any]: + captured["process_kwargs"] = kwargs return {"meshes": {"heart": heart_mesh}} monkeypatch.setattr( @@ -119,7 +119,7 @@ def run_workflow(self, **kwargs: Any) -> dict[str, Any]: captured["init_kwargs"]["segmentation_method"], SegmentHeartSimplewareTrimmedBranches, ) - assert captured["run_kwargs"]["anatomy_groups"] == ["heart"] + assert captured["process_kwargs"]["anatomy_groups"] == ["heart"] assert workflow.patient_models == [heart_mesh] diff --git a/tutorials/tutorial_02_ct_to_vtk.py b/tutorials/tutorial_02_ct_to_vtk.py index 895edde..6226942 100644 --- a/tutorials/tutorial_02_ct_to_vtk.py +++ b/tutorials/tutorial_02_ct_to_vtk.py @@ -4,8 +4,8 @@ Purpose ------- Segment one 3D CT frame into anatomical groups and save combined VTK surface -and voxel mesh files. The output can be inspected directly in PyVista or used -as input for Tutorial 5. +and tetrahedral volume mesh files. The output can be inspected directly in +PyVista or used as input for Tutorial 5. Data Required ------------- @@ -23,6 +23,7 @@ import itk import pyvista as pv +from physiomotion4d.contour_tools import ContourTools from physiomotion4d.segment_chest_total_segmentator_with_contrast import ( SegmentChestTotalSegmentatorWithContrast, ) @@ -79,19 +80,27 @@ # %% # Workflow execution - result = workflow.run_workflow(input_image=ct_image) + # + # surface_target_reduction decimates each exported VTP surface; + # mesh_target_reduction decimates the surface netgen tetrahedralizes into + # each VTU volume mesh (a coarser input surface yields a coarser mesh). + result = workflow.process( + input_image=ct_image, + surface_target_reduction=0.5, + mesh_target_reduction=0.7, + ) # %% # Result saving surface_file = Path( - WorkflowConvertImageToVTK.save_combined_surface( + ContourTools.save_combined_surface( result["surfaces"], str(output_dir), prefix="patient", ) ) mesh_file = Path( - WorkflowConvertImageToVTK.save_combined_mesh( + ContourTools.save_combined_mesh( result["meshes"], str(output_dir), prefix="patient",