From f787d7eaba560cbc1244d5774e5840367eacc80e Mon Sep 17 00:00:00 2001 From: Stephen Aylward Date: Thu, 13 Aug 2026 17:58:45 -0400 Subject: [PATCH 1/4] ENH: Movement evaluation workflow and end-to-end tutorials 11-13 Add WorkflowEvaluateMovement, which scores an inferred moving anatomy per anatomical structure (Dice, volume, surface RMSE) and writes a CSV, volume plots, and a markdown report with provenance. Add WorkflowInferMovement.process_time_series to predict one subject across a whole time series, plus TransformTools helpers to smooth scalar arrays and wrap them as ITK images so per-timepoint deformation magnitude and RMSE fields can be written out. Tutorials: - 11 (heart, lung): evaluate PhysicsNeMo inference against ground truth - 12 (heart, lung): end-to-end inference from a raw 4D image - 13: combined heart and lung motion - rename tutorial_10_duke_heart_infer_physicsnemo to *_mgn for consistency with the lung tutorial Fix ContourTools.split_labeled_surface to pass an integer array to extract_cells so a label with no cells is reported instead of raising. Tests: test_workflow_evaluate_movement, test_workflow_train_physicsnemo, and added transform_tools coverage. Repository hygiene: - stop tracking tutorials/network_weights in LFS; weights are produced by running Tutorials 2 and 9, not distributed - MANIFEST.in drops setup.py, scripts/, and bundled weights; ships AGENTS.md - CI test paths and comments follow the renamed test modules - nightly status.json is published through the nightly-status branch Docs: tutorials, architecture, quickstart, installation, testing, and API pages updated for the new workflow; new assets for tutorials 8-13; drop stale brain_vessel_modeling and lung_gated_ct CLI pages. --- .gitattributes | 8 - .github/scripts/build_dashboard.py | 17 +- .github/workflows/README.md | 22 +- .github/workflows/ci.yml | 10 +- .gitignore | 14 +- MANIFEST.in | 21 +- data/Duke-Heart-4DLabelmaps/README.md | 15 +- docs/api/index.rst | 7 + docs/api/physicsnemo/evaluate.rst | 73 ++ docs/api/physicsnemo/index.rst | 8 +- docs/api/physicsnemo/manifest.rst | 7 + docs/api/utilities/index.rst | 5 + docs/api/workflows.rst | 6 +- docs/architecture.rst | 27 +- docs/assets/example.gif | 3 - ...nal.gif => tutorial_03_heart_original.gif} | 0 ..._recon.gif => tutorial_03_heart_recon.gif} | 0 ...heart-2png.png => tutorial_04_heart-2.png} | 0 .../assets/tutorial_08_duke_heart_def_mag.gif | 3 + docs/assets/tutorial_08_lung.gif | 4 +- ...al_09_duke_heart_deformation_magnitude.gif | 3 + docs/assets/tutorial_09_duke_heart_motion.gif | 3 + docs/assets/tutorial_09_duke_heart_rmse.gif | 3 + ...tutorial_09_lung_deformation_magnitude.gif | 3 + docs/assets/tutorial_09_lung_motion.gif | 3 + docs/assets/tutorial_09_lung_rmse.gif | 3 + .../tutorial_10_duke_heart_motion_usd.gif | 3 + docs/assets/tutorial_10_lung_motion_usd.gif | 3 + docs/assets/tutorial_11_duke_heart_stats.png | 3 + .../assets/tutorial_11_duke_heart_volumes.png | 3 + docs/assets/tutorial_11_lung_stats.png | 3 + docs/assets/tutorial_11_lung_volumes.png | 3 + docs/assets/tutorial_12_duke_heart.gif | 3 + docs/assets/tutorial_12_lung.gif | 3 + docs/assets/tutorial_13_combined_motion.gif | 3 + docs/cli_scripts/brain_vessel_modeling.rst | 16 - docs/cli_scripts/download_data.rst | 11 +- docs/cli_scripts/lung_gated_ct.rst | 15 - docs/developer/workflows.rst | 14 + docs/faq.rst | 12 +- docs/index.rst | 15 + docs/installation.rst | 22 +- docs/quickstart.rst | 24 +- docs/testing.rst | 12 + docs/tutorials.rst | 426 ++++++++- docs/viewing_usd.rst | 8 +- pyproject.toml | 15 +- src/physiotwin4d/__init__.py | 2 + src/physiotwin4d/cli/__init__.py | 4 + src/physiotwin4d/contour_tools.py | 6 +- src/physiotwin4d/transform_tools.py | 145 ++- .../workflow_evaluate_movement.py | 550 ++++++++++++ src/physiotwin4d/workflow_infer_movement.py | 231 ++++- .../workflow_infer_physicsnemo.py | 2 + statistics.md | 80 +- tests/test_transform_tools.py | 113 +++ tests/test_workflow_evaluate_movement.py | 216 +++++ tests/test_workflow_train_physicsnemo.py | 160 ++++ tutorials/README.md | 44 +- ...ial_09_duke_heart_train_physicsnemo_mgn.py | 2 +- ...al_10_duke_heart_infer_physicsnemo_mgn.py} | 127 +-- .../tutorial_10_lung_infer_physicsnemo_mgn.py | 118 +-- ...rial_11_duke_heart_evaluate_physicsnemo.py | 220 +++++ .../tutorial_11_lung_evaluate_physicsnemo.py | 230 +++++ ...rial_12_duke_heart_end_to_end_inference.py | 351 ++++++++ .../tutorial_12_lung_end_to_end_inference.py | 348 +++++++ .../tutorial_13_heart_and_lung_motion.py | 849 ++++++++++++++++++ 67 files changed, 4345 insertions(+), 338 deletions(-) create mode 100644 docs/api/physicsnemo/evaluate.rst delete mode 100644 docs/assets/example.gif rename docs/assets/{Tutorial_03_heart_original.gif => tutorial_03_heart_original.gif} (100%) rename docs/assets/{Tutorial_03_heart_recon.gif => tutorial_03_heart_recon.gif} (100%) rename docs/assets/{tutorial_04_heart-2png.png => tutorial_04_heart-2.png} (100%) create mode 100644 docs/assets/tutorial_08_duke_heart_def_mag.gif create mode 100644 docs/assets/tutorial_09_duke_heart_deformation_magnitude.gif create mode 100644 docs/assets/tutorial_09_duke_heart_motion.gif create mode 100644 docs/assets/tutorial_09_duke_heart_rmse.gif create mode 100644 docs/assets/tutorial_09_lung_deformation_magnitude.gif create mode 100644 docs/assets/tutorial_09_lung_motion.gif create mode 100644 docs/assets/tutorial_09_lung_rmse.gif create mode 100644 docs/assets/tutorial_10_duke_heart_motion_usd.gif create mode 100644 docs/assets/tutorial_10_lung_motion_usd.gif create mode 100644 docs/assets/tutorial_11_duke_heart_stats.png create mode 100644 docs/assets/tutorial_11_duke_heart_volumes.png create mode 100644 docs/assets/tutorial_11_lung_stats.png create mode 100644 docs/assets/tutorial_11_lung_volumes.png create mode 100644 docs/assets/tutorial_12_duke_heart.gif create mode 100644 docs/assets/tutorial_12_lung.gif create mode 100644 docs/assets/tutorial_13_combined_motion.gif delete mode 100644 docs/cli_scripts/brain_vessel_modeling.rst delete mode 100644 docs/cli_scripts/lung_gated_ct.rst create mode 100644 src/physiotwin4d/workflow_evaluate_movement.py create mode 100644 tests/test_workflow_evaluate_movement.py create mode 100644 tests/test_workflow_train_physicsnemo.py rename tutorials/{tutorial_10_duke_heart_infer_physicsnemo.py => tutorial_10_duke_heart_infer_physicsnemo_mgn.py} (59%) create mode 100644 tutorials/tutorial_11_duke_heart_evaluate_physicsnemo.py create mode 100644 tutorials/tutorial_11_lung_evaluate_physicsnemo.py create mode 100644 tutorials/tutorial_12_duke_heart_end_to_end_inference.py create mode 100644 tutorials/tutorial_12_lung_end_to_end_inference.py create mode 100644 tutorials/tutorial_13_heart_and_lung_motion.py diff --git a/.gitattributes b/.gitattributes index fd851d0f..a75aeae1 100644 --- a/.gitattributes +++ b/.gitattributes @@ -2,11 +2,3 @@ *.mha filter=lfs diff=lfs merge=lfs -text docs/assets/*.gif filter=lfs diff=lfs merge=lfs -text docs/assets/*.png filter=lfs diff=lfs merge=lfs -text -tutorials/network_weights/physicsnemo_mgn_heart/*.pt filter=lfs diff=lfs merge=lfs -text -tutorials/network_weights/physicsnemo_mgn_heart/*.vtu filter=lfs diff=lfs merge=lfs -text -tutorials/network_weights/physicsnemo_mgn_heart/*.vtp filter=lfs diff=lfs merge=lfs -text -tutorials/network_weights/physicsnemo_mgn_heart/pca_model.json filter=lfs diff=lfs merge=lfs -text -tutorials/network_weights/physicsnemo_mgn_lung_motion/*.pt filter=lfs diff=lfs merge=lfs -text -tutorials/network_weights/physicsnemo_mgn_lung_motion/*.vtu filter=lfs diff=lfs merge=lfs -text -tutorials/network_weights/physicsnemo_mgn_lung_motion/*.vtp filter=lfs diff=lfs merge=lfs -text -tutorials/network_weights/physicsnemo_mgn_lung_motion/pca_model.json filter=lfs diff=lfs merge=lfs -text diff --git a/.github/scripts/build_dashboard.py b/.github/scripts/build_dashboard.py index 3ec64782..1bef8133 100644 --- a/.github/scripts/build_dashboard.py +++ b/.github/scripts/build_dashboard.py @@ -14,18 +14,17 @@ --timestamp "2026-03-31T07:05:42Z" \\ --health-outcome "success" -Artifact publishing: - ``status.json`` is uploaded by ``nightly-health.yml`` as a standalone - artifact named ``nightly-status-json`` (90-day retention). ``docs.yml`` - downloads that artifact during its ``deploy`` job and copies - ``status.json`` into the Pages output directory so that the file is - served at the live URL: +Status publishing: + ``nightly-health.yml`` uploads the whole output directory as a + ``health-dashboard`` artifact, then force-pushes ``status.json`` alone to + the orphan ``nightly-status`` branch. ``docs.yml`` fetches it from that + branch through the GitHub contents API during its ``deploy`` job and copies + it into the Pages output, so the file is served at the live URL: https:///status.json - The copy step uses ``continue-on-error: true`` so the first docs deploy - (before any nightly run has produced the artifact) succeeds without - ``status.json`` being present. + The fetch tolerates a missing branch, so the first docs deploy (before any + nightly run has pushed ``status.json``) succeeds without it. """ from __future__ import annotations diff --git a/.github/workflows/README.md b/.github/workflows/README.md index ed81f499..e73f7af2 100644 --- a/.github/workflows/README.md +++ b/.github/workflows/README.md @@ -11,7 +11,7 @@ Runs on every push and pull request to main branches. Includes: - **unit-tests**: Cross-platform unit tests - Runs on Ubuntu and Windows - Python 3.11 and 3.12 - - Uses PyTorch CPU version to avoid GPU dependencies + - Installs `.[test]` only, so no CUDA toolchain is pulled in - Excludes slow tests and tests requiring external data - Generates coverage reports @@ -35,7 +35,8 @@ Runs on every push and pull request to main branches. Includes: ### `test-slow.yml` - Long-Running Tests -Runs nightly at 2 AM UTC or on manual trigger. Includes: +Manual trigger only (`workflow_dispatch`); the nightly schedule lives in +`nightly-health.yml`. Includes: - **test-slow-gpu**: Slow tests requiring GPU - Tests marked with `slow` marker @@ -61,6 +62,23 @@ Two-job workflow for building and deploying Sphinx documentation: This separation ensures PRs can build and validate docs without triggering environment protection rules. +The `deploy` job also fetches `status.json` from the orphan `nightly-status` +branch (via the GitHub contents API) and copies it into the Pages output, so the +nightly-health badge in the top-level `README.md` resolves. + +### `nightly-health.yml` - Nightly Full-Suite Health Check + +Runs at 07:00 UTC daily, or on manual trigger with a `reason` input. On the +self-hosted Windows GPU runner it installs +`.[test,docs,cuda13,dev,physicsnemo]` and runs the entire suite with +`--run-all`, which enables every opt-in bucket. The run itself is +`continue-on-error`, so a failing test records a red status rather than +failing the workflow. + +A second `build-dashboard` job turns the JUnit XML and coverage JSON into an +HTML dashboard via `.github/scripts/build_dashboard.py`, then force-pushes +`status.json` to the orphan `nightly-status` branch for `docs.yml` to pick up. + ### `release.yml` - Build and Publish Distributions Builds the wheel and source distribution, validates them with Twine, and diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ba25d1e5..d867dcd3 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -235,7 +235,7 @@ jobs: - name: Run USD conversion tests run: | xvfb-run -a --server-args="-screen 0 1024x768x24" \ - pytest tests/test_convert_vtk_to_usd_polymesh.py -v --cov=physiotwin4d --cov-append --cov-report=xml + pytest tests/test_convert_vtk_to_usd.py -v --cov=physiotwin4d --cov-append --cov-report=xml continue-on-error: true - name: Run USD utility tests @@ -408,13 +408,13 @@ jobs: # The following tests are excluded from CI and should be run locally: # # Slow/GPU-intensive tests: -# - tests/test_register_images_ANTS.py (slow, computationally intensive) -# - tests/test_register_images_ICON.py (requires CUDA for ICON) +# - tests/test_register_images_ants.py (slow, computationally intensive) +# - tests/test_register_images_icon.py (requires CUDA for ICON) # - tests/test_transform_tools.py (depends on slow registration tests) # - tests/test_segment_chest_total_segmentator.py (requires CUDA for TotalSegmentator) # # Tutorial tests (SLOW - hours to complete): -# - tests/test_tutorials.py (runs every script in tutorials/ end-to-end) +# - tests/test_tutorials.py (runs 9 of the 29 tutorial scripts end-to-end) # These tests are NEVER run in the PR CI and must be opted into # They execute end-to-end workflows that may take multiple hours # @@ -422,7 +422,7 @@ jobs: # pytest tests/ -v --run-slow # Run all slow tests # pytest tests/ -v --run-gpu --run-slow # GPU + slow (typical local dev profile) # pytest tests/ -v --run-simpleware --run-gpu --run-slow # Full Simpleware coverage -# pytest tests/test_register_images_ANTS.py -v --run-slow +# pytest tests/test_register_images_ants.py -v --run-slow # # Self-hosted GPU runner enables ALL buckets via --run-all # (--run-gpu --run-slow --run-simpleware --run-physicsnemo --run-tutorials). diff --git a/.gitignore b/.gitignore index 6627d71f..b9cf66f6 100644 --- a/.gitignore +++ b/.gitignore @@ -28,17 +28,15 @@ pr*review_summary.md docs/_build/ docs/_static/.buildinfo -# Tutorial videos (GIF versions are tracked in LFS instead) +# Tutorial videos (GIF versions are tracked in LFS instead). The MP4 sources +# are kept beside the GIFs in docs/assets/mp4/ and are never tracked. docs/assets/*.mp4 +docs/assets/mp4/ -# Network weights +# Network weights. Weights are produced by running Tutorials 2 and 9, not +# distributed with the repository, so nothing under network_weights is tracked. network_weights - -# Track the MGN heart and lung model weights (other network_weights stay ignored) -!tutorials/network_weights/ -tutorials/network_weights/* -!tutorials/network_weights/physicsnemo_mgn_heart/ -!tutorials/network_weights/physicsnemo_mgn_lung_motion/ +tutorials/network_weights/ # Data files *.gz diff --git a/MANIFEST.in b/MANIFEST.in index 0498a90e..11fd8ead 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,24 +1,21 @@ -# Include essential documentation +# Include essential documentation. AGENTS.md is the cross-tool sibling of +# CLAUDE.md; both describe repository conventions and ship together. include README.md include LICENSE include CLAUDE.md -include CHANGELOG.md +include AGENTS.md -# Include configuration files +# Include configuration files. The build is pure PEP 621, so there is no +# setup.py to ship. include pyproject.toml -include setup.py - -# Include example scripts (but not in the package itself) -recursive-include scripts *.py -recursive-include scripts *.md # Include all source code recursive-include src *.py -# Include network weights and models -recursive-include src/physiotwin4d/network_weights * - -# Include data files +# Include data files. src/ currently holds no non-Python package data; these +# rules are kept so that any future bundled config ships automatically. Network +# weights are deliberately excluded - they are produced by running Tutorials 2 +# and 9, not distributed with the package. recursive-include src/physiotwin4d *.json recursive-include src/physiotwin4d *.yaml recursive-include src/physiotwin4d *.yml diff --git a/data/Duke-Heart-4DLabelmaps/README.md b/data/Duke-Heart-4DLabelmaps/README.md index 1af9a9db..fadbb6e0 100644 --- a/data/Duke-Heart-4DLabelmaps/README.md +++ b/data/Duke-Heart-4DLabelmaps/README.md @@ -4,9 +4,11 @@ Gated 4D cardiac labelmaps acquired at Duke University by Dr. Paul Segars. ## Availability -This dataset is **not currently available**. It is being considered for public -release; until that happens it cannot be downloaded, and it is not distributed -with this repository. +This dataset is **scheduled for public release soon**. It is not distributed +with this repository and has no automatic downloader yet. + +In the meantime, contact Stephen Aylward () to request +access. ## Effect on the tutorials @@ -15,9 +17,10 @@ their organ field, for example: - `tutorials/tutorial_02_duke_heart_distancemap_finetune_icon.py` -These `duke_heart` tutorials will not run without the data. Every other -tutorial uses a publicly available dataset and is unaffected — see -[../README.md](../README.md) for download instructions. +There are ten of them, forming their own chain: Tutorial 4 (duke heart) -> 5 -> +6 -> 7 -> 8 -> 9 -> 10 -> 11 -> 12. They will not run until the data is +available. The other 19 tutorial scripts use publicly available datasets and are +unaffected — see [../README.md](../README.md) for download instructions. Downstream tutorials that consume `duke_heart` outputs (such as the finetuned distance-map ICON weights used by diff --git a/docs/api/index.rst b/docs/api/index.rst index 1bcf1a16..70d37c39 100644 --- a/docs/api/index.rst +++ b/docs/api/index.rst @@ -43,16 +43,23 @@ By Category * :class:`~physiotwin4d.WorkflowTrainPhysicsNeMo` - Train a mesh-stage model * :class:`~physiotwin4d.WorkflowInferPhysicsNeMo` - Predict per-point targets * :class:`~physiotwin4d.WorkflowInferMovement` - Turn predictions back into geometry + * :class:`~physiotwin4d.WorkflowEvaluateMovement` - Score predictions against the acquired frames **Segmentation** * :class:`~physiotwin4d.SegmentAnatomyBase` - Base segmentation class * :class:`~physiotwin4d.SegmentChestTotalSegmentator` - TotalSegmentator + * :class:`~physiotwin4d.SegmentChestTotalSegmentatorWithContrast` - TotalSegmentator for contrast-enhanced CT * :class:`~physiotwin4d.SegmentHeartSimpleware` - Simpleware cardiac segmentation + * :class:`~physiotwin4d.SegmentHeartSimplewareTrimmedBranches` - Simpleware with trimmed great vessels + * :class:`~physiotwin4d.SegmentNVSegmentCTMRI` - NV-Segment-CTMR, CT *and* MRI **Image Registration** * :class:`~physiotwin4d.RegisterImagesBase` - Base registration class * :class:`~physiotwin4d.RegisterImagesANTS` - ANTs registration + * :class:`~physiotwin4d.RegisterImagesGreedy` - Greedy classical deformable registration * :class:`~physiotwin4d.RegisterImagesICON` - Icon deep learning registration + * :class:`~physiotwin4d.RegisterImagesChain` - Run registrations back to back + * :class:`~physiotwin4d.RegisterImagesGreedyICON` - Greedy then ICON, as a preset chain * :class:`~physiotwin4d.RegisterTimeSeriesImages` - 4D time series registration **Model Registration** diff --git a/docs/api/physicsnemo/evaluate.rst b/docs/api/physicsnemo/evaluate.rst new file mode 100644 index 00000000..7c27a3ea --- /dev/null +++ b/docs/api/physicsnemo/evaluate.rst @@ -0,0 +1,73 @@ +========================================== +Scoring a Mesh-Stage Model Against Images +========================================== + +.. module:: physiotwin4d.workflow_evaluate_movement +.. currentmodule:: physiotwin4d + +A mm error against the surfaces a registration produced says how well the +network reproduces that registration. :class:`WorkflowEvaluateMovement` asks the +other question: how close are the size and shape of the inferred anatomy to the +anatomy that was actually imaged, structure by structure. See Tutorial 11 in +:doc:`../../tutorials`. + +For every gated time point it carries the reference frame's labelmap into that +time point with the network's own deformation, and compares the result to the +labelmap of the frame that was acquired: volume difference, Dice and surface +RMSE per lung lobe or per heart chamber. + +Per-structure scoring +===================== + +.. autoclass:: WorkflowEvaluateMovement + :members: + :undoc-members: + :show-inheritance: + +Example +======= + +.. code-block:: python + + from physiotwin4d import ( + WorkflowEvaluateMovement, + WorkflowInferMovement, + WorkflowInferPhysicsNeMo, + ) + + evaluate = WorkflowEvaluateMovement( + movement_workflow=WorkflowInferMovement( + WorkflowInferPhysicsNeMo(model_directory=model_dir) + ), + label_names={28: "lung_upper_lobe_left", 29: "lung_lower_lobe_left"}, + ) + result = evaluate.process( + case_id="Case1Pack", + shape_parameters=pca_coefficients_file, + reference_mesh=ssm_surface_file, + reference_labelmap=reference_labelmap, + ground_truth_labelmaps={0.0: frame_00, 0.1: frame_10}, + output_directory=out_dir, + ) + print(result["report_file"], result["csv_file"]) + +Notes +===== + +**Why labelmaps rather than the model's surface.** The lung shape model carries +its five lobes as per-cell labels, but the heart model is a single structure --- +the whole heart minus its chamber cavities --- so its chambers exist only in the +acquired labelmaps. Warping those labelmaps scores every structure the +acquisition contains, whether or not the shape model represents it separately. + +**The evaluation grid.** Everything is measured on one isotropic grid built +around the reference anatomy, so a case whose gated frames carry different slice +pitches is still scored on a single, stated voxel volume. Its pitch sets both +that voxel volume and the memory the per-stage deformation fields take, which +grows with its cube. + +See Also +======== + +* :doc:`infer` +* :doc:`train` diff --git a/docs/api/physicsnemo/index.rst b/docs/api/physicsnemo/index.rst index aa7152d5..ce07948e 100644 --- a/docs/api/physicsnemo/index.rst +++ b/docs/api/physicsnemo/index.rst @@ -6,7 +6,7 @@ PhysioTwin4D trains and runs PhysicsNeMo mesh-stage models: given a subject's shape parameters and a stage (a point in the cardiac or respiratory cycle), predict a per-vertex target on the shared template mesh. When that target is a displacement, the prediction replaces a per-phase registration solve with one -forward pass — see Tutorials 9 and 10 in :doc:`../../tutorials`. +forward pass — see Tutorials 9 through 13 in :doc:`../../tutorials`. The layer follows the same has-a shape as the rest of the workflow tier: a workflow owns the data and the artifacts, and a *method* object owns the @@ -24,7 +24,10 @@ network. - Loads a trained model and predicts raw per-point targets * - :class:`~physiotwin4d.WorkflowInferMovement` - Interprets 3-component targets as displacements: deformed meshes, mm - error statistics, rasterized deformation fields + error statistics, rasterized deformation fields, warped images and USD + * - :class:`~physiotwin4d.WorkflowEvaluateMovement` + - Scores those predictions per structure against the acquired frames: + volume difference, Dice and surface RMSE * - :class:`~physiotwin4d.TrainPhysicsNeMoMGN` / :class:`~physiotwin4d.TrainPhysicsNeMoMLP` - The networks to train: MeshGraphNet or fully connected @@ -46,3 +49,4 @@ imports happen lazily inside the methods that need them. manifest train infer + evaluate diff --git a/docs/api/physicsnemo/manifest.rst b/docs/api/physicsnemo/manifest.rst index dc02bd44..a97b9489 100644 --- a/docs/api/physicsnemo/manifest.rst +++ b/docs/api/physicsnemo/manifest.rst @@ -45,6 +45,13 @@ decides which domain the model lives on. Reference ========= +These live in :mod:`physiotwin4d.physicsnemo_tools`, which is not re-exported +from the top-level package — import it by module: + +.. code-block:: python + + from physiotwin4d.physicsnemo_tools import SubjectManifest, parse_manifest + .. autoclass:: SubjectManifest :exclude-members: subject_id, reference_mesh, pca_coefficients, target_array, phases diff --git a/docs/api/utilities/index.rst b/docs/api/utilities/index.rst index fae50ec2..eef82ce5 100644 --- a/docs/api/utilities/index.rst +++ b/docs/api/utilities/index.rst @@ -12,9 +12,13 @@ Overview Utility modules provide low-level operations: * **Image Tools**: Image I/O, preprocessing, manipulation +* **Labelmap Tools**: Labelmap to registration-mask conversion * **Transform Tools**: Coordinate transforms and warping +* **Landmark Tools**: Landmark-based registration validation metrics * **Contour Tools**: Contour extraction and processing * **4D Image Conversion**: 4D image to 3D time-series conversion utilities +* **Test Tools**: Baseline and result comparison helpers +* **Data Download Tools**: Optional dataset download helpers Quick Links =========== @@ -23,6 +27,7 @@ Quick Links * :doc:`image_tools` - Image processing utilities * :doc:`labelmap_tools` - Labelmap to registration-mask conversion * :doc:`transform_tools` - Transform operations + * :doc:`landmark_tools` - Landmark-based registration validation * :doc:`contour_tools` - Contour processing * :doc:`image_conversion` - 4D image to 3D time-series utilities * :doc:`test_tools` - Baseline / result comparison helpers diff --git a/docs/api/workflows.rst b/docs/api/workflows.rst index 4a4bf3e6..d9c4e6a5 100644 --- a/docs/api/workflows.rst +++ b/docs/api/workflows.rst @@ -46,8 +46,10 @@ Available Workflows - Finetune uniGradICON on your own cohort and return the weights :class:`RegisterImagesICON` can load. -The PhysicsNeMo AI-surrogate workflows have their own section — see -:doc:`physicsnemo/index`. +The PhysicsNeMo AI-surrogate workflows — :class:`WorkflowTrainPhysicsNeMo`, +:class:`WorkflowInferPhysicsNeMo`, :class:`WorkflowInferMovement` and +:class:`WorkflowEvaluateMovement` — have their own section, since they need the +optional ``[physicsnemo]`` extra. See :doc:`physicsnemo/index`. Convert Image to USD ==================== diff --git a/docs/architecture.rst b/docs/architecture.rst index e7226ae6..ff54e452 100644 --- a/docs/architecture.rst +++ b/docs/architecture.rst @@ -87,10 +87,16 @@ Primary Workflows :mod:`physiotwin4d.vtk_to_usd` package exposes advanced file conversion primitives. +``WorkflowEvaluateMovement`` + Scores predicted motion per anatomic structure against the frames that were + actually acquired, reporting volume difference, Dice and surface RMSE on a + single isotropic evaluation grid. It wraps ``WorkflowInferMovement``, so it + measures whatever that produces. + AI Surrogate Workflows (PhysicsNeMo) ===================================== -The final tier of tutorials (``tutorial_08`` through ``tutorial_10``) turns a +The final tier of tutorials (``tutorial_08`` through ``tutorial_13``) turns a fitted statistical shape model into a trained AI physiological surrogate, replacing the explicit per-phase registration solve with a learned model at inference time: @@ -121,6 +127,25 @@ inference time: produced the training data, and able to predict stages that were never acquired. +``tutorial_11_lung_evaluate_physicsnemo.py`` + Scores the same prediction with ``WorkflowEvaluateMovement``, against the + *images* rather than against the registration Tutorial 10 compares to. It + carries the reference frame's labelmap into each gated time point through the + network's own deformation and compares it to the labelmap of the frame that + was acquired, reporting volume difference and surface RMSE per structure + (plus Dice per chamber in the ``duke_heart`` variant). + +``tutorial_12_lung_end_to_end_inference.py`` + Collapses the chain into one script: segment the reference frame, fit the + shape model to that patient, and infer every stage — no registration + anywhere, and nothing read from Tutorial 8. This is the shape the deployed + pipeline takes, and why it runs in minutes where Tutorial 8 runs in hours. + +``tutorial_13_heart_and_lung_motion.py`` + Drives *two* trained networks over a single static clinical CT, animating + respiratory and cardiac motion together on a scan that has no 4D acquisition + behind it at all. + These tutorials are thin drivers over the ``WorkflowTrainPhysicsNeMo`` / ``WorkflowInferPhysicsNeMo`` workflow classes; each workflow owns the data side (manifests, normalization, datasets, saving) and delegates the network to a diff --git a/docs/assets/example.gif b/docs/assets/example.gif deleted file mode 100644 index 3b1e2e99..00000000 --- a/docs/assets/example.gif +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:8b434d88c3681856ef855b6f9d9b9fc9d3096855e955b543f25494db36f48329 -size 1196946 diff --git a/docs/assets/Tutorial_03_heart_original.gif b/docs/assets/tutorial_03_heart_original.gif similarity index 100% rename from docs/assets/Tutorial_03_heart_original.gif rename to docs/assets/tutorial_03_heart_original.gif diff --git a/docs/assets/Tutorial_03_heart_recon.gif b/docs/assets/tutorial_03_heart_recon.gif similarity index 100% rename from docs/assets/Tutorial_03_heart_recon.gif rename to docs/assets/tutorial_03_heart_recon.gif diff --git a/docs/assets/tutorial_04_heart-2png.png b/docs/assets/tutorial_04_heart-2.png similarity index 100% rename from docs/assets/tutorial_04_heart-2png.png rename to docs/assets/tutorial_04_heart-2.png diff --git a/docs/assets/tutorial_08_duke_heart_def_mag.gif b/docs/assets/tutorial_08_duke_heart_def_mag.gif new file mode 100644 index 00000000..0d0df0f3 --- /dev/null +++ b/docs/assets/tutorial_08_duke_heart_def_mag.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4b8190d3770083d6ed1ef020466af6359697fb90c868f014bc7cb38d8ad22e0e +size 1442112 diff --git a/docs/assets/tutorial_08_lung.gif b/docs/assets/tutorial_08_lung.gif index 0811432a..f46a1d6c 100644 --- a/docs/assets/tutorial_08_lung.gif +++ b/docs/assets/tutorial_08_lung.gif @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:a990d8ef57e463c57c5a8b2f5fed9c61c2bbf2cfde61e3e450f713fb89711141 -size 4443886 +oid sha256:10a88e62663f8980a46208fdc687de5cbbd0e8f17aff01456d707912f290d675 +size 1179606 diff --git a/docs/assets/tutorial_09_duke_heart_deformation_magnitude.gif b/docs/assets/tutorial_09_duke_heart_deformation_magnitude.gif new file mode 100644 index 00000000..59d819bd --- /dev/null +++ b/docs/assets/tutorial_09_duke_heart_deformation_magnitude.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c421a14678114059b889a852d2b857d0415fbc97cb3015d1d182c436d1f6009f +size 869859 diff --git a/docs/assets/tutorial_09_duke_heart_motion.gif b/docs/assets/tutorial_09_duke_heart_motion.gif new file mode 100644 index 00000000..e713edd4 --- /dev/null +++ b/docs/assets/tutorial_09_duke_heart_motion.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b13fdf5b0b558074dacf2abd58de75f77a0ec76334d5c336c0f31c60df27a049 +size 6240228 diff --git a/docs/assets/tutorial_09_duke_heart_rmse.gif b/docs/assets/tutorial_09_duke_heart_rmse.gif new file mode 100644 index 00000000..700698b7 --- /dev/null +++ b/docs/assets/tutorial_09_duke_heart_rmse.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:96928a24b5e8fb6d444a166025de40c3e46b46b7d88e6f28d12054b822bfb4f1 +size 742650 diff --git a/docs/assets/tutorial_09_lung_deformation_magnitude.gif b/docs/assets/tutorial_09_lung_deformation_magnitude.gif new file mode 100644 index 00000000..edea1797 --- /dev/null +++ b/docs/assets/tutorial_09_lung_deformation_magnitude.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4b436bd00c1d8c73f28fd85a77b22a19dc870738ac36dcbb59728d4009c22bb0 +size 1641656 diff --git a/docs/assets/tutorial_09_lung_motion.gif b/docs/assets/tutorial_09_lung_motion.gif new file mode 100644 index 00000000..c7c70b92 --- /dev/null +++ b/docs/assets/tutorial_09_lung_motion.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f009e967c6889f2782b507ab2091add7445f5da01996a4135e2fec6edd2f03b8 +size 7102792 diff --git a/docs/assets/tutorial_09_lung_rmse.gif b/docs/assets/tutorial_09_lung_rmse.gif new file mode 100644 index 00000000..6324ac20 --- /dev/null +++ b/docs/assets/tutorial_09_lung_rmse.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:471b9db1e7cc53bd1618725fdb79575ee96a4ad5808be6ba4bc17c7a8096c3a0 +size 1598406 diff --git a/docs/assets/tutorial_10_duke_heart_motion_usd.gif b/docs/assets/tutorial_10_duke_heart_motion_usd.gif new file mode 100644 index 00000000..3f86bfef --- /dev/null +++ b/docs/assets/tutorial_10_duke_heart_motion_usd.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:def9da175690bbaba6e32b4ee4e5fd2c194d12ce2603747d165b3790dd818609 +size 5110071 diff --git a/docs/assets/tutorial_10_lung_motion_usd.gif b/docs/assets/tutorial_10_lung_motion_usd.gif new file mode 100644 index 00000000..5d9ba430 --- /dev/null +++ b/docs/assets/tutorial_10_lung_motion_usd.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:a8dde5b6b430b8ac83a8772591b7fef4d031039ddb4c3c8082f04382ab23eb94 +size 4979517 diff --git a/docs/assets/tutorial_11_duke_heart_stats.png b/docs/assets/tutorial_11_duke_heart_stats.png new file mode 100644 index 00000000..38ea189a --- /dev/null +++ b/docs/assets/tutorial_11_duke_heart_stats.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:224162182d91e167a350509ae674058912cd1a691ed35440e12ecd4c6e0e5a17 +size 20875 diff --git a/docs/assets/tutorial_11_duke_heart_volumes.png b/docs/assets/tutorial_11_duke_heart_volumes.png new file mode 100644 index 00000000..dfc4d71a --- /dev/null +++ b/docs/assets/tutorial_11_duke_heart_volumes.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:be856b93a4f42ef2847af9497e2b25e379d77427cd3654c0fce073ff9ae6a934 +size 136090 diff --git a/docs/assets/tutorial_11_lung_stats.png b/docs/assets/tutorial_11_lung_stats.png new file mode 100644 index 00000000..b7f96c4f --- /dev/null +++ b/docs/assets/tutorial_11_lung_stats.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3f4be84ccd994b1a84c696021f07cb7525a589525e6649b95159518279156278 +size 19174 diff --git a/docs/assets/tutorial_11_lung_volumes.png b/docs/assets/tutorial_11_lung_volumes.png new file mode 100644 index 00000000..81ba7a30 --- /dev/null +++ b/docs/assets/tutorial_11_lung_volumes.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8974bb522536d55ba95e5271e15f2870f64dd9cc8bb89f9754ba1ba6994d236d +size 87122 diff --git a/docs/assets/tutorial_12_duke_heart.gif b/docs/assets/tutorial_12_duke_heart.gif new file mode 100644 index 00000000..eff21abb --- /dev/null +++ b/docs/assets/tutorial_12_duke_heart.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f8506c5ebc3f604f35e0af9f5997a29f9787e372b99e219276ef2b3d222eb6c4 +size 4672178 diff --git a/docs/assets/tutorial_12_lung.gif b/docs/assets/tutorial_12_lung.gif new file mode 100644 index 00000000..9aee3b91 --- /dev/null +++ b/docs/assets/tutorial_12_lung.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ac377b18fa94a08d2e7dbee4461eeca1c6d99030d02ef2305e7c1d06d70d034b +size 3687280 diff --git a/docs/assets/tutorial_13_combined_motion.gif b/docs/assets/tutorial_13_combined_motion.gif new file mode 100644 index 00000000..27a64582 --- /dev/null +++ b/docs/assets/tutorial_13_combined_motion.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:fdd50797c4e63c62bcaa8687ffde257d1fed79a8c886738be751f27313eab203 +size 2169333 diff --git a/docs/cli_scripts/brain_vessel_modeling.rst b/docs/cli_scripts/brain_vessel_modeling.rst deleted file mode 100644 index a9869495..00000000 --- a/docs/cli_scripts/brain_vessel_modeling.rst +++ /dev/null @@ -1,16 +0,0 @@ -:orphan: - -============================ -Brain Vessel Modeling Status -============================ - -PhysioTwin4D does not currently install a brain-vessel modeling CLI command. -Earlier research experiments helped shape the toolkit, but they are not -supported usage examples for users or developers. - -For supported open-source release workflows, start with: - -* :doc:`overview` -* :doc:`heart_gated_ct` -* :doc:`vtk_to_usd` -* :doc:`../tutorials` diff --git a/docs/cli_scripts/download_data.rst b/docs/cli_scripts/download_data.rst index f37346f0..76b42f95 100644 --- a/docs/cli_scripts/download_data.rst +++ b/docs/cli_scripts/download_data.rst @@ -113,10 +113,13 @@ existing non-empty file, so re-running resumes an interrupted download. See Also ======== -* :doc:`../tutorials` — Tutorials 1 and 4-7 use ``Slicer-Heart-CT``, - ``KCL-Heart-Model``, and ``Chest-CT`` (the lung variant of Tutorial 7); - ``DirLab-4DCT`` (Tutorials 2 and 3) is manual-only, see - ``data/DirLab-4DCT/README.md``. +* :doc:`../tutorials` — ``Slicer-Heart-CT`` drives Heart Tutorials 1, 3 and 4; + ``KCL-Heart-Model`` drives Heart Tutorial 6; ``Chest-CT`` drives Lung + Tutorial 7 and Tutorial 13. ``DirLab-4DCT`` — Lung Tutorials 1, 2, 3, 4, 6, 8 + and 10-12, plus Heart Tutorial 7 — is manual-only, see + ``data/DirLab-4DCT/README.md``. ``Duke-Heart-4DLabelmaps``, which drives the + ten ``duke_heart`` variants, is being released soon; see + ``data/Duke-Heart-4DLabelmaps/README.md``. * :doc:`byod_tutorials` * :doc:`heart_gated_ct` * :doc:`overview` diff --git a/docs/cli_scripts/lung_gated_ct.rst b/docs/cli_scripts/lung_gated_ct.rst deleted file mode 100644 index 33c9b966..00000000 --- a/docs/cli_scripts/lung_gated_ct.rst +++ /dev/null @@ -1,15 +0,0 @@ -:orphan: - -===================== -Lung-Gated CT Status -===================== - -PhysioTwin4D does not currently install a dedicated lung-gated CT CLI command. -Respiratory 4D CT work is represented by the high-resolution reconstruction -workflow and Tutorial 3, which requires manually prepared DirLab-4DCT data. - -Use these supported resources instead: - -* :doc:`4dct_reconstruction` -* :doc:`../tutorials` -* :class:`physiotwin4d.WorkflowReconstructHighres4DCT` diff --git a/docs/developer/workflows.rst b/docs/developer/workflows.rst index 8d9bbe78..29cb18dc 100644 --- a/docs/developer/workflows.rst +++ b/docs/developer/workflows.rst @@ -26,6 +26,20 @@ Current Workflow Mapping - :class:`physiotwin4d.WorkflowFitStatisticalModelToPatient` * - ``physiotwin4d-reconstruct-highres-4d-ct`` - :class:`physiotwin4d.WorkflowReconstructHighres4DCT` + * - ``physiotwin4d-train-physicsnemo`` + - :class:`physiotwin4d.WorkflowTrainPhysicsNeMo` + * - ``physiotwin4d-infer-physicsnemo`` + - :class:`physiotwin4d.WorkflowInferPhysicsNeMo` + * - ``physiotwin4d-convert-image-4d-to-3d`` + - :class:`physiotwin4d.ConvertImage4DTo3D` (a converter, not a workflow) + * - ``physiotwin4d-download-data`` + - :class:`physiotwin4d.DataDownloadTools` (a utility, not a workflow) + * - ``physiotwin4d-visualize-pca-modes`` + - Reads a ``pca_model.json`` directly; no workflow class + +That is all eleven installed commands. Two workflow classes have no CLI +wrapper: :class:`physiotwin4d.WorkflowFinetuneICONRegistration` and +:class:`physiotwin4d.WorkflowEvaluateMovement`. Workflow Example ================ diff --git a/docs/faq.rst b/docs/faq.rst index f606d904..d8335e2f 100644 --- a/docs/faq.rst +++ b/docs/faq.rst @@ -68,7 +68,10 @@ torchvision, and torchaudio are sourced from What Python version is required? --------------------------------- -Python 3.11 or 3.12 are supported. +Python 3.10, 3.11 and 3.12 are supported. + +The one exception is the optional ``[physicsnemo]`` extra: ``nvidia-physicsnemo`` +requires Python >= 3.11, so the AI-surrogate tutorials need 3.11 or 3.12. Usage Questions =============== @@ -97,8 +100,13 @@ See :doc:`api/segmentation/index` for comparison. Which registration method should I use? ---------------------------------------- -* **ICON**: Recommended for cardiac/lung (fast, GPU) +* **Greedy**: CPU-capable classical deformable registration; what Tutorials 1 + and 3 use by default +* **ICON**: Recommended for cardiac/lung (fast, GPU), and finetunable on your + own cohort — see Tutorial 2 * **ANTs**: Best for brain imaging and general purpose +* **Greedy+ICON** (``RegisterImagesGreedyICON``, a ``RegisterImagesChain`` + preset): Greedy for the coarse alignment, ICON for the refinement See :doc:`api/registration/index` for comparison. diff --git a/docs/index.rst b/docs/index.rst index 2ab48998..3fd90f36 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -80,6 +80,21 @@

Predict Motion With the Surrogate

Replace the registration solve with one forward pass, then export to USD.

+ + 11 +

Score the Surrogate Against the Images

+

Volume and surface RMSE per lobe, plus Dice per chamber, on the held-out case.

+ + + 12 +

The Whole Inference Pipeline in One Script

+

Go from a gated series to an animated prediction without registering a single phase.

+ + + 13 +

Breathe and Beat a Static Clinical CT

+

Animate one routine breath-hold scan with both rhythms, from two networks at once.

+
diff --git a/docs/installation.rst b/docs/installation.rst index 58353c1c..0b461a19 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -206,14 +206,18 @@ Expected output: Command-Line Tools ================== -PhysioTwin4D provides command-line interfaces that should be available after installation: +PhysioTwin4D installs eleven command-line tools, each prefixed +``physiotwin4d-``. There is no bare ``physiotwin4d`` command; check the install +with any one of them: .. code-block:: bash # Check CLI is available - physiotwin4d --help + physiotwin4d-download-data --help physiotwin4d-convert-image-to-usd --help +See :doc:`cli_scripts/overview` for the full list. + GPU Setup ========= @@ -229,6 +233,20 @@ A plain ``pip install physiotwin4d`` installs a CPU-only build. It runs without error but emits a ``UserWarning`` at import time and will be significantly slower than a GPU-enabled install. +Optional External Software +-------------------------- + +One segmentation backend is not a Python dependency and cannot be installed +with pip: + +* **Synopsys Simpleware Medical** — required by + :class:`~physiotwin4d.SegmentHeartSimpleware` and + :class:`~physiotwin4d.SegmentHeartSimplewareTrimmedBranches`, and therefore by + Tutorial 13, which segments the heart it fits. It needs a local licensed + installation; see :doc:`api/segmentation/simpleware`. Everything else in the + toolkit runs without it, and the ``requires_simpleware`` tests skip cleanly + when it is absent. + If CUDA is not yet installed, download the CUDA Toolkit from `NVIDIA's website `_, then verify: diff --git a/docs/quickstart.rst b/docs/quickstart.rst index f40a3965..25d401fb 100644 --- a/docs/quickstart.rst +++ b/docs/quickstart.rst @@ -84,19 +84,27 @@ Which dataset each tutorial needs: - Heart Tutorials 1, 3, 4 * - ``DirLab-4DCT`` - Manual - - Lung Tutorials 1, 2, 3, 4, 6, 8, and Heart Tutorial 7 + - Lung Tutorials 1, 2, 3, 4, 6, 8, 10, 11, 12, and Heart Tutorial 7 * - ``KCL-Heart-Model`` - CLI - Heart Tutorial 6 * - ``Chest-CT`` - CLI - - Lung Tutorial 7 + - Lung Tutorial 7, and Tutorial 13 + * - ``Duke-Heart-4DLabelmaps`` + - Releasing soon + - The ten ``duke_heart`` variants, Tutorials 2 and 4 through 12 * - ``CHOP-Valve4D`` - CLI - No tutorial - used by the valve experiments under ``experiments/`` -Tutorials 5, 9 and 10 need no dataset of their own: they consume the outputs of -Tutorials 4, 8 and 9 respectively. +Tutorials 5 and 9 need no dataset of their own: they consume the outputs of +Tutorials 4 and 8 respectively. + +``Duke-Heart-4DLabelmaps`` is scheduled for public release soon. Until then the +ten ``duke_heart`` variants cannot be run; contact Stephen Aylward +(saylward@nvidia.com) to request access, and see +``data/Duke-Heart-4DLabelmaps/README.md``. ``DirLab-4DCT`` is the one dataset with no automatic downloader: DIR-Lab distributes each case individually and may require registration, so download it @@ -110,9 +118,9 @@ as a plain script: python tutorials/tutorial_01_heart_gated_ct_to_usd.py -:doc:`tutorials` is the full guide — ten tutorials with previews of what each -one produces, the run order, and per-tutorial notes on pointing them at your -own data. The rest of this page is the same functionality as a CLI call and as +:doc:`tutorials` is the full guide — thirteen numbered stages with previews of +what each one produces, the run order, and per-tutorial notes on pointing them +at your own data. The rest of this page is the same functionality as a CLI call and as a Python API call, for when you would rather not start from a script. Basic Workflow @@ -301,7 +309,7 @@ layouts for every dataset. DirLab-4DCT data is manual-only; see ``data/DirLab-4DCT/README.md``. It drives the whole lung pipeline — Lung Tutorials 1, 2, 3, 4, 6 and 8, plus Heart -Tutorial 7 — which then feeds the AI-surrogate Tutorials 9 and 10. Those two +Tutorial 7 — which then feeds the AI-surrogate Tutorials 9 through 12. Those additionally require the optional ``physicsnemo`` extra (``pip install "physiotwin4d[physicsnemo]"``, plus ``torch-geometric`` for the MeshGraphNet); PhysicsNeMo itself requires diff --git a/docs/testing.rst b/docs/testing.rst index 33890fe9..cb275f2c 100644 --- a/docs/testing.rst +++ b/docs/testing.rst @@ -39,6 +39,16 @@ once. The self-hosted CI GPU runner uses it (after installing pytest tests/ -v --run-all +Three further markers are declared but not gated by a flag: ``unit`` and +``integration`` are descriptive, and ``xdist_group`` is a ``pytest-xdist`` +builtin used to keep related tests on one worker. Tests carrying only these +markers always run. + +``tests/test_tutorials.py`` holds the ``tutorial`` bucket. It is not +parametrized over the tutorials directory — it is one hand-written class per +covered script, currently 9 of the 29 tutorial scripts, so adding a tutorial +does not automatically add a test. + Test Categories =============== @@ -60,6 +70,8 @@ Specific Areas pytest tests/test_contour_tools.py -v pytest tests/test_transform_tools.py -v pytest tests/test_image_tools.py -v + pytest tests/test_workflow_train_physicsnemo.py -v + pytest tests/test_workflow_evaluate_movement.py -v Real Data and GPU Tests ======================= diff --git a/docs/tutorials.rst b/docs/tutorials.rst index 68e797b3..bd68b3db 100644 --- a/docs/tutorials.rst +++ b/docs/tutorials.rst @@ -11,9 +11,9 @@ Tutorials

PhysioTwin4D tutorials

From a CT scan to an animated digital twin

- Ten numbered stages across 24 Python scripts, 16 of them runnable today: - the eight duke_heart variants wait on a dataset that is not - public yet. + Thirteen numbered stages across 29 Python scripts, 19 of them runnable + today: the ten duke_heart variants wait on a dataset that + is being released soon. Each one drives the real workflow classes end-to-end on downloadable data, shows what it produced, and ends with the handful of constants to change so it runs on your own scans. @@ -46,11 +46,17 @@ relative to the current working directory: physiotwin4d-download-data Chest-CT --directory data/Chest-CT That covers Heart Tutorials 1, 3, 4 and 6 (``Slicer-Heart-CT`` and -``KCL-Heart-Model``) and Lung Tutorial 7 (``Chest-CT``). ``DirLab-4DCT`` — used -by Lung Tutorials 1, 2, 3, 4, 6 and 8, and by Heart Tutorial 7 — is **not** -auto-downloaded: DIR-Lab distributes each case individually and may require -registration. Tutorials 5, 9 and 10 need no dataset of their own; they consume -the outputs of Tutorials 4, 8 and 9. See ``data/DirLab-4DCT/README.md``, and +``KCL-Heart-Model``) and Lung Tutorial 7 (``Chest-CT``), which Tutorial 13 also +animates. ``DirLab-4DCT`` — used by Lung Tutorials 1, 2, 3, 4, 6, 8, 10, 11 and +12, and by Heart Tutorial 7 — is **not** auto-downloaded: DIR-Lab distributes +each case individually and may require registration. + +Tutorials 5 and 9 need no dataset of their own; they consume the outputs of +Tutorials 4 and 8. ``Duke-Heart-4DLabelmaps`` drives the ten ``duke_heart`` +variants, which form their own chain from Tutorial 4 through Tutorial 12; it is +being released soon, and until then access can be requested from Stephen Aylward +(saylward@nvidia.com). See ``data/DirLab-4DCT/README.md``, +``data/Duke-Heart-4DLabelmaps/README.md``, and :doc:`cli_scripts/download_data` for every dataset's size and source. **3. Know where output lands.** Every tutorial writes to @@ -120,6 +126,24 @@ second run is cheap and later tutorials pick up earlier results automatically.

Replace the registration solve with one forward pass, then export to USD.

Tutorials 8 and 9 output + + 11 +

Score the Surrogate Against the Images

+

Volume and surface RMSE per lobe, plus Dice per chamber, on the held-out case.

+ Tutorials 8, 9 and 10 output + + + 12 +

The Whole Inference Pipeline in One Script

+

Go from a gated series to an animated prediction without registering a single phase.

+ Tutorials 6 and 9 output + + + 13 +

Breathe and Beat a Static Clinical CT

+

Animate one routine breath-hold scan with both rhythms, from two networks at once.

+ Chest-CT · Tutorials 7 and 9 output +
Recommended Run Order @@ -129,7 +153,7 @@ Tutorials are straightforward Python scripts: run one with ``python tutorials/tutorial_01_heart_gated_ct_to_usd.py``, or open it in your editor and read it top to bottom. Numbers 1, 4 and 5 are the fastest way to see the toolkit -work end-to-end; 6 through 10 build the statistical-model and AI-surrogate +work end-to-end; 6 through 13 build the statistical-model and AI-surrogate pipeline on top. 1. **Tutorial 1** — after downloading Slicer-Heart-CT. @@ -145,6 +169,12 @@ pipeline on top. 8. **Tutorial 8** — after Tutorial 6 (lung); Tutorial 2 optional. 9. **Tutorial 9** — after Tutorial 8, whose fitted meshes it trains on. 10. **Tutorial 10** — after Tutorial 9, whose checkpoint it loads. +11. **Tutorial 11** — after Tutorial 9. The lung variant segments every gated + frame of the held-out case, so it needs a GPU and the segmentation weights. +12. **Tutorial 12** — after Tutorial 6 and Tutorial 9 for its anatomy; it fits + the model to the patient itself, so nothing is read from Tutorial 8. +13. **Tutorial 13** — after Tutorial 7 (lung) and Tutorial 9 for both anatomies. + It also needs Simpleware Medical, which segments the heart it fits. Tutorial 1: Gated 4D CT to Animated USD ======================================= @@ -234,8 +264,10 @@ Script with lung ones. The per-organ values live in ``tutorials/parameters_lung_ct_dirlab.py`` for the lung variant and ``tutorials/parameters_duke_heart_labelmaps.py`` for this one. This is a - ``duke_heart`` tutorial: Duke-Heart-4DLabelmaps is not publicly available - yet, so it cannot be run — see ``data/Duke-Heart-4DLabelmaps/README.md``. + ``duke_heart`` tutorial: Duke-Heart-4DLabelmaps is being released soon (see + `Before You Start`_), and until then access can be requested from Stephen + Aylward (saylward@nvidia.com) — see + ``data/Duke-Heart-4DLabelmaps/README.md``. Workflow :class:`~physiotwin4d.WorkflowFinetuneICONRegistration`, then @@ -325,13 +357,13 @@ Requirements ``[30, 15, 7, 3]``. Preview - .. figure:: assets/Tutorial_03_heart_original.gif + .. figure:: assets/tutorial_03_heart_original.gif :alt: Acquired cardiac phases :width: 90% The acquired cardiac phases. - .. figure:: assets/Tutorial_03_heart_recon.gif + .. figure:: assets/tutorial_03_heart_recon.gif :alt: Cardiac phases reconstructed at the reference resolution :width: 90% @@ -386,6 +418,10 @@ Script ``tutorials/tutorial_04_lung_ct_to_vtk.py`` + ``tutorials/tutorial_04_duke_heart_labelmap_to_vtk.py`` — starts from gated + labelmaps rather than CT, and also extracts tetrahedral meshes. Needs + Duke-Heart-4DLabelmaps (see `Before You Start`_). + Workflow :class:`~physiotwin4d.WorkflowConvertImageToVTK` with :class:`~physiotwin4d.SegmentChestTotalSegmentatorWithContrast` (heart) or @@ -411,6 +447,13 @@ Preview The same workflow on a DIR-Lab respiratory case. + .. figure:: assets/tutorial_04_duke_heart.png + :alt: Heart surfaces extracted from a gated Duke labelmap + :width: 90% + + The ``duke_heart`` variant, which starts from a gated labelmap rather than + a CT and also writes tetrahedral meshes. + Inner API usage .. code-block:: python @@ -451,6 +494,10 @@ Tutorial 5: VTK Surfaces to Animated USD Script ``tutorials/tutorial_05_heart_vtk_to_usd.py`` + ``tutorials/tutorial_05_duke_heart_vtk_to_usd.py`` — the 4D counterpart, + animating Tutorial 4 (duke heart)'s per-phase surfaces. Needs + Duke-Heart-4DLabelmaps (see `Before You Start`_). + Workflow :class:`~physiotwin4d.WorkflowConvertVTKToUSD`. @@ -516,6 +563,10 @@ Script ``tutorials/tutorial_06_lung_create_statistical_model.py`` + ``tutorials/tutorial_06_duke_heart_create_statistical_model.py`` — builds the + cardiac model the ``duke_heart`` surrogate chain trains against. Needs + Duke-Heart-4DLabelmaps (see `Before You Start`_). + Workflow :class:`~physiotwin4d.WorkflowCreateStatisticalModel`; the lung variant first builds an unbiased atlas with @@ -587,6 +638,10 @@ Script ``tutorials/tutorial_07_lung_fit_statistical_model_to_patient.py`` + ``tutorials/tutorial_07_duke_heart_fit_statistical_model_to_patient.py`` — + fits the Tutorial 6 (duke heart) model. Needs Duke-Heart-4DLabelmaps (see + `Before You Start`_). + Workflow :class:`~physiotwin4d.WorkflowFitStatisticalModelToPatient`. @@ -654,6 +709,11 @@ Tutorial 8: Propagate the Shape Model Through 4D Script ``tutorials/tutorial_08_lung_fit_model_to_4d_patients.py`` + ``tutorials/tutorial_08_duke_heart_fit_model_to_4d_patients.py`` — the same + fit-then-propagate pass over cardiac phases, using + :class:`~physiotwin4d.RegisterModelsDistanceMaps` in place of the image + registration. Needs Duke-Heart-4DLabelmaps (see `Before You Start`_). + Workflow :class:`~physiotwin4d.WorkflowFitStatisticalModelToPatient` at the reference phase, then :class:`~physiotwin4d.WorkflowReconstructHighres4DCT` to carry @@ -676,6 +736,13 @@ Preview The fitted shape-model surface propagated across the phases of a DIR-Lab case. + .. figure:: assets/tutorial_08_duke_heart_def_mag.gif + :alt: Deformation magnitude over the propagated heart surface + :width: 90% + + The ``duke_heart`` variant, coloured by deformation magnitude across the + cardiac phases. + Inner API usage .. code-block:: python @@ -718,6 +785,10 @@ Tutorial 9: Train a PhysicsNeMo Surrogate Script ``tutorials/tutorial_09_lung_train_physicsnemo_mgn.py`` + ``tutorials/tutorial_09_duke_heart_train_physicsnemo_mgn.py`` — trains the + cardiac network Tutorial 13 uses for its heartbeat. Needs + Duke-Heart-4DLabelmaps (see `Before You Start`_). + Workflow :class:`~physiotwin4d.WorkflowTrainPhysicsNeMo` driving :class:`~physiotwin4d.TrainPhysicsNeMoMGN`, then @@ -740,12 +811,34 @@ Requirements Python >= 3.11. 1500 epochs by default. Preview - .. figure:: assets/example.gif - :alt: Tutorial 9 output preview (capture pending) - :width: 60% + .. figure:: assets/tutorial_09_lung_motion.gif + :alt: Predicted lung motion across the respiratory cycle + :width: 90% - Capture pending — the tutorial writes ``predicted_surface.png`` and - ``rmse_surface.png`` when it runs. + The held-out lung case, predicted at every stage by the trained network. + + .. figure:: assets/tutorial_09_lung_rmse.gif + :alt: Per-vertex RMSE of the predicted lung surface + :width: 90% + + The same surface coloured by per-vertex error against the registration + that produced the training data. + + .. figure:: assets/tutorial_09_lung_deformation_magnitude.gif + :alt: Deformation magnitude over the lung surface + :width: 90% + + Deformation magnitude, which is what the error above should be read + against — the largest errors sit where the motion is largest. + + .. figure:: assets/tutorial_09_duke_heart_motion.gif + :alt: Predicted heart motion across the cardiac cycle + :width: 90% + + The ``duke_heart`` variant over a cardiac cycle, with its own RMSE and + deformation-magnitude captures in + ``tutorial_09_duke_heart_rmse.gif`` and + ``tutorial_09_duke_heart_deformation_magnitude.gif``. Inner API usage .. code-block:: python @@ -790,6 +883,10 @@ Tutorial 10: Predict Motion With the Surrogate Script ``tutorials/tutorial_10_lung_infer_physicsnemo_mgn.py`` + ``tutorials/tutorial_10_duke_heart_infer_physicsnemo_mgn.py`` — the same + prediction over a cardiac cycle. Needs Duke-Heart-4DLabelmaps (see + `Before You Start`_). + Workflow :class:`~physiotwin4d.WorkflowInferPhysicsNeMo` for the raw prediction, :class:`~physiotwin4d.WorkflowInferMovement` to turn it back into geometry, @@ -799,16 +896,22 @@ Dataset Tutorial 8's fitted surfaces for one case, and Tutorial 9's checkpoint. Requirements - The ``[physicsnemo]`` extra; otherwise trivial — one forward pass replaces - the per-phase registration solve that produced the training data. + The ``[physicsnemo]`` extra; otherwise trivial — one forward pass per stage + replaces the per-phase registration solve that produced the training data. Preview - .. figure:: assets/example.gif - :alt: Tutorial 10 output preview (capture pending) - :width: 60% + .. figure:: assets/tutorial_10_lung_motion_usd.gif + :alt: Animated USD of the predicted lung motion + :width: 90% - Capture pending — the tutorial writes ``predicted_surface.png`` and - ``ground_truth_surface.png`` when it runs. + The exported USD scene, played back over the respiratory cycle — every + frame a forward pass rather than a registration solve. + + .. figure:: assets/tutorial_10_duke_heart_motion_usd.gif + :alt: Animated USD of the predicted heart motion + :width: 90% + + The ``duke_heart`` variant over a cardiac cycle. Inner API usage .. code-block:: python @@ -817,12 +920,17 @@ Inner API usage model_directory=model_dir, epoch=epoch, ) - infer_result = WorkflowInferMovement(infer_workflow).predict_single( + infer_result = WorkflowInferMovement(infer_workflow).process_time_series( shape_parameters=pca_file, - stage=test_stage, - reference_mesh=reference_file, - ground_truth=ground_truth_file, + stages=stages, output_directory=output_dir, + reference_mesh=reference_file, + ground_truth=phase_files, + reference_image=itk.imread(str(reference_ct_file)), + warp_interpolation="linear", + warp_background_value=-1000.0, + usd_project_name=f"{case_id}_mgn_motion", + anatomy_type="lung", ) Run @@ -831,18 +939,266 @@ Run python tutorials/tutorial_10_lung_infer_physicsnemo_mgn.py Outputs - The predicted surface, its error statistics against the ground-truth phase - in millimetres, and a USD scene, under - ``tutorials/output/tutorial_10_lung_mgn//``. + One predicted surface and one warped CT per stage, one animated USD across + all of them, and ``statistics_per_stage.csv`` with the mm error against each + acquired phase, under ``tutorials/output/tutorial_10_lung_mgn//``. Adapt to your data - Change ``case_id`` and ``stage_fraction`` to predict a different subject, or - a stage that was never acquired — which is the point of the surrogate. Omit + Change ``case_id`` to predict a different subject, or pass ``stages`` that + were never acquired — which is the point of the surrogate. Omit ``reference_mesh`` to displace the mesh reconstructed from the PCA - coefficients alone, needing no per-subject geometry at all. Use + coefficients alone, needing no per-subject geometry at all. Omit + ``reference_image`` to write meshes without warping anything. Use :class:`~physiotwin4d.WorkflowInferPhysicsNeMo` on its own to get the raw target array when your model predicts something other than displacement. +Tutorial 11: Score the Surrogate Against the Images +=================================================== + +Script + ``tutorials/tutorial_11_lung_evaluate_physicsnemo.py`` + + ``tutorials/tutorial_11_duke_heart_evaluate_physicsnemo.py`` + +Workflow + :class:`~physiotwin4d.WorkflowEvaluateMovement`, driving + :class:`~physiotwin4d.WorkflowInferMovement` and, for the lung variant, + :class:`~physiotwin4d.SegmentNVSegmentCTMRI`. + +Dataset + The gated sequence itself — DIR-Lab for the lung, Duke-Heart-4DLabelmaps for + the heart — plus Tutorial 8's fitted surface and Tutorial 9's checkpoint for + the held-out case. + +Requirements + The ``[physicsnemo]`` extra. The lung variant also segments every gated frame + on first run, so it needs a GPU and the segmentation weights; the labelmaps + are cached, and a re-run skips them. + +Preview + .. figure:: assets/tutorial_11_lung_volumes.png + :alt: Acquired and predicted lobe volumes across the respiratory cycle + :width: 90% + + ``volume_vs_stage.png`` for the held-out lung case: acquired volume solid, + predicted dashed, one pair per lobe across every gated stage. + + .. figure:: assets/tutorial_11_lung_stats.png + :alt: Per-lobe volume difference and surface RMSE for the lung case + :width: 90% + + The same run summarised per lobe. No Dice column — see the note below. + + .. figure:: assets/tutorial_11_duke_heart_stats.png + :alt: Per-chamber Dice, volume difference and surface RMSE for the heart + :width: 90% + + The ``duke_heart`` variant, which does report Dice per chamber, alongside + its own ``tutorial_11_duke_heart_volumes.png``. + +Inner API usage + .. code-block:: python + + evaluate = WorkflowEvaluateMovement( + movement_workflow=WorkflowInferMovement(infer_workflow), + label_names=lobe_names, + ) + result = evaluate.process( + case_id=case_id, + shape_parameters=pca_file, + reference_mesh=reference_mesh_file, + reference_labelmap=itk.imread(str(reference_labelmap_file)), + ground_truth_labelmaps=ground_truth_labelmaps, + output_directory=output_dir, + evaluation_spacing_mm=2.0, + include_dice=False, + ) + +Run + .. code-block:: bash + + python tutorials/tutorial_11_lung_evaluate_physicsnemo.py + +Outputs + ``evaluation_report.md``, ``evaluation_metrics.csv`` and + ``volume_vs_stage.png`` under ``tutorials/output/tutorial_11_lung//``, + carrying volume difference and surface RMSE per lobe at every gated stage; + the duke variant adds Dice per chamber. The plot traces each structure's + acquired and predicted volume across the stages. Report and CSV both record + the hold-out case name, its shape parameters, and the network weights path + with its dates, so a number can be traced back to the run that produced it. + + The lung variant passes ``include_dice=False``. Dice is an overlap fraction, + so a lobe that moves a few millimeters against its own bulk scores over 0.96 + however well or badly the motion is predicted; the column would describe the + lobe rather than the model. Chambers change shape enough over a heartbeat for + it to discriminate, so the duke variant keeps it. + +Adapt to your data + Change ``LOBE_LABEL_IDS`` (or ``HEART_LABEL_IDS``) to score a different set + of structures — any label your segmenter writes and your reference frame + contains. Raise ``evaluation_spacing_mm`` if the deformation fields do not + fit in memory; lower it to resolve a thin wall, at the cost of its cube. + +Tutorial 12: The Whole Inference Pipeline in One Script +======================================================= + +Script + ``tutorials/tutorial_12_lung_end_to_end_inference.py`` + + ``tutorials/tutorial_12_duke_heart_end_to_end_inference.py`` + +Workflow + :class:`~physiotwin4d.WorkflowConvertImageToVTK` (lung) or + :class:`~physiotwin4d.ContourTools` (heart), + :class:`~physiotwin4d.WorkflowFitStatisticalModelToPatient`, then + :meth:`~physiotwin4d.WorkflowInferMovement.process_time_series`. + +Dataset + The gated sequence alone — DIR-Lab for the lung, Duke-Heart-4DLabelmaps for + the heart — plus the Tutorial 6 shape model and the Tutorial 9 checkpoint. + Unlike Tutorial 10, nothing is read from Tutorial 8: this script fits the + model to the patient itself, so the chain from image to animation runs in one + place. + +Requirements + The ``[physicsnemo]`` extra. The output directory is emptied at the start of + every run, so nothing is reused and the reported runtimes are the whole + pipeline's. Neither variant registers a phase — that is what the network + replaces, and it is why this runs in minutes where Tutorial 8 runs in hours. + +Preview + .. figure:: assets/tutorial_12_lung.gif + :alt: Lung motion predicted end-to-end from a gated series + :width: 90% + + The whole chain on one DIR-Lab case: segment, fit, infer, animate — no + phase registered anywhere in it. + + .. figure:: assets/tutorial_12_duke_heart.gif + :alt: Heart motion predicted end-to-end from gated labelmaps + :width: 90% + + The ``duke_heart`` variant, starting from gated labelmaps instead of CT. + +Inner API usage + .. code-block:: python + + # The fit puts the model in this patient: coefficients condition the + # network, and the fitted surface is what its displacements move. + fit = WorkflowFitStatisticalModelToPatient( + template_model=pca_mean_surface, + patient_models=[lung_surface], + patient_image=reference_image, + patient_labelmap=lung_labelmap, + ) + fit.set_use_pca_registration( + use_pca_registration=True, + pca_model=pca_model, + number_of_pca_components=6, + use_surface=False, + ) + fit_result = fit.process() + + infer_result = WorkflowInferMovement(infer_workflow).process_time_series( + shape_parameters=pca_coefficients_file, + stages=stages, + output_directory=output_dir, + reference_mesh=reference_mesh_file, + reference_image=reference_image, + usd_project_name=f"{case_id}_mgn_motion", + anatomy_type="lung", + ) + +Run + .. code-block:: bash + + python tutorials/tutorial_12_lung_end_to_end_inference.py + +Outputs + Under ``tutorials/output/tutorial_12_lung//`` (or + ``tutorial_12_duke_heart``): the patient's fitted + ``_ssm_surface.vtp`` and ``_ssm_pca_coefficients.json``, one + predicted ``*_pred.vtp`` surface and one ``*_warped.mha`` volume per stage, + ``_mgn_motion.usd`` animating the whole cycle, and + ``_runtimes.csv`` timing each step of the run. + +Adapt to your data + Point the script at any case of the same cohort by changing ``case_id``; the + stages come from the filenames, so a sequence with a different number of + phases needs no other change. To predict stages the acquisition never + sampled, pass your own ``stages`` list — the network is continuous in stage, + and nothing downstream requires a matching image. + +Tutorial 13: Breathe and Beat a Static Clinical CT +================================================== + +Script + ``tutorials/tutorial_13_heart_and_lung_motion.py`` + +Workflow + :class:`~physiotwin4d.WorkflowInferMovement` over both Tutorial 9 networks, + :class:`~physiotwin4d.WorkflowFitStatisticalModelToPatient` for the heart fit, + and :class:`~physiotwin4d.ConvertVTKToUSD` with + :class:`~physiotwin4d.USDAnatomyTools` for the animation. + +Dataset + ``data/Chest-CT/Chest-CT.mha``, one routine breath-hold scan, plus Tutorial 7 + (lung)'s fit of it and both Tutorial 9 checkpoints. No 4D acquisition is + involved: every deformation comes from a network, none from a registration. + +Requirements + The ``[physicsnemo]`` extra, and Simpleware Medical for the heart + segmentation. Both segmentations and the heart fit are cached, so a re-run + goes straight to inference. Budget disk: 100 combined frames, each with its + own warped CT and labelmap, come to roughly 43 GB. + +Preview + .. figure:: assets/tutorial_13_combined_motion.gif + :alt: Combined heart and lung motion on a static clinical CT + :width: 90% + + ``heart_and_lung_motion.usd``: one routine breath-hold scan, breathing and + beating at once, with every deformation coming from a network and none + from a registration. + +Inner API usage + .. code-block:: python + + infer = WorkflowInferMovement( + WorkflowInferPhysicsNeMo(model_directory=lung_model_dir) + ) + # "forward" moves mesh vertices; "inverse" is what resampling an image + # into the stage's frame needs. + field = infer.create_deformation_field( + shape_parameters=lung_coefficients_file, + stage=0.0, + reference_image=patient_image, + reference_mesh=lung_reference_mesh_file, + direction="forward", + ) + transform = TransformTools().smooth_deformation_field_transform( + field["deformation_field"], 15.0, field["weight_image"] + ) + +Run + .. code-block:: bash + + python tutorials/tutorial_13_heart_and_lung_motion.py + +Outputs + Under ``tutorials/output/tutorial_13_heart_and_lung/``: one 4D USD per rhythm + (``breathing_lungs.usd``, ``beating_heart.usd``), 100 combined frames as VTP + plus ``heart_and_lung_motion.usd`` split by anatomy and painted with organ + materials, and the CT and labelmap warped by the same per-frame deformation. + +Adapt to your data + Point ``patient_image_file`` at your own chest CT and rerun Tutorial 7 (lung) + on it to get the lung fit; the heart fit happens inside this script. Change + ``cardiac_cycles_per_phase`` to re-time the heartbeat against the breath, and + the two ``*_sigma_mm`` values to change how far each rhythm's surface motion + is carried into the surrounding tissue. + Where to Go Next ================ diff --git a/docs/viewing_usd.rst b/docs/viewing_usd.rst index b5662e4f..2d9c78dc 100644 --- a/docs/viewing_usd.rst +++ b/docs/viewing_usd.rst @@ -2,7 +2,8 @@ Viewing USD Files ================== -Every USD-producing workflow in PhysioTwin4D — Tutorials 1, 5 and 10, and the +Every USD-producing workflow in PhysioTwin4D — Tutorials 1, 5, 10, 12 and 13, +and the ``physiotwin4d-convert-image-to-usd`` and ``physiotwin4d-convert-vtk-to-usd`` commands — writes an OpenUSD scene: anatomy split into per-organ prims, painted with OmniSurface materials, and time-sampled when the input was a series. To @@ -83,8 +84,9 @@ coordinate and unit details. Before USD: viewing the meshes directly ======================================= -The intermediate ``.vtp`` and ``.vtu`` files that Tutorials 4, 6, 7, 8 and 9 -write need no USD tooling at all — PyVista, already a dependency, opens them: +The intermediate ``.vtp`` and ``.vtu`` files that Tutorials 4, 6, 7, 8, 9, 10, +11 and 12 write need no USD tooling at all — PyVista, already a dependency, +opens them: .. code-block:: python diff --git a/pyproject.toml b/pyproject.toml index a4d94979..dd3a51a2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,7 +5,7 @@ build-backend = "setuptools.build_meta" [project] name = "physiotwin4d" version = "2026.07.3" -description = "Methods, workflows, tutorials, and CLI for creating personalized physiological digital twins from 3D medical images" +description = "Methods, workflows, tutorials, and CLI for creating personalized physiological digital twins from 3D/4D medical images" authors = [ {name = "Stephen R. Aylward", email = "saylward@nvidia.com"} ] @@ -50,7 +50,11 @@ keywords = [ "ICON", "ANTS", "physiological-motion", - "4d-visualization" + "4d-visualization", + "digital-twin", + "statistical-shape-model", + "physicsnemo", + "meshgraphnet" ] dependencies = [ # Core medical imaging @@ -353,8 +357,13 @@ module = [ "tutorial_08_lung_fit_model_to_4d_patients", "tutorial_09_duke_heart_train_physicsnemo_mgn", "tutorial_09_lung_train_physicsnemo_mgn", - "tutorial_10_duke_heart_infer_physicsnemo", + "tutorial_10_duke_heart_infer_physicsnemo_mgn", "tutorial_10_lung_infer_physicsnemo_mgn", + "tutorial_11_duke_heart_evaluate_physicsnemo", + "tutorial_11_lung_evaluate_physicsnemo", + "tutorial_12_duke_heart_end_to_end_inference", + "tutorial_12_lung_end_to_end_inference", + "tutorial_13_heart_and_lung_motion", ] disable_error_code = ["import-not-found", "import-untyped"] diff --git a/src/physiotwin4d/__init__.py b/src/physiotwin4d/__init__.py index 5e3108f9..86431ad4 100644 --- a/src/physiotwin4d/__init__.py +++ b/src/physiotwin4d/__init__.py @@ -93,6 +93,7 @@ ) from .workflow_infer_physicsnemo import WorkflowInferPhysicsNeMo from .workflow_infer_movement import WorkflowInferMovement +from .workflow_evaluate_movement import WorkflowEvaluateMovement from .infer_physicsnemo_base import InferPhysicsNeMoBase from .infer_physicsnemo_mgn import InferPhysicsNeMoMGN from .infer_physicsnemo_mlp import InferPhysicsNeMoMLP @@ -114,6 +115,7 @@ "WorkflowTrainPhysicsNeMo", "WorkflowInferPhysicsNeMo", "WorkflowInferMovement", + "WorkflowEvaluateMovement", # Training method classes "TrainPhysicsNeMoBase", "TrainPhysicsNeMoMGN", diff --git a/src/physiotwin4d/cli/__init__.py b/src/physiotwin4d/cli/__init__.py index 35f17949..06837878 100644 --- a/src/physiotwin4d/cli/__init__.py +++ b/src/physiotwin4d/cli/__init__.py @@ -2,10 +2,14 @@ __all__ = [ "convert_image_to_usd", + "convert_image_to_vtk", "convert_image_4d_to_3d", "convert_vtk_to_usd", "create_statistical_model", "download_data", "fit_statistical_model_to_patient", + "infer_physicsnemo", + "reconstruct_highres_4d_ct", + "train_physicsnemo", "visualize_pca_modes", ] diff --git a/src/physiotwin4d/contour_tools.py b/src/physiotwin4d/contour_tools.py index dc35c4ac..0584982c 100644 --- a/src/physiotwin4d/contour_tools.py +++ b/src/physiotwin4d/contour_tools.py @@ -369,7 +369,11 @@ def extract_label_surfaces( merged_ids = np.asarray(merged.cell_data["LabelId"]) surfaces: dict[int, pv.PolyData] = {} for label_id in label_ids: - cell_ids: list[int] = np.flatnonzero(merged_ids == label_id).tolist() + # Kept as an array rather than a list: a label with no cells gives an + # empty selection, and an empty list has no integer dtype for + # extract_cells to recognize it by. Empty here means no surface, + # which is what the next branch reports. + cell_ids = np.flatnonzero(merged_ids == label_id) surface = self.extract_surface(merged.extract_cells(cell_ids)).triangulate() if surface.n_cells == 0: # A label smaller than the isotropic grid loses its vote to its diff --git a/src/physiotwin4d/transform_tools.py b/src/physiotwin4d/transform_tools.py index 49e3aeb4..e3244172 100644 --- a/src/physiotwin4d/transform_tools.py +++ b/src/physiotwin4d/transform_tools.py @@ -12,7 +12,7 @@ """ import logging -from typing import Type, Union, cast +from typing import Optional, Type, Union, cast import itk import numpy as np @@ -641,32 +641,147 @@ def smooth_transform( return tfm_smooth def smooth_deformation_field_transform( - self, field: itk.Image, sigma: float + self, + field: itk.Image, + sigma: float, + weight_image: Optional[itk.Image] = None, + normal_image: Optional[itk.Image] = None, + interior_mask: Optional[itk.Image] = None, ) -> itk.DisplacementFieldTransform: - """Wrap a deformation field as a Gaussian-smoothed field transform. - - The float vector ``field`` is converted to a double-precision vector - field, wrapped as a :class:`itk.DisplacementFieldTransform` and - Gaussian-smoothed by ``sigma`` (physical millimeters). Smoothing spreads - a thin surface-shell field into a continuous deformation (and attenuates - its peak magnitude). + """Spread a sparsely sampled deformation field into a continuous one. + + ``field`` is treated as a weighted set of displacement *samples* rather + than as an image: the weighted samples and their weights are each + Gaussian-smoothed by ``sigma`` (physical millimeters) and then divided, + which is a Gaussian-weighted average of the nearby samples. A thin + surface shell therefore becomes a continuous deformation that keeps the + displacement magnitude the samples carried, instead of being diluted by + the empty voxels a plain blur would average in. Far from every sample + the smoothed weight vanishes and the field decays to zero. + + That spread is otherwise isotropic, and carries the whole displacement + vector outward. Giving ``normal_image`` and ``interior_mask`` splits each + sample into the component along the surface normal, which expansion and + contraction live in, and the tangential remainder, which sliding lives + in, and spreads only the normal component outside the mask. Tissue + beyond an organ is then pushed and pulled by it without being dragged + along it, which is how a slip interface such as the pleura or the + pericardium behaves. Inside the mask the full vector is spread, so the + organ's own contents still follow its surface. Args: - field (itk.Image): Input vector deformation field. + field (itk.Image): Input vector deformation field, sampled where + ``weight_image`` is non-zero. sigma (float): Standard deviation of the Gaussian smoothing kernel in physical units (millimeters). + weight_image (Optional[itk.Image]): Per-voxel sample weight, such as + the vertex count + :meth:`WorkflowInferMovement.create_deformation_field` returns. + Omit to weight every voxel holding a non-zero displacement + equally, which cannot tell an empty voxel from a genuinely + zero-displacement one. + normal_image (Optional[itk.Image]): Per-voxel unit surface normal on + ``field``'s grid, as + :meth:`WorkflowInferMovement.create_deformation_field` returns + alongside the field. Samples whose normal is zero are spread + whole, having no direction to project onto. + interior_mask (Optional[itk.Image]): Scalar image on ``field``'s + grid, 1 where the full displacement should be spread and 0 where + only its normal component should be. Soften its edge to set the + width of the band the tangential motion dies out over; a binary + mask makes the boundary a discontinuity. Returns: itk.DisplacementFieldTransform: Smoothed field transform. + + Raises: + ValueError: If only one of ``normal_image`` and ``interior_mask`` is + given, if either does not lie on ``field``'s grid, or if the + field holds no non-zero samples to spread. """ - field_double = ImageTools().convert_array_to_image_of_vectors( - itk.array_from_image(field), reference_image=field, ptype=itk.D + if (normal_image is None) != (interior_mask is None): + raise ValueError( + "normal_image and interior_mask must be given together: the " + "normals say what to project onto, the mask says where to." + ) + + field_arr = itk.array_from_image(field).astype(np.float64) + if weight_image is not None: + weights = itk.array_from_image(weight_image).astype(np.float64) + else: + weights = (np.linalg.norm(field_arr, axis=3) > 0.0).astype(np.float64) + + # Outside the mask only the normal component of each sample is spread. + # Both component sets share one denominator, so the weights are smoothed + # once however many fields are being spread through them. + sample_sets = [field_arr] + mask: Optional[np.ndarray] = None + if normal_image is not None and interior_mask is not None: + normals = itk.array_from_image(normal_image).astype(np.float64) + mask = itk.array_from_image(interior_mask).astype(np.float64) + if normals.shape != field_arr.shape or mask.shape != field_arr.shape[:3]: + raise ValueError( + f"normal_image {normals.shape} and interior_mask " + f"{mask.shape} must lie on the field's grid " + f"{field_arr.shape}." + ) + projected = (field_arr * normals).sum(axis=3, keepdims=True) * normals + # A vertex interior to a volumetric template carries a zero normal. + # Projecting it would delete a sample the weights still count in the + # denominator, biasing the result toward zero rather than leaving the + # sample unprojected, so those keep their full displacement. + unoriented = np.linalg.norm(normals, axis=3) == 0.0 + projected[unoriented] = field_arr[unoriented] + sample_sets.append(projected) + + smoothed_sets = [np.zeros_like(field_arr) for _ in sample_sets] + for samples, into in zip(sample_sets, smoothed_sets): + for dim in range(field_arr.shape[3]): + into[:, :, :, dim] = self._smooth_scalar_array( + samples[:, :, :, dim] * weights, sigma, field + ) + smoothed = smoothed_sets[0] + smoothed_weights = self._smooth_scalar_array(weights, sigma, field) + + # Add a floor to the denominator rather than clamping to it. ITK's + # recursive Gaussian is an IIR approximation, so far from every sample + # both smoothed arrays ring around zero; clamping a denominator that + # small turns that ringing into displacements several times larger than + # any the samples carried, while adding to it lets the quotient fall off + # to zero there, which is what a field with no nearby sample should do. + weight_floor = 1.0e-3 * float(smoothed_weights.max()) + if weight_floor <= 0.0: + raise ValueError("Deformation field has no non-zero samples to spread.") + denominator = (np.maximum(smoothed_weights, 0.0) + weight_floor)[..., None] + for spread in smoothed_sets: + spread /= denominator + + if mask is not None: + inside = np.clip(mask, 0.0, 1.0)[..., None] + smoothed = inside * smoothed_sets[0] + (1.0 - inside) * smoothed_sets[1] + + smoothed_field = ImageTools().convert_array_to_image_of_vectors( + smoothed, reference_image=field, ptype=itk.D ) field_transform = itk.DisplacementFieldTransform[itk.D, 3].New() - field_transform.SetDisplacementField(field_double) - return self.smooth_transform( - field_transform, sigma=sigma, reference_image=field + field_transform.SetDisplacementField(smoothed_field) + return field_transform + + @staticmethod + def _smooth_scalar_array( + array: np.ndarray, sigma: float, reference_image: itk.Image + ) -> np.ndarray: + """Gaussian-smooth a scalar array on ``reference_image``'s grid. + + The array is wrapped with the reference geometry before filtering, so + ``sigma`` is in millimeters rather than in voxels. + """ + image = itk.image_from_array(np.ascontiguousarray(array)) + image.CopyInformation(reference_image) + smoothed: np.ndarray = itk.array_from_image( + itk.smoothing_recursive_gaussian_image_filter(image, Sigma=sigma) ) + return smoothed def combine_transforms_with_masks( self, diff --git a/src/physiotwin4d/workflow_evaluate_movement.py b/src/physiotwin4d/workflow_evaluate_movement.py new file mode 100644 index 00000000..911f0c83 --- /dev/null +++ b/src/physiotwin4d/workflow_evaluate_movement.py @@ -0,0 +1,550 @@ +"""Accuracy of an inferred moving anatomy, per anatomical structure. + +:class:`WorkflowEvaluateMovement` scores a +:class:`physiotwin4d.WorkflowInferMovement` against geometry extracted from a +gated image sequence. For every gated time point it carries the reference +frame's labelmap into that time point with the network's own deformation and +compares the result, structure by structure, to the labelmap of the frame that +was actually acquired: volume difference, Dice, and surface RMSE per lung lobe +or per heart chamber. + +Going through labelmaps rather than through the model's own surface is what lets +one workflow serve both anatomies. The lung shape model carries its five lobes +as per-cell labels, but the heart model is a single structure -- the whole heart +minus its chamber cavities -- so its chambers exist only in the acquired +labelmaps. Warping those labelmaps scores every structure the acquisition +contains, whether or not the shape model represents it separately. + +Everything is measured on one isotropic evaluation grid built around the +reference anatomy, so a case whose gated frames carry different slice pitches is +still scored on a single, stated voxel volume. +""" + +from __future__ import annotations + +import csv +import json +import logging +from datetime import datetime, timezone +from pathlib import Path +from typing import Any, cast + +import itk +import numpy as np +import pyvista as pv + +from . import physicsnemo_tools as pnt +from .contour_tools import ContourTools +from .physiotwin4d_base import PhysioTwin4DBase +from .workflow_infer_movement import WorkflowInferMovement + + +class WorkflowEvaluateMovement(PhysioTwin4DBase): + """Score inferred motion per anatomical structure against acquired frames. + + Args: + movement_workflow: The displacement decoder whose predictions are + scored. + label_names: Structures to score, ``{label_id: name}``. Ids the + reference frame does not contain are dropped with a warning; ids a + single acquired frame does not contain are skipped for that frame + alone, since a structure can leave the field of view. + log_level: Logging level. Default: ``logging.INFO``. + """ + + # Volume-plot series colors, assigned in this order and never cycled: eight + # hues whose neighbors stay apart under the common color-vision deficiencies. + _SERIES_COLORS = ( + "#2a78d6", + "#eb6834", + "#1baf7a", + "#eda100", + "#e87ba4", + "#008300", + "#4a3aa7", + "#e34948", + ) + + def __init__( + self, + movement_workflow: WorkflowInferMovement, + label_names: dict[int, str], + log_level: int | str = logging.INFO, + ) -> None: + super().__init__(class_name=self.__class__.__name__, log_level=log_level) + self.movement_workflow = movement_workflow + self.label_names = dict(label_names) + self.contour_tools = ContourTools(log_level=log_level) + + # ─────────────────────────── Public API ──────────────────────────────── + def process( + self, + case_id: str, + shape_parameters: Path, + reference_mesh: Path, + reference_labelmap: itk.Image, + ground_truth_labelmaps: dict[float, itk.Image], + output_directory: Path, + smoothing_sigma_mm: float = 10.0, + evaluation_spacing_mm: float = 1.0, + include_dice: bool = True, + ) -> dict[str, Any]: + """Score every gated time point of one case. + + Args: + case_id: Name of the case being scored, recorded in every output. + shape_parameters: JSON file with the case's PCA coefficient vector. + reference_mesh: The case's fitted reference-frame SSM surface. The + predicted displacements are added to its points, and its extent + defines the evaluation grid. + reference_labelmap: Labelmap of the reference frame, the anatomy + carried into every other time point. + ground_truth_labelmaps: Acquired labelmap per stage, keyed by the + normalized stage in ``[0, 1]``. + output_directory: Directory the report, the CSV and the per-stage + geometry are written to. + smoothing_sigma_mm: Gaussian sigma, in millimeters, that turns the + network's surface-shell deformation into a continuous field. + evaluation_spacing_mm: Isotropic pitch every metric is measured on. + It sets both the voxel volume the Dice and volume figures are + quantized to and the resolution of the deformation fields, whose + memory grows with its cube. + include_dice: Report the Dice overlap. Turn it off for a structure + whose motion is small against its own size: Dice is an overlap + fraction, so a lung lobe scores over 0.96 undeformed and the + column says more about the organ's bulk than about the motion. + The volume and surface figures still resolve it. + + Returns: + Dict with ``rows`` (every metric row), ``csv_file``, + ``report_file``, ``volume_plot_file``, ``predicted_surfaces`` and + ``warped_labelmaps``. + + Raises: + ValueError: If ``ground_truth_labelmaps`` is empty, or none of the + requested labels are in the reference frame. + """ + if not ground_truth_labelmaps: + raise ValueError("No ground-truth labelmaps to evaluate against.") + + out_dir = Path(output_directory) + out_dir.mkdir(parents=True, exist_ok=True) + stages = sorted(ground_truth_labelmaps) + self.log_section("EVALUATE MOVEMENT [%s]: %d stages", case_id, len(stages)) + + grid = self.contour_tools.create_reference_image( + mesh=cast(pv.DataSet, pv.read(str(reference_mesh))), + spatial_resolution=evaluation_spacing_mm, + buffer_factor=0.25, + ptype=itk.template(reference_labelmap)[1][0], + ) + self.log_info( + "Evaluation grid: %s voxels at %.2f mm", + list(itk.size(grid)), + evaluation_spacing_mm, + ) + reference_on_grid = self._resample_labelmap(reference_labelmap, grid) + scored_labels = self._labels_present(reference_on_grid) + provenance = self._provenance(case_id, shape_parameters) + + # One deformation per stage, from the network's own predictions. Each + # stage's warped image is the reference labelmap carried into that + # stage, which is exactly what the metrics below compare. + series = self.movement_workflow.process_time_series( + shape_parameters=shape_parameters, + stages=stages, + output_directory=out_dir, + reference_mesh=reference_mesh, + reference_image=reference_on_grid, + warp_interpolation="nearest", + warp_background_value=0.0, + smoothing_sigma_mm=smoothing_sigma_mm, + ) + + rows: list[dict[str, Any]] = [] + for index, stage in enumerate(stages): + truth = self._resample_labelmap(ground_truth_labelmaps[stage], grid) + truth_surfaces = self._label_surfaces(truth) + predicted = itk.imread(str(series["warped_images"][index])) + rows.extend( + self._score( + case_id, + stage, + truth, + truth_surfaces, + predicted, + self._label_surfaces(predicted), + scored_labels, + provenance, + include_dice, + ) + ) + + csv_file = self._write_csv(rows, out_dir) + plot_file = self._write_volume_plot(rows, out_dir) + report_file = self._write_report( + rows, + provenance, + stages, + smoothing_sigma_mm, + evaluation_spacing_mm, + plot_file, + out_dir, + ) + return { + "rows": rows, + "csv_file": csv_file, + "report_file": report_file, + "volume_plot_file": plot_file, + "predicted_surfaces": series["predicted_surfaces"], + "warped_labelmaps": series["warped_images"], + } + + # ──────────────────────────── Metrics ────────────────────────────────── + @staticmethod + def dice(truth: np.ndarray, predicted: np.ndarray, label: int) -> float: + """Dice overlap of one label. ``nan`` when neither volume contains it.""" + truth_mask = truth == label + predicted_mask = predicted == label + denominator = np.count_nonzero(truth_mask) + np.count_nonzero(predicted_mask) + if denominator == 0: + return float("nan") + return float(2.0 * np.count_nonzero(truth_mask & predicted_mask) / denominator) + + @staticmethod + def volume_mm3(labels: np.ndarray, label: int, voxel_volume_mm3: float) -> float: + """Volume of one label, in cubic millimeters.""" + return float(np.count_nonzero(labels == label) * voxel_volume_mm3) + + @staticmethod + def surface_rmse_mm(truth: pv.PolyData, predicted: pv.PolyData) -> float: + """Symmetric point-to-surface RMSE, in millimeters. + + Both directions are pooled before the root-mean-square. A one-sided RMSE + misses a prediction that covers the truth everywhere but also bulges + somewhere the truth does not reach. + """ + forward = predicted.copy().compute_implicit_distance(truth) + reverse = truth.copy().compute_implicit_distance(predicted) + distances = np.concatenate( + [ + np.asarray(forward["implicit_distance"], dtype=np.float64), + np.asarray(reverse["implicit_distance"], dtype=np.float64), + ] + ) + return float(np.sqrt(np.mean(distances**2))) + + # ──────────────────────────── Internals ──────────────────────────────── + @staticmethod + def _resample_labelmap(labelmap: itk.Image, grid: itk.Image) -> itk.Image: + """Resample a labelmap onto ``grid``, preserving its discrete values.""" + return itk.resample_image_filter( + labelmap, + use_reference_image=True, + reference_image=grid, + interpolator=itk.NearestNeighborInterpolateImageFunction.New(labelmap), + default_pixel_value=0, + ) + + def _labels_present(self, reference_labelmap: itk.Image) -> dict[int, str]: + """Drop the requested labels the reference frame does not contain.""" + present = set(np.unique(itk.GetArrayViewFromImage(reference_labelmap)).tolist()) + scored = { + label: name for label, name in self.label_names.items() if label in present + } + missing = sorted(set(self.label_names) - set(scored)) + if missing: + self.log_warning( + "Reference frame has no voxels for label(s) %s; not scored.", missing + ) + if not scored: + raise ValueError( + "None of the requested labels are present in the reference frame." + ) + self.log_info( + "Scoring %d structure(s): %s", + len(scored), + ", ".join(scored[label] for label in sorted(scored)), + ) + return scored + + def _label_surfaces(self, labelmap: itk.Image) -> dict[int, pv.PolyData]: + """Contour every label of one labelmap on the evaluation grid's pitch.""" + return self.contour_tools.extract_label_surfaces(labelmap) + + def _provenance(self, case_id: str, shape_parameters: Path) -> dict[str, Any]: + """Case name, shape parameters and network weights, with their dates.""" + inference = self.movement_workflow.inference_workflow + checkpoint = Path(inference.checkpoint_file) + info = checkpoint.stat() + coefficients = pnt.load_pca_coefficients(shape_parameters) + provenance: dict[str, Any] = { + "case_id": case_id, + "shape_parameters_file": str(shape_parameters), + "network_weights_file": str(checkpoint), + "network_weights_created": self._timestamp(info.st_ctime), + "network_weights_modified": self._timestamp(info.st_mtime), + "network_epoch": "final" if inference.epoch is None else inference.epoch, + } + for index, coefficient in enumerate(coefficients, start=1): + provenance[f"pca_c{index:02d}"] = float(coefficient) + return provenance + + @staticmethod + def _timestamp(seconds: float) -> str: + """Format a filesystem timestamp as an ISO-8601 UTC string.""" + return datetime.fromtimestamp(seconds, tz=timezone.utc).isoformat( + timespec="seconds" + ) + + def _score( + self, + case_id: str, + stage: float, + truth: itk.Image, + truth_surfaces: dict[int, pv.PolyData], + predicted: itk.Image, + predicted_surfaces: dict[int, pv.PolyData], + scored_labels: dict[int, str], + provenance: dict[str, Any], + include_dice: bool = True, + ) -> list[dict[str, Any]]: + """One metric row per scored label of one stage.""" + truth_array = itk.GetArrayViewFromImage(truth) + predicted_array = itk.GetArrayViewFromImage(predicted) + voxel_volume_mm3 = float(np.prod(np.asarray(truth.GetSpacing()))) + + rows: list[dict[str, Any]] = [] + for label in sorted(scored_labels): + truth_volume = self.volume_mm3(truth_array, label, voxel_volume_mm3) + if truth_volume == 0.0: + self.log_info( + "stage %.3f: %s absent from the acquired frame; skipped.", + stage, + scored_labels[label], + ) + continue + predicted_volume = self.volume_mm3(predicted_array, label, voxel_volume_mm3) + rmse = ( + self.surface_rmse_mm(truth_surfaces[label], predicted_surfaces[label]) + if label in truth_surfaces and label in predicted_surfaces + else float("nan") + ) + row: dict[str, Any] = { + "case_id": case_id, + "stage": stage, + "label_id": label, + "label_name": scored_labels[label], + } + if include_dice: + row["dice"] = self.dice(truth_array, predicted_array, label) + row.update( + { + "volume_truth_mm3": truth_volume, + "volume_predicted_mm3": predicted_volume, + "volume_difference_mm3": predicted_volume - truth_volume, + "volume_difference_percent": ( + 100.0 * (predicted_volume - truth_volume) / truth_volume + ), + "surface_rmse_mm": rmse, + } + ) + row.update( + {key: value for key, value in provenance.items() if key != "case_id"} + ) + rows.append(row) + self.log_info( + "stage %.3f %-24s %sdV=%+.2f%% rmse=%.3f mm", + stage, + scored_labels[label], + f"dice={row['dice']:.4f} " if include_dice else "", + row["volume_difference_percent"], + rmse, + ) + return rows + + @staticmethod + def _write_csv(rows: list[dict[str, Any]], out_dir: Path) -> Path: + """Write every metric row, provenance included, to one CSV.""" + csv_file = out_dir / "evaluation_metrics.csv" + with csv_file.open("w", newline="", encoding="utf-8") as fh: + writer = csv.DictWriter(fh, fieldnames=list(rows[0].keys())) + writer.writeheader() + writer.writerows(rows) + return csv_file + + def _write_volume_plot(self, rows: list[dict[str, Any]], out_dir: Path) -> Path: + """Plot the acquired and predicted volume of every structure against stage. + + One color per structure, taken in a fixed order from a hue set separable + under color-vision deficiency; the acquired volume is solid and the + predicted volume dashed, so the two never rest on color alone. + """ + import matplotlib.pyplot as plt + from matplotlib.lines import Line2D + + plot_file = out_dir / "volume_vs_stage.png" + fig, ax = plt.subplots(figsize=(8.0, 4.5)) + try: + ends: list[tuple[float, float, str]] = [] + for index, label in enumerate( + sorted({int(row["label_id"]) for row in rows}) + ): + matching = sorted( + (row for row in rows if row["label_id"] == label), + key=lambda row: float(row["stage"]), + ) + stages = [float(row["stage"]) for row in matching] + truth = [float(row["volume_truth_mm3"]) / 1000.0 for row in matching] + predicted = [ + float(row["volume_predicted_mm3"]) / 1000.0 for row in matching + ] + color = self._SERIES_COLORS[index] + ax.plot(stages, truth, color=color, linewidth=2.0, marker="o", ms=5) + ax.plot(stages, predicted, color=color, linewidth=2.0, linestyle="--") + ends.append((stages[-1], truth[-1], str(matching[0]["label_name"]))) + + # Three of the hues fall below 3:1 against a white page, so each line + # is named where it ends rather than in a color key alone. Structures + # of similar size end on top of each other, so the names are pushed + # apart, largest first, before they are drawn. + span = float(np.ptp(ax.get_ylim())) + previous = float("inf") + for x_end, y_end, name in sorted(ends, key=lambda end: -end[1]): + text_y = min(y_end, previous - 0.05 * span) + ax.annotate( + name, + xy=(x_end, text_y), + xytext=(6, 0), + textcoords="offset points", + color="#52514e", + fontsize=9, + va="center", + ) + previous = text_y + + ax.set_xlabel("Stage", color="#52514e") + ax.set_ylabel("Volume (mL)", color="#52514e") + ax.grid(True, color="#e1e0d9", linewidth=0.8) + ax.set_axisbelow(True) + ax.spines["top"].set_visible(False) + ax.spines["right"].set_visible(False) + ax.spines["left"].set_color("#c3c2b7") + ax.spines["bottom"].set_color("#c3c2b7") + ax.tick_params(colors="#898781", labelsize=9) + ax.legend( + handles=[ + Line2D([], [], color="#898781", linewidth=2.0, label="acquired"), + Line2D( + [], + [], + color="#898781", + linewidth=2.0, + linestyle="--", + label="predicted", + ), + ], + frameon=False, + loc="best", + fontsize=9, + labelcolor="#52514e", + ) + fig.savefig(str(plot_file), bbox_inches="tight", dpi=150) + finally: + plt.close(fig) + self.log_info("Volume plot: %s", plot_file) + return plot_file + + def _write_report( + self, + rows: list[dict[str, Any]], + provenance: dict[str, Any], + stages: list[float], + smoothing_sigma_mm: float, + evaluation_spacing_mm: float, + plot_file: Path, + out_dir: Path, + ) -> Path: + """Write the markdown report beside the CSV.""" + coefficients = [ + provenance[key] for key in sorted(provenance) if key.startswith("pca_c") + ] + # Both tables carry whichever metrics the rows were scored with. + has_dice = "dice" in rows[0] + lines = [ + f"# Movement accuracy: {provenance['case_id']}", + "", + "## Run", + "", + f"- Hold-out case: `{provenance['case_id']}`", + f"- Stages evaluated: {len(stages)} " + f"({', '.join(f'{stage:.2f}' for stage in stages)})", + f"- Shape parameters: `{provenance['shape_parameters_file']}`", + "- Shape parameters (standard deviations): " + + json.dumps([round(value, 4) for value in coefficients]), + f"- Network weights: `{provenance['network_weights_file']}`", + f"- Network weights created: {provenance['network_weights_created']}", + f"- Network weights modified: {provenance['network_weights_modified']}", + f"- Network epoch: {provenance['network_epoch']}", + f"- Deformation smoothing sigma: {smoothing_sigma_mm:.1f} mm", + f"- Evaluation grid pitch: {evaluation_spacing_mm:.2f} mm isotropic", + "", + "Every score compares the reference frame carried into that stage " + "with the inferred deformation against the frame acquired there.", + "", + "## Volume over the stages", + "", + f"![Structure volume against stage]({plot_file.name})", + "", + "Solid: the volume acquired at that stage. Dashed: the volume the " + "prediction carries there.", + "", + "## Per structure, averaged over stages", + "", + ] + metrics = (["Dice"] if has_dice else []) + [ + "Volume difference (%)", + "Surface RMSE (mm)", + ] + lines += self._table_header(["Structure"], metrics) + for label in sorted({int(row["label_id"]) for row in rows}): + matching = [row for row in rows if row["label_id"] == label] + cells = [str(matching[0]["label_name"])] + if has_dice: + cells.append(f"{self._mean(matching, 'dice'):.4f}") + cells += [ + f"{self._mean(matching, 'volume_difference_percent'):+.2f}", + f"{self._mean(matching, 'surface_rmse_mm'):.3f}", + ] + lines.append("| " + " | ".join(cells) + " |") + + lines += ["", "## Per stage", ""] + lines += self._table_header(["Stage", "Structure"], metrics) + for row in rows: + cells = [f"{row['stage']:.2f}", str(row["label_name"])] + if has_dice: + cells.append(f"{row['dice']:.4f}") + cells += [ + f"{row['volume_difference_percent']:+.2f}", + f"{row['surface_rmse_mm']:.3f}", + ] + lines.append("| " + " | ".join(cells) + " |") + + report_file = out_dir / "evaluation_report.md" + report_file.write_text("\n".join(lines) + "\n", encoding="utf-8") + self.log_info("Report: %s", report_file) + return report_file + + @staticmethod + def _table_header(keys: list[str], metrics: list[str]) -> list[str]: + """Markdown header and alignment rows: keys left, metrics right.""" + return [ + "| " + " | ".join(keys + metrics) + " |", + "| " + " | ".join(["---"] * len(keys) + ["---:"] * len(metrics)) + " |", + ] + + @staticmethod + def _mean(rows: list[dict[str, Any]], key: str) -> float: + """Mean of one column, ignoring the rows where it could not be measured.""" + values = [float(row[key]) for row in rows if not np.isnan(float(row[key]))] + return float(np.mean(values)) if values else float("nan") diff --git a/src/physiotwin4d/workflow_infer_movement.py b/src/physiotwin4d/workflow_infer_movement.py index 9bbabbf7..98723639 100644 --- a/src/physiotwin4d/workflow_infer_movement.py +++ b/src/physiotwin4d/workflow_infer_movement.py @@ -15,8 +15,9 @@ import csv import logging +from collections.abc import Sequence from pathlib import Path -from typing import Any, Optional, cast +from typing import Any, Literal, Optional, cast import itk import numpy as np @@ -24,6 +25,8 @@ from . import physicsnemo_tools as pnt from .physiotwin4d_base import PhysioTwin4DBase +from .transform_tools import TransformTools +from .workflow_convert_vtk_to_usd import WorkflowConvertVTKToUSD from .workflow_infer_physicsnemo import WorkflowInferPhysicsNeMo @@ -249,6 +252,178 @@ def predict_single( ) return result + def process_time_series( + self, + shape_parameters: Path, + stages: Sequence[float], + output_directory: Path, + reference_mesh: Optional[Path] = None, + ground_truth: Optional[Sequence[Path]] = None, + reference_image: Optional[itk.Image] = None, + warp_interpolation: str = "linear", + warp_background_value: float = 0.0, + smoothing_sigma_mm: float = 10.0, + usd_project_name: Optional[str] = None, + anatomy_type: Optional[str] = None, + separate_by_connectivity: bool = False, + ) -> dict[str, Any]: + """Predict one subject across a whole time series and write its geometry. + + One prediction per entry of ``stages``, each written as a mesh. When + ``reference_image`` is supplied, each stage also gets a deformation + field, which is Gaussian-smoothed into a continuous + :class:`itk.DisplacementFieldTransform` and used to carry + ``reference_image`` into that stage's frame. The smoothing spreads a + surface-shell field into the volume, so the warped image is an + interpolation of the surface motion, not an independent registration. + + Args: + shape_parameters: JSON file with the subject PCA coefficient vector. + stages: Stages to predict, in the order they are to be animated. + output_directory: Directory every artifact is written to. + reference_mesh: The subject's reference mesh; omit to displace the + PCA reconstruction instead. + ground_truth: One mesh per stage whose points are the true stage + positions, for error reporting. Must align with ``stages``. + reference_image: Image carried through each stage's deformation, and + the grid the deformation field is rasterized on. Omit to write + meshes only. + warp_interpolation: Interpolation used to resample + ``reference_image``: ``"linear"`` for intensity images, + ``"nearest"`` for labelmaps and masks. + warp_background_value: Value written where a stage's grid samples + outside ``reference_image``. ``0.0`` suits labelmaps; CT needs + ``-1000.0``, which is air in Hounsfield units. + smoothing_sigma_mm: Gaussian sigma, in millimeters, that turns the + sparse surface-shell field into a continuous deformation. + usd_project_name: When given, the stage meshes are also written as + one animated USD under this name, one time sample per stage. + anatomy_type: Anatomy whose materials color that USD. + separate_by_connectivity: Whether that USD splits each frame into + separate objects by connectivity. + + Returns: + Dict with ``stages``, ``predicted_surfaces``, ``warped_images``, + ``transforms``, ``usd_file``, ``statistics`` and + ``statistics_file``. Entries that were not requested are empty + lists or ``None``. + + Raises: + ValueError: If ``stages`` is empty, or ``ground_truth`` is given + with a different length. + """ + if not stages: + raise ValueError("process_time_series needs at least one stage.") + if ground_truth is not None and len(ground_truth) != len(stages): + raise ValueError( + f"ground_truth has {len(ground_truth)} entries but there are " + f"{len(stages)} stages." + ) + + workflow = self.inference_workflow + coeffs = pnt.load_pca_coefficients(shape_parameters) + ref_mesh = ( + cast(pv.DataSet, pv.read(str(reference_mesh))) + if reference_mesh is not None + else None + ) + ref_points = self._reference_points(coeffs, ref_mesh) + template = ref_mesh if ref_mesh is not None else workflow.template_mesh + + out_dir = Path(output_directory) + out_dir.mkdir(parents=True, exist_ok=True) + stem = Path(shape_parameters).stem + suffix = ".vtp" if isinstance(template, pv.PolyData) else ".vtu" + self.log_section("INFER MOVEMENT TIME SERIES [%s]", stem) + + transform_tools = TransformTools(log_level=self.log_level) + stage_meshes: list[pv.DataSet] = [] + surfaces: list[Path] = [] + warped_images: list[Path] = [] + transforms: list[itk.Transform] = [] + stats: list[dict] = [] + + for index, stage in enumerate(stages): + tag = f"s{int(stage * 100):03d}" + pred_points = ref_points + workflow.predict(coeffs, stage) + pred_mesh = template.copy(deep=True) + pred_mesh.points = pred_points + surface_file = out_dir / f"{stem}_{tag}_pred{suffix}" + pred_mesh.save(str(surface_file)) + stage_meshes.append(pred_mesh) + surfaces.append(surface_file) + + if reference_image is not None: + field = self.create_deformation_field( + shape_parameters=shape_parameters, + stage=stage, + reference_image=reference_image, + reference_mesh=reference_mesh, + direction="inverse", + ) + transform = transform_tools.smooth_deformation_field_transform( + field["deformation_field"], + sigma=smoothing_sigma_mm, + weight_image=field["weight_image"], + ) + transforms.append(transform) + warped = transform_tools.transform_image( + reference_image, + transform, + reference_image=reference_image, + interpolation_method=warp_interpolation, + background_value=warp_background_value, + ) + warped_file = out_dir / f"{stem}_{tag}_warped.mha" + itk.imwrite(warped, str(warped_file), compression=True) + warped_images.append(warped_file) + + if ground_truth is not None: + actual = np.asarray( + pv.read(str(ground_truth[index])).points, dtype=np.float32 + ) + stats.append(self._error_row(stem, stage, pred_points, actual)) + self.log_info( + "stage %.3f: mean=%.3f mm max=%.3f mm", + stage, + stats[-1]["mean_error_mm"], + stats[-1]["max_error_mm"], + ) + else: + self.log_info("stage %.3f -> %s", stage, surface_file.name) + + statistics_file: Optional[Path] = None + if stats: + statistics_file = out_dir / "statistics_per_stage.csv" + with statistics_file.open("w", newline="", encoding="utf-8") as fh: + writer = csv.DictWriter(fh, fieldnames=list(stats[0].keys())) + writer.writeheader() + writer.writerows(stats) + + usd_file: Optional[Path] = None + if usd_project_name is not None: + usd_workflow = WorkflowConvertVTKToUSD( + input_meshes=stage_meshes, + usd_project_name=usd_project_name, + output_directory=out_dir, + appearance="anatomy" if anatomy_type is not None else "solid", + anatomy_type=anatomy_type, + separate_by_connectivity=separate_by_connectivity, + frames_per_second=float(len(stage_meshes)), + log_level=self.log_level, + ) + usd_file = Path(usd_workflow.process()["usd_file"]) + + return { + "stages": list(stages), + "predicted_surfaces": surfaces, + "warped_images": warped_images, + "transforms": transforms, + "usd_file": usd_file, + "statistics": stats, + "statistics_file": statistics_file, + } + def create_deformation_field( self, shape_parameters: Path, @@ -256,6 +431,7 @@ def create_deformation_field( reference_image: itk.Image, output_directory: Optional[Path] = None, reference_mesh: Optional[Path] = None, + direction: Literal["forward", "inverse"] = "forward", ) -> dict[str, Any]: """Rasterize the inferred deformation onto a reference image grid. @@ -266,6 +442,15 @@ def create_deformation_field( (renormalized) reference-surface normal of those vertices. Empty voxels are zero. + That is the ``"forward"`` field, which maps reference positions to stage + positions and is what transforming a *mesh* needs. Resampling an + *image*, though, maps each output point through the transform to find + where to sample the input, so carrying the reference image into the + stage frame needs the opposite mapping. ``direction="inverse"`` builds + it exactly rather than by negating the forward field: each vertex is + binned by its **deformed** position ``reference + displacement`` and + contributes ``-displacement``. + The binning positions come from ``reference_mesh``, so a patient scan whose statistical-model fit applied a pose transform not captured by the shape coefficients is binned where it actually aligns with @@ -278,16 +463,23 @@ def create_deformation_field( stage: Target stage for the deformation. reference_image: The frame's image; defines the output grid geometry (size, spacing, origin, direction). - output_directory: If given, the two images are written there as + output_directory: If given, the three images are written there as compressed ``.mha`` files. reference_mesh: Mesh whose points supply the binning positions and normals; omit to use the PCA reconstruction. Must share the template topology (same point count and ordering). + direction: ``"forward"`` for the reference-to-stage field that + deforms meshes, ``"inverse"`` for the stage-to-reference field + that resamples images into the stage frame. Returns: Dict with ``deformation_field`` and ``normal_image`` (ITK vector - images), ``deformed_surface`` (the stage mesh as ``pv.DataSet``) - and, when written, their paths. + images), ``weight_image`` (the vertex count per voxel, which + distinguishes an empty voxel from one whose displacement happens to + be zero and is what + :meth:`TransformTools.smooth_deformation_field_transform` normalizes + by), ``deformed_surface`` (the stage mesh as ``pv.DataSet``) and, + when written, their paths. """ workflow = self.inference_workflow template = workflow.template_mesh @@ -325,12 +517,19 @@ def create_deformation_field( normal_sum = np.zeros((sz, sy, sx, 3), dtype=np.float64) count = np.zeros((sz, sy, sx), dtype=np.float64) - for i in range(ref_points.shape[0]): - point = [float(c) for c in ref_points[i]] + if direction == "inverse": + bin_points = ref_points + disps + bin_disps = -disps + else: + bin_points = ref_points + bin_disps = disps + + for i in range(bin_points.shape[0]): + point = [float(c) for c in bin_points[i]] idx = reference_image.TransformPhysicalPointToIndex(point) ix, iy, iz = int(idx[0]), int(idx[1]), int(idx[2]) if 0 <= ix < sx and 0 <= iy < sy and 0 <= iz < sz: - disp_sum[iz, iy, ix] += disps[i] + disp_sum[iz, iy, ix] += bin_disps[i] normal_sum[iz, iy, ix] += normals[i] count[iz, iy, ix] += 1.0 @@ -347,11 +546,14 @@ def create_deformation_field( deformation_image = self._vector_image_like(disp_field, reference_image) normal_image = self._vector_image_like(normal_field, reference_image) + weight_image = self._scalar_image_like( + count.astype(np.float32), reference_image + ) self.log_info( "Deformation field: %d/%d voxels populated by %d vertices", int(occupied.sum()), sx * sy * sz, - ref_points.shape[0], + bin_points.shape[0], ) # Deformed (stage) mesh: reference positions displaced by the network, @@ -362,6 +564,7 @@ def create_deformation_field( result: dict[str, Any] = { "deformation_field": deformation_image, "normal_image": normal_image, + "weight_image": weight_image, "deformed_surface": deformed_surface, } if output_directory is not None: @@ -370,12 +573,15 @@ def create_deformation_field( suffix = ".vtp" if isinstance(template, pv.PolyData) else ".vtu" field_path = out_dir / "deformation_field.mha" normal_path = out_dir / "surface_normal_field.mha" + weight_path = out_dir / "deformation_weight.mha" surface_path = out_dir / f"deformed_surface{suffix}" itk.imwrite(deformation_image, str(field_path), compression=True) itk.imwrite(normal_image, str(normal_path), compression=True) + itk.imwrite(weight_image, str(weight_path), compression=True) deformed_surface.save(str(surface_path)) result["deformation_field_file"] = field_path result["normal_image_file"] = normal_path + result["weight_image_file"] = weight_path result["deformed_surface_file"] = surface_path return result @@ -388,6 +594,15 @@ def _vector_image_like(array: np.ndarray, reference_image: itk.Image) -> itk.Ima image.SetDirection(reference_image.GetDirection()) return image + @staticmethod + def _scalar_image_like(array: np.ndarray, reference_image: itk.Image) -> itk.Image: + """Wrap a ``(z, y, x)`` array as an ITK scalar image on ``reference``'s grid.""" + image = itk.image_from_array(np.ascontiguousarray(array)) + image.SetSpacing(reference_image.GetSpacing()) + image.SetOrigin(reference_image.GetOrigin()) + image.SetDirection(reference_image.GetDirection()) + return image + @staticmethod def _error_row( subject_id: str, stage: float, pred: np.ndarray, actual: np.ndarray diff --git a/src/physiotwin4d/workflow_infer_physicsnemo.py b/src/physiotwin4d/workflow_infer_physicsnemo.py index c153f5c5..737a27e0 100644 --- a/src/physiotwin4d/workflow_infer_physicsnemo.py +++ b/src/physiotwin4d/workflow_infer_physicsnemo.py @@ -95,6 +95,8 @@ def __init__( checkpoint_file = self.model_directory / f"{tag}_stage_model.pt" if not checkpoint_file.exists(): raise FileNotFoundError(f"Model checkpoint not found: {checkpoint_file}") + self.epoch = epoch + self.checkpoint_file = checkpoint_file self._device = torch.device("cuda" if torch.cuda.is_available() else "cpu") self.log_info("Loading %s model from %s", tag.upper(), checkpoint_file) diff --git a/statistics.md b/statistics.md index a5a85941..177443af 100644 --- a/statistics.md +++ b/statistics.md @@ -1,24 +1,28 @@ # PhysioTwin4D - Software Development Statistics -**Report Generated:** August 5, 2026 +**Report Generated:** August 13, 2026 **Project Version:** 2026.07.3 **Status:** Beta (Development Status: 4 - Beta) +Line counts below are total lines per file (`cat | wc -l`), including blanks and +comments. Stating the method matters: the previous revision of this report used +an unstated one, so its figures are not directly comparable to these. + --- ## Executive Summary PhysioTwin4D is a collection of methods, workflows, tutorials, and CLI tools -for creating personalized physiological digital twins from 3D medical images. +for creating personalized physiological digital twins from 3D/4D medical images. This report summarizes development effort, code quality, and project maturity. ### Key Metrics at a Glance | Metric | Value | | ------------------------------ | ---------------------------------------------- | -| **Total Lines of Code** | ~53,400 | -| **Development Period** | December 5, 2025 - August 5, 2026 (~8 months) | -| **Total Commits** | 114 | +| **Total Lines of Code** | ~74,100 | +| **Development Period** | December 5, 2025 - August 13, 2026 (~8 months) | +| **Total Commits** | 122 | | **Primary Developer** | 1 (Stephen Aylward), plus 1 outside contributor | --- @@ -29,14 +33,19 @@ This report summarizes development effort, code quality, and project maturity. | Category | Files | Lines of Code | Percentage | | ---------------------------------------- | -------------- | -------------- | ---------- | -| **Core Python Source (`src/`)** | 73 files | 22,663 | 42.4% | -| **Test Suite (`tests/`)** | 36 files | 8,225 | 15.4% | -| **Experiment Scripts (`experiments/`)** | 46 files | 7,897 | 14.8% | -| **Tutorial Scripts (`tutorials/`)** | 20 files | 4,153 | 7.8% | -| **Utility Scripts (`utils/`)** | 3 files | 1,460 | 2.7% | -| **Documentation (`docs/*.rst`)** | 85 files | 6,069 | 11.4% | -| **Markdown (repo-wide READMEs, guides)** | 35 files | 2,958 | 5.5% | -| **TOTAL** | **298 files** | **~53,400** | **100%** | +| **Core Python Source (`src/`)** | 74 files | 29,316 | 39.6% | +| **Test Suite (`tests/`)** | 40 files | 11,601 | 15.7% | +| **Experiment Scripts (`experiments/`)** | 46 files | 9,219 | 12.4% | +| **Tutorial Scripts (`tutorials/`)** | 32 files | 9,107 | 12.3% | +| **Utility Scripts (`utils/`)** | 3 files | 1,739 | 2.3% | +| **Documentation (`docs/*.rst`)** | 85 files | 8,871 | 12.0% | +| **Markdown (repo-wide READMEs, guides)** | 38 files | 4,204 | 5.7% | +| **TOTAL** | **318 files** | **~74,100** | **100%** | + +The 32 files under `tutorials/` are 29 numbered tutorial scripts plus 3 +per-organ parameter modules (`parameters_heart_ct_kcl.py`, +`parameters_lung_ct_dirlab.py`, `parameters_duke_heart_labelmaps.py`) that +carry the constants the tutorials share. All experiment and tutorial sources are plain `.py` files run with `python