diff --git a/.github/scripts/build_dashboard.py b/.github/scripts/build_dashboard.py
index 601e6d2..3ec6478 100644
--- a/.github/scripts/build_dashboard.py
+++ b/.github/scripts/build_dashboard.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python3
 """
 Build the nightly health dashboard for PhysioTwin4D.
 
diff --git a/README.md b/README.md
index da68eaf..8a1ff9c 100644
--- a/README.md
+++ b/README.md
@@ -25,7 +25,7 @@ for users and contributors. Key sections:
 
 - **Personalized digital twins**: build subject-specific anatomic models and physiological AI surrogates from 3D/4D medical images
 - **Statistical shape models**: capture subject-specific anatomy and establish cross-subject correspondence, aiding AI surrogate generalization and simplifying traditional solver setup
-- **Simplified workflows on industry-leading open-source tools**: ICON and Greedy for registration; MONAI with TotalSegmentator and Simpleware for segmentation; Netgen for meshing (LGPL license); scikit-learn for statistical shape modeling; ITK for image processing; PyVista and OpenUSD/Omniverse for geometry manipulation; CuPy for accelerated computing; and PhysicsNeMo for AI surrogates
+- **Simplified workflows on industry-leading open-source tools**: ICON and Greedy for registration; MONAI with TotalSegmentator and Simpleware for segmentation; scikit-learn for statistical shape modeling; ITK for image processing; PyVista and OpenUSD/Omniverse for geometry manipulation; CuPy for accelerated computing; and PhysicsNeMo for AI surrogates
 - **Extensible class hierarchy**: add new segmentation and registration methods, and extend to new data types, organs, and physiological processes, without reworking the workflow layer
 - **Physiological motion**: cardiac and respiratory motion today, expanding to electrophysiology, blood flow, and organ perfusion
 - **NVIDIA Omniverse as the simulation hub**: the end goal for simulation — a simulation-information hub and gateway to other engines (e.g., Ansys solvers), interactive simulations for treatment planning (e.g., Isaac Sim, Newton), visualization systems (e.g., AR/VR devices), and physical systems (e.g., robots via ROS)
@@ -71,9 +71,8 @@ See the [contributing guide](https://project-monai.github.io/physiotwin4d/contri
 
 This project is licensed under the Apache 2.0 License - see the LICENSE file for details.
 
-Note that some optional dependencies carry their own license terms: meshing
-via Netgen is LGPL-licensed, and NVIDIA Omniverse is distributed under its
-own customer license, which is free for academic and commercial use. To our
-knowledge, none of the dependencies carry commercially prohibitive licenses
-such as GPL, but no guarantees are provided — review the license of each
-dependency for your own use case.
+Note that some optional dependencies carry their own license terms: NVIDIA
+Omniverse is distributed under its own customer license, which is free for
+academic and commercial use. To our knowledge, none of the dependencies
+carry commercially prohibitive licenses such as GPL, but no guarantees are
+provided — review the license of each dependency for your own use case.
diff --git a/data/DirLab-4DCT/Convert4DCTToMHD.py b/data/DirLab-4DCT/Convert4DCTToMHD.py
index a6c77ee..025cfc5 100644
--- a/data/DirLab-4DCT/Convert4DCTToMHD.py
+++ b/data/DirLab-4DCT/Convert4DCTToMHD.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %%
 import os
 
diff --git a/data/Slicer-Heart-CT/README.md b/data/Slicer-Heart-CT/README.md
index e184120..f65b569 100644
--- a/data/Slicer-Heart-CT/README.md
+++ b/data/Slicer-Heart-CT/README.md
@@ -18,14 +18,17 @@ assert DataDownloadTools.VerifySlicerHeartCTData("data/Slicer-Heart-CT")
 ```
 
 This fetches a single ~1.2 GB file from
-[github.com/SlicerHeart/SlicerHeart](https://github.com/SlicerHeart/SlicerHeart/releases/download/TestingData/TruncalValve_4DCT.seq.nrrd).
-An existing non-empty file is reused, so re-running the command resumes an
-interrupted download.
+[github.com/SlicerHeart/SlicerHeart](https://github.com/SlicerHeart/SlicerHeart/releases/download/TestingData/TruncalValve_4DCT.seq.nrrd),
+then splits it into per-phase 3D `slice_???.mha` volumes in the same
+directory via `ConvertImage4DTo3D`. An existing non-empty `.seq.nrrd` is
+reused, and the split is skipped once the `slice_???.mha` files are present —
+so re-running the command resumes an interrupted download or conversion.
 
 **Directory structure after download:**
-```
+```text
 data/Slicer-Heart-CT/
 ├── TruncalValve_4DCT.seq.nrrd
+├── slice_000.mha ... slice_020.mha
 └── README.md (this file)
 ```
 
@@ -54,11 +57,14 @@ Philadelphia):
 - Primary dataset for tutorials
 - Primary dataset for `experiments/Heart-GatedCT_To_USD/` and
   `experiments/Heart-VTKSeries_To_USD/`, whose
-  `0-download_and_convert_4d_to_3d.py` scripts download this sequence and
-  split it into per-phase 3D `.mha` slices via `ConvertImage4DTo3D`
+  `0-download_and_convert_4d_to_3d.py` scripts call
+  `DataDownloadTools.DownloadSlicerHeartCTData`, which downloads this
+  sequence and splits it into per-phase 3D `.mha` slices
 - Used in the test suite (`tests/test_download_heart_data.py`)
 - Example data for cardiac motion visualization in NVIDIA Omniverse
 
 ### Files in This Directory
 
 - `TruncalValve_4DCT.seq.nrrd` — the downloaded 4D CT sequence
+- `slice_000.mha` ... `slice_020.mha` — the 21 cardiac phases split out of
+  the sequence above
diff --git a/docs/api/workflows.rst b/docs/api/workflows.rst
index 477b7b0..4798d43 100644
--- a/docs/api/workflows.rst
+++ b/docs/api/workflows.rst
@@ -27,7 +27,7 @@ Available Workflows
    * - :class:`WorkflowConvertImageToUSD`
      - Convert a 4D cardiac CT sequence into animated USD anatomy.
    * - :class:`WorkflowConvertImageToVTK`
-     - Segment one CT image and export anatomy-group VTK surfaces and meshes.
+     - Segment one CT image and export anatomy-group VTK surfaces.
    * - :class:`WorkflowConvertVTKToUSD`
      - Convert VTK/VTP/VTU meshes or time series into USD.
    * - :class:`WorkflowCreateStatisticalModel`
diff --git a/docs/architecture.rst b/docs/architecture.rst
index 1ccf967..92065dc 100644
--- a/docs/architecture.rst
+++ b/docs/architecture.rst
@@ -61,8 +61,7 @@ Primary Workflows
    contours and painted animated USD files.
 
 ``WorkflowConvertImageToVTK``
-   Segments a 3D CT image and exports anatomy groups as VTK surfaces and voxel
-   meshes.
+   Segments a 3D CT image and exports anatomy groups as VTK surfaces.
 
 ``WorkflowCreateStatisticalModel``
    Aligns a population of meshes to a reference and builds a PCA statistical
diff --git a/docs/cli_scripts/download_data.rst b/docs/cli_scripts/download_data.rst
index e759ca9..16660a8 100644
--- a/docs/cli_scripts/download_data.rst
+++ b/docs/cli_scripts/download_data.rst
@@ -54,15 +54,18 @@ Options
 Output
 ======
 
-For ``Slicer-Heart-CT``, the command downloads or reuses:
+For ``Slicer-Heart-CT``, the command downloads the 4-D sequence and splits it
+into per-phase 3-D volumes:
 
 .. code-block:: text
 
    data/Slicer-Heart-CT/TruncalValve_4DCT.seq.nrrd
+   data/Slicer-Heart-CT/slice_000.mha ... slice_020.mha
 
 The command uses
 :meth:`physiotwin4d.data_download_tools.DataDownloadTools.DownloadSlicerHeartCTData`,
-so repeated runs reuse the existing non-empty file.
+so repeated runs reuse the existing non-empty file and skip the split once
+the ``slice_???.mha`` files are present.
 
 For ``KCL-Heart-Model``, the command downloads, extracts, and reuses:
 
diff --git a/docs/cli_scripts/overview.rst b/docs/cli_scripts/overview.rst
index ab098ad..238d2ff 100644
--- a/docs/cli_scripts/overview.rst
+++ b/docs/cli_scripts/overview.rst
@@ -54,7 +54,7 @@ Current Scripts
    * - :doc:`heart_gated_ct`
      - Process cardiac gated CT to animated heart models with physiological motion
    * - ``physiotwin4d-convert-image-to-vtk``
-     - Segment one 3D image and export anatomy-group VTK surfaces and meshes
+     - Segment one 3D image and export anatomy-group VTK surfaces
    * - ``physiotwin4d-convert-image-4d-to-3d``
      - Split a 4D medical image into a 3D time series using ITK readers
    * - :doc:`create_statistical_model`
diff --git a/docs/developer/architecture.rst b/docs/developer/architecture.rst
index d1c0bc4..8cd6885 100644
--- a/docs/developer/architecture.rst
+++ b/docs/developer/architecture.rst
@@ -40,7 +40,7 @@ Workflow Classes
    transformation, and animated USD generation.
 
 ``WorkflowConvertImageToVTK``
-   Converts one 3D image into labeled VTK surface and voxel-mesh outputs.
+   Converts one 3D image into labeled VTK surface outputs.
 
 ``WorkflowCreateStatisticalModel``
    Builds a PCA statistical shape model from aligned population meshes.
diff --git a/docs/tutorials.rst b/docs/tutorials.rst
index df036f0..dae135e 100644
--- a/docs/tutorials.rst
+++ b/docs/tutorials.rst
@@ -201,7 +201,7 @@ Preview
 
 Inner API usage
    The workflow owns a segmentation method and turns each anatomy group into
-   decimated VTK surfaces and volume meshes:
+   a decimated VTK surface:
 
    .. code-block:: python
 
@@ -213,11 +213,17 @@ Inner API usage
       result = workflow.process(
           input_image=ct_image,
           surface_target_reduction=0.5,
-          mesh_target_reduction=0.7,
       )
 
    Use ``SegmentChestTotalSegmentator`` instead for non-contrast studies.
 
+   By default the script saves one combined ``patient_surfaces.vtp``.  Set
+   ``SAVE_GROUP_SURFACES = True`` and/or ``SAVE_LABEL_SURFACES = True`` near
+   the top of the script to additionally save one VTP per anatomy group
+   (e.g. ``patient_heart.vtp``) and/or one VTP per individual anatomical
+   structure (e.g. ``patient_left_ventricle.vtp``); the latter passes
+   ``extract_label_surfaces=True`` to :meth:`WorkflowConvertImageToVTK.process`.
+
 Run
    .. code-block:: bash
 
diff --git a/experiments/Colormap-VTK_To_USD/colormap_vtk_to_usd.py b/experiments/Colormap-VTK_To_USD/colormap_vtk_to_usd.py
index 8098be6..22f1422 100644
--- a/experiments/Colormap-VTK_To_USD/colormap_vtk_to_usd.py
+++ b/experiments/Colormap-VTK_To_USD/colormap_vtk_to_usd.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %% [markdown]
 # # Colormap Features for VTK to USD Conversion
 #
diff --git a/experiments/Convert_VTK_To_USD/convert_chop_alterra_valve_to_usd.py b/experiments/Convert_VTK_To_USD/convert_chop_alterra_valve_to_usd.py
index dd033fc..f66b87d 100644
--- a/experiments/Convert_VTK_To_USD/convert_chop_alterra_valve_to_usd.py
+++ b/experiments/Convert_VTK_To_USD/convert_chop_alterra_valve_to_usd.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %% [markdown]
 # # Cardiac Valve 4D Time-Series Conversion to USD
 #
diff --git a/experiments/Convert_VTK_To_USD/convert_chop_heart_vtk_to_usd.py b/experiments/Convert_VTK_To_USD/convert_chop_heart_vtk_to_usd.py
index 5bf0b4a..2a7f6f2 100644
--- a/experiments/Convert_VTK_To_USD/convert_chop_heart_vtk_to_usd.py
+++ b/experiments/Convert_VTK_To_USD/convert_chop_heart_vtk_to_usd.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %%
 from pathlib import Path
 
diff --git a/experiments/Convert_VTK_To_USD/convert_chop_tpv25_valve_to_usd.py b/experiments/Convert_VTK_To_USD/convert_chop_tpv25_valve_to_usd.py
index 5b006ae..ccd8dff 100644
--- a/experiments/Convert_VTK_To_USD/convert_chop_tpv25_valve_to_usd.py
+++ b/experiments/Convert_VTK_To_USD/convert_chop_tpv25_valve_to_usd.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %% [markdown]
 # # Cardiac Valve 4D Time-Series Conversion to USD
 #
diff --git a/experiments/Convert_VTK_To_USD/convert_vtk_to_usd_using_class.py b/experiments/Convert_VTK_To_USD/convert_vtk_to_usd_using_class.py
index 0b47e18..44c4463 100644
--- a/experiments/Convert_VTK_To_USD/convert_vtk_to_usd_using_class.py
+++ b/experiments/Convert_VTK_To_USD/convert_vtk_to_usd_using_class.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %% [markdown]
 # # VTK to USD Converter Test Notebook
 #
diff --git a/experiments/DisplacementField_To_USD/displacement_field_to_usd.py b/experiments/DisplacementField_To_USD/displacement_field_to_usd.py
index e9c5f1f..b7b053f 100644
--- a/experiments/DisplacementField_To_USD/displacement_field_to_usd.py
+++ b/experiments/DisplacementField_To_USD/displacement_field_to_usd.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %% [markdown]
 # # Displacement Field to USD for Omniverse Visualization
 #
diff --git a/experiments/Heart-Create_Statistical_Model/1-input_meshes_to_input_surfaces.py b/experiments/Heart-Create_Statistical_Model/1-input_meshes_to_input_surfaces.py
index 3762872..c98fba5 100644
--- a/experiments/Heart-Create_Statistical_Model/1-input_meshes_to_input_surfaces.py
+++ b/experiments/Heart-Create_Statistical_Model/1-input_meshes_to_input_surfaces.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %% [markdown]
 # # Convert VTK to VTP Surface Files
 #
diff --git a/experiments/Heart-Create_Statistical_Model/2-input_surfaces_to_surfaces_aligned.py b/experiments/Heart-Create_Statistical_Model/2-input_surfaces_to_surfaces_aligned.py
index 2cbe175..0e3d4ae 100644
--- a/experiments/Heart-Create_Statistical_Model/2-input_surfaces_to_surfaces_aligned.py
+++ b/experiments/Heart-Create_Statistical_Model/2-input_surfaces_to_surfaces_aligned.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %% [markdown]
 # # ICP Affine Registration: Align Heart Models to Average
 #
diff --git a/experiments/Heart-Create_Statistical_Model/3-registration_based_correspondence.py b/experiments/Heart-Create_Statistical_Model/3-registration_based_correspondence.py
index a27d628..ffd495b 100644
--- a/experiments/Heart-Create_Statistical_Model/3-registration_based_correspondence.py
+++ b/experiments/Heart-Create_Statistical_Model/3-registration_based_correspondence.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %% [markdown]
 # # Registration-Based Correspondence
 #
diff --git a/experiments/Heart-Create_Statistical_Model/4-surfaces_aligned_correspond_to_pca_inputs.py b/experiments/Heart-Create_Statistical_Model/4-surfaces_aligned_correspond_to_pca_inputs.py
index 92fde1b..26e125e 100644
--- a/experiments/Heart-Create_Statistical_Model/4-surfaces_aligned_correspond_to_pca_inputs.py
+++ b/experiments/Heart-Create_Statistical_Model/4-surfaces_aligned_correspond_to_pca_inputs.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %%
 from pathlib import Path
 
diff --git a/experiments/Heart-Create_Statistical_Model/5-compute_pca_model.py b/experiments/Heart-Create_Statistical_Model/5-compute_pca_model.py
index 2ddaf44..365309c 100644
--- a/experiments/Heart-Create_Statistical_Model/5-compute_pca_model.py
+++ b/experiments/Heart-Create_Statistical_Model/5-compute_pca_model.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %% [markdown]
 # # Compute PCA Statistical Model
 #
diff --git a/experiments/Heart-GatedCT_To_USD/0-download_and_convert_4d_to_3d.py b/experiments/Heart-GatedCT_To_USD/0-download_and_convert_4d_to_3d.py
index 25f56ca..6de2dec 100644
--- a/experiments/Heart-GatedCT_To_USD/0-download_and_convert_4d_to_3d.py
+++ b/experiments/Heart-GatedCT_To_USD/0-download_and_convert_4d_to_3d.py
@@ -1,9 +1,7 @@
-#!/usr/bin/env python
 # %%
 import shutil
 from pathlib import Path
 
-from physiotwin4d.convert_image_4d_to_3d import ConvertImage4DTo3D
 from physiotwin4d.data_download_tools import DataDownloadTools
 
 _HERE = Path(__file__).resolve().parent
@@ -12,15 +10,11 @@
 data_dir = _HERE.parent.parent / "data" / "Slicer-Heart-CT"
 output_dir = _HERE / "results"
 
-data_dir.mkdir(parents=True, exist_ok=True)
 output_dir.mkdir(parents=True, exist_ok=True)
 
-input_image_filename = DataDownloadTools.DownloadSlicerHeartCTData(data_dir)
+# Downloads TruncalValve_4DCT.seq.nrrd and splits it into slice_???.mha.
+DataDownloadTools.DownloadSlicerHeartCTData(data_dir)
 
 # %%
-conv = ConvertImage4DTo3D()
-conv.load_image_4d(str(input_image_filename))
-conv.save_3d_images(data_dir, "slice")
-
 # Save the mid-stroke slice as the fixed/reference image
 shutil.copyfile(data_dir / "slice_007.mha", output_dir / "slice_fixed.mha")
diff --git a/experiments/Heart-GatedCT_To_USD/1-register_images.py b/experiments/Heart-GatedCT_To_USD/1-register_images.py
index a665aec..4cf632e 100644
--- a/experiments/Heart-GatedCT_To_USD/1-register_images.py
+++ b/experiments/Heart-GatedCT_To_USD/1-register_images.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %%
 from pathlib import Path
 
diff --git a/experiments/Heart-GatedCT_To_USD/2-generate_segmentation.py b/experiments/Heart-GatedCT_To_USD/2-generate_segmentation.py
index 0df8de6..f8859c1 100644
--- a/experiments/Heart-GatedCT_To_USD/2-generate_segmentation.py
+++ b/experiments/Heart-GatedCT_To_USD/2-generate_segmentation.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %%
 import os
 
diff --git a/experiments/Heart-GatedCT_To_USD/3-transform_dynamic_and_static_contours.py b/experiments/Heart-GatedCT_To_USD/3-transform_dynamic_and_static_contours.py
index 0e7a660..166ddde 100644
--- a/experiments/Heart-GatedCT_To_USD/3-transform_dynamic_and_static_contours.py
+++ b/experiments/Heart-GatedCT_To_USD/3-transform_dynamic_and_static_contours.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %%
 import os
 
diff --git a/experiments/Heart-GatedCT_To_USD/4-merge_dynamic_and_static_usd.py b/experiments/Heart-GatedCT_To_USD/4-merge_dynamic_and_static_usd.py
index 5276e8e..9a66ef6 100644
--- a/experiments/Heart-GatedCT_To_USD/4-merge_dynamic_and_static_usd.py
+++ b/experiments/Heart-GatedCT_To_USD/4-merge_dynamic_and_static_usd.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %%
 import os
 
diff --git a/experiments/Heart-Simpleware_Segmentation/simpleware_heart_segmentation.py b/experiments/Heart-Simpleware_Segmentation/simpleware_heart_segmentation.py
index a040431..37a66df 100644
--- a/experiments/Heart-Simpleware_Segmentation/simpleware_heart_segmentation.py
+++ b/experiments/Heart-Simpleware_Segmentation/simpleware_heart_segmentation.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %% [markdown]
 # # Heart Segmentation using Simpleware Medical ASCardio
 #
diff --git a/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_model_icp_itk.py b/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_model_icp_itk.py
index ee82914..73c29d2 100644
--- a/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_model_icp_itk.py
+++ b/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_model_icp_itk.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %% [markdown]
 # # ICP via ITK Heart Model to Image Registration Experiment
 #
diff --git a/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_model_registration_pca.py b/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_model_registration_pca.py
index b210f66..19dc87c 100644
--- a/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_model_registration_pca.py
+++ b/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_model_registration_pca.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %% [markdown]
 # # PCA-based Heart Model to Image Registration Experiment
 #
diff --git a/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_patient-CHOPValve.py b/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_patient-CHOPValve.py
index 5678ecf..59ff4f8 100644
--- a/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_patient-CHOPValve.py
+++ b/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_patient-CHOPValve.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %% [markdown]
 # ## Setup and Imports
 
diff --git a/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_patient.py b/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_patient.py
index 17bfed3..66a846e 100644
--- a/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_patient.py
+++ b/experiments/Heart-Statistical_Model_To_Patient/heart_model_to_patient.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %% [markdown]
 # ## Setup and Imports
 
diff --git a/experiments/Heart-VTKSeries_To_USD/0-download_and_convert_4d_to_3d.py b/experiments/Heart-VTKSeries_To_USD/0-download_and_convert_4d_to_3d.py
index d129e3b..2851f6f 100644
--- a/experiments/Heart-VTKSeries_To_USD/0-download_and_convert_4d_to_3d.py
+++ b/experiments/Heart-VTKSeries_To_USD/0-download_and_convert_4d_to_3d.py
@@ -1,8 +1,6 @@
-#!/usr/bin/env python
 # %%
 import os
 
-from physiotwin4d.convert_image_4d_to_3d import ConvertImage4DTo3D
 from physiotwin4d.data_download_tools import DataDownloadTools
 
 _HERE = os.path.dirname(os.path.abspath(__file__))
@@ -11,16 +9,8 @@
 data_dir = os.path.join(_HERE, "..", "..", "data", "Slicer-Heart-CT")
 output_dir = os.path.join(_HERE, "results")
 
-if not os.path.exists(data_dir):
-    os.makedirs(data_dir)
-
 if not os.path.exists(output_dir):
     os.makedirs(output_dir)
 
-input_image_filename = DataDownloadTools.DownloadSlicerHeartCTData(data_dir)
-
-# %%
-if not os.path.exists(f"{data_dir}/slice_000.mha"):
-    conv = ConvertImage4DTo3D()
-    conv.load_image_4d(str(input_image_filename))
-    conv.save_3d_images(data_dir, "slice")
+# Downloads TruncalValve_4DCT.seq.nrrd and splits it into slice_???.mha.
+DataDownloadTools.DownloadSlicerHeartCTData(data_dir)
diff --git a/experiments/Heart-VTKSeries_To_USD/1-heart_vtkseries_to_usd.py b/experiments/Heart-VTKSeries_To_USD/1-heart_vtkseries_to_usd.py
index 7a239b4..8137abf 100644
--- a/experiments/Heart-VTKSeries_To_USD/1-heart_vtkseries_to_usd.py
+++ b/experiments/Heart-VTKSeries_To_USD/1-heart_vtkseries_to_usd.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %%
 import glob
 import os
diff --git a/experiments/Lung-GatedCT_To_USD/0-register_dirlab_4dct.py b/experiments/Lung-GatedCT_To_USD/0-register_dirlab_4dct.py
index ba8c21b..b509c09 100644
--- a/experiments/Lung-GatedCT_To_USD/0-register_dirlab_4dct.py
+++ b/experiments/Lung-GatedCT_To_USD/0-register_dirlab_4dct.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %%
 import os
 from typing import Optional
diff --git a/experiments/Lung-GatedCT_To_USD/1-make_dirlab_models.py b/experiments/Lung-GatedCT_To_USD/1-make_dirlab_models.py
index fb32134..290d069 100644
--- a/experiments/Lung-GatedCT_To_USD/1-make_dirlab_models.py
+++ b/experiments/Lung-GatedCT_To_USD/1-make_dirlab_models.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %%
 import itk
 import numpy as np
diff --git a/experiments/Lung-GatedCT_To_USD/2-paint_dirlab_models.py b/experiments/Lung-GatedCT_To_USD/2-paint_dirlab_models.py
index e53087c..7448b59 100644
--- a/experiments/Lung-GatedCT_To_USD/2-paint_dirlab_models.py
+++ b/experiments/Lung-GatedCT_To_USD/2-paint_dirlab_models.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %%
 from data_dirlab_4d_ct import DataDirLab4DCT
 from pxr import Usd
diff --git a/experiments/Lung-GatedCT_To_USD/Experiment_ArrangeOnStage.py b/experiments/Lung-GatedCT_To_USD/Experiment_ArrangeOnStage.py
index b9974d8..5923342 100644
--- a/experiments/Lung-GatedCT_To_USD/Experiment_ArrangeOnStage.py
+++ b/experiments/Lung-GatedCT_To_USD/Experiment_ArrangeOnStage.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %%
 import os
 
diff --git a/experiments/Lung-GatedCT_To_USD/Experiment_CombineModels.py b/experiments/Lung-GatedCT_To_USD/Experiment_CombineModels.py
index 8fafe5b..3f7cf54 100644
--- a/experiments/Lung-GatedCT_To_USD/Experiment_CombineModels.py
+++ b/experiments/Lung-GatedCT_To_USD/Experiment_CombineModels.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %%
 import os
 
diff --git a/experiments/Lung-GatedCT_To_USD/Experiment_SegReg.py b/experiments/Lung-GatedCT_To_USD/Experiment_SegReg.py
index cb996c4..ddb88fa 100644
--- a/experiments/Lung-GatedCT_To_USD/Experiment_SegReg.py
+++ b/experiments/Lung-GatedCT_To_USD/Experiment_SegReg.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %%
 import os
 
diff --git a/experiments/Lung-GatedCT_To_USD/Experiment_SubSurfaceScatter.py b/experiments/Lung-GatedCT_To_USD/Experiment_SubSurfaceScatter.py
index 8dd68e7..b41aa47 100644
--- a/experiments/Lung-GatedCT_To_USD/Experiment_SubSurfaceScatter.py
+++ b/experiments/Lung-GatedCT_To_USD/Experiment_SubSurfaceScatter.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %%
 import os
 
diff --git a/experiments/Lung-VesselsAirways/0-GenData.py b/experiments/Lung-VesselsAirways/0-GenData.py
index d150366..c0771b4 100644
--- a/experiments/Lung-VesselsAirways/0-GenData.py
+++ b/experiments/Lung-VesselsAirways/0-GenData.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %%
 from datetime import datetime
 
diff --git a/experiments/Reconstruct4DCT/reconstruct_4d_ct.py b/experiments/Reconstruct4DCT/reconstruct_4d_ct.py
index c9b1439..a01c3e6 100644
--- a/experiments/Reconstruct4DCT/reconstruct_4d_ct.py
+++ b/experiments/Reconstruct4DCT/reconstruct_4d_ct.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %%
 import os
 
diff --git a/experiments/Reconstruct4DCT/reconstruct_4d_ct_class.py b/experiments/Reconstruct4DCT/reconstruct_4d_ct_class.py
index ad1b320..931ca99 100644
--- a/experiments/Reconstruct4DCT/reconstruct_4d_ct_class.py
+++ b/experiments/Reconstruct4DCT/reconstruct_4d_ct_class.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 # %% [markdown]
 # # 4D CT Reconstruction Using RegisterTimeSeriesImages Class
 #
diff --git a/pyproject.toml b/pyproject.toml
index 945fb45..ec945e1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -78,7 +78,6 @@ dependencies = [
     "pyvista[all]>=0.47.0",
     "usd-core>=23.11",
     "trimesh>=4.0.0",
-    "netgen-mesher>=6.2.2606",
 
     # Utilities
     "ipykernel>=6.0.0",
@@ -229,8 +228,6 @@ module = [
     "itk.*",
     "matplotlib",
     "matplotlib.*",
-    "netgen",
-    "netgen.*",
     "nibabel",
     "nibabel.*",
     "nrrd",
diff --git a/src/physiotwin4d/cli/convert_image_4d_to_3d.py b/src/physiotwin4d/cli/convert_image_4d_to_3d.py
index 7de5a4f..f4fe923 100644
--- a/src/physiotwin4d/cli/convert_image_4d_to_3d.py
+++ b/src/physiotwin4d/cli/convert_image_4d_to_3d.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """Command-line interface for splitting a 3D/4D image into a 3D time series.
 
 Reads a 3D or 4D medical image and writes one ``.mha`` file per temporal frame
