ENH: Update fit-to-patient for non-contrast heart models

data/DirLab-4DCT/README.md

@@ -8,10 +8,16 @@ DIR-Lab distributes each case individually and may require registration.
 
 1. Visit the DIR-Lab 4D-CT page and request/download the case archives:
    https://med.emory.edu/departments/radiation-oncology/research-laboratories/deformable-image-registration/downloads-and-reference-data/4dct.html
-2. Extract each case's raw images into `data/DirLab-4DCT/` (see layout
-   below). The `.mhd` headers already committed in this directory point at
-   those raw files, so no conversion step is required — see
-   "About the Committed `.mhd` Files" below.
+2. Extract each case's raw images into `data/DirLab-4DCT/downloaded_data/`
+   (see layout below). The `.mhd` headers already committed in that
+   directory point at those raw files, so no header-writing step is
+   required — see "About the Committed `.mhd` Files" below.
+3. Run `data/DirLab-4DCT/fix_downloaded_data.py`. DIR-Lab's raw volumes are
+   not in Hounsfield units (see "Fixing Raw Intensities" below); this
+   script reads every `.mhd` header in `downloaded_data/` with its backing
+   pixel data present, corrects the intensities, and writes one compressed
+   `.mha` volume per phase directly into `data/DirLab-4DCT/`. Tutorials and
+   experiments read these top-level `.mha` files, not the raw `.mhd` ones.
 
 Once populated, check the layout with:
 
@@ -21,26 +27,31 @@ from physiotwin4d import DataDownloadTools
 assert DataDownloadTools.VerifyDirLab4DCTData("data/DirLab-4DCT")
 ```
 
-**Directory structure after download:**
+**Directory structure after download and fixing:**
 ```
 data/DirLab-4DCT/
-├── Case1Pack/
-│   ├── Images/              # T00-T90 phase images
-│   ├── ExtremePhases/       # T00 and T50 (max inhale/exhale)
-│   └── Sampled4D/           # Sampled time points
-├── Case1Pack_T00.mhd        # Already-committed headers, per phase
-├── Case1Pack_T10.mhd
+├── downloaded_data/
+│   ├── Case1Pack/
+│   │   ├── Images/              # T00-T90 phase images
+│   │   ├── ExtremePhases/       # T00 and T50 (max inhale/exhale)
+│   │   └── Sampled4D/           # Sampled time points
+│   ├── Case1Pack_T00.mhd        # Already-committed headers, per phase
+│   ├── Case1Pack_T10.mhd
+│   ...
+│   ├── Case10Pack/
+│   └── Convert4DCTToMHD.py      # Documents how the .mhd headers were generated
+├── Case1Pack_T00.mha             # Written by fix_downloaded_data.py
+├── Case1Pack_T10.mha
 ...
-├── Case10Pack/
-├── Convert4DCTToMHD.py      # Documents how the .mhd headers were generated
+├── fix_downloaded_data.py
 └── README.md (this file)
 ```
 
 ### About the Committed `.mhd` Files
 
-The `Case*Pack_T*.mhd` files in this directory are already committed to the
-repository, but they are only MetaImage **headers** (a few hundred bytes
-each) — for example:
+The `Case*Pack_T*.mhd` files in `downloaded_data/` are already committed to
+the repository, but they are only MetaImage **headers** (a few hundred
+bytes each) — for example:
 
 ```
 ObjectType = Image
@@ -60,6 +71,16 @@ until you complete the manual download above and the referenced
 raw DIR-Lab archives. It is included for provenance/documentation only —
 you do not need to run it; the `.mhd` files are already committed.
 
+### Fixing Raw Intensities
+
+DIR-Lab's raw `.mhd`/`.img` volumes store intensities offset by +1024 from
+Hounsfield units, not real HU. `fix_downloaded_data.py` calls
+`DataDownloadTools.FixDirLab4DCTData`, which subtracts 1024 and clips the
+result to `[-1024, 1024]` for every downloaded phase, writing one
+compressed `.mha` file per phase into `data/DirLab-4DCT/` (`.gitignore`d,
+like the raw data). Run it once after downloading; tutorials and
+experiments assume the `.mha` files it produces already exist.
+
 ## Overview
 
 Benchmark dataset for 4D CT respiratory motion analysis. Contains 10 cases
