ENH: Per-organ tutorial parameters and distance-map registration fixes - #119
ENH: Per-organ tutorial parameters and distance-map registration fixes#119aylward wants to merge 2 commits into
Conversation
Introduce tutorials/parameters_lung_ct_dirlab.py and tutorials/parameters_heart_ct_kcl.py as the single source for each use case's mask dilation, distance-map saturation radius, PCA component counts, Greedy iteration schedule, segmenter class, and (heart only) the interior chamber label ids. Every tutorial that rasterizes or registers a distance map now reads the same values, so the maps a network is finetuned on match the maps it later infers over. No paths live in these modules; each tutorial keeps its own inputs and outputs. Add tutorials/tutorial_02_heart_distancemap_finetune_icon.py, which finetunes uniGradICON on heart distance maps built from the Duke-Heart-4DLabelmaps labelmaps with the chambers excluded. The heart needs its own run rather than reusing the lung weights: its registration mask is much tighter, so its distance maps saturate over a shorter radius and do not share an intensity distribution with the lung ones. Library fixes: - transform_tools.transform_image gains an explicit background_value. Resampling previously fell back to ITK's default of 0, which for CT is water, not air, so pre-warped moving images carried a false soft-tissue shell wherever they had no data. register_images_base now fills with -1000 HU for CT (exactly uniGradICON's window floor) on the image warp only; masks and labelmaps keep 0. - register_models_distance_maps composed the Greedy and ICON transforms in the wrong order. ITK's CompositeTransform applies back to front, so the residual must be added last in the forward and first in the inverse. Also drops two unconditional debug_*.nii.gz writes that crashed when mask_dilation_mm was 0. - register_from restores moving_image and clears moving_image_registered, so a registrar can be reused after an initialized run. - workflow_fit_statistical_model_to_patient grids the PCA field on a template-frame reference image rather than the patient image, and pads physically from spacing. - Default registrar for intensity registration switches from RegisterImagesGreedyICON/ICON to RegisterImagesGreedy in register_time_series_images and workflow_convert_image_to_usd. Distance-map registration keeps ICON, now with the finetuned weights. - convert_vtk_to_usd validates object names as USD identifiers and rejects duplicates; workflow_convert_vtk_to_usd wraps raw vtkDataSet. Rename number_of_components / number_of_modes to number_of_pca_components throughout the workflows and tutorials. tutorial_02_lung_finetune_icon now writes difference images (fixed minus registered) instead of the resampled volumes, and reports the chain's Greedy-stage-only score as its own row. The chain remains unconditional: on DIR-Lab, ICON's 175^3 residual grid is about 1.4 mm over the FOV, coarser than the 1.10 mm Greedy already achieves, so it cannot refine and the tutorial reports that honestly. tutorial_02_lung_distancemap_finetune_icon restricts its cached labelmaps to the lung labels. They previously held all 97 whole-body classes, and uniGradICON's Dice loss one-hots every shared class at 175^3 by batch 4, which saturated GPU memory. Baselines for the slow and GPU buckets will need refreshing: the composition-order fix and the PCA field grid change alter registration output.
WalkthroughThe PR changes registration defaults, adds modality-aware resampling and mesh conversion, centralizes tutorial parameters, and adds Duke heart labelmap, statistical-model, and PhysicsNeMo workflows. ChangesRegistration and mesh processing
Tutorial configuration and workflows
Estimated code review effort: 5 (Critical) | ~90 minutes Sequence Diagram(s)sequenceDiagram
participant DukeHeartLabelmaps
participant ContourTools
participant StatisticalModel
participant PhysicsNeMo
DukeHeartLabelmaps->>ContourTools: provide labeled frames
ContourTools->>StatisticalModel: provide reference surfaces
StatisticalModel->>PhysicsNeMo: provide fitted surfaces and displacement targets
Possibly related PRs
✨ Finishing Touches🧪 Generate unit tests (beta)
|
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## main #119 +/- ##
==========================================
+ Coverage 42.18% 44.56% +2.37%
==========================================
Files 72 72
Lines 8742 9039 +297
==========================================
+ Hits 3688 4028 +340
+ Misses 5054 5011 -43
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 10
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (1)
src/physiotwin4d/transform_tools.py (1)
452-475: 🗄️ Data Integrity & Integration | 🟠 Major | ⚡ Quick winUse modality-aware fill values for off-grid intensity resampling.
Keep
background_value=0.0for masks and labelmaps. Use_prewarp_background_value()where intensity data is resampled to avoid zero-filling CT data as water. UpdateRegisterImagesBase.get_registered_image()and the time-series reconstruction path, then add a regression test for off-grid CT fill values.🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@src/physiotwin4d/transform_tools.py` around lines 452 - 475, Update RegisterImagesBase.get_registered_image() and the time-series reconstruction path to pass _prewarp_background_value() for intensity-image resampling, while preserving background_value=0.0 for masks and labelmaps. Add a regression test verifying off-grid CT voxels receive the modality-appropriate fill value rather than zero.Source: Coding guidelines
🧹 Nitpick comments (3)
tutorials/parameters_heart_ct_kcl.py (1)
64-69: 📐 Maintainability & Code Quality | 🔵 Trivial | 💤 Low valueMake the two interior-id fields use the same type.
interior_object_ids_totalsegmentatorisOptional[list[int]]andinterior_object_ids_simplewareislist[int]. Both always hold a list. The consumer parameterlabelmap_interior_object_idsacceptsOptional[list], solist[int]works for both. Drop theOptionalfor consistency, and then theOptionalimport becomes unused.♻️ Proposed change
- interior_object_ids_totalsegmentator: Optional[list[int]] = field( + interior_object_ids_totalsegmentator: list[int] = field( default_factory=lambda: [141, 142, 143, 144] )🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@tutorials/parameters_heart_ct_kcl.py` around lines 64 - 69, Update the interior_object_ids_totalsegmentator field to use list[int], matching interior_object_ids_simpleware and the labelmap_interior_object_ids consumer, and remove the now-unused Optional import.src/physiotwin4d/workflow_fit_statistical_model_to_patient.py (1)
248-248: 📐 Maintainability & Code Quality | 🔵 Trivial | 💤 Low valueAdd
distancemap_squared_maxto the classAttributesdocstring.The class docstring at lines 68-110 lists the configurable state, including
mask_dilation_mm. The new attribute is missing there. Add one line so readers find the knob and its default derivation.🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@src/physiotwin4d/workflow_fit_statistical_model_to_patient.py` at line 248, Update the Attributes class docstring to document distancemap_squared_max alongside the existing configurable state, including that its default is derived as appropriate from the implementation. Leave the attribute declaration and surrounding behavior unchanged.tutorials/tutorial_02_duke_heart_distancemap_finetune_icon.py (1)
95-97: 📐 Maintainability & Code Quality | 🔵 Trivial | ⚡ Quick winConsider defaulting
run_finetuningtoFalse.
run_finetuning = Truecombined with lines 225-230 deletesexperiment_diron every run.tutorials/tutorial_07_heart_fit_statistical_model_to_patient.pyreads the checkpoint from that tree, so a re-run of this tutorial discards the weights Tutorial 7 depends on before it retrains them. The sibling tutorialtutorials/tutorial_02_lung_finetune_icon.pydocuments the opposite default: finetuning is off so runs reuse the checkpoint. Align the two defaults, or state in the module docstring why the heart run always retrains.🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@tutorials/tutorial_02_duke_heart_distancemap_finetune_icon.py` around lines 95 - 97, Default the run_finetuning setting in the tutorial to False, matching tutorials/tutorial_02_lung_finetune_icon.py so existing checkpoints are reused and experiment_dir is not deleted on reruns; retain the True option for explicitly requested retraining.
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Inline comments:
In `@docs/api/registration/chained.rst`:
- Around line 12-13: Update the documentation for RegisterImagesGreedyICON to
remove the incorrect claim that Tutorial 2 uses it. Either name a verified
consumer of the chained registrar or describe it only as the Greedy-then-ICON
pairing, while preserving the accurate statistical-model fit reference.
In `@docs/tutorials.rst`:
- Line 14: Update the tutorial summary near “Ten numbered stages across 17
runnable Python scripts” to avoid claiming all scripts are runnable: use “17
Python scripts” or explicitly distinguish the 16 runnable scripts from the
dataset-gated tutorial_02_duke_heart_distancemap_finetune_icon.py.
In `@src/physiotwin4d/register_images_base.py`:
- Around line 137-169: Add regression tests covering all branches of
_prewarp_background_value: an explicit override, CT returning -1000.0, and
non-CT using the moving image minimum. Add nearest-neighbor mask and labelmap
cases that verify zero fill, and assert the actual off-grid voxel value rather
than only output shape.
In `@src/physiotwin4d/transform_tools.py`:
- Around line 520-535: Update the background-value conversion in the resampling
flow before itk.resample_image_filter: for integer and discrete types, round to
an integer, validate it against np.iinfo(dtype), and reject out-of-range values
before invoking ITK; include boolean types via np.issubdtype(dtype, np.bool_).
Preserve floating-point conversion for non-discrete pixel types.
In `@src/physiotwin4d/workflow_reconstruct_highres_4d_ct.py`:
- Line 58: Update the class and constructor docstrings in
src/physiotwin4d/workflow_reconstruct_highres_4d_ct.py to consistently document
the standalone RegisterImagesGreedy backend and remove references to the retired
combined backend. In docs/tutorials.rst lines 313-321, remove Tutorial 3 from
Tutorial 2’s weight dependency; in lines 658-660, state that Tutorial 2 is
optional when stock weights are acceptable.
In `@tests/test_workflow_convert_image_to_usd.py`:
- Around line 121-122: Extend the migration test around workflow.registrar after
setting iterations to validate registration semantics, not just the
RegisterImagesGreedy type and artifact creation. Compare the registered image or
contours against an appropriate TestTools baseline, or assert a near-identity
result for the same-frame input, ensuring incorrect transform direction or
composition fails the test.
In `@tutorials/tutorial_02_duke_heart_distancemap_finetune_icon.py`:
- Around line 276-289: Update read_landmarks to validate the loaded markups
metadata before using control points: assert that the file’s coordinateSystem
field is exactly “LPS”, and fail clearly otherwise. Preserve the existing point
extraction only after this validation.
- Around line 198-217: After the training-case loop, validate that the collected
cohort is non-empty before the existing “Finetuning cohort” log and before
constructing WorkflowFinetuneICONRegistration. Mirror the held-out validation
behavior and fail immediately when subject_distance_map_files (and corresponding
subject IDs/labels) contains no surviving cases.
- Around line 139-151: The distance-map cache in distance_map_for currently keys
files only by frame stem, causing collisions across case directories. Include
the case directory name when constructing distance_map_file, while preserving
the existing derived_dir location and cache lookup behavior.
In `@tutorials/tutorial_07_heart_fit_statistical_model_to_patient.py`:
- Around line 164-169: Update the set_use_pca_registration call in the pca_model
block to pass the configured HEART_CT_KCL.pca_components(test_mode) value as
number_of_pca_components, matching the heart PCA builder and lung fit tutorial
instead of relying on the default.
---
Outside diff comments:
In `@src/physiotwin4d/transform_tools.py`:
- Around line 452-475: Update RegisterImagesBase.get_registered_image() and the
time-series reconstruction path to pass _prewarp_background_value() for
intensity-image resampling, while preserving background_value=0.0 for masks and
labelmaps. Add a regression test verifying off-grid CT voxels receive the
modality-appropriate fill value rather than zero.
---
Nitpick comments:
In `@src/physiotwin4d/workflow_fit_statistical_model_to_patient.py`:
- Line 248: Update the Attributes class docstring to document
distancemap_squared_max alongside the existing configurable state, including
that its default is derived as appropriate from the implementation. Leave the
attribute declaration and surrounding behavior unchanged.
In `@tutorials/parameters_heart_ct_kcl.py`:
- Around line 64-69: Update the interior_object_ids_totalsegmentator field to
use list[int], matching interior_object_ids_simpleware and the
labelmap_interior_object_ids consumer, and remove the now-unused Optional
import.
In `@tutorials/tutorial_02_duke_heart_distancemap_finetune_icon.py`:
- Around line 95-97: Default the run_finetuning setting in the tutorial to
False, matching tutorials/tutorial_02_lung_finetune_icon.py so existing
checkpoints are reused and experiment_dir is not deleted on reruns; retain the
True option for explicitly requested retraining.
🪄 Autofix
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: Organization UI
Review profile: CHILL
Plan: Pro Plus
Run ID: e3e3d426-db8d-4bc2-a2e9-01e6df24641a
📒 Files selected for processing (30)
data/Duke-Heart-4DLabelmaps/.gitignoredata/Duke-Heart-4DLabelmaps/README.mddocs/api/registration/chained.rstdocs/tutorials.rstpyproject.tomlsrc/physiotwin4d/register_images_base.pysrc/physiotwin4d/register_images_chain.pysrc/physiotwin4d/register_models_distance_maps.pysrc/physiotwin4d/register_time_series_images.pysrc/physiotwin4d/transform_tools.pysrc/physiotwin4d/workflow_convert_image_to_usd.pysrc/physiotwin4d/workflow_fit_statistical_model_to_patient.pysrc/physiotwin4d/workflow_reconstruct_highres_4d_ct.pytests/test_workflow_convert_image_to_usd.pytests/test_workflow_reconstruct_highres_4d_ct.pytutorials/README.mdtutorials/parameters_heart_ct_kcl.pytutorials/parameters_lung_ct_dirlab.pytutorials/tutorial_01_heart_gated_ct_to_usd.pytutorials/tutorial_01_lung_gated_ct_to_usd.pytutorials/tutorial_02_duke_heart_distancemap_finetune_icon.pytutorials/tutorial_02_lung_distancemap_finetune_icon.pytutorials/tutorial_02_lung_finetune_icon.pytutorials/tutorial_03_heart_reconstruct_highres_4d_ct.pytutorials/tutorial_03_lung_reconstruct_highres_4d_ct.pytutorials/tutorial_06_heart_create_statistical_model.pytutorials/tutorial_06_lung_create_statistical_model.pytutorials/tutorial_07_heart_fit_statistical_model_to_patient.pytutorials/tutorial_07_lung_fit_statistical_model_to_patient.pytutorials/tutorial_08_lung_fit_model_to_4d_patients.py
| ``RegisterImagesGreedyICON`` is the Greedy-then-ICON pairing, used by Tutorial 2 | ||
| and by the distance-map stage of the statistical-model fit. |
There was a problem hiding this comment.
📐 Maintainability & Code Quality | 🟡 Minor | ⚡ Quick win
🧩 Analysis chain
🏁 Script executed:
#!/bin/bash
set -euo pipefail
echo "Files matching chained.rst:"
fd -a 'chained\.rst$' . || true
echo
echo "Target file excerpt:"
if [ -f docs/api/registration/chained.rst ]; then
nl -ba docs/api/registration/chained.rst | sed -n '1,80p'
fi
echo
echo "Search for RegisterImagesGreedyICON/RegisterModelsDistanceMaps/Tutorial 2:"
rg -n "RegisterImagesGreedyICON|RegisterModelsDistanceMaps|RegisterImagesGreedy|RegisterImagesICON|Tutorial 2|Tutorial.*2|distance-map|distance map" . || trueRepository: Project-MONAI/physiotwin4d
Length of output: 308
🏁 Script executed:
#!/bin/bash
set -euo pipefail
echo "Repository files relevant to registration:"
git ls-files | rg 'registration|tutorial|models|fit|distance' || true
echo
echo "All definitions/usages with some context:"
rg -n "class RegisterImagesGreedyICON|RegisterImagesGreedyICON\s*=|class RegisterModelsDistanceMaps|RegisterModelsDistanceMaps\s*=|RegisterImagesGreedy|RegisterImagesICON" . || true
echo
echo "Candidate files by name:"
git ls-files | rg 'registration|tutorial|statistical|model|distance' | head -200 || trueRepository: Project-MONAI/physiotwin4d
Length of output: 40451
Correct the documented consumers of RegisterImagesGreedyICON.
Tutorial 2 distance-map stages use separate RegisterImagesGreedy and RegisterImagesICON instances; they do not use RegisterImagesGreedyICON. Replace the Tutorial 2 claim with a consumer that actually uses the chained registrar, or keep the wording as the documented Greedy-then-ICON pattern.
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@docs/api/registration/chained.rst` around lines 12 - 13, Update the
documentation for RegisterImagesGreedyICON to remove the incorrect claim that
Tutorial 2 uses it. Either name a verified consumer of the chained registrar or
describe it only as the Greedy-then-ICON pairing, while preserving the accurate
statistical-model fit reference.
| assert isinstance(workflow.registrar, RegisterImagesGreedy) | ||
| workflow.registrar.set_number_of_iterations([2]) |
There was a problem hiding this comment.
🗄️ Data Integrity & Integration | 🟠 Major | 🏗️ Heavy lift
Add a semantic registration assertion to this migration test.
The test verifies the Greedy type and artifact creation, but it does not verify registration correctness. Because the backend and transform composition changed, an incorrect transform direction can still pass. Compare the registered image or contours with a TestTools baseline, or assert a near-identity result for the same-frame input.
As per coding guidelines: “When a test produces an image or surface, compare it with a baseline using utilities such as TestTools.”
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@tests/test_workflow_convert_image_to_usd.py` around lines 121 - 122, Extend
the migration test around workflow.registrar after setting iterations to
validate registration semantics, not just the RegisterImagesGreedy type and
artifact creation. Compare the registered image or contours against an
appropriate TestTools baseline, or assert a near-identity result for the
same-frame input, ensuring incorrect transform direction or composition fails
the test.
Source: Coding guidelines
ContourTools gains the mesh generation the labelmap tutorials need:
extract_watertight_surface (pad + blur + flying edges, outward normals),
extract_label_surfaces (per-label surfaces that conform at shared walls via
signed distance maps), extract_tetrahedra + trim_tetrahedra_to_surface,
remesh_and_smooth_surface (ACVD clustering + Taubin), and
apply_anatomy_color. extract_contours now resamples anisotropic labelmaps
onto an isotropic grid and Taubin-smooths, so contours no longer terrace on
the voxel pitch. Adds pyacvd dependency.
USD export: WorkflowConvertVTKToUSD takes label_names/segmenter and splits
each mesh on its per-cell label array, giving one prim per structure at
/World/{project}/{group}/{structure} so identity survives a time series.
Labels now take precedence over static_merge (combining both raises).
UsdMeshConverter time-samples faceVertexCounts/Indices when frames disagree
on topology, instead of silently reusing frame 0's triangulation.
USDTools.get_mesh_paths descends the full subtree.
Registration: out-of-FOV voxels fill with the modality background value
rather than 0 in get_registered_image and the 4D reconstruction;
transform_tools validates background_value against the pixel type range and
handles bool images.
ImageTools.resample_image_by_scale resamples by voxel-count multiplier while
preserving physical extent.
Breaking: surface_target_reduction renamed to surface_reduction_rate
(WorkflowConvertImageToVTK.process, convert-image-to-vtk CLI flag);
WorkflowConvertImageToUSD gains surface_reduction_rate.
Tutorials: new Duke heart labelmap track (parameters module plus tutorials
04-10), per-organ parameter modules updated, docs/tutorials.rst and
tutorials/README.md refreshed, LFS + gitignore rules for the lung MGN
weights.
Tests: test_contour_mesh_extraction.py, test_register_images_base.py, plus
coverage for labeled and topology-varying USD export.
There was a problem hiding this comment.
Actionable comments posted: 16
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (1)
docs/tutorials.rst (1)
232-236: 📐 Maintainability & Code Quality | 🟡 Minor | ⚡ Quick winName
parameters_duke_heart_labelmaps.pyfor the Duke heart tutorials. The Duke distance-map tutorial imports this module, not the lung or KCL heart parameter modules.🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@docs/tutorials.rst` around lines 232 - 236, Update the Duke heart tutorial documentation near the referenced parameter-module descriptions to name parameters_duke_heart_labelmaps.py as the module used by the Duke distance-map tutorial, replacing the unrelated lung and KCL heart module references while preserving the existing availability note.
🧹 Nitpick comments (4)
tests/test_image_tools.py (1)
521-538: 🎯 Functional Correctness | 🔵 Trivial | ⚡ Quick winStrengthen the nearest-neighbor assertion and cover the upsampling border.
set(np.unique(...)) <= {0.0, 7.0}also passes when the output is entirely zeros, so it does not prove that the block survived. Add an assertion that7.0is present.No test covers
scale > 1.0output values. That is the path whereresample_image_by_scalesamples before the first input voxel center andResampleImageFilterwrites its default 0 (see the comment onsrc/physiotwin4d/image_tools.pylines 316-326).💚 Proposed test changes
out = image_tools.resample_image_by_scale(itk_image, 0.5, interpolate=False) - assert set(np.unique(itk.array_from_image(out))) <= {0.0, 7.0} + values = set(np.unique(itk.array_from_image(out))) + assert values <= {0.0, 7.0} + assert 7.0 in values, "nearest-neighbor resampling dropped the block" + + def test_upsampling_border_is_not_zero_filled( + self, image_tools: ImageTools + ) -> None: + """Upsampling must not write a background shell into a constant image.""" + arr = np.full((4, 4, 4), -1000.0, dtype=np.float32) + itk_image = _make_synthetic_itk_image((4, 4, 4), arr=arr) + + out = image_tools.resample_image_by_scale(itk_image, 2.0) + + assert np.allclose(itk.array_from_image(out), -1000.0)🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@tests/test_image_tools.py` around lines 521 - 538, Strengthen test_nearest_neighbor_keeps_input_values by asserting that 7.0 is present in the unique output values in addition to restricting values to {0.0, 7.0}. Add coverage for a scale greater than 1.0 that verifies resample_image_by_scale preserves the expected input intensities while allowing the documented default zero values at the upsampling border.src/physiotwin4d/workflow_convert_image_to_vtk.py (1)
316-334: 📐 Maintainability & Code Quality | 🔵 Trivial | ⚡ Quick winTwo reduction algorithms now share the name
surface_reduction_rate.
ContourTools.remesh_and_smooth_surfacereplaceddecimate_prowith ACVD remeshing, and its docstring states thatdecimate_proleaves a watertight input non-watertight. This workflow still callsdecimate_prodirectly for both the group surface and the label surfaces. The same parameter therefore selects a different algorithm depending on the entry point.Route both reductions through
self._contour_tools.remesh_and_smooth_surfaceso the parameter has one meaning.♻️ Proposed refactor
export_surface = base_surface if surface_reduction_rate > 0.0: - export_surface = export_surface.decimate_pro( - surface_reduction_rate, preserve_topology=True - ) + export_surface = self._contour_tools.remesh_and_smooth_surface( + export_surface, surface_reduction_rate + )if surface_reduction_rate > 0.0: - label_surface = label_surface.decimate_pro( - surface_reduction_rate, preserve_topology=True - ) + label_surface = self._contour_tools.remesh_and_smooth_surface( + label_surface, surface_reduction_rate + )Update the
processdocstring at lines 236-239 accordingly.🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@src/physiotwin4d/workflow_convert_image_to_vtk.py` around lines 316 - 334, In the process flow, replace both direct decimate_pro calls for export_surface and label_surface with self._contour_tools.remesh_and_smooth_surface, preserving the existing surface_reduction_rate control and topology behavior through that helper. Update the process docstring to describe the parameter as ACVD remeshing rather than decimation.src/physiotwin4d/contour_tools.py (2)
147-157: 🚀 Performance & Scalability | 🔵 Trivial | ⚖️ Poor tradeoffConsider making the isotropic resample opt-out or bounded.
The resample always targets the finest spacing. On a 0.7 x 0.7 x 2.5 mm CT the voxel count grows by about 3.6x on the slice axis, and the triangle count grows with it. Existing callers such as
WorkflowConvertImageToVTK._extract_surfaceandworkflow_convert_image_to_usd.py(line 363) get this cost without any new argument.A parameter (target spacing or a disable flag) would let a caller cap the cost on large chest volumes.
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@src/physiotwin4d/contour_tools.py` around lines 147 - 157, Add an opt-out or configurable target-spacing parameter to the contouring flow around the anisotropy check and _resample_labelmap_isotropic, allowing callers such as _extract_surface to avoid or cap expensive finest-spacing resampling. Preserve the current finest-spacing behavior by default, while ensuring the selected target spacing is isotropic and does not increase resolution beyond the configured limit.
667-667: 🎯 Functional Correctness | 🔵 Trivial | ⚡ Quick winGuard the tetrahedra-only assumption on the clipped mesh.
cells_dict[np.uint8(pv.CellType.TETRA)]raisesKeyErrorwhen the clipped mesh holds no tetrahedra, and it silently ignores non-tetrahedral cells when the mesh is mixed.qualityat line 716 covers every cell, soconnectivity[below_bound]then indexes the tetrahedra array with cell positions that belong to other cell types.The method is public and documents only
tetrahedraas the input, so a mixed-cell caller gets a wrong damping set rather than an error. Reject a non-tetrahedral input explicitly.♻️ Proposed guard
+ if set(relaxed.celltypes) != {pv.CellType.TETRA}: + raise ValueError( + "trim_tetrahedra_to_surface requires a tetrahedral mesh; got cell " + f"types {sorted(set(relaxed.celltypes))}" + ) connectivity = relaxed.cells_dict[np.uint8(pv.CellType.TETRA)]🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@src/physiotwin4d/contour_tools.py` at line 667, Update the method containing connectivity and quality so it validates that the clipped mesh contains only tetrahedral cells before accessing cells_dict[np.uint8(pv.CellType.TETRA)]. Reject empty or mixed-cell meshes explicitly with an appropriate error, preserving the existing damping logic only for tetrahedra-only input.
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Inline comments:
In `@src/physiotwin4d/contour_tools.py`:
- Around line 491-499: Update the is_watertight method to explicitly return
False when the triangulated surface has no faces, before computing edge counts;
preserve the existing exactly-two-shared-edges check for non-empty surfaces.
- Around line 530-542: Update extract_tetrahedra around the mask_arr
bounding-box calculation to detect when mask_arr contains no non-zero voxels
before calling axis_extent.min(). Return the same empty pv.UnstructuredGrid used
by the existing analogous empty path, while preserving the current cropping and
tetrahedra generation for non-empty masks.
- Around line 369-388: In the label-processing loop around extract_surface,
check whether the extracted and triangulated surface has zero cells before
removing data or computing normals. For empty surfaces, skip adding the label to
surfaces and emit a warning; retain the existing cleanup, normal computation,
and debug logging for labels with extracted cells.
In `@src/physiotwin4d/convert_vtk_to_usd.py`:
- Around line 662-669: Update the labeled-mesh dispatch around
_convert_with_labels so a static merge with exactly one labeled mesh uses
create_mesh() rather than create_time_varying_mesh(), while preserving the
existing labeled conversion for non-static or multi-mesh cases. Add a workflow
test assertion confirming the labeled static output contains no point or
topology time samples.
In `@src/physiotwin4d/image_tools.py`:
- Around line 316-326: The resample_image_by_scale path must support
configurable border fill values: add a background_value parameter following the
transform_image convention and apply it with SetDefaultPixelValue before
Update(). In src/physiotwin4d/image_tools.py lines 316-326, make this change; in
tests/test_image_tools.py lines 521-538, add a scale greater than 1 test using a
constant-valued image that asserts every output voxel retains that value, and
ensure the nearest-neighbor test still asserts the 7.0 block value.
In `@tests/test_contour_mesh_extraction.py`:
- Around line 47-70: Update the docstrings of _left_handed_box and
_touching_boxes to explicitly state their synthetic volume dimensions: 6 × 6 × 6
and 60 × 48 × 14 voxels respectively, while preserving the existing fixture
descriptions.
In `@tests/test_register_images_base.py`:
- Around line 1-5: Update the module docstring for
tests/test_register_images_base.py to explicitly state that the synthetic images
are 10×10×10 voxels, while preserving its existing description of the off-grid
resampling scenario.
In `@tutorials/README.md`:
- Around line 37-51: Update the tutorial index’s recommended-order documentation
to include the Duke heart dependency chain from Tutorial 4 through Tutorial 10,
or explicitly mark the Duke sequence as unavailable until Duke-Heart-4DLabelmaps
is released. Ensure the descriptions for the Duke Tutorials 8–10 do not imply
they use DIR-Lab inputs.
In `@tutorials/tutorial_04_duke_heart_labelmap_to_vtk.py`:
- Around line 189-198: Guard the volume percentage calculation in the
reporter.log_debug call within the mesh refinement loop so surface.volume equal
to zero cannot cause division by zero. Use a safe fallback for the logged ratio
while preserving the existing percentage calculation for positive surface
volumes.
- Around line 263-272: Update the whole_heart extraction flow around
contour_tools.extract_label_surfaces to store its returned mapping, check
whether label id 1 exists before indexing it, and skip the current frame with a
warning when the mapping is empty or lacks that label. Preserve the existing
coloring and annotation behavior for valid whole-heart surfaces.
- Around line 234-244: Update the per-case output handling in the loop over
case_dir so case_output_dir is created as output_dir / case_dir.name rather than
output_dir itself. Keep directory creation and surface/whole-heart file
generation using this case-specific directory, producing the pm????/ layout
expected by Tutorial 5.
In `@tutorials/tutorial_05_duke_heart_vtk_to_usd.py`:
- Around line 69-77: Restore dynamic case discovery by assigning case_dirs to
the sorted directories matching input_dir.glob("pm[0-9][0-9][0-9][0-9]"),
filtering to directories, and remove the hard-coded pm0004 assignment and
commented-out code. Keep the existing FileNotFoundError guard so missing cohorts
fail as documented.
In `@tutorials/tutorial_06_duke_heart_create_statistical_model.py`:
- Around line 88-92: The Tutorial 4 surface consumers use the wrong output
layout. In tutorials/tutorial_06_duke_heart_create_statistical_model.py lines
88-92, update the sample_files glob to search
*_ref_heart_minus_interior_chambers.vtp within pm???? case directories while
preserving the hold-out filter; in
tutorials/tutorial_08_duke_heart_fit_model_to_4d_patients.py lines 189-192,
include labelmap_file.parent.name between tutorial_04_dir and the surface
filename when constructing the surface path.
- Around line 117-120: Update the template_surface argument in
WorkflowCreateMeanSurface to use a valid entry from sample_surfaces, such as the
middle element, while preserving the existing minimum-count check and the
remaining workflow configuration.
In `@tutorials/tutorial_07_duke_heart_fit_statistical_model_to_patient.py`:
- Around line 114-118: Update the PCA model loading logic around pca_model and
pca_json so that pca_json is required whenever pca_mean_file exists; when it is
absent, raise the same prerequisite error used for the missing mean surface
instead of continuing with pca_model unset. Preserve JSON loading when the file
exists and ensure fitting cannot proceed in a non-PCA mode under these
conditions.
In `@tutorials/tutorial_10_duke_heart_infer_physicsnemo.py`:
- Around line 130-135: Clamp test_index in the frame-selection logic before
indexing stages and phase_files, ensuring values such as stage_fraction=1.0
resolve to the final frame rather than len(stages). Preserve the existing
fractional selection behavior for in-range values and use the last valid index
as the upper bound.
---
Outside diff comments:
In `@docs/tutorials.rst`:
- Around line 232-236: Update the Duke heart tutorial documentation near the
referenced parameter-module descriptions to name
parameters_duke_heart_labelmaps.py as the module used by the Duke distance-map
tutorial, replacing the unrelated lung and KCL heart module references while
preserving the existing availability note.
---
Nitpick comments:
In `@src/physiotwin4d/contour_tools.py`:
- Around line 147-157: Add an opt-out or configurable target-spacing parameter
to the contouring flow around the anisotropy check and
_resample_labelmap_isotropic, allowing callers such as _extract_surface to avoid
or cap expensive finest-spacing resampling. Preserve the current finest-spacing
behavior by default, while ensuring the selected target spacing is isotropic and
does not increase resolution beyond the configured limit.
- Line 667: Update the method containing connectivity and quality so it
validates that the clipped mesh contains only tetrahedral cells before accessing
cells_dict[np.uint8(pv.CellType.TETRA)]. Reject empty or mixed-cell meshes
explicitly with an appropriate error, preserving the existing damping logic only
for tetrahedra-only input.
In `@src/physiotwin4d/workflow_convert_image_to_vtk.py`:
- Around line 316-334: In the process flow, replace both direct decimate_pro
calls for export_surface and label_surface with
self._contour_tools.remesh_and_smooth_surface, preserving the existing
surface_reduction_rate control and topology behavior through that helper. Update
the process docstring to describe the parameter as ACVD remeshing rather than
decimation.
In `@tests/test_image_tools.py`:
- Around line 521-538: Strengthen test_nearest_neighbor_keeps_input_values by
asserting that 7.0 is present in the unique output values in addition to
restricting values to {0.0, 7.0}. Add coverage for a scale greater than 1.0 that
verifies resample_image_by_scale preserves the expected input intensities while
allowing the documented default zero values at the upsampling border.
🪄 Autofix
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: Organization UI
Review profile: CHILL
Plan: Pro Plus
Run ID: 80e9d1ca-cbfb-4d02-9e8f-fe45d0f90ca8
⛔ Files ignored due to path filters (3)
docs/assets/tutorial_04_duke_heart.pngis excluded by!**/*.pngdocs/assets/tutorial_04_heart-2png.pngis excluded by!**/*.pngdocs/assets/tutorial_04_lung-2.pngis excluded by!**/*.png
📒 Files selected for processing (47)
.gitattributes.gitignoredocs/tutorials.rstpyproject.tomlsrc/physiotwin4d/cli/convert_image_to_vtk.pysrc/physiotwin4d/contour_tools.pysrc/physiotwin4d/convert_vtk_to_usd.pysrc/physiotwin4d/image_tools.pysrc/physiotwin4d/register_images_base.pysrc/physiotwin4d/register_time_series_images.pysrc/physiotwin4d/transform_tools.pysrc/physiotwin4d/usd_tools.pysrc/physiotwin4d/vtk_to_usd/usd_mesh_converter.pysrc/physiotwin4d/workflow_convert_image_to_usd.pysrc/physiotwin4d/workflow_convert_image_to_vtk.pysrc/physiotwin4d/workflow_convert_vtk_to_usd.pysrc/physiotwin4d/workflow_fit_statistical_model_to_patient.pysrc/physiotwin4d/workflow_reconstruct_highres_4d_ct.pytests/test_contour_mesh_extraction.pytests/test_convert_vtk_to_usd.pytests/test_image_tools.pytests/test_register_images_base.pytests/test_tutorials.pytests/test_workflow_convert_vtk_to_usd.pytutorials/README.mdtutorials/parameters_duke_heart_labelmaps.pytutorials/parameters_heart_ct_kcl.pytutorials/parameters_lung_ct_dirlab.pytutorials/tutorial_01_heart_gated_ct_to_usd.pytutorials/tutorial_01_lung_gated_ct_to_usd.pytutorials/tutorial_02_duke_heart_distancemap_finetune_icon.pytutorials/tutorial_04_duke_heart_labelmap_to_vtk.pytutorials/tutorial_04_heart_ct_to_vtk.pytutorials/tutorial_04_lung_ct_to_vtk.pytutorials/tutorial_05_duke_heart_vtk_to_usd.pytutorials/tutorial_06_duke_heart_create_statistical_model.pytutorials/tutorial_06_heart_create_statistical_model.pytutorials/tutorial_06_lung_create_statistical_model.pytutorials/tutorial_07_duke_heart_fit_statistical_model_to_patient.pytutorials/tutorial_07_heart_fit_statistical_model_to_patient.pytutorials/tutorial_07_lung_fit_statistical_model_to_patient.pytutorials/tutorial_08_duke_heart_fit_model_to_4d_patients.pytutorials/tutorial_08_lung_fit_model_to_4d_patients.pytutorials/tutorial_09_duke_heart_train_physicsnemo_mgn.pytutorials/tutorial_09_lung_train_physicsnemo_mgn.pytutorials/tutorial_10_duke_heart_infer_physicsnemo.pytutorials/tutorial_10_lung_infer_physicsnemo_mgn.py
🚧 Files skipped from review as they are similar to previous changes (10)
- tutorials/tutorial_07_heart_fit_statistical_model_to_patient.py
- tutorials/tutorial_01_heart_gated_ct_to_usd.py
- tutorials/tutorial_01_lung_gated_ct_to_usd.py
- tutorials/tutorial_08_lung_fit_model_to_4d_patients.py
- src/physiotwin4d/transform_tools.py
- src/physiotwin4d/register_time_series_images.py
- src/physiotwin4d/workflow_fit_statistical_model_to_patient.py
- tutorials/tutorial_02_duke_heart_distancemap_finetune_icon.py
- src/physiotwin4d/register_images_base.py
- src/physiotwin4d/workflow_reconstruct_highres_4d_ct.py
| for label_id in label_ids: | ||
| cell_ids: list[int] = np.flatnonzero(merged_ids == label_id).tolist() | ||
| surface = self.extract_surface(merged.extract_cells(cell_ids)).triangulate() | ||
| # The bookkeeping arrays of the merge and the split; the label is | ||
| # the key of the returned mapping, so it is not data on the mesh. | ||
| for array_name in ("LabelId", "vtkOriginalCellIds", "vtkOriginalPointIds"): | ||
| surface.cell_data.pop(array_name, None) | ||
| surface.point_data.pop(array_name, None) | ||
| # VTK winds faces for a right-handed direction matrix, so an LPS | ||
| # image with a negative-determinant direction comes out inward. | ||
| surfaces[label_id] = surface.compute_normals( | ||
| auto_orient_normals=True, consistent_normals=True | ||
| ) | ||
| self.log_debug( | ||
| "Label %d: %d triangles, %.3f of its voxel volume", | ||
| label_id, | ||
| surfaces[label_id].n_cells, | ||
| float(surfaces[label_id].volume) | ||
| / (int(np.count_nonzero(labels == label_id)) * voxel_volume), | ||
| ) |
There was a problem hiding this comment.
🩺 Stability & Availability | 🟡 Minor | ⚡ Quick win
🧩 Analysis chain
🏁 Script executed:
#!/bin/bash
# Description: Check for existing empty-mesh guards and tests around extract_label_surfaces.
rg -n -C4 'extract_label_surfaces' --type=py
rg -n -C3 'n_cells == 0|n_cells > 0' src/physiotwin4d/contour_tools.pyRepository: Project-MONAI/physiotwin4d
Length of output: 164
🏁 Script executed:
#!/bin/bash
set -e
printf '%s\n' '--- contour_tools.py outline ---'
ast-grep outline src/physiotwin4d/contour_tools.py --view expanded
printf '%s\n' '--- relevant implementation ---'
sed -n '260,405p' src/physiotwin4d/contour_tools.py
printf '%s\n' '--- surface extraction and dependency references ---'
rg -n -C5 'def extract_surface|extract_surface\(|pyvista|PyVista|pyacvd' src pyproject.toml tests 2>/dev/null || trueRepository: Project-MONAI/physiotwin4d
Length of output: 50383
🏁 Script executed:
#!/bin/bash
set -u
python3 - <<'PY'
import numpy as np
try:
import pyvista as pv
except Exception as exc:
print(f"IMPORT_ERROR: {type(exc).__name__}: {exc}")
raise SystemExit(0)
print(f"PYVISTA_VERSION: {getattr(pv, '__version__', 'unknown')}")
# Reproduce the empty-label branch with a standalone PyVista mesh.
empty = pv.PolyData()
for name, operation in [
("extract_cells", lambda: empty.extract_cells([])),
("triangulate", lambda: empty.triangulate()),
("compute_normals", lambda: empty.compute_normals(
auto_orient_normals=True, consistent_normals=True
)),
]:
try:
result = operation()
print(
f"{name}: OK type={type(result).__name__} "
f"n_points={result.n_points} n_cells={result.n_cells}"
)
except Exception as exc:
print(f"{name}: ERROR {type(exc).__name__}: {exc}")
# Reproduce the complete selection sequence used by the review target.
try:
points = np.array(
[[0, 0, 0], [1, 0, 0], [0, 1, 0]], dtype=float
)
nonempty = pv.PolyData(points, faces=np.array([3, 0, 1, 2]))
nonempty.cell_data["LabelId"] = np.array([1], dtype=np.int32)
empty_label = pv.PolyData()
empty_label.cell_data["LabelId"] = np.array([], dtype=np.int32)
merged = pv.merge([nonempty, empty_label], merge_points=True)
print(
f"merged: n_points={merged.n_points} n_cells={merged.n_cells} "
f"LabelId={np.asarray(merged.cell_data['LabelId']).tolist()}"
)
merged_ids = np.asarray(merged.cell_data["LabelId"])
cell_ids = np.flatnonzero(merged_ids == 2).tolist()
selected = merged.extract_cells(cell_ids)
surface = selected.triangulate()
result = surface.compute_normals(
auto_orient_normals=True, consistent_normals=True
)
print(
f"missing_label_sequence: OK cell_ids={cell_ids} "
f"type={type(result).__name__} n_points={result.n_points} "
f"n_cells={result.n_cells}"
)
except Exception as exc:
print(
f"missing_label_sequence: ERROR {type(exc).__name__}: {exc}"
)
PYRepository: Project-MONAI/physiotwin4d
Length of output: 282
🏁 Script executed:
#!/bin/bash
set -u
printf '%s\n' '--- available repository-local Python environments ---'
find .. -maxdepth 3 -type f \( -path '*/.venv/bin/python*' -o -path '*/venv/bin/python*' -o -path '*/.venv/Scripts/python.exe' -o -path '*/venv/Scripts/python.exe' \) -print 2>/dev/null
printf '%s\n' '--- available Python executables ---'
command -v py || true
command -v python || true
command -v python3 || true
printf '%s\n' '--- installed PyVista candidates ---'
find /usr /opt /root -type f \( -name 'pyvista.py' -o -path '*/pyvista/__init__.py' \) -print 2>/dev/null | head -20Repository: Project-MONAI/physiotwin4d
Length of output: 326
🌐 Web query:
PyVista empty PolyData extract_cells([]) triangulate compute_normals behavior
💡 Result:
When working with PyVista meshes, passing an empty list (or an empty array) to extract_cells([]) returns an empty mesh of the same type as the input [1][2]. PyVista is designed to handle empty meshes gracefully to maintain pipeline continuity, consistent with VTK's behavior [2]. Following this, the behavior of triangulate and compute_normals on such an empty mesh is as follows: 1. triangulate(): Calling this method on an empty PolyData will return an empty PolyData object. It does not perform operations on non-existent cells [3][4]. 2. compute_normals(): Attempting to compute normals on an empty mesh (or a mesh containing only vertex/line cells) generally does not produce a "Normals" array [5][6]. - Historically, this could lead to a KeyError if the user attempted to access the "Normals" array immediately after the call [5]. - Recent versions of PyVista include protections and have introduced warnings or TypeError exceptions when attempting to compute normals on meshes that contain only vertex or line cells, as surface normals cannot be mathematically defined for such geometry [5][7]. In summary, empty inputs propagate empty outputs, and subsequent filters generally detect the lack of valid surface geometry and will not create the requested normal data arrays [5][6].
Citations:
- 1: https://docs.pyvista.org/api/core/_autosummary/pyvista.DataSetFilters.extract_cells.html
- 2: patched edge case for empty mesh pyvista/pyvista#1030
- 3: https://docs.pyvista.org/api/core/_autosummary/pyvista.PolyDataFilters.triangulate.html
- 4: https://github.com/pyvista/pyvista/blob/main/pyvista/core/filters/poly_data.py
- 5:
compute_normalsdoes not throw error if normals cannot be computed pyvista/pyvista#4709 - 6: https://docs.pyvista.org/api/core/_autosummary/pyvista.PolyDataFilters.compute_normals.html
- 7: Raise
TypeErrorforcompute_normalswith vertex or line cells pyvista/pyvista#5754
🌐 Web query:
PyVista 0.47.0 compute_normals empty PolyData source vtkPolyDataNormals empty mesh
💡 Result:
In PyVista 0.47.0, when calling compute_normals on an empty PolyData source or a mesh that contains only points/lines (no polygons or triangle strips), the compute_normals filter does not produce a Normals array [1][2]. The underlying VTK filter, vtkPolyDataNormals, is designed to calculate normals exclusively for polygons and triangle strips [3][4]. If a mesh lacks these cell types, the filter performs no computation and does not add the Normals data array to the output [3][5]. Because compute_normals does not raise an error in these cases, attempts to subsequently access the mesh["Normals"] array will result in a KeyError [1]. While there have been community discussions and pull requests regarding adding explicit error or warning handling for these scenarios [1][6][7], PyVista maintains its reliance on the standard VTK behavior where the output simply lacks the requested data if it cannot be computed [3][5]. To avoid this, you should check for the presence of polygons or triangle strips in your mesh before invoking compute_normals, or handle the potential KeyError when accessing the normals data [1].
Citations:
- 1:
compute_normalsdoes not throw error if normals cannot be computed pyvista/pyvista#4709 - 2: https://docs.pyvista.org/api/core/_autosummary/pyvista.PolyDataFilters.compute_normals.html
- 3: https://vtk.org/doc/nightly/html/classvtkPolyDataNormals.html
- 4: https://vtk.org/doc/nightly/release/9.1/html/classvtkPolyDataNormals.html
- 5: https://github.com/Kitware/VTK/blob/master/Filters/Core/vtkPolyDataNormals.cxx
- 6: Crash when
compute_normalspyvista/pyvista#2124 - 7: Raise
TypeErrorforcompute_normalswith vertex or line cells pyvista/pyvista#5754
🏁 Script executed:
#!/bin/bash
set -u
printf '%s\n' '--- call sites ---'
rg -n -C4 'extract_label_surfaces' src tests || true
printf '%s\n' '--- consumers of returned surfaces ---'
rg -n -C3 'surfaces\[|for .*surface|surface\.n_cells|compute_normals\(' src/physiotwin4d tests/test_contour_mesh_extraction.py | head -300Repository: Project-MONAI/physiotwin4d
Length of output: 36781
Skip labels with no extracted cells. The empty mesh path returns an empty pv.PolyData without raising, but it has no surface cells or normals. Omit that label from surfaces and emit a warning.
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@src/physiotwin4d/contour_tools.py` around lines 369 - 388, In the
label-processing loop around extract_surface, check whether the extracted and
triangulated surface has zero cells before removing data or computing normals.
For empty surfaces, skip adding the label to surfaces and emit a warning; retain
the existing cleanup, normal computation, and debug logging for labels with
extracted cells.
| @staticmethod | ||
| def smooth_and_decimate_surface( | ||
| def is_watertight(surface: pv.PolyData) -> bool: | ||
| """Report whether every edge of *surface* is shared by exactly two faces.""" | ||
| faces = surface.triangulate().faces.reshape(-1, 4)[:, 1:] | ||
| edges = np.sort( | ||
| np.vstack([faces[:, [0, 1]], faces[:, [1, 2]], faces[:, [2, 0]]]), axis=1 | ||
| ) | ||
| _, counts = np.unique(edges, axis=0, return_counts=True) | ||
| return bool(np.all(counts == 2)) |
There was a problem hiding this comment.
🎯 Functional Correctness | 🟡 Minor | ⚡ Quick win
is_watertight returns True for an empty surface.
With no faces, counts is empty and np.all(counts == 2) evaluates to True. extract_watertight_surface at line 242 then logs no warning for a surface that holds no geometry.
Reject an empty surface explicitly.
🐛 Proposed fix
`@staticmethod`
def is_watertight(surface: pv.PolyData) -> bool:
"""Report whether every edge of *surface* is shared by exactly two faces."""
faces = surface.triangulate().faces.reshape(-1, 4)[:, 1:]
+ if faces.size == 0:
+ return False
edges = np.sort(
np.vstack([faces[:, [0, 1]], faces[:, [1, 2]], faces[:, [2, 0]]]), axis=1
)
_, counts = np.unique(edges, axis=0, return_counts=True)
return bool(np.all(counts == 2))📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| @staticmethod | |
| def smooth_and_decimate_surface( | |
| def is_watertight(surface: pv.PolyData) -> bool: | |
| """Report whether every edge of *surface* is shared by exactly two faces.""" | |
| faces = surface.triangulate().faces.reshape(-1, 4)[:, 1:] | |
| edges = np.sort( | |
| np.vstack([faces[:, [0, 1]], faces[:, [1, 2]], faces[:, [2, 0]]]), axis=1 | |
| ) | |
| _, counts = np.unique(edges, axis=0, return_counts=True) | |
| return bool(np.all(counts == 2)) | |
| `@staticmethod` | |
| def is_watertight(surface: pv.PolyData) -> bool: | |
| """Report whether every edge of *surface* is shared by exactly two faces.""" | |
| faces = surface.triangulate().faces.reshape(-1, 4)[:, 1:] | |
| if faces.size == 0: | |
| return False | |
| edges = np.sort( | |
| np.vstack([faces[:, [0, 1]], faces[:, [1, 2]], faces[:, [2, 0]]]), axis=1 | |
| ) | |
| _, counts = np.unique(edges, axis=0, return_counts=True) | |
| return bool(np.all(counts == 2)) |
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@src/physiotwin4d/contour_tools.py` around lines 491 - 499, Update the
is_watertight method to explicitly return False when the triangulated surface
has no faces, before computing edge counts; preserve the existing
exactly-two-shared-edges check for non-empty surfaces.
| mask_arr = itk.GetArrayViewFromImage(mask_image) != 0 | ||
| # mask_arr axes are reversed relative to the ITK image, so the per-axis | ||
| # extents come back as (z, y, x) and are flipped to (x, y, z). | ||
| starts, stops = [], [] | ||
| for axis_extent in np.nonzero(mask_arr): | ||
| starts.append(int(axis_extent.min())) | ||
| stops.append(int(axis_extent.max()) + 1) | ||
| start_zyx, stop_zyx = np.array(starts), np.array(stops) | ||
| cropped_arr = mask_arr[ | ||
| start_zyx[0] : stop_zyx[0], | ||
| start_zyx[1] : stop_zyx[1], | ||
| start_zyx[2] : stop_zyx[2], | ||
| ] |
There was a problem hiding this comment.
🩺 Stability & Availability | 🟠 Major | ⚡ Quick win
An empty mask crashes extract_tetrahedra.
If mask_image holds no non-zero voxel, np.nonzero(mask_arr) returns empty index arrays and axis_extent.min() raises ValueError: zero-size array to reduction operation minimum. The method already returns an empty pv.UnstructuredGrid for the analogous case at line 569, so the two empty paths behave differently.
Return an empty grid for an empty mask as well.
🐛 Proposed fix
mask_arr = itk.GetArrayViewFromImage(mask_image) != 0
+ if not mask_arr.any():
+ self.log_warning("Mask holds no voxel; its tetrahedral mesh is empty")
+ return pv.UnstructuredGrid()
# mask_arr axes are reversed relative to the ITK image, so the per-axis
# extents come back as (z, y, x) and are flipped to (x, y, z).
starts, stops = [], []📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| mask_arr = itk.GetArrayViewFromImage(mask_image) != 0 | |
| # mask_arr axes are reversed relative to the ITK image, so the per-axis | |
| # extents come back as (z, y, x) and are flipped to (x, y, z). | |
| starts, stops = [], [] | |
| for axis_extent in np.nonzero(mask_arr): | |
| starts.append(int(axis_extent.min())) | |
| stops.append(int(axis_extent.max()) + 1) | |
| start_zyx, stop_zyx = np.array(starts), np.array(stops) | |
| cropped_arr = mask_arr[ | |
| start_zyx[0] : stop_zyx[0], | |
| start_zyx[1] : stop_zyx[1], | |
| start_zyx[2] : stop_zyx[2], | |
| ] | |
| mask_arr = itk.GetArrayViewFromImage(mask_image) != 0 | |
| if not mask_arr.any(): | |
| self.log_warning("Mask holds no voxel; its tetrahedral mesh is empty") | |
| return pv.UnstructuredGrid() | |
| # mask_arr axes are reversed relative to the ITK image, so the per-axis | |
| # extents come back as (z, y, x) and are flipped to (x, y, z). | |
| starts, stops = [], [] | |
| for axis_extent in np.nonzero(mask_arr): | |
| starts.append(int(axis_extent.min())) | |
| stops.append(int(axis_extent.max()) + 1) | |
| start_zyx, stop_zyx = np.array(starts), np.array(stops) | |
| cropped_arr = mask_arr[ | |
| start_zyx[0] : stop_zyx[0], | |
| start_zyx[1] : stop_zyx[1], | |
| start_zyx[2] : stop_zyx[2], | |
| ] |
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@src/physiotwin4d/contour_tools.py` around lines 530 - 542, Update
extract_tetrahedra around the mask_arr bounding-box calculation to detect when
mask_arr contains no non-zero voxels before calling axis_extent.min(). Return
the same empty pv.UnstructuredGrid used by the existing analogous empty path,
while preserving the current cropping and tetrahedra generation for non-empty
masks.
| # Process meshes. Labels win over the static layout: a per-cell label | ||
| # array names the structures outright, which the static layout's one | ||
| # prim per input mesh cannot. | ||
| if self.mask_ids: | ||
| # Split by anatomical regions | ||
| self._convert_with_labels(stage, root_path, material_mgr, mesh_converter) | ||
| elif self._is_static_merge: | ||
| self._convert_static_merge(stage, root_path, material_mgr, mesh_converter) |
There was a problem hiding this comment.
🎯 Functional Correctness | 🟠 Major | ⚡ Quick win
Keep labeled static merges static.
When static_merge=True and one labeled mesh is provided, this dispatch calls
_convert_with_labels(). That method always calls create_time_varying_mesh().
The output then has time samples even though the constructor documents a static
layout and the workflow suppresses stage time metadata for static output.
Call create_mesh() for the single labeled mesh when _is_static_merge is true.
Add an assertion in tests/test_workflow_convert_vtk_to_usd.py that the labeled
static mesh has no point or topology time samples.
Proposed fix
- mesh_converter.create_time_varying_mesh(
- label_mesh_sequence, mesh_path, label_time_codes, bind_material=True
- )
+ if self._is_static_merge:
+ mesh_converter.create_mesh(
+ label_mesh_sequence[0], mesh_path, bind_material=True
+ )
+ else:
+ mesh_converter.create_time_varying_mesh(
+ label_mesh_sequence,
+ mesh_path,
+ label_time_codes,
+ bind_material=True,
+ )📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| # Process meshes. Labels win over the static layout: a per-cell label | |
| # array names the structures outright, which the static layout's one | |
| # prim per input mesh cannot. | |
| if self.mask_ids: | |
| # Split by anatomical regions | |
| self._convert_with_labels(stage, root_path, material_mgr, mesh_converter) | |
| elif self._is_static_merge: | |
| self._convert_static_merge(stage, root_path, material_mgr, mesh_converter) | |
| if self._is_static_merge: | |
| mesh_converter.create_mesh( | |
| label_mesh_sequence[0], mesh_path, bind_material=True | |
| ) | |
| else: | |
| mesh_converter.create_time_varying_mesh( | |
| label_mesh_sequence, | |
| mesh_path, | |
| label_time_codes, | |
| bind_material=True, | |
| ) |
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@src/physiotwin4d/convert_vtk_to_usd.py` around lines 662 - 669, Update the
labeled-mesh dispatch around _convert_with_labels so a static merge with exactly
one labeled mesh uses create_mesh() rather than create_time_varying_mesh(),
while preserving the existing labeled conversion for non-static or multi-mesh
cases. Add a workflow test assertion confirming the labeled static output
contains no point or topology time samples.
| direction = itk.array_from_matrix(image.GetDirection()) | ||
| resampler = itk.ResampleImageFilter[ImageType, ImageType].New() | ||
| resampler.SetInput(image) | ||
| resampler.SetInterpolator(interpolator) | ||
| resampler.SetOutputSpacing([float(v) for v in new_spacing]) | ||
| resampler.SetSize([int(n) for n in new_size]) | ||
| resampler.SetOutputOrigin( | ||
| np.asarray(image.GetOrigin(), dtype=np.float64) | ||
| + direction @ ((new_spacing - spacing) / 2.0) | ||
| ) | ||
| resampler.SetOutputDirection(image.GetDirection()) |
There was a problem hiding this comment.
🎯 Functional Correctness | 🟠 Major | ⚡ Quick win
resample_image_by_scale never sets a default pixel value, and no test covers the border it fills. The origin shift at line 322-325 places the first output voxel center before the first input voxel center whenever scale > 1.0. ITK reports that continuous index as outside the buffer and writes DefaultPixelValue, which is 0. For CT, 0 HU is water rather than air, the exact problem TransformTools.transform_image already documents.
src/physiotwin4d/image_tools.py#L316-L326: add abackground_valueparameter and callresampler.SetDefaultPixelValue(background_value)beforeUpdate(), following thetransform_imageconvention.tests/test_image_tools.py#L521-L538: add ascale > 1.0test that resamples a constant-valued image and asserts every output voxel keeps that value, and assert that the nearest-neighbor test's block value7.0is still present.
📍 Affects 2 files
src/physiotwin4d/image_tools.py#L316-L326(this comment)tests/test_image_tools.py#L521-L538
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@src/physiotwin4d/image_tools.py` around lines 316 - 326, The
resample_image_by_scale path must support configurable border fill values: add a
background_value parameter following the transform_image convention and apply it
with SetDefaultPixelValue before Update(). In src/physiotwin4d/image_tools.py
lines 316-326, make this change; in tests/test_image_tools.py lines 521-538, add
a scale greater than 1 test using a constant-valued image that asserts every
output voxel retains that value, and ensure the nearest-neighbor test still
asserts the 7.0 block value.
| case_dirs = [input_dir / "pm0004"] | ||
| # sorted( | ||
| # path for path in input_dir.glob("pm[0-9][0-9][0-9][0-9]") if path.is_dir() | ||
| # ) | ||
| if not case_dirs: | ||
| raise FileNotFoundError( | ||
| f"No pm???? case directories found under {input_dir}.\n" | ||
| "Run tutorial_04_duke_heart_labelmap_to_vtk.py first." | ||
| ) |
There was a problem hiding this comment.
📐 Maintainability & Code Quality | 🟠 Major | ⚡ Quick win
Restore the cohort glob; the hard-coded case is a leftover debug artifact.
Line 69 pins the run to pm0004 and comments out the case discovery. The docstring (Line 29) documents pm????/*_surfaces.vtp, so the script no longer matches its own contract. The if not case_dirs guard at Line 73 is also dead, because the list literal is never empty. If pm0004 is absent, the run only logs a warning and writes nothing.
♻️ Proposed fix
- case_dirs = [input_dir / "pm0004"]
- # sorted(
- # path for path in input_dir.glob("pm[0-9][0-9][0-9][0-9]") if path.is_dir()
- # )
+ case_dirs = sorted(
+ path for path in input_dir.glob("pm[0-9][0-9][0-9][0-9]") if path.is_dir()
+ )
if not case_dirs:📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| case_dirs = [input_dir / "pm0004"] | |
| # sorted( | |
| # path for path in input_dir.glob("pm[0-9][0-9][0-9][0-9]") if path.is_dir() | |
| # ) | |
| if not case_dirs: | |
| raise FileNotFoundError( | |
| f"No pm???? case directories found under {input_dir}.\n" | |
| "Run tutorial_04_duke_heart_labelmap_to_vtk.py first." | |
| ) | |
| case_dirs = sorted( | |
| path for path in input_dir.glob("pm[0-9][0-9][0-9][0-9]") if path.is_dir() | |
| ) | |
| if not case_dirs: | |
| raise FileNotFoundError( | |
| f"No pm???? case directories found under {input_dir}.\n" | |
| "Run tutorial_04_duke_heart_labelmap_to_vtk.py first." | |
| ) |
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@tutorials/tutorial_05_duke_heart_vtk_to_usd.py` around lines 69 - 77, Restore
dynamic case discovery by assigning case_dirs to the sorted directories matching
input_dir.glob("pm[0-9][0-9][0-9][0-9]"), filtering to directories, and remove
the hard-coded pm0004 assignment and commented-out code. Keep the existing
FileNotFoundError guard so missing cohorts fail as documented.
| sample_files = [ | ||
| path | ||
| for path in sorted(input_dir.glob("*_ref_heart_minus_interior_chambers.vtp")) | ||
| if not path.name.startswith(DUKE_HEART.hold_out_case) | ||
| ] |
There was a problem hiding this comment.
🗄️ Data Integrity & Integration | 🟠 Major | ⚡ Quick win
Read Duke Tutorial 4 outputs from their case directories.
Both consumers assume that Tutorial 4 writes surfaces directly under its output root. The documented output layout stores each surface under a pm???? case directory. Tutorial 6 finds no samples, and Tutorial 8 never reuses Tutorial 4 surfaces.
tutorials/tutorial_06_duke_heart_create_statistical_model.py#L88-L92: searchpm????/*_ref_heart_minus_interior_chambers.vtp.tutorials/tutorial_08_duke_heart_fit_model_to_4d_patients.py#L189-L192: includelabelmap_file.parent.namebetweentutorial_04_dirand the surface filename.
📍 Affects 2 files
tutorials/tutorial_06_duke_heart_create_statistical_model.py#L88-L92(this comment)tutorials/tutorial_08_duke_heart_fit_model_to_4d_patients.py#L189-L192
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@tutorials/tutorial_06_duke_heart_create_statistical_model.py` around lines 88
- 92, The Tutorial 4 surface consumers use the wrong output layout. In
tutorials/tutorial_06_duke_heart_create_statistical_model.py lines 88-92, update
the sample_files glob to search *_ref_heart_minus_interior_chambers.vtp within
pm???? case directories while preserving the hold-out filter; in
tutorials/tutorial_08_duke_heart_fit_model_to_4d_patients.py lines 189-192,
include labelmap_file.parent.name between tutorial_04_dir and the surface
filename when constructing the surface path.
| mean_workflow = WorkflowCreateMeanSurface( | ||
| surfaces=sample_surfaces, | ||
| template_surface=sample_surfaces[4], | ||
| log_level=log_level, |
There was a problem hiding this comment.
🎯 Functional Correctness | 🟠 Major | ⚡ Quick win
Select a template surface that exists.
Test mode limits sample_files to three entries. Line 119 then selects the fifth surface and raises IndexError. Select a valid sample, such as the middle surface, after the existing minimum-count check.
Proposed fix
- template_surface=sample_surfaces[4],
+ template_surface=sample_surfaces[len(sample_surfaces) // 2],📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| mean_workflow = WorkflowCreateMeanSurface( | |
| surfaces=sample_surfaces, | |
| template_surface=sample_surfaces[4], | |
| log_level=log_level, | |
| mean_workflow = WorkflowCreateMeanSurface( | |
| surfaces=sample_surfaces, | |
| template_surface=sample_surfaces[len(sample_surfaces) // 2], | |
| log_level=log_level, |
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@tutorials/tutorial_06_duke_heart_create_statistical_model.py` around lines
117 - 120, Update the template_surface argument in WorkflowCreateMeanSurface to
use a valid entry from sample_surfaces, such as the middle element, while
preserving the existing minimum-count check and the remaining workflow
configuration.
| pca_model: Optional[dict[str, Any]] = None | ||
| if pca_json.exists(): | ||
| with pca_json.open(encoding="utf-8") as f: | ||
| pca_model = json.load(f) | ||
|
|
There was a problem hiding this comment.
🎯 Functional Correctness | 🟠 Major | ⚡ Quick win
Require the PCA model JSON before fitting.
If pca_mean_file exists but pca_json is absent, this code disables PCA registration and continues. The tutorial then performs a non-PCA fit while its inputs and outputs claim PCA behavior. Reject a missing pca_json with the same prerequisite error used for the mean surface.
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@tutorials/tutorial_07_duke_heart_fit_statistical_model_to_patient.py` around
lines 114 - 118, Update the PCA model loading logic around pca_model and
pca_json so that pca_json is required whenever pca_mean_file exists; when it is
absent, raise the same prerequisite error used for the missing mean surface
instead of continuing with pca_model unset. Preserve JSON loading when the file
exists and ensure fitting cannot proceed in a non-PCA mode under these
conditions.
| # Step 1: pick the test frame - the one 70% of the way through the case's | ||
| # ordered gated frames - and read its stage and ground-truth surface. | ||
| stages = [_cardiac_stage_from_filename(f) for f in phase_files] | ||
| test_index = int(stage_fraction * len(stages)) | ||
| test_stage = stages[test_index] | ||
| ground_truth_file = phase_files[test_index] |
There was a problem hiding this comment.
🩺 Stability & Availability | 🟡 Minor | ⚡ Quick win
Clamp the frame index.
Line 133 computes test_index = int(stage_fraction * len(stages)). The docs invite users to change stage_fraction; a value of 1.0 produces an index equal to len(stages) and raises IndexError at Line 134. Clamp the index to the last frame.
🛡️ Proposed fix
- test_index = int(stage_fraction * len(stages))
+ test_index = min(int(stage_fraction * len(stages)), len(stages) - 1)📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| # Step 1: pick the test frame - the one 70% of the way through the case's | |
| # ordered gated frames - and read its stage and ground-truth surface. | |
| stages = [_cardiac_stage_from_filename(f) for f in phase_files] | |
| test_index = int(stage_fraction * len(stages)) | |
| test_stage = stages[test_index] | |
| ground_truth_file = phase_files[test_index] | |
| # Step 1: pick the test frame - the one 70% of the way through the case's | |
| # ordered gated frames - and read its stage and ground-truth surface. | |
| stages = [_cardiac_stage_from_filename(f) for f in phase_files] | |
| test_index = min(int(stage_fraction * len(stages)), len(stages) - 1) | |
| test_stage = stages[test_index] | |
| ground_truth_file = phase_files[test_index] |
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@tutorials/tutorial_10_duke_heart_infer_physicsnemo.py` around lines 130 -
135, Clamp test_index in the frame-selection logic before indexing stages and
phase_files, ensuring values such as stage_fraction=1.0 resolve to the final
frame rather than len(stages). Preserve the existing fractional selection
behavior for in-range values and use the last valid index as the upper bound.
1 participant
Introduce tutorials/parameters_lung_ct_dirlab.py and tutorials/parameters_heart_ct_kcl.py as the single source for each use case's mask dilation, distance-map saturation radius, PCA component counts, Greedy iteration schedule, segmenter class, and (heart only) the interior chamber label ids. Every tutorial that rasterizes or registers a distance map now reads the same values, so the maps a network is finetuned on match the maps it later infers over. No paths live in these modules; each tutorial keeps its own inputs and outputs.
Add tutorials/tutorial_02_heart_distancemap_finetune_icon.py, which finetunes uniGradICON on heart distance maps built from the Duke-Heart-4DLabelmaps labelmaps with the chambers excluded. The heart needs its own run rather than reusing the lung weights: its registration mask is much tighter, so its distance maps saturate over a shorter radius and do not share an intensity distribution with the lung ones.
Library fixes:
Rename number_of_components / number_of_modes to
number_of_pca_components throughout the workflows and tutorials.
tutorial_02_lung_finetune_icon now writes difference images (fixed minus registered) instead of the resampled volumes, and reports the chain's Greedy-stage-only score as its own row. The chain remains unconditional: on DIR-Lab, ICON's 175^3 residual grid is about 1.4 mm over the FOV, coarser than the 1.10 mm Greedy already achieves, so it cannot refine and the tutorial reports that honestly.
tutorial_02_lung_distancemap_finetune_icon restricts its cached labelmaps to the lung labels. They previously held all 97 whole-body classes, and uniGradICON's Dice loss one-hots every shared class at 175^3 by batch 4, which saturated GPU memory.
Baselines for the slow and GPU buckets will need refreshing: the composition-order fix and the PCA field grid change alter registration output.
Summary by CodeRabbit
New Features
Documentation
Updates