Eliminate tubetk

README.md

@@ -159,7 +159,7 @@ print(WorkflowConvertImageToUSD.__name__)
 
 ### Key Dependencies
 
-- **Medical Imaging**: ITK, TubeTK, MONAI, nibabel, PyVista
+- **Medical Imaging**: ITK, MONAI, nibabel, PyVista
 - **AI/ML**: PyTorch, CuPy (CUDA 13), transformers, MONAI
 - **Registration**: icon-registration, unigradicon
 - **Visualization**: USD-core, PyVista

docs/api/workflows.rst

@@ -166,7 +166,8 @@ High-Resolution 4D CT Reconstruction
        registration_method=RegisterImagesGreedyICON(),
    )
 
-   result = workflow.run_workflow(upsample_to_fixed_resolution=True)
+   workflow.set_upsample_to_fixed_resolution(True)
+   result = workflow.run_workflow()
 
 See Also
 ========

docs/cli_scripts/4dct_reconstruction.rst

@@ -27,16 +27,15 @@ Basic Usage
        --fixed-image highres_reference.mha \
        --output-dir ./results
 
-With Upsampling
-===============
+Choosing a Reference Frame
+===========================
 
 .. code-block:: bash
 
    physiomotion4d-reconstruct-highres-4d-ct \
        --time-series-images frame_000.mha frame_001.mha frame_002.mha \
        --fixed-image highres_reference.mha \
        --reference-frame 0 \
-       --upsample \
        --output-dir ./results
 
 Registration Options

docs/conf.py

@@ -26,7 +26,6 @@ def __getattr__(cls, name):
 
 
 # Mock modules that need special handling
-sys.modules["itk.TubeTK"] = Mock()
 sys.modules["icon_registration.losses"] = Mock()
 sys.modules["icon_registration.network_wrappers"] = Mock()
 

docs/examples.rst

@@ -69,7 +69,8 @@ phase images:
        registration_method=RegisterImagesGreedyICON(),
    )
 
-   result = workflow.run_workflow(upsample_to_fixed_resolution=True)
+   workflow.set_upsample_to_fixed_resolution(True)
+   result = workflow.run_workflow()
    reconstructed_images = result["reconstructed_images"]
 
 Segmentation Examples

docs/installation.rst

@@ -21,7 +21,7 @@ Software Dependencies
 
 PhysioMotion4D relies on several key packages:
 
-* **Medical Imaging**: ITK, TubeTK, MONAI, nibabel, PyVista
+* **Medical Imaging**: ITK, MONAI, nibabel, PyVista
 * **AI/ML**: PyTorch, CuPy (CUDA 13), transformers, MONAI
 * **Registration**: icon-registration, unigradicon
 * **Visualization**: USD-core, PyVista

experiments/Colormap-VTK_To_USD/colormap_vtk_to_usd.py

@@ -10,13 +10,6 @@
 # 2. **Custom intensity ranges**: Control value-to-color mapping
 # 3. **Point data visualization**: Map scalar data to colors on 3D meshes
 # 4. **Time-varying data**: Create animated USD files with colored meshes
-#
-# ## Requirements
-#
-# ```bash
-# pip install physiomotion4d pyvista numpy
-# ```
-
 # %% [markdown]
 # ## Setup and Imports
 
@@ -114,7 +107,7 @@ def create_example_mesh_with_data(time_step):
 # Convert to USD
 output_file = output_dir / "example1_plasma_auto.usd"
 converter.convert(str(output_file))
-print(f"\n✓ Created: {output_file}")
+print(f"\nCreated: {output_file}")
 
 # %% [markdown]
 # ## Example 2: Rainbow Colormap with Custom Range
@@ -141,7 +134,7 @@ def create_example_mesh_with_data(time_step):
 output_file = output_dir / "example2_rainbow_custom.usd"
 print("Creating file:", output_file)
 stage = converter.convert(str(output_file))
-print(f"✓ Created: {output_file}")
+print(f"Created: {output_file}")
 
 # %% [markdown]
 # ## Example 3: Heat Colormap for Temperature Data
@@ -166,7 +159,7 @@ def create_example_mesh_with_data(time_step):
 
 output_file = output_dir / "example3_heat_temperature.usd"
 stage = converter.convert(str(output_file))