diff --git a/src/physiotwin4d/cli/convert_image_to_usd.py b/src/physiotwin4d/cli/convert_image_to_usd.py
index 309ce33..4a3c691 100644
--- a/src/physiotwin4d/cli/convert_image_to_usd.py
+++ b/src/physiotwin4d/cli/convert_image_to_usd.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Command-line interface for the Image-to-USD workflow.
 
diff --git a/src/physiotwin4d/cli/convert_image_to_vtk.py b/src/physiotwin4d/cli/convert_image_to_vtk.py
index bf21ee0..c38e4cb 100644
--- a/src/physiotwin4d/cli/convert_image_to_vtk.py
+++ b/src/physiotwin4d/cli/convert_image_to_vtk.py
@@ -1,9 +1,7 @@
-#!/usr/bin/env python
 """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 tetrahedral volume meshes annotated with anatomy labels and
-colors.
+surfaces annotated with anatomy labels and colors.
 """
 
 import argparse
@@ -27,7 +25,7 @@
 def main() -> int:
     """CLI entry point for image to VTK conversion."""
     parser = argparse.ArgumentParser(
-        description="Segment a 3D image and export anatomy groups as VTK surfaces and meshes.",
+        description="Segment a 3D image and export anatomy groups as VTK surfaces.",
         formatter_class=argparse.RawDescriptionHelpFormatter,
         epilog="""
 Anatomy groups
@@ -38,12 +36,14 @@ def main() -> int:
 Output files — combined mode (default)
 ---------------------------------------
   {prefix}_surfaces.vtp   all surfaces merged into one file
-  {prefix}_meshes.vtu     all tetrahedral volume meshes merged into one file
 
-Output files — split mode (--split-files)
-------------------------------------------
+Output files — group mode (--output-mode group)
+---------------------------------------------------
   {prefix}_{group}.vtp    one surface per anatomy group
-  {prefix}_{group}.vtu    one tetrahedral volume mesh per anatomy group
+
+Output files — label mode (--output-mode label)
+---------------------------------------------------
+  {prefix}_{label}.vtp    one surface per individual anatomical structure
 
 Examples
 --------
@@ -118,17 +118,6 @@ def main() -> int:
             "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(
@@ -137,12 +126,13 @@ def main() -> int:
         help="Filename prefix for output files (default: no prefix).",
     )
     parser.add_argument(
-        "--split-files",
-        action="store_true",
-        default=False,
+        "--output-mode",
+        choices=("combined", "group", "label"),
+        default="combined",
         help=(
-            "Write one VTP and one VTU file per anatomy group instead of "
-            "merging all groups into a single VTP and VTU (default: combined)."
+            "combined (default): merge all surfaces into one VTP.  "
+            "group: one VTP per anatomy group.  "
+            "label: one VTP per individual anatomical structure."
         ),
     )
     parser.add_argument(
@@ -188,17 +178,18 @@ def main() -> int:
             input_image=input_image,
             anatomy_groups=args.anatomy_groups,
             surface_target_reduction=args.surface_target_reduction,
-            mesh_target_reduction=args.mesh_target_reduction,
+            extract_label_surfaces=(args.output_mode == "label"),
         )
     except (ValueError, RuntimeError, OSError) as exc:
         print(f"Error during workflow: {exc}")
         traceback.print_exc()
         return 1
 
-    surfaces = result["surfaces"]
-    meshes = result["meshes"]
+    surfaces = (
+        result["label_surfaces"] if args.output_mode == "label" else result["surfaces"]
+    )
 
-    if not surfaces and not meshes:
+    if not surfaces:
         print("No anatomy groups produced any output.  Check the input image.")
         return 1
 
@@ -209,32 +200,18 @@ def main() -> int:
     prefix = args.output_prefix
 
     try:
-        if args.split_files:
-            # One file per anatomy group
-            if 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 = ContourTools.save_meshes(
-                    meshes, args.output_dir, prefix=prefix
-                )
-                for group, path in saved_meshes.items():
-                    print(f"  Mesh     [{group:15s}] -> {path}")
+        if args.output_mode == "combined":
+            surface_file = ContourTools.save_combined_surface(
+                surfaces, args.output_dir, prefix=prefix
+            )
+            print(f"  Combined surface -> {surface_file}")
         else:
-            # Combined single-file output
-            if surfaces:
-                surface_file = ContourTools.save_combined_surface(
-                    surfaces, args.output_dir, prefix=prefix
-                )
-                print(f"  Combined surface -> {surface_file}")
-            if meshes:
-                mesh_file = ContourTools.save_combined_mesh(
-                    meshes, args.output_dir, prefix=prefix
-                )
-                print(f"  Combined mesh    -> {mesh_file}")
+            # One file per anatomy group or per individual label
+            saved_surfaces = ContourTools.save_surfaces(
+                surfaces, args.output_dir, prefix=prefix
+            )
+            for name, path in saved_surfaces.items():
+                print(f"  Surface  [{name:20s}] -> {path}")
 
         if args.save_labelmap:
             labelmap = result["labelmap"]
diff --git a/src/physiotwin4d/cli/convert_vtk_to_usd.py b/src/physiotwin4d/cli/convert_vtk_to_usd.py
index daef1a9..56feaea 100644
--- a/src/physiotwin4d/cli/convert_vtk_to_usd.py
+++ b/src/physiotwin4d/cli/convert_vtk_to_usd.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Command-line interface for VTK to USD conversion workflow.
 
diff --git a/src/physiotwin4d/cli/create_statistical_model.py b/src/physiotwin4d/cli/create_statistical_model.py
index 5f1b139..62bc248 100644
--- a/src/physiotwin4d/cli/create_statistical_model.py
+++ b/src/physiotwin4d/cli/create_statistical_model.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Command-line interface for Create Statistical Model workflow.
 
diff --git a/src/physiotwin4d/cli/download_data.py b/src/physiotwin4d/cli/download_data.py
index 990e4ec..5399ec6 100644
--- a/src/physiotwin4d/cli/download_data.py
+++ b/src/physiotwin4d/cli/download_data.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """Command-line interface for downloading PhysioTwin4D example data."""
 
 from __future__ import annotations
diff --git a/src/physiotwin4d/cli/fit_statistical_model_to_patient.py b/src/physiotwin4d/cli/fit_statistical_model_to_patient.py
index 1bdd7f0..85976a9 100644
--- a/src/physiotwin4d/cli/fit_statistical_model_to_patient.py
+++ b/src/physiotwin4d/cli/fit_statistical_model_to_patient.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Command-line interface for Heart Model to Patient Registration workflow.
 
diff --git a/src/physiotwin4d/cli/reconstruct_highres_4d_ct.py b/src/physiotwin4d/cli/reconstruct_highres_4d_ct.py
index 32d44a2..3e0fc7d 100644
--- a/src/physiotwin4d/cli/reconstruct_highres_4d_ct.py
+++ b/src/physiotwin4d/cli/reconstruct_highres_4d_ct.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Command-line interface for high-resolution 4D CT reconstruction workflow.
 
diff --git a/src/physiotwin4d/cli/visualize_pca_modes.py b/src/physiotwin4d/cli/visualize_pca_modes.py
index 2ab8dbf..51d1c04 100644
--- a/src/physiotwin4d/cli/visualize_pca_modes.py
+++ b/src/physiotwin4d/cli/visualize_pca_modes.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Command-line interface to visualize PCA modes of variation.
 
diff --git a/src/physiotwin4d/contour_tools.py b/src/physiotwin4d/contour_tools.py
index 533e7e9..2c46267 100644
--- a/src/physiotwin4d/contour_tools.py
+++ b/src/physiotwin4d/contour_tools.py
@@ -6,7 +6,7 @@
 
 import logging
 import os
-from typing import Optional, cast
+from typing import cast
 
 import itk
 import numpy as np
@@ -71,84 +71,6 @@ 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,
@@ -571,33 +493,6 @@ def save_surfaces(
             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],
@@ -633,38 +528,3 @@ def save_combined_surface(
         )
         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/physiotwin4d/data_download_tools.py b/src/physiotwin4d/data_download_tools.py
index 8e699a4..879352e 100644
--- a/src/physiotwin4d/data_download_tools.py
+++ b/src/physiotwin4d/data_download_tools.py
@@ -18,6 +18,8 @@
 from pathlib import Path
 from typing import Union
 
+from .convert_image_4d_to_3d import ConvertImage4DTo3D
+
 _DOWNLOAD_TIMEOUT_SECONDS = 60.0
 _logger = logging.getLogger(__name__)
 
@@ -31,12 +33,20 @@ class DataDownloadTools:
     )
     SLICER_HEART_CT_FILENAME = "TruncalValve_4DCT.seq.nrrd"
 
+    SLICER_HEART_CT_SLICE_BASENAME = "slice"
+
     @staticmethod
     def DownloadSlicerHeartCTData(dirname: Union[str, Path]) -> Path:  # noqa: N802
         """Download the Slicer-Heart-CT 4-D CT sample into ``dirname``.
 
+        Also splits the downloaded 4-D sequence into per-frame
+        ``slice_???.mha`` volumes in the same directory (via
+        ``ConvertImage4DTo3D``), skipping the split if those files already
+        exist.
+
         Args:
-            dirname: Directory where ``TruncalValve_4DCT.seq.nrrd`` should live.
+            dirname: Directory where ``TruncalValve_4DCT.seq.nrrd`` and the
+                converted ``slice_???.mha`` files should live.
 
         Returns:
             Path to the downloaded or already-cached ``.seq.nrrd`` file.
@@ -45,41 +55,60 @@ def DownloadSlicerHeartCTData(dirname: Union[str, Path]) -> Path:  # noqa: N802
         data_dir.mkdir(parents=True, exist_ok=True)
 
         data_file = data_dir / DataDownloadTools.SLICER_HEART_CT_FILENAME
-        if data_file.exists() and data_file.stat().st_size > 0:
-            return data_file
-
-        # Stream to a unique temp file in the same directory with an explicit
-        # timeout, then atomically replace the target on success. The temp
-        # name is unique so concurrent callers do not clobber each other.
-        # Avoids partial files on interrupt and the indefinite hang that
-        # urlretrieve has without a timeout.
-        tmp_handle = tempfile.NamedTemporaryFile(
-            dir=str(data_dir),
-            prefix=f".{DataDownloadTools.SLICER_HEART_CT_FILENAME}.",
-            suffix=".tmp",
-            delete=False,
-        )
-        tmp_file = Path(tmp_handle.name)
-        try:
-            with (
-                urllib.request.urlopen(  # noqa: S310
-                    DataDownloadTools.SLICER_HEART_CT_URL,
-                    timeout=_DOWNLOAD_TIMEOUT_SECONDS,
-                ) as response,
-                tmp_handle as out,
-            ):
-                shutil.copyfileobj(response, out)
-            if tmp_file.stat().st_size == 0:
-                raise RuntimeError(
-                    f"Downloaded file is empty: {DataDownloadTools.SLICER_HEART_CT_URL}"
+        if not (data_file.exists() and data_file.stat().st_size > 0):
+            # Stream to a unique temp file in the same directory with an
+            # explicit timeout, then atomically replace the target on
+            # success. The temp name is unique so concurrent callers do not
+            # clobber each other. Avoids partial files on interrupt and the
+            # indefinite hang that urlretrieve has without a timeout.
+            tmp_handle = tempfile.NamedTemporaryFile(
+                dir=str(data_dir),
+                prefix=f".{DataDownloadTools.SLICER_HEART_CT_FILENAME}.",
+                suffix=".tmp",
+                delete=False,
+            )
+            tmp_file = Path(tmp_handle.name)
+            try:
+                with (
+                    urllib.request.urlopen(  # noqa: S310
+                        DataDownloadTools.SLICER_HEART_CT_URL,
+                        timeout=_DOWNLOAD_TIMEOUT_SECONDS,
+                    ) as response,
+                    tmp_handle as out,
+                ):
+                    shutil.copyfileobj(response, out)
+                if tmp_file.stat().st_size == 0:
+                    raise RuntimeError(
+                        f"Downloaded file is empty: {DataDownloadTools.SLICER_HEART_CT_URL}"
+                    )
+                tmp_file.replace(data_file)
+                _logger.info(
+                    "Downloaded %s", DataDownloadTools.SLICER_HEART_CT_FILENAME
                 )
-            tmp_file.replace(data_file)
-            _logger.info("Downloaded %s", DataDownloadTools.SLICER_HEART_CT_FILENAME)
-        except BaseException:
-            tmp_handle.close()
-            if tmp_file.exists():
-                tmp_file.unlink()
-            raise
+            except BaseException:
+                tmp_handle.close()
+                if tmp_file.exists():
+                    tmp_file.unlink()
+                raise
+
+        slice_basename = DataDownloadTools.SLICER_HEART_CT_SLICE_BASENAME
+        if not any(data_dir.glob(f"{slice_basename}_???.mha")):
+            # Convert into a temp directory first, then atomically move each
+            # finished frame into data_dir, so a failed or interrupted
+            # conversion leaves no partial slice_???.mha files behind and
+            # concurrent callers never observe a half-written frame.
+            with tempfile.TemporaryDirectory(dir=str(data_dir)) as tmp_dir_name:
+                tmp_dir = Path(tmp_dir_name)
+                conv = ConvertImage4DTo3D()
+                conv.load_image_4d(str(data_file))
+                conv.save_3d_images(tmp_dir, slice_basename)
+                for tmp_slice_file in sorted(tmp_dir.glob(f"{slice_basename}_???.mha")):
+                    tmp_slice_file.replace(data_dir / tmp_slice_file.name)
+            _logger.info(
+                "Converted %s to %s_???.mha frames",
+                DataDownloadTools.SLICER_HEART_CT_FILENAME,
+                slice_basename,
+            )
         return data_file
 
     @staticmethod
diff --git a/src/physiotwin4d/workflow_convert_image_to_vtk.py b/src/physiotwin4d/workflow_convert_image_to_vtk.py
index 2770a93..66a35b6 100644
--- a/src/physiotwin4d/workflow_convert_image_to_vtk.py
+++ b/src/physiotwin4d/workflow_convert_image_to_vtk.py
@@ -1,8 +1,7 @@
-"""Workflow for segmenting a CT image and converting anatomy groups to VTK surfaces and meshes.
+"""Workflow for segmenting a CT image and converting anatomy groups to VTK surfaces.
 
-The workflow segments a 3D CT image using a chosen backend, then extracts one VTP
-(surface) and one VTU (tetrahedral volume mesh, generated by netgen from that
-surface) per non-empty anatomy group.  Each output object carries anatomy
+The workflow segments a 3D CT image using a chosen backend, then extracts one
+VTP surface per non-empty anatomy group.  Each output surface 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.
@@ -18,17 +17,17 @@
     ct = itk.imread('chest_ct.nii.gz')
     segmenter = SegmentChestTotalSegmentatorWithContrast()
     workflow = WorkflowConvertImageToVTK(segmentation_method=segmenter)
-    result = workflow.process(
-        ct, surface_target_reduction=0.5, mesh_target_reduction=0.7
-    )
+    result = workflow.process(ct, surface_target_reduction=0.5)
 
     # Combined single-file output (default)
     ContourTools.save_combined_surface(result['surfaces'], './out', prefix='patient')
-    ContourTools.save_combined_mesh(result['meshes'], './out', prefix='patient')
 
     # Per-group split output
     ContourTools.save_surfaces(result['surfaces'], './out', prefix='patient')
-    ContourTools.save_meshes(result['meshes'], './out', prefix='patient')
+
+    # Per-label split output (one VTP per individual anatomical structure)
+    result = workflow.process(ct, extract_label_surfaces=True)
+    ContourTools.save_surfaces(result['label_surfaces'], './out', prefix='patient')
 """
 
 import logging
@@ -48,7 +47,7 @@
 
 
 class WorkflowConvertImageToVTK(PhysioTwin4DBase):
-    """Segment a CT image and produce per-anatomy-group VTK surfaces and meshes.
+    """Segment a CT image and produce per-anatomy-group VTK surfaces.
 
     ``segmentation_method`` accepts a pre-configured
     :class:`SegmentAnatomyBase` instance (e.g. :class:`SegmentChestTotalSegmentator`,
@@ -62,12 +61,14 @@ class WorkflowConvertImageToVTK(PhysioTwin4DBase):
 
     Determined by the active segmenter's :attr:`SegmentAnatomyBase.taxonomy`
     (see :attr:`ANATOMY_GROUPS`).  Groups that are empty after segmentation
-    are silently skipped.
+    are silently skipped.  Pass ``extract_label_surfaces=True`` to
+    :meth:`process` to additionally extract one surface per individual
+    structure (label) within each group, e.g. ``left_ventricle`` separately
+    from ``right_ventricle`` within ``heart``.
 
     **VTK object annotation**
 
-    Each :class:`pyvista.PolyData` surface and :class:`pyvista.UnstructuredGrid` mesh
-    returned by :meth:`process` carries:
+    Each :class:`pyvista.PolyData` surface returned by :meth:`process` carries:
 
     - ``field_data['AnatomyGroup']`` — anatomy group name, e.g. ``'heart'``.
     - ``field_data['SegmentationLabelNames']`` — individual structure names within the
@@ -80,9 +81,8 @@ class WorkflowConvertImageToVTK(PhysioTwin4DBase):
 
     :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
+    :meth:`ContourTools.save_surfaces` and
+    :meth:`ContourTools.save_combined_surface` — or the CLI
     ``physiotwin4d-convert-image-to-vtk`` — to write results to disk.
     """
 
@@ -187,25 +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, surface: pv.PolyData, mesh_target_reduction: float = 0.0
-    ) -> Optional[pv.UnstructuredGrid]:
-        """Generate a tetrahedral volume mesh (VTU) from a closed surface.
+    def _extract_label_surface(
+        self, labelmap_image: Any, labelmap_arr: np.ndarray, label_id: int
+    ) -> Optional[pv.PolyData]:
+        """Extract a smoothed triangulated surface for one individual label.
 
-        Delegates to :meth:`ContourTools.extract_mesh`.
+        Isolates *label_id* out of *labelmap_arr* into its own binary mask
+        before delegating to :meth:`ContourTools.extract_contours`.
 
         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.
+            labelmap_image: Source image, used only for ``CopyInformation``
+                (origin/spacing/direction) on the isolated label mask.
+            labelmap_arr: ``labelmap_image``'s voxel array. Callers extracting
+                multiple labels from the same labelmap (as :meth:`process`
+                does) should compute this once — e.g. via
+                ``itk.GetArrayViewFromImage`` — and pass the same array to
+                every call, rather than re-deriving it per label.
+            label_id: Integer label id to isolate.
 
         Returns:
-            :class:`pyvista.UnstructuredGrid` of tetrahedral cells, or
-            ``None`` if the surface is empty or no volume mesh could be
-            generated from it.
+            Smoothed :class:`pyvista.PolyData`, or ``None`` if *label_id* has
+            no voxels in *labelmap_arr*.
         """
-        return self._contour_tools.extract_mesh(surface, mesh_target_reduction)
+        label_mask_arr = (labelmap_arr == label_id).astype(np.uint8)
+        if int(label_mask_arr.sum()) == 0:
+            return None
+        label_mask = itk.GetImageFromArray(label_mask_arr)
+        label_mask.CopyInformation(labelmap_image)
+        return self._contour_tools.extract_contours(label_mask)
 
     # ─────────────────────────── Main workflow ─────────────────────────────
 
@@ -214,9 +223,9 @@ def process(
         input_image: Any,
         anatomy_groups: Optional[list[str]] = None,
         surface_target_reduction: float = 0.0,
-        mesh_target_reduction: float = 0.0,
+        extract_label_surfaces: bool = False,
     ) -> dict[str, Any]:
-        """Segment the CT image and extract per-anatomy-group VTK objects.
+        """Segment the CT image and extract per-anatomy-group VTK surfaces.
 
         Args:
             input_image: Input 3D CT image (``itk.Image``).
@@ -227,19 +236,20 @@ def process(
             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.
+                Applied to both group and (when requested) label surfaces.
+            extract_label_surfaces: When ``True``, also extract one surface per
+                individual anatomical structure (label) within each processed
+                group — e.g. ``left_ventricle`` and ``right_ventricle``
+                separately within the ``heart`` group — in addition to the
+                per-group surfaces.  ``False`` (default) skips this and leaves
+                ``'label_surfaces'`` empty.
 
         Returns:
             ``dict`` with the following keys:
 
             - ``'surfaces'`` — ``dict[str, pv.PolyData]``: smoothed surface per group.
-            - ``'meshes'`` — ``dict[str, pv.UnstructuredGrid]``: tetrahedral volume
-              mesh per group, generated by netgen from that group's surface.
+            - ``'label_surfaces'`` — ``dict[str, pv.PolyData]``: smoothed surface per
+              individual label, populated only when *extract_label_surfaces* is True.
             - ``'labelmap'`` — ``itk.Image``: detailed per-structure segmentation
               labelmap from the segmenter.
             - ``'segmentation_masks'`` — ``dict[str, itk.Image]``: per-group binary
@@ -268,10 +278,15 @@ def process(
         self.log_section("Running segmentation")
         seg_result: dict[str, Any] = self._segmenter.segment(input_image)
 
+        # A zero-copy view, computed once and shared across every label
+        # surface extracted below so a multi-group/multi-label run doesn't
+        # repeatedly re-copy the (often large) labelmap volume.
+        labelmap_arr = itk.GetArrayViewFromImage(seg_result["labelmap"])
+
         # Extract VTK objects per anatomy group
         self.log_section("Extracting VTK objects")
         surfaces: dict[str, pv.PolyData] = {}
-        meshes: dict[str, pv.UnstructuredGrid] = {}
+        label_surfaces: dict[str, pv.PolyData] = {}
         seg_masks: dict[str, Any] = {}
 
         for group in groups_to_process:
@@ -305,19 +320,30 @@ def process(
             self._annotate(export_surface, group, label_names, label_ids, color)
             surfaces[group] = export_surface
 
-            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
+            if extract_label_surfaces:
+                self.log_info("  Extracting label surfaces for: %s", group)
+                for label_id, label_name in zip(label_ids, label_names, strict=True):
+                    label_surface = self._extract_label_surface(
+                        seg_result["labelmap"], labelmap_arr, label_id
+                    )
+                    if label_surface is None:
+                        continue
+                    if surface_target_reduction > 0.0:
+                        label_surface = label_surface.decimate_pro(
+                            surface_target_reduction, preserve_topology=True
+                        )
+                    self._annotate(
+                        label_surface, group, [label_name], [label_id], color
+                    )
+                    label_surfaces[label_name] = label_surface
 
         self.log_section("IMAGE TO VTK WORKFLOW COMPLETE")
-        self.log_info("Surfaces extracted: %d", len(surfaces))
-        self.log_info("Meshes extracted:   %d", len(meshes))
+        self.log_info("Surfaces extracted:       %d", len(surfaces))
+        self.log_info("Label surfaces extracted: %d", len(label_surfaces))
 
         return {
             "surfaces": surfaces,
-            "meshes": meshes,
+            "label_surfaces": label_surfaces,
             "labelmap": seg_result["labelmap"],
             "segmentation_masks": seg_masks,
         }
diff --git a/src/physiotwin4d/workflow_fit_statistical_model_to_patient.py b/src/physiotwin4d/workflow_fit_statistical_model_to_patient.py
index dab5bcb..94c25ab 100644
--- a/src/physiotwin4d/workflow_fit_statistical_model_to_patient.py
+++ b/src/physiotwin4d/workflow_fit_statistical_model_to_patient.py
@@ -194,7 +194,7 @@ def __init__(
                 input_image=patient_image,
                 anatomy_groups=["heart"],
             )
-            patient_models = [patient_models_data["meshes"]["heart"]]
+            patient_models = [patient_models_data["surfaces"]["heart"]]
         elif patient_models is None:
             raise ValueError("Either patient_models or patient_image must be provided.")
         self.patient_models = patient_models
diff --git a/statistics.md b/statistics.md
index 7b42604..679c151 100644
--- a/statistics.md
+++ b/statistics.md
@@ -112,7 +112,7 @@ PhysioTwin4D operates across several technically demanding domains:
 | **Medical Imaging**    | ITK, MONAI, nibabel, pydicom, pynrrd                 |
 | **Deep Learning**      | PyTorch, CuPy (CUDA 13), transformers                |
 | **Registration**       | ANTs (antspyx), picsl-greedy, icon-registration, UniGradICON |
-| **3D Graphics / USD**  | VTK, PyVista, USD-core, trimesh, netgen-mesher        |
+| **3D Graphics / USD**  | VTK, PyVista, USD-core, trimesh                       |
 | **AI Segmentation**    | TotalSegmentator                                      |
 | **Development Tools**  | pytest, pytest-cov, pytest-xdist, ruff, mypy, sphinx, uv |
 
diff --git a/tests/conftest.py b/tests/conftest.py
index 3e13c68..af220d8 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Shared pytest fixtures for PhysioTwin4D tests.
 
@@ -15,7 +14,6 @@
 import pytest
 
 from physiotwin4d.contour_tools import ContourTools
-from physiotwin4d.convert_image_4d_to_3d import ConvertImage4DTo3D
 from physiotwin4d.data_download_tools import DataDownloadTools
 from physiotwin4d.register_images_ants import RegisterImagesANTS
 from physiotwin4d.register_images_greedy import RegisterImagesGreedy
@@ -456,21 +454,12 @@ def test_images(
     download_test_data: Path,
     test_directories: dict[str, Path],
 ) -> list[Any]:
-    """Convert and resample 4D NRRD data; return pre-resampled time points."""
+    """Resample DownloadSlicerHeartCTData's 3D time series; return time points."""
     data_dir = test_directories["slicer_heart_data"]
     small_data_dir = test_directories["slicer_heart_small_data"]
 
-    # Convert 4D NRRD to 3D time series if not already done
-    slice_000 = data_dir / "slice_000.mha"
-    slice_007 = data_dir / "slice_007.mha"
-    if not slice_000.exists() or not slice_007.exists():
-        print("\nConverting 4D image to 3D time series...")
-        conv = ConvertImage4DTo3D()
-        conv.load_image_4d(str(download_test_data))
-        conv.save_3d_images(data_dir, "slice")
-    else:
-        print("\n3D slice files already exist")
-
+    # DownloadSlicerHeartCTData() already split the 4D NRRD into
+    # slice_???.mha 3D time-series volumes.
     # Resample each slice_???.mha to 1.5x1.5x1.5 mm into slicer_heart_small.
     target_spacing = [1.5, 1.5, 1.5]
     for slice_file in sorted(data_dir.glob("slice_???.mha")):
diff --git a/tests/test_contour_tools.py b/tests/test_contour_tools.py
index 010a5f8..d01bdbd 100644
--- a/tests/test_contour_tools.py
+++ b/tests/test_contour_tools.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Test for contour tools functionality.
 
diff --git a/tests/test_convert_image_4d_to_3d.py b/tests/test_convert_image_4d_to_3d.py
index aa8eff8..718174e 100644
--- a/tests/test_convert_image_4d_to_3d.py
+++ b/tests/test_convert_image_4d_to_3d.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Test for converting a 4D image to a 3D time series using ITK readers.
 
diff --git a/tests/test_convert_vtk_to_usd.py b/tests/test_convert_vtk_to_usd.py
index 3bf626d..35c1ec6 100644
--- a/tests/test_convert_vtk_to_usd.py
+++ b/tests/test_convert_vtk_to_usd.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Test for VTK to USD conversion.
 
diff --git a/tests/test_download_heart_data.py b/tests/test_download_heart_data.py
index 9d15e8f..0143e6c 100644
--- a/tests/test_download_heart_data.py
+++ b/tests/test_download_heart_data.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Test for downloading and converting Slicer-Heart-CT data.
 
diff --git a/tests/test_experiments.py b/tests/test_experiments.py
index 5060f9d..9d73368 100644
--- a/tests/test_experiments.py
+++ b/tests/test_experiments.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Test suite for running experiment scripts.
 
diff --git a/tests/test_image_tools.py b/tests/test_image_tools.py
index a1edbc0..13af859 100644
--- a/tests/test_image_tools.py
+++ b/tests/test_image_tools.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Tests for ImageTools functionality.
 
diff --git a/tests/test_labelmap_tools.py b/tests/test_labelmap_tools.py
index ed0df03..35636df 100644
--- a/tests/test_labelmap_tools.py
+++ b/tests/test_labelmap_tools.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Tests for LabelmapTools functionality.
 
diff --git a/tests/test_register_images_ants.py b/tests/test_register_images_ants.py
index a36fe41..a40a523 100644
--- a/tests/test_register_images_ants.py
+++ b/tests/test_register_images_ants.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Test for ANTs-based image registration.
 
diff --git a/tests/test_register_images_greedy.py b/tests/test_register_images_greedy.py
index b832060..d079c70 100644
--- a/tests/test_register_images_greedy.py
+++ b/tests/test_register_images_greedy.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Tests for Greedy-based image registration.
 
diff --git a/tests/test_register_images_icon.py b/tests/test_register_images_icon.py
index c700987..8b25431 100644
--- a/tests/test_register_images_icon.py
+++ b/tests/test_register_images_icon.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Test for ICON-based image registration.
 
diff --git a/tests/test_register_time_series_images.py b/tests/test_register_time_series_images.py
index 9db3c9a..e7268d9 100644
--- a/tests/test_register_time_series_images.py
+++ b/tests/test_register_time_series_images.py
@@ -1,4 +1,3 @@
-#! /usr/bin/env python
 """
 Test for time series image registration.
 
diff --git a/tests/test_segment_chest_total_segmentator.py b/tests/test_segment_chest_total_segmentator.py
index cb0b962..e11c332 100644
--- a/tests/test_segment_chest_total_segmentator.py
+++ b/tests/test_segment_chest_total_segmentator.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Test for chest CT segmentation using TotalSegmentator.
 
diff --git a/tests/test_segment_heart_simpleware.py b/tests/test_segment_heart_simpleware.py
index 3d5cb9b..f7c3f06 100644
--- a/tests/test_segment_heart_simpleware.py
+++ b/tests/test_segment_heart_simpleware.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Tests for heart CT segmentation using SegmentHeartSimpleware (Simpleware Medical ASCardio).
 
diff --git a/tests/test_segment_heart_simpleware_trimmed_branches.py b/tests/test_segment_heart_simpleware_trimmed_branches.py
index 96847a9..a05d6c2 100644
--- a/tests/test_segment_heart_simpleware_trimmed_branches.py
+++ b/tests/test_segment_heart_simpleware_trimmed_branches.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Tests for SegmentHeartSimplewareTrimmedBranches.
 
diff --git a/tests/test_transform_tools.py b/tests/test_transform_tools.py
index b90b06c..5ed9938 100644
--- a/tests/test_transform_tools.py
+++ b/tests/test_transform_tools.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Test for transform tools functionality.
 
diff --git a/tests/test_tutorials.py b/tests/test_tutorials.py
index 5abb201..12fd88c 100644
--- a/tests/test_tutorials.py
+++ b/tests/test_tutorials.py
@@ -113,7 +113,6 @@ def test_run(self, test_directories: dict[str, Path]) -> None:
         out_dir = _REPO_ROOT / "tutorials" / "output" / "tutorial_02"
         results = _run_tutorial_script("tutorial_02_ct_to_vtk.py")
         assert results["surface_file"].exists(), "Combined VTP surface should exist"
-        assert results["mesh_file"].exists(), "Combined VTU mesh should exist"
 
         tt = TestTools(
             class_name=self._class_name,
@@ -137,7 +136,7 @@ class TestTutorial03CreateStatisticalModel:
 
     def test_run(self, test_directories: dict[str, Path]) -> None:
         kcl_dir = test_directories["data"] / "KCL-Heart-Model"
-        if not (kcl_dir / "pca_mean.vtu").exists():
+        if not (kcl_dir / "average_mesh.vtk").exists():
             pytest.skip(
                 "KCL-Heart-Model not downloaded. See data/README.md for instructions."
             )
@@ -164,7 +163,7 @@ class TestTutorial04FitStatisticalModelToPatient:
 
     def test_run(self, test_directories: dict[str, Path]) -> None:
         kcl_dir = test_directories["data"] / "KCL-Heart-Model"
-        if not (kcl_dir / "pca_mean.vtu").exists():
+        if not (kcl_dir / "average_mesh.vtk").exists():
             pytest.skip(
                 "KCL-Heart-Model not downloaded. See data/README.md for instructions."
             )
diff --git a/tests/test_vtk_to_usd_library.py b/tests/test_vtk_to_usd_library.py
index ff30a3d..4ea02c9 100644
--- a/tests/test_vtk_to_usd_library.py
+++ b/tests/test_vtk_to_usd_library.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Tests for VTK-to-USD conversion through ConvertVTKToUSD.
 
diff --git a/tests/test_workflow_convert_image_to_vtk.py b/tests/test_workflow_convert_image_to_vtk.py
index 83d2455..f3fbee3 100644
--- a/tests/test_workflow_convert_image_to_vtk.py
+++ b/tests/test_workflow_convert_image_to_vtk.py
@@ -4,6 +4,8 @@
 
 from typing import Any
 
+import itk
+import numpy as np
 import pytest
 
 from physiotwin4d.segment_chest_total_segmentator import SegmentChestTotalSegmentator
@@ -30,3 +32,27 @@ def test_caller_supplied_instance_is_used_as_is() -> None:
     segmenter = SegmentHeartSimpleware()
     workflow = WorkflowConvertImageToVTK(segmentation_method=segmenter)
     assert workflow._segmenter is segmenter
+
+
+def test_extract_label_surface_isolates_single_label() -> None:
+    """_extract_label_surface must isolate one label id from a multi-value labelmap.
+
+    Regression coverage for the extract_label_surfaces=True path added to
+    process(): each label's surface should come from only that label's
+    voxels, not the whole labelmap, and an absent label id should yield
+    None (mirrors the empty-mask skip used for whole-group surfaces).
+    """
+    arr = np.zeros((10, 10, 10), dtype=np.uint8)
+    arr[2:5, 2:5, 2:5] = 1
+    arr[6:9, 6:9, 6:9] = 2
+    labelmap = itk.image_from_array(arr)
+
+    workflow = WorkflowConvertImageToVTK()
+
+    surface_1 = workflow._extract_label_surface(labelmap, arr, 1)
+    surface_2 = workflow._extract_label_surface(labelmap, arr, 2)
+    surface_absent = workflow._extract_label_surface(labelmap, arr, 99)
+
+    assert surface_1 is not None and surface_1.n_points > 0
+    assert surface_2 is not None and surface_2.n_points > 0
+    assert surface_absent is None
diff --git a/tests/test_workflow_fit_statistical_model_to_patient.py b/tests/test_workflow_fit_statistical_model_to_patient.py
index 6bceea7..d65ae47 100644
--- a/tests/test_workflow_fit_statistical_model_to_patient.py
+++ b/tests/test_workflow_fit_statistical_model_to_patient.py
@@ -82,7 +82,7 @@ def test_fit_workflow_routes_default_to_image_to_vtk_with_trimmed_branches(
             dtype=np.float64,
         )
     )
-    heart_mesh = pv.PolyData(
+    heart_surface = pv.PolyData(
         np.array(
             [
                 [0.0, 0.0, 0.0],
@@ -102,7 +102,7 @@ def __init__(self, **kwargs: Any) -> None:
 
         def process(self, **kwargs: Any) -> dict[str, Any]:
             captured["process_kwargs"] = kwargs
-            return {"meshes": {"heart": heart_mesh}}
+            return {"surfaces": {"heart": heart_surface}}
 
     monkeypatch.setattr(
         "physiotwin4d.workflow_fit_statistical_model_to_patient."
@@ -120,7 +120,7 @@ def process(self, **kwargs: Any) -> dict[str, Any]:
         SegmentHeartSimplewareTrimmedBranches,
     )
     assert captured["process_kwargs"]["anatomy_groups"] == ["heart"]
-    assert workflow.patient_models == [heart_mesh]
+    assert workflow.patient_models == [heart_surface]
 
 
 def test_image_to_vtk_segmenter_uses_supplied_instance() -> None:
diff --git a/tutorials/tutorial_02_ct_to_vtk.py b/tutorials/tutorial_02_ct_to_vtk.py
index 571efdf..30369ec 100644
--- a/tutorials/tutorial_02_ct_to_vtk.py
+++ b/tutorials/tutorial_02_ct_to_vtk.py
@@ -3,9 +3,9 @@
 
 Purpose
 -------
-Segment one 3D CT frame into anatomical groups and save combined VTK surface
-and tetrahedral volume mesh files. The output can be inspected directly in
-PyVista or used as input for Tutorial 5.
+Segment one 3D CT frame into anatomical groups and save a combined VTK
+surface file. The output can be inspected directly in PyVista or used as
+input for Tutorial 5.
 
 Data Required
 -------------
@@ -48,6 +48,12 @@
     BASELINES_DIR = REPO_ROOT / "tests" / "baselines"
     LOG_LEVEL = logging.INFO
 
+    # In addition to the combined surface file always saved below, also
+    # save one VTP per anatomy group (e.g. heart.vtp, lung.vtp) and/or one
+    # VTP per individual anatomical structure (e.g. left_ventricle.vtp).
+    SAVE_GROUP_SURFACES = True
+    SAVE_LABEL_SURFACES = True
+
     # %%
     # Data reading
     test_mode = TestTools.running_as_test()
@@ -81,13 +87,11 @@
     # %%
     # Workflow execution
     #
-    # 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).
+    # surface_target_reduction decimates each exported VTP surface.
     result = workflow.process(
         input_image=ct_image,
         surface_target_reduction=0.5,
-        mesh_target_reduction=0.7,
+        extract_label_surfaces=SAVE_LABEL_SURFACES,
     )
 
     # %%
@@ -99,13 +103,14 @@
             prefix="patient",
         )
     )
-    mesh_file = Path(
-        ContourTools.save_combined_mesh(
-            result["meshes"],
-            str(output_dir),
-            prefix="patient",
+    if SAVE_GROUP_SURFACES:
+        ContourTools.save_surfaces(
+            result["surfaces"], str(output_dir), prefix="patient"
+        )
+    if SAVE_LABEL_SURFACES:
+        ContourTools.save_surfaces(
+            result["label_surfaces"], str(output_dir), prefix="patient"
         )
-    )
     labelmap_file = output_dir / "patient_labelmap.mha"
     itk.imwrite(result["labelmap"], str(labelmap_file), compression=True)
 
@@ -148,7 +153,6 @@
     tutorial_results = {
         "result": result,
         "surface_file": surface_file,
-        "mesh_file": mesh_file,
         "labelmap_file": labelmap_file,
         "screenshots": screenshots,
     }
diff --git a/tutorials/tutorial_03_create_statistical_model.py b/tutorials/tutorial_03_create_statistical_model.py
index 9726f80..a411ddf 100644
--- a/tutorials/tutorial_03_create_statistical_model.py
+++ b/tutorials/tutorial_03_create_statistical_model.py
@@ -46,7 +46,6 @@
     OUTPUT_DIR = TUTORIALS_DIR / "output" / "tutorial_03"
     BASELINES_DIR = REPO_ROOT / "tests" / "baselines"
     PCA_COMPONENTS = 10
-    MAX_SAMPLES = 20
     LOG_LEVEL = logging.INFO
 
     # %%
@@ -59,26 +58,25 @@
 
     if test_mode:
         pca_components = min(PCA_COMPONENTS, 5)
-        max_samples = min(MAX_SAMPLES, 10)
     else:
         pca_components = PCA_COMPONENTS
-        max_samples = MAX_SAMPLES
 
     output_dir.mkdir(parents=True, exist_ok=True)
 
-    reference_file = data_dir / "pca_mean.vtu"
+    reference_file = data_dir / "average_mesh.vtk"
     if not reference_file.exists():
         raise FileNotFoundError(
             f"KCL-Heart-Model reference mesh not found: {reference_file}\n"
             "See data/README.md for download instructions."
         )
 
-    sample_dir = data_dir / "sample_meshes"
-    sample_files = sorted(sample_dir.glob("*.vtu"))
+    sample_dir = data_dir / "input_meshes"
+    sample_files = sorted(sample_dir.glob("*.vtk"))
     if not sample_files:
-        sample_files = sorted(data_dir.glob("*.vtu"))
-    sample_files = [path for path in sample_files if path.name != "pca_mean.vtu"]
-    sample_files = sample_files[:max_samples]
+        sample_files = sorted(data_dir.glob("*.vtk"))
+        sample_files = [
+            path for path in sample_files if path.name != reference_file.name
+        ]
     if len(sample_files) < 3:
         raise FileNotFoundError(
             f"Need at least 3 sample meshes under {sample_dir} or {data_dir}.\n"
diff --git a/utils/ai_agent_github_reviews.py b/utils/ai_agent_github_reviews.py
index 7e760f8..2d3d530 100644
--- a/utils/ai_agent_github_reviews.py
+++ b/utils/ai_agent_github_reviews.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 ai_agent_github_reviews.py — Screen GitHub PR review comments with an AI agent.
 
diff --git a/utils/generate_api_map.py b/utils/generate_api_map.py
index 2b83380..e34e91e 100644
--- a/utils/generate_api_map.py
+++ b/utils/generate_api_map.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 Generate a structured API index of public Python symbols in the repository.
 
diff --git a/utils/setup_feature_worktree.py b/utils/setup_feature_worktree.py
index c7a38bc..5676b69 100644
--- a/utils/setup_feature_worktree.py
+++ b/utils/setup_feature_worktree.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 setup_feature_worktree.py — Automate creation of a Git feature worktree on Windows.