@@ -104,9 +125,16 @@ If you use this dataset, please cite:
 
 ### Files in This Directory
 
-- `Convert4DCTToMHD.py` — documents how the committed `.mhd` headers were
-  generated from the raw DIR-Lab archives; not needed to run
-- `Case*Pack_T*.mhd` — MetaImage headers for each case/phase (see above)
+- `fix_downloaded_data.py` — converts `downloaded_data/*.mhd` to
+  HU-corrected `Case*Pack_T*.mha` files in this directory (see "Fixing Raw
+  Intensities" above)
+- `downloaded_data/Convert4DCTToMHD.py` — documents how the committed
+  `.mhd` headers were generated from the raw DIR-Lab archives; not needed
+  to run
+- `downloaded_data/Case*Pack_T*.mhd` — MetaImage headers for each
+  case/phase (see above)
+- `Case*Pack_T*.mha` — HU-corrected volumes written by
+  `fix_downloaded_data.py`; what tutorials and experiments actually read
 
 ### Additional Resources
 

data/DirLab-4DCT/fix_downloaded_data.py

@@ -0,0 +1,14 @@
+# %%
+from pathlib import Path
+
+from physiotwin4d.data_download_tools import DataDownloadTools
+
+_HERE = Path(__file__).resolve().parent
+
+# %%
+input_dir = _HERE / "downloaded_data"
+output_dir = _HERE
+
+# Converts raw DirLab-4DCT intensities (.mhd) to clipped Hounsfield units and
+# writes each result as a compressed .mha file in output_dir.
+DataDownloadTools.FixDirLab4DCTData(input_dir, output_dir)

docs/api/model_registration/index.rst

@@ -42,7 +42,7 @@ Most users should access model registration through the workflow:
        patient_image=itk.imread("patient_ct.nii.gz"),
    )
 
-   result = workflow.run_workflow()
+   result = workflow.process()
 
 See Also
 ========

docs/api/workflows.rst

@@ -147,7 +147,7 @@ Statistical Shape Modeling
        patient_image=itk.imread("patient_ct.nii.gz"),
    )
 
-   result = workflow.run_workflow()
+   result = workflow.process()
 
 High-Resolution 4D CT Reconstruction
 ====================================
@@ -172,7 +172,7 @@ High-Resolution 4D CT Reconstruction
    )
 
    workflow.set_upsample_to_fixed_resolution(True)
-   result = workflow.run_workflow()
+   result = workflow.process()
 
 See Also
 ========

docs/developer/registration_models.rst

@@ -22,7 +22,7 @@ Recommended Entry Point
        patient_image=itk.imread("patient_ct.nii.gz"),
    )
 
-   result = workflow.run_workflow()
+   result = workflow.process()
 
 Lower-Level Classes
 ===================

docs/developer/segmentation.rst

@@ -72,8 +72,9 @@ produces by populating ``self.taxonomy``. The base class owns the
            ...
 
 The base class then provides the shared pre/post-processing, contrast-agent
-fusion, and ``create_anatomy_group_masks()`` that walks the taxonomy and
-emits one mask per registered group plus ``"other"``.
+fusion, and ``create_anatomy_group_labelmaps()`` that walks the taxonomy and
+emits one labelmap (preserving the original label ids) per registered group
+plus ``"other"``.
 
 Custom Anatomy Looks
 ====================

docs/tutorials.rst

@@ -315,7 +315,7 @@ Inner API usage
       )
       if pca_model is not None:
           workflow.set_use_pca_registration(True, pca_model=pca_model)
-      result = workflow.run_workflow()
+      result = workflow.process()
       registered_surface = result["registered_template_model_surface"]
 
 Run
@@ -425,7 +425,7 @@ Inner API usage
           registration_method=registration_method,
       )
       workflow.set_modality("ct")
-      result = workflow.run_workflow()
+      result = workflow.process()
       reconstructed_images = result["reconstructed_images"]
 
 Run