-print(f"✓ Created: {output_file}")
+print(f"Created: {output_file}")
 
 # %% [markdown]
 # ## Example 4: Coolwarm (Diverging) Colormap
@@ -191,7 +184,7 @@ def create_example_mesh_with_data(time_step):
 
 output_file = output_dir / "example4_coolwarm_diverging.usd"
 stage = converter.convert(str(output_file))
-print(f"✓ Created: {output_file}")
+print(f"Created: {output_file}")
 
 # %% [markdown]
 # ## Example 5: Grayscale Colormap
@@ -214,7 +207,7 @@ def create_example_mesh_with_data(time_step):
 
 output_file = output_dir / "example5_grayscale.usd"
 stage = converter.convert(str(output_file))
-print(f"✓ Created: {output_file}")
+print(f"Created: {output_file}")
 
 # %% [markdown]
 # ## Example 6: Random Colormap for Categorical Data
@@ -243,7 +236,7 @@ def create_example_mesh_with_data(time_step):
 
 output_file = output_dir / "example6_random_categorical.usd"
 stage = converter.convert(str(output_file))
-print(f"✓ Created: {output_file}")
+print(f"Created: {output_file}")
 
 # %% [markdown]
 # ## Example 7: Method Chaining for Concise API Usage
@@ -269,7 +262,7 @@ def create_example_mesh_with_data(time_step):
     .convert(str(output_file))
 )
 
-print(f"✓ Created: {output_file}")
+print(f"Created: {output_file}")
 
 # %% [markdown]
 # ## Summary: Available Colormaps and Features

experiments/Convert_VTK_To_USD/convert_vtk_to_usd_using_class.py

@@ -300,7 +300,7 @@ def create_deformed_pv_mesh(
             )
 
 print("\n" + "=" * 60)
-print("✓ All conversions completed successfully!")
+print("All conversions completed successfully!")
 print("=" * 60)
 
 
@@ -318,15 +318,15 @@ def verify_usd_file(usd_path):
 
     stage = Usd.Stage.Open(str(usd_path))
     if not stage:
-        print("  ✗ Failed to open stage")
+        print("  FAILED: could not open stage")
         return False
 
     # Check default prim
     default_prim = stage.GetDefaultPrim()
     if not default_prim:
-        print("  ✗ No default prim")
+        print("  FAILED: no default prim")
         return False
-    print(f"  ✓ Default prim: {default_prim.GetPath()}")
+    print(f"  Default prim: {default_prim.GetPath()}")
 
     # Find mesh prims
     mesh_count = 0
@@ -336,13 +336,13 @@ def verify_usd_file(usd_path):
             mesh = UsdGeom.Mesh(prim)
             points = mesh.GetPointsAttr().Get()
             if points:
-                print(f"  ✓ Mesh '{prim.GetName()}': {len(points)} points")
+                print(f"  Mesh '{prim.GetName()}': {len(points)} points")
 
     if mesh_count == 0:
-        print("  ✗ No meshes found")
+        print("  FAILED: no meshes found")
         return False
 
-    print(f"  ✓ Total meshes: {mesh_count}")
+    print(f"  Total meshes: {mesh_count}")
     return True
 
 
@@ -358,9 +358,9 @@ def verify_usd_file(usd_path):
 
 print("\n" + "=" * 60)
 if all_valid:
-    print("✓ All USD files are valid!")
+    print("All USD files are valid!")
 else:
-    print("✗ Some USD files have issues")
+    print("FAILED: some USD files have issues")
 print("=" * 60)
 
 # %% [markdown]

experiments/LongitudinalRegistration/0-cardiacGatedCT_segment_and_landmark.py → experiments/Heart-GatedCT-OptimizedLongitudinalRegistration/0-cardiacGatedCT_segment_and_landmark.py

@@ -14,7 +14,7 @@
 #
 # Output files follow the input stem:
 # - `<stem>_labelmap.nii.gz`
-# - `<stem>_landmark.csv`
+# - `<stem>_landmark.mrk.json`
 #
 
 # %%
@@ -93,7 +93,7 @@ def segment_images(
             if not os.path.exists(labelmap_path) or not os.path.exists(landmark_path):
                 image_path = os.path.join(src_dir, f)
                 input_image = itk.imread(image_path, pixel_type=itk.F)
-                results = segmenter.segment(input_image, contrast_enhanced_study=False)
+                results = segmenter.segment(input_image)
                 labelmap = results["labelmap"]
                 itk.imwrite(labelmap, str(labelmap_path), compression=True)
 

experiments/LongitudinalRegistration/1-initial_registration.py → experiments/Heart-GatedCT-OptimizedLongitudinalRegistration/1-initial_registration.py

@@ -178,13 +178,19 @@ def crop_image_to_mask(
         print("No mask provided, using image as mask")
     else:
         mask_arr = itk.array_from_image(mask)
+    shape_z, shape_y, shape_x = mask_arr.shape
     bounding_box = np.where(mask_arr > 0)
-    min_x = np.min(bounding_box[2])
-    max_x = np.max(bounding_box[2])
-    min_y = np.min(bounding_box[1])
-    max_y = np.max(bounding_box[1])
-    min_z = np.min(bounding_box[0])
-    max_z = np.max(bounding_box[0])
+    if not mask_arr.any() or bounding_box[0].size == 0:
+        min_x, max_x = 0, shape_x - 1
+        min_y, max_y = 0, shape_y - 1
+        min_z, max_z = 0, shape_z - 1
+    else:
+        min_x = int(np.min(bounding_box[2]))
+        max_x = int(np.max(bounding_box[2]))
+        min_y = int(np.min(bounding_box[1]))
+        max_y = int(np.max(bounding_box[1]))
+        min_z = int(np.min(bounding_box[0]))
+        max_z = int(np.max(bounding_box[0]))
     margin_x = int((max_x - min_x) * margin_fraction)
     margin_y = int((max_y - min_y) * margin_fraction)
     margin_z = int((max_z - min_z) * margin_fraction)
@@ -194,25 +200,20 @@ def crop_image_to_mask(
     max_y += margin_y
     min_z -= margin_z
     max_z += margin_z
-    if min_x < 0:
-        min_x = 0
-    if min_y < 0:
-        min_y = 0
-    if min_z < 0:
-        min_z = 0
-    max_size = image.GetLargestPossibleRegion().GetSize()
-    if max_x >= max_size[0]:
-        max_x = max_size[0] - 1
-    if max_y >= max_size[1]:
-        max_y = max_size[1] - 1
-    if max_z >= max_size[2]:
-        max_z = max_size[2] - 1
+    min_x = max(0, min(min_x, shape_x - 1))
+    min_y = max(0, min(min_y, shape_y - 1))
+    min_z = max(0, min(min_z, shape_z - 1))
+    max_x = max(0, min(max_x, shape_x - 1))
+    max_y = max(0, min(max_y, shape_y - 1))
+    max_z = max(0, min(max_z, shape_z - 1))
     print(f"array shape: {mask_arr.shape}")
     print(
         f"min_x: {min_x}, max_x: {max_x}, min_y: {min_y}, max_y: {max_y}, min_z: {min_z}, max_z: {max_z}"
     )
     new_image_arr = itk.array_from_image(image)
-    new_image_arr = new_image_arr[min_z:max_z, min_y:max_y, min_x:max_x]
+    new_image_arr = new_image_arr[
+        min_z : max_z + 1, min_y : max_y + 1, min_x : max_x + 1
+    ]
     new_origin = image.TransformIndexToPhysicalPoint(
         [int(min_x), int(min_y), int(min_z)]
     )
@@ -223,15 +224,19 @@ def crop_image_to_mask(
 
     if labelmap is not None:
         new_labelmap_arr = itk.array_from_image(labelmap)
-        new_labelmap_arr = new_labelmap_arr[min_z:max_z, min_y:max_y, min_x:max_x]
+        new_labelmap_arr = new_labelmap_arr[
+            min_z : max_z + 1, min_y : max_y + 1, min_x : max_x + 1
+        ]
         new_labelmap = itk.image_from_array(new_labelmap_arr)
         new_labelmap.CopyInformation(new_image)
     else:
         new_labelmap = None
 
     if mask is not None:
         new_mask_arr = itk.array_from_image(mask)
-        new_mask_arr = new_mask_arr[min_z:max_z, min_y:max_y, min_x:max_x]
+        new_mask_arr = new_mask_arr[
+            min_z : max_z + 1, min_y : max_y + 1, min_x : max_x + 1
+        ]
         new_mask = itk.image_from_array(new_mask_arr)
         new_mask.CopyInformation(new_image)
     else:

experiments/LongitudinalRegistration/3-eval_icon.py → experiments/Heart-GatedCT-OptimizedLongitudinalRegistration/3-eval_icon.py

@@ -121,21 +121,16 @@ def _build_registrar(
     },
 ]
 
-output_dir.mkdir(parents=True, exist_ok=True)
-detail_file = output_dir / "landmark_errors_by_point.csv"
-summary_file = output_dir / "registration_summary.csv"
-warped_ref_detail_file = output_dir / "warped_ref_landmark_errors_by_point.csv"
-if detail_file.exists():
-    detail_file.unlink()
-if warped_ref_detail_file.exists():
-    warped_ref_detail_file.unlink()
-
-# %%
 all_patient_ids = sorted(
     p.name
     for p in timepoint_base_dir.iterdir()
     if p.is_dir() and p.name.startswith("pm00")
 )
+if len(all_patient_ids) < 2:
+    raise ValueError(
+        f"Expected at least 2 subjects under {timepoint_base_dir}, "
+        f"found {len(all_patient_ids)}."
+    )
 n_train = max(
     1, min(len(all_patient_ids) - 1, round(train_fraction * len(all_patient_ids)))
 )
@@ -146,6 +141,15 @@ def _build_registrar(
 )
 print(f"Held-out test subjects: {test_subjects}")
 
+output_dir.mkdir(parents=True, exist_ok=True)
+detail_file = output_dir / "landmark_errors_by_point.csv"
+summary_file = output_dir / "registration_summary.csv"
+warped_ref_detail_file = output_dir / "warped_ref_landmark_errors_by_point.csv"
+if detail_file.exists():
+    detail_file.unlink()
+if warped_ref_detail_file.exists():
+    warped_ref_detail_file.unlink()
+
 # %% [markdown]
 # ## 2. Validate that every test subject has exactly one reference file
 

experiments/LongitudinalRegistration/composite_time_series_mid_slice.py → experiments/Heart-GatedCT-OptimizedLongitudinalRegistration/composite_time_series_mid_slice.py

@@ -2,9 +2,7 @@
 
 import argparse
 import re
-import tkinter as tk
 from pathlib import Path
-from tkinter import filedialog
 from typing import Optional
 
 import itk
@@ -16,10 +14,21 @@
 
 def select_directory() -> Optional[Path]:
     """Open a directory chooser and return the selected directory."""
+    try:
+        import tkinter as tk
+        from tkinter import filedialog
+    except ImportError as exc:
+        raise RuntimeError(
+            "tkinter is unavailable; pass the directory path on the command line."
+        ) from exc
+
     try:
         root = tk.Tk()
-    except tk.TclError:
-        return None
+    except tk.TclError as exc:
+        raise RuntimeError(
+            "tkinter could not open a directory picker; pass the directory path "
+            "on the command line."
+        ) from exc
     try:
         root.withdraw()
         root.update()
@@ -149,7 +158,11 @@ def main(argv: Optional[list[str]] = None) -> int:
 
     input_dir = args.directory
     if input_dir is None:
-        input_dir = select_directory()
+        try:
+            input_dir = select_directory()
+        except RuntimeError as exc:
+            print(f"Error: {exc}")
+            return 1
     if input_dir is None:
         print("No directory selected")
         return 1

experiments/LongitudinalRegistration/registration_results_analysis.py → experiments/Heart-GatedCT-OptimizedLongitudinalRegistration/registration_results_analysis.py

@@ -64,7 +64,8 @@ def _split_by_occurrence(
         Ordered dict mapping group name → sub-data-frame.
     """
     occurrence = frame.groupby(key_cols, sort=False).cumcount() + 1  # 1-indexed
-    max_occ = int(occurrence.max())
+    max_occ_value = occurrence.max()
+    max_occ = 0 if frame.empty or pd.isna(max_occ_value) else int(max_occ_value)
     result: dict[str, pd.DataFrame] = {}
     for n in range(1, max_occ + 1):
         sub = frame[occurrence == n]