@@ -490,9 +490,9 @@ Inner API usage
           labelmap_interior_object_ids=LABELMAP_INTERIOR_OBJECT_IDS,
       )
       ssm_fit_workflow.set_use_pca_registration(
-          use_pca_registration=True, pca_model=ssm_model, pca_uses_surface=False,
+          use_pca_registration=True, pca_model=ssm_model, use_surface=False,
       )
-      ssm_fit_workflow_result = ssm_fit_workflow.run_workflow()
+      ssm_fit_workflow_result = ssm_fit_workflow.process()
       ssm_pca_coefficients = ssm_fit_workflow.pca_coefficients
 
       icon_registration_method = RegisterImagesICON()
@@ -502,7 +502,7 @@ Inner API usage
           fixed_image=ref_image,
           registration_method=icon_registration_method,
       )
-      reg_result = reg_workflow.run_workflow()
+      reg_result = reg_workflow.process()
       reconstructed_images = reg_result["reconstructed_images"]
 
 Run

experiments/Heart-Create_Statistical_Model/README.md

@@ -161,7 +161,7 @@ workflow = WorkflowFitStatisticalModelToPatient(
     pca_number_of_modes=10
 )
 
-registered_mesh = workflow.run_workflow()
+registered_mesh = workflow.process()["registered_template_model_surface"]
 ```
 
 ## Requirements

experiments/Heart-GatedCT_To_USD/2-generate_segmentation.py

@@ -81,7 +81,7 @@
         )
     else:
         labelmap_image = itk.imread(os.path.join(output_dir, f"{outname}.all_mask.mha"))
-        result = seg.create_anatomy_group_masks(labelmap_image)
+        result = seg.create_anatomy_group_labelmaps(labelmap_image)
 
     lung_mask = result["lung"]
     heart_mask = result["heart"]

experiments/Heart-Statistical_Model_To_Patient/heart_model_to_patient-CHOPValve.py

@@ -67,7 +67,7 @@
 )
 
 # %%
-results = registrar.run_workflow()
+results = registrar.process()
 
 # %%
 registered_model = results["registered_template_model"]

experiments/Lung-GatedCT_To_USD/0-register_dirlab_4dct.py

@@ -98,11 +98,12 @@ def register_image(
     os.makedirs(output_dir, exist_ok=True)
 
     for case_name in case_names:
+        # .mha files are DirLab-4DCT data already converted to HU by
+        # data/DirLab-4DCT/fix_downloaded_data.py.
         fixed_image_filename = (
-            f"{input_dir}/{case_name}_T{fixed_image_num * 10:02d}.mhd"
+            f"{input_dir}/{case_name}_T{fixed_image_num * 10:02d}.mha"
         )
         fixed_image = itk.imread(fixed_image_filename)
-        fixed_image = DataDirLab4DCT().fix_image(fixed_image)
 
         print("Segmenting fixed image...")
         fixed_result = seg_image.segment(fixed_image)
@@ -151,9 +152,8 @@ def register_image(
         for image_num in images:
             if image_num != fixed_image_num:
                 moving_image = itk.imread(
-                    os.path.join(input_dir, f"{case_name}_T{image_num * 10:02d}.mhd")
+                    os.path.join(input_dir, f"{case_name}_T{image_num * 10:02d}.mha")
                 )
-                moving_image = DataDirLab4DCT().fix_image(moving_image)
 
                 print("***")
                 print(

experiments/Lung-GatedCT_To_USD/Experiment_SegReg.py

@@ -3,8 +3,6 @@
 
 import itk
 
-from data_dirlab_4d_ct import DataDirLab4DCT
-
 from physiotwin4d import RegisterImagesICON
 from physiotwin4d import SegmentChestTotalSegmentator
 
@@ -18,12 +16,10 @@
     _RESULTS_DIR = os.path.join(_HERE, "results_SegReg")
 
     # %%
-    fixed_image = DataDirLab4DCT().fix_image(
-        itk.imread(os.path.join(_DATA_DIR, "Case1Pack_T30.mhd"))
-    )
-    moving_image = DataDirLab4DCT().fix_image(
-        itk.imread(os.path.join(_DATA_DIR, "Case1Pack_T00.mhd"))
-    )
+    # .mha files are DirLab-4DCT data already converted to HU by
+    # data/DirLab-4DCT/fix_downloaded_data.py.
+    fixed_image = itk.imread(os.path.join(_DATA_DIR, "Case1Pack_T30.mha"))
+    moving_image = itk.imread(os.path.join(_DATA_DIR, "Case1Pack_T00.mha"))
 
     # %%
     # Register images

experiments/Lung-GatedCT_To_USD/data_dirlab_4d_ct.py

@@ -2,9 +2,13 @@
 class pmDataDirLab4dCT:
 This module contains the pmDataDirLab4DCT class, which is used to store the
 data for the DirLab 4DCT dataset.
-"""
 
-import itk
+DirLab-4DCT's raw ``.mhd``/``.img`` volumes are not in Hounsfield units; run
+``data/DirLab-4DCT/fix_downloaded_data.py`` (backed by
+``DataDownloadTools.FixDirLab4DCTData``) once to write corrected ``.mha``
+volumes before using this dataset. That script's output is what this
+directory's ``case_names`` are meant to be read from.
+"""
 
 
 class DataDirLab4DCT:
@@ -30,13 +34,3 @@ def __init__(self):
     def get_case_names(self) -> list[str]:
         """Get the case names"""
         return self.case_names
-
-    def fix_image(self, input_image: itk.image) -> itk.image:
-        """Fix DirLab_4DCT intensities to conform to HU"""
-        input_image_arr = itk.GetArrayViewFromImage(input_image)
-        input_image_arr -= 1024
-        input_image_arr = input_image_arr.clip(-1024, 3071)
-        new_image = itk.GetImageFromArray(input_image_arr)
-        new_image.CopyInformation(input_image)
-
-        return new_image

src/physiotwin4d/cli/fit_statistical_model_to_patient.py

@@ -262,7 +262,7 @@ def main() -> int:
         # Execute registration workflow
         print("\nStarting registration pipeline...")
         print("=" * 70)
-        result = workflow.run_workflow(
+        result = workflow.process(
             use_ICON_registration_refinement=args.use_ICON_refinement,
         )
 

src/physiotwin4d/cli/reconstruct_highres_4d_ct.py

@@ -312,7 +312,7 @@ def main() -> int:
         print("\nStarting reconstruction pipeline...")
         print("=" * 70)
         workflow.set_upsample_to_fixed_resolution(True)
-        result = workflow.run_workflow()
+        result = workflow.process()
 
         # Save results
         print("\n" + "=" * 70)

src/physiotwin4d/contour_tools.py

@@ -31,8 +31,8 @@ def __init__(self, log_level: int | str = logging.INFO):
         """
         super().__init__(class_name=self.__class__.__name__, log_level=log_level)
 
+    @staticmethod
     def extract_contours(
-        self,
         labelmap_image: itk.image,
     ) -> pv.PolyData:
         """
@@ -71,8 +71,8 @@ def extract_contours(
 
         return contours
 
+    @staticmethod
     def transform_contours(
-        self,
         contours: pv.PolyData,
         tfm: itk.Transform,
         with_deformation_magnitude: bool = False,
@@ -135,8 +135,8 @@ def merge_meshes(
 
         return merged_mesh, pv_meshes
 
+    @staticmethod
     def create_reference_image(
-        self,
         mesh: pv.DataSet,
         spatial_resolution: float = 0.5,
         buffer_factor: float = 0.25,
@@ -164,8 +164,8 @@ def create_reference_image(
         reference_image.Allocate()
         return reference_image
 
+    @staticmethod
     def create_mask_from_mesh(
-        self,
         mesh: pv.DataSet | pv.UnstructuredGrid,
         reference_image: itk.Image,
     ) -> itk.Image:
@@ -360,8 +360,8 @@ def create_distance_map(
 
         return distance_image
 
+    @staticmethod
     def create_deformation_field(
-        self,
         points: np.ndarray,
         point_displacements: np.ndarray,
         reference_image: itk.Image,