diff --git a/.github/workflows/README.md b/.github/workflows/README.md
index 23493c2..3c41b65 100644
--- a/.github/workflows/README.md
+++ b/.github/workflows/README.md
@@ -25,7 +25,7 @@ Runs on every push and pull request to main branches. Includes:
     - Manually triggered via workflow_dispatch, OR
     - PR has the `run-gpu-tests` label
   - Requires self-hosted runner with `[self-hosted, linux, gpu]` labels
-  - Uses PyTorch with CUDA 12.6 support
+  - Uses PyTorch with CUDA 13.0 support
   - Timeout: 30 minutes
 
 - **code-quality**: Static code analysis
@@ -87,7 +87,7 @@ To run GPU tests, you must either:
 
 GPU tests require self-hosted runners with:
 - Linux OS
-- NVIDIA GPU with CUDA 12.6+ support
+- NVIDIA GPU with CUDA 13.0 support
 - Runner labels: `[self-hosted, linux, gpu]`
 
 **Why are GPU tests disabled by default?**
@@ -107,13 +107,18 @@ GPU tests require self-hosted runners with:
    # Install NVIDIA drivers
    sudo apt-get install nvidia-driver-535
    
-   # Install CUDA toolkit 12.6
-   wget https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2204/x86_64/cuda-keyring_1.0-1_all.deb
-   sudo dpkg -i cuda-keyring_1.0-1_all.deb
+   # Install CUDA toolkit 13.0
+   # CUDA 13.0 requires Ubuntu 22.04 LTS or later. These cuda-keyring_1.1-1_all.deb
+   # and cuda-toolkit-13 commands fail on Ubuntu 20.04 runners.
+   wget https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2204/x86_64/cuda-keyring_1.1-1_all.deb
+   sudo dpkg -i cuda-keyring_1.1-1_all.deb
    sudo apt-get update
-   sudo apt-get install cuda-toolkit-12-6
+   sudo apt-get install cuda-toolkit-13
    ```
 
+   Self-hosted GPU runners should be upgraded to Ubuntu 22.04 LTS or later
+   before installing CUDA 13.0.
+
 3. **Configure Runner Labels**:
    - Add labels: `self-hosted`, `linux`, `gpu`
    - Verify GPU is accessible: `nvidia-smi`
@@ -140,8 +145,7 @@ GPU tests require self-hosted runners with:
 **Option 3: Run Locally**
 ```bash
 # Install with CUDA support
-pip install torch torchvision torchaudio --index-url https://download.pytorch.org/whl/cu126
-pip install -e ".[test]"
+uv pip install -e ".[test,cuda13]"
 
 # Run GPU tests
 pytest tests/ -v -m "not slow"
@@ -214,8 +218,7 @@ pytest tests/ -m "unit and not requires_gpu" --cov=physiomotion4d
 ### GPU Tests
 ```bash
 # Install with CUDA support
-pip install torch torchvision torchaudio --index-url https://download.pytorch.org/whl/cu126
-pip install -e ".[test]"
+uv pip install -e ".[test,cuda13]"
 
 # Run all tests (including GPU)
 pytest tests/ -m "not slow"
diff --git a/README.md b/README.md
index 09e5e3c..56f6ab9 100644
--- a/README.md
+++ b/README.md
@@ -26,7 +26,7 @@ PhysioMotion4D is a comprehensive medical imaging package that converts 4D CT sc
 ### Prerequisites
 
 - Python 3.10+ (Python 3.10, 3.11, or 3.12 recommended)
-- NVIDIA GPU with CUDA 13 (default) or CUDA 12 — recommended for production use; CPU-only installation is supported but slow
+- NVIDIA GPU with CUDA 13 — recommended for production use; CPU-only installation is supported but slow
 - 16GB+ RAM (32GB+ recommended for large datasets)
 - NVIDIA Omniverse (for USD visualization)
 - **Git LFS** (required for running tests: baseline files in `tests/baselines/` are stored with Git LFS; install from [git-lfs.github.com](https://git-lfs.github.com), then run `git lfs install` and `git lfs pull` after cloning)
@@ -34,20 +34,16 @@ PhysioMotion4D is a comprehensive medical imaging package that converts 4D CT sc
 ### Installation from PyPI
 
 ```bash
-# CPU-only install — works out of the box; a runtime warning points to the GPU extras
+# CPU-only PyPI install — works out of the box; a runtime warning points to the GPU extra
 pip install physiomotion4d
 
 # CUDA 13 install (recommended for production)
 uv pip install "physiomotion4d[cuda13]"
-
-# CUDA 12 install
-uv pip install "physiomotion4d[cuda12]"
 ```
 
-The `[cuda13]` and `[cuda12]` extras install both CuPy and the matching
-CUDA-built PyTorch wheel in one step — there is no need to install PyTorch
-separately. PyTorch is listed in the extras so that uv's dependency resolver
-fetches the GPU wheel from the PyTorch index instead of the CPU wheel from PyPI.
+The `[cuda13]` extra installs CuPy. In uv-managed source environments, PyTorch,
+torchvision, and torchaudio resolve from the CUDA 13.0 PyTorch wheel index.
+There is no need to install PyTorch separately.
 
 For development with NVIDIA NIM cloud services:
 ```bash
@@ -77,14 +73,11 @@ pip install physiomotion4d[nim]
 
 4. **Install PhysioMotion4D**:
    ```bash
-   # CPU-only (evaluation / no GPU)
+   # CUDA 13 PyTorch is the default for uv-managed source environments
    uv pip install -e "."
 
-   # CUDA 13 (recommended for production)
+   # Add CuPy for CUDA 13 GPU acceleration
    uv pip install -e ".[cuda13]"
-
-   # CUDA 12
-   uv pip install -e ".[cuda12]"
    ```
 
 ### Verify Installation
@@ -136,11 +129,40 @@ print(f"PhysioMotion4D version: {physiomotion4d.__version__}")
 ### Key Dependencies
 
 - **Medical Imaging**: ITK, TubeTK, MONAI, nibabel, PyVista
-- **AI/ML**: PyTorch, CuPy (CUDA 13 default; CUDA 12 via `[cuda12]` extra), transformers, MONAI
+- **AI/ML**: PyTorch, CuPy (CUDA 13), transformers, MONAI
 - **Registration**: icon-registration, unigradicon
 - **Visualization**: USD-core, PyVista
 - **Segmentation**: TotalSegmentator
 
+## Getting Started: Tutorials
+
+The `tutorials/` directory contains six end-to-end Python scripts, one for each
+major workflow. They are the recommended starting point for new users.
+
+| # | Script | Workflow | Dataset |
+|---|--------|----------|---------|
+| 1 | `tutorials/tutorial_01_heart_gated_ct_to_usd.py` | Heart-gated CT to animated USD | Slicer-Heart-CT (prepare first) |
+| 2 | `tutorials/tutorial_02_ct_to_vtk.py` | CT to VTK surfaces | Slicer-Heart-CT (prepare first) |
+| 3 | `tutorials/tutorial_03_fit_statistical_model_to_patient.py` | Fit statistical model to patient | KCL-Heart-Model (manual) |
+| 4 | `tutorials/tutorial_04_create_statistical_model.py` | Build PCA shape model | KCL-Heart-Model (manual) |
+| 5 | `tutorials/tutorial_05_vtk_to_usd.py` | VTK surfaces to animated USD | output of tutorial 2 |
+| 6 | `tutorials/tutorial_06_reconstruct_highres_4d_ct.py` | Reconstruct high-res 4D CT | DirLab-4DCT (manual) |
+
+Each script is runnable directly:
+
+```bash
+# Tutorial 1 (CPU-safe ANTs registration; requires Slicer-Heart-CT data)
+python tutorials/tutorial_01_heart_gated_ct_to_usd.py \
+    --data-dir ./data --output-dir ./output/tutorial_01
+
+# Tutorial 2 (CT to VTK)
+python tutorials/tutorial_02_ct_to_vtk.py \
+    --data-dir ./data --output-dir ./output/tutorial_02
+```
+
+See `tutorials/README.md` for the full tutorial index, dataset preparation
+instructions, recommended run order, and experiment-test instructions.
+
 ## 🎯 Quick Start
 
 ### Command-Line Interface
diff --git a/docs/API_MAP.md b/docs/API_MAP.md
index 58bfad0..9c261cb 100644
--- a/docs/API_MAP.md
+++ b/docs/API_MAP.md
@@ -6,8 +6,8 @@ _Re-run `py utils/generate_api_map.py` whenever public APIs change._
 ## docs/conf.py
 
 - **class Mock** (line 20)
-- `def autodoc_skip_member(app, what, name, obj, skip, options)` (line 215): Custom function to skip certain members during autodoc processing.
-- `def setup(app)` (line 223): Custom setup function for Sphinx.
+- `def autodoc_skip_member(app, what, name, obj, skip, options)` (line 216): Custom function to skip certain members during autodoc processing.
+- `def setup(app)` (line 224): Custom setup function for Sphinx.
 
 ## experiments/Colormap-VTK_To_USD/colormap_vtk_to_usd.py
 
@@ -286,18 +286,20 @@ _Re-run `py utils/generate_api_map.py` whenever public APIs change._
 ## src/physiomotion4d/test_tools.py
 
 - `def set_create_baseline_if_missing(value)` (line 28): Set whether to create baseline files when missing (used by pytest conftest).
-- **class TestTools** (line 34): Utilities for pytest image comparison: baseline directory, results directory,
-  - `def __init__(self, results_dir, baselines_dir, class_name, *, log_level=logging.INFO)` (line 44)
-  - `def image_pass_fail_and_pixels_above_tolerance(self)` (line 78): Return (pass, value) for number of pixels above tolerance from the most
-  - `def image_pass_fail_and_total_absolute_error(self)` (line 93): Return (pass, value) for total absolute error from the most recent
-  - `def image_difference(self)` (line 108): Return the difference image (itk.Image) from the most recent
-  - `def transform_pass_fail_and_number_of_values_above_tolerance(self)` (line 115): Return (pass, value) for number of values above tolerance from the most recent compare_result_to_baseline_transform call.
-  - `def transform_pass_fail_and_total_absolute_error(self)` (line 131): Return (pass, value) for total absolute error from the most recent compare_result_to_baseline_transform call.
-  - `def transform_difference(self)` (line 145): Return the difference transform (itk.Transform) from the most recent compare_result_to_baseline_transform call.
-  - `def write_result_image(self, image, filename)` (line 151): Write the image to the results directory.
-  - `def write_result_transform(self, transform, filename)` (line 155): Write the transform to the results directory.
-  - `def compare_result_to_baseline_transform(self, filename, *, per_value_absolute_error_tol=0.0, max_number_of_values_above_tol=0, total_absolute_error_tol=0.0)` (line 161): Compare the transform to the baseline transform.
-  - `def compare_result_to_baseline_image(self, filename, *, per_pixel_absolute_error_tol=0.0, max_number_of_pixels_above_tol=0, total_absolute_error_tol=0.0)` (line 239): Load a 3D result image and a 3D baseline image (.mha), compare the full
+- **class TestTools** (line 34): Utilities for pytest image comparison: baseline directory, result directory,
+  - `def __init__(self, results_dir, baselines_dir, class_name, *, results_output_dir=None, log_level=logging.INFO)` (line 44): Initialize test helpers.
+  - `def image_pass_fail_and_pixels_above_tolerance(self)` (line 95): Return (pass, value) for number of pixels above tolerance from the most
+  - `def image_pass_fail_and_total_absolute_error(self)` (line 110): Return (pass, value) for total absolute error from the most recent
+  - `def image_difference(self)` (line 125): Return the difference image (itk.Image) from the most recent
+  - `def transform_pass_fail_and_number_of_values_above_tolerance(self)` (line 132): Return (pass, value) for number of values above tolerance from the most recent compare_result_to_baseline_transform call.
+  - `def transform_pass_fail_and_total_absolute_error(self)` (line 148): Return (pass, value) for total absolute error from the most recent compare_result_to_baseline_transform call.
+  - `def transform_difference(self)` (line 162): Return the difference transform (itk.Transform) from the most recent compare_result_to_baseline_transform call.
+  - `def write_result_image(self, image, filename)` (line 168): Write the image to the configured result artifact directory.
+  - `def write_result_transform(self, transform, filename)` (line 172): Write the transform to the configured result artifact directory.
+  - `def compare_result_to_baseline_transform(self, filename, *, per_value_absolute_error_tol=0.0, max_number_of_values_above_tol=0, total_absolute_error_tol=0.0)` (line 178): Compare the transform to the baseline transform.
+  - `def compare_result_to_baseline_image(self, filename, *, per_pixel_absolute_error_tol=0.0, max_number_of_pixels_above_tol=0, total_absolute_error_tol=0.0)` (line 256): Load a 3D result image and a 3D baseline image (.mha), compare the full
+  - `def save_screenshot_mesh(self, mesh, filename, *, camera_position='iso', window_size=(800, 600), color='pink', opacity=0.9)` (line 363): Render a PyVista mesh off-screen and save a PNG.
+  - `def save_screenshot_image_slice(self, image, filename, *, axis=0, slice_fraction=0.5, colormap='gray', vmin=None, vmax=None, overlay_mask=None, overlay_alpha=0.4)` (line 412): Extract one slice from an ITK image and save a PNG via matplotlib.
 
 ## src/physiomotion4d/transform_tools.py
 
@@ -378,7 +380,7 @@ _Re-run `py utils/generate_api_map.py` whenever public APIs change._
 - **class UsdMeshConverter** (line 25): Converts MeshData to UsdGeomMesh with full feature support.
   - `def __init__(self, stage, settings, material_mgr)` (line 36): Initialize mesh converter.
   - `def create_mesh(self, mesh_data, mesh_path, time_code=None, bind_material=True)` (line 53): Create a UsdGeomMesh from MeshData.
-  - `def create_time_varying_mesh(self, mesh_data_sequence, mesh_path, time_codes, bind_material=True)` (line 288): Create a mesh with time-varying attributes.
+  - `def create_time_varying_mesh(self, mesh_data_sequence, mesh_path, time_codes, bind_material=True)` (line 289): Create a mesh with time-varying attributes.
 
 ## src/physiomotion4d/vtk_to_usd/usd_utils.py
 
@@ -684,6 +686,27 @@ _Re-run `py utils/generate_api_map.py` whenever public APIs change._
   - `def test_multiple_transform_applications(self, transform_tools, test_transforms, test_images)` (line 628): Test applying multiple transforms in sequence.
   - `def test_identity_transform(self, transform_tools, test_images)` (line 656): Test that identity transform doesn't change the image.
 
+## tests/test_tutorials.py
+
+- `def test_testtools_results_output_dir_override(tmp_path)` (line 78): Store result artifacts in an explicit directory when requested.
+- `def test_tutorial_01_contour_png_mesh_uses_current_run_results()` (line 103): Select current in-memory contours instead of disk VTP outputs.
+- `def test_tutorial_01_reference_png_uses_workflow_fixed_image(tmp_path)` (line 119): Select the actual workflow reference image over cached slice images.
+- `def test_tutorial_01_overlay_uses_workflow_fixed_segmentation(tmp_path)` (line 134): Select the current fixed labelmap over cached slice labelmaps.
+- `def test_tutorial_01_overlay_falls_back_to_fixed_image_mask(tmp_path)` (line 149): Read fixed_image_mask.mha before stale slice labelmap files.
+- **class TestTutorial01HeartGatedCTToUSD** (line 176): End-to-end test for tutorial_01_heart_gated_ct_to_usd.py.
+  - `def test_run(self, test_directories)` (line 181)
+- `def test_tutorial_03_extract_surface_uses_dataset_surface()` (line 207): Use the robust dataset_surface algorithm for VTK surface extraction.
+- **class TestTutorial02CTToVTK** (line 227): End-to-end test for tutorial_02_ct_to_vtk.py.
+  - `def test_run(self, test_directories)` (line 232)
+- **class TestTutorial03FitStatisticalModelToPatient** (line 260): End-to-end test for tutorial_03_fit_statistical_model_to_patient.py.
+  - `def test_run(self, test_directories)` (line 265)
+- **class TestTutorial04CreateStatisticalModel** (line 298): End-to-end test for tutorial_04_create_statistical_model.py.
+  - `def test_run(self, test_directories)` (line 303)
+- **class TestTutorial05VTKToUSD** (line 339): End-to-end test for tutorial_05_vtk_to_usd.py.
+  - `def test_run(self, test_directories)` (line 344)
+- **class TestTutorial06ReconstructHighres4DCT** (line 389): End-to-end test for tutorial_06_reconstruct_highres_4d_ct.py.
+  - `def test_run(self, test_directories)` (line 394)
+
 ## tests/test_usd_merge.py
 
 - `def analyze_usd_file(filepath)` (line 17): Analyze a USD file for materials and time samples.
@@ -738,6 +761,30 @@ _Re-run `py utils/generate_api_map.py` whenever public APIs change._
   - `def test_normals_remain_unit_length(self, tmp_path)` (line 342): Normal vectors must not be scaled.
   - `def test_stage_meters_per_unit(self, tmp_path)` (line 362): Stage metersPerUnit metadata must be 1.0.
 
+## tutorials/tutorial_01_heart_gated_ct_to_usd.py
+
+- `def run_tutorial(data_dir, output_dir, *, registration_method='ants', log_level=logging.INFO)` (line 164): Run Tutorial 1: Heart-Gated CT to Animated USD.
+
+## tutorials/tutorial_02_ct_to_vtk.py
+
+- `def run_tutorial(data_dir, output_dir, *, log_level=logging.INFO)` (line 90): Run Tutorial 2: CT Segmentation to VTK Surfaces.
+
+## tutorials/tutorial_03_fit_statistical_model_to_patient.py
+
+- `def run_tutorial(data_dir, output_dir, *, log_level=logging.INFO)` (line 99): Run Tutorial 3: Fit Statistical Shape Model to Patient Data.
+
+## tutorials/tutorial_04_create_statistical_model.py
+
+- `def run_tutorial(data_dir, output_dir, *, pca_components=10, max_samples=20, log_level=logging.INFO)` (line 93): Run Tutorial 4: Create a PCA Statistical Shape Model.
+
+## tutorials/tutorial_05_vtk_to_usd.py
+
+- `def run_tutorial(data_dir, output_dir, *, vtk_file=None, log_level=logging.INFO)` (line 80): Run Tutorial 5: VTK Surface Series to Animated USD.
+
+## tutorials/tutorial_06_reconstruct_highres_4d_ct.py
+
+- `def run_tutorial(data_dir, output_dir, *, case=1, max_frames=4, registration_method='ants', log_level=logging.INFO)` (line 97): Run Tutorial 6: Reconstruct High-Resolution 4D CT.
+
 ## utils/claude_github_reviews.py
 
 - `def git_fetch(repo_root, remote, branch)` (line 67): Run ``git fetch <remote> <branch>``, printing progress.
diff --git a/docs/_static/custom.css b/docs/_static/custom.css
index b9e4525..b35459f 100644
--- a/docs/_static/custom.css
+++ b/docs/_static/custom.css
@@ -12,6 +12,13 @@
     /* Minimal border radius */
     --radius-sm: 0.25rem;
     --radius-md: 0.375rem;
+    --nvidia-green: #76b900;
+    --pm4d-ink: #111111;
+    --pm4d-charcoal: #1a1a1a;
+    --pm4d-muted: #5f6670;
+    --pm4d-line: #d9dee6;
+    --pm4d-surface: #ffffff;
+    --pm4d-page: #f4f6f8;
 }
 
 /* ==================== GLOBAL STYLES ==================== */
@@ -19,6 +26,8 @@
 /* Improved Typography - minimal override */
 body {
     line-height: 1.6;
+    background: var(--pm4d-page);
+    color: var(--pm4d-ink);
 }
 
 /* Smooth Scrolling */
@@ -36,6 +45,7 @@ h4,
 h5,
 h6 {
     line-height: 1.3;
+    color: var(--pm4d-ink);
 }
 
 /* ==================== LINKS ==================== */
@@ -80,15 +90,47 @@ table {
     border-radius: var(--radius-md);
 }
 
+.wy-side-nav-search {
+    background: var(--pm4d-ink);
+    border-bottom: 3px solid var(--nvidia-green);
+}
+
+.wy-side-nav-search > a,
+.wy-side-nav-search .wy-dropdown > a {
+    color: #ffffff;
+}
+
+.wy-side-nav-search img {
+    width: 148px;
+    border-radius: 0;
+    background: var(--nvidia-green);
+    padding: 8px 12px;
+}
+
+.wy-nav-side {
+    background: var(--pm4d-charcoal);
+}
+
+.wy-menu-vertical header,
+.wy-menu-vertical p.caption {
+    color: var(--nvidia-green);
+}
+
 .wy-menu-vertical a {
     border-radius: var(--radius-sm);
 }
 
+.wy-menu-vertical a:hover {
+    background: #2b2b2b;
+    color: #ffffff;
+}
+
 /* ==================== CONTENT AREA ==================== */
 
 /* Use default content area styling */
 .wy-nav-content {
     max-width: 1200px;
+    background: var(--pm4d-surface);
 }
 
 /* ==================== BUTTONS ==================== */
@@ -193,4 +235,194 @@ dt.sig {
     border-radius: var(--radius-md);
 }
 
-/* ==================== END ==================== */
\ No newline at end of file
+/* ==================== NVIDIA-STYLE TUTORIAL LANDING ==================== */
+
+.pm4d-hero {
+    margin: -1.618em -3.236em 2rem;
+    padding: 4rem 3.236em 4.5rem;
+    color: #ffffff;
+    background: #111111;
+    border-bottom: 6px solid var(--nvidia-green);
+}
+
+.pm4d-hero__brand img {
+    width: 150px;
+    border-radius: 0;
+    margin-bottom: 1.5rem;
+}
+
+.pm4d-kicker {
+    margin: 0 0 0.75rem;
+    color: var(--nvidia-green);
+    font-size: 0.85rem;
+    font-weight: 700;
+    letter-spacing: 0.08em;
+    text-transform: uppercase;
+}
+
+.pm4d-hero h1 {
+    max-width: 780px;
+    margin: 0 0 1rem;
+    color: #ffffff;
+    font-size: 2.55rem;
+    line-height: 1.08;
+}
+
+.pm4d-hero p:not(.pm4d-kicker) {
+    max-width: 760px;
+    margin: 0;
+    color: #d8dde6;
+    font-size: 1.1rem;
+}
+
+.pm4d-card-grid {
+    display: grid;
+    grid-template-columns: repeat(auto-fit, minmax(250px, 1fr));
+    gap: 1.15rem;
+    margin: -4rem 0 2.25rem;
+    position: relative;
+    z-index: 2;
+}
+
+.pm4d-card,
+.pm4d-card:visited {
+    display: flex;
+    min-height: 214px;
+    padding: 1.25rem;
+    color: var(--pm4d-ink);
+    text-decoration: none;
+    background: #ffffff;
+    border: 1px solid var(--pm4d-line);
+    border-top: 4px solid var(--nvidia-green);
+    border-radius: 8px;
+    box-shadow: 0 16px 34px rgba(17, 17, 17, 0.16);
+    flex-direction: column;
+    transition: transform 160ms ease, box-shadow 160ms ease, border-color 160ms ease;
+}
+
+.pm4d-card:hover {
+    color: var(--pm4d-ink);
+    text-decoration: none;
+    border-color: var(--nvidia-green);
+    box-shadow: 0 20px 42px rgba(17, 17, 17, 0.22);
+    transform: translateY(-4px);
+}
+
+.pm4d-card:focus-visible {
+    color: var(--pm4d-ink);
+    text-decoration: none;
+    border-color: var(--nvidia-green);
+    outline: 3px solid #111111;
+    outline-offset: 3px;
+    box-shadow: 0 0 0 5px rgba(118, 185, 0, 0.35);
+}
+
+.pm4d-card__number {
+    color: var(--nvidia-green);
+    font-size: 0.82rem;
+    font-weight: 800;
+    letter-spacing: 0.12em;
+}
+
+.pm4d-card h2 {
+    margin: 0.55rem 0 0.65rem;
+    color: var(--pm4d-ink);
+    font-size: 1.18rem;
+    line-height: 1.2;
+}
+
+.pm4d-card p {
+    margin: 0 0 1rem;
+    color: var(--pm4d-muted);
+    font-size: 0.95rem;
+    line-height: 1.5;
+}
+
+.pm4d-card__meta {
+    margin-top: auto;
+    color: var(--pm4d-ink);
+    font-size: 0.78rem;
+    font-weight: 700;
+    letter-spacing: 0.04em;
+    text-transform: uppercase;
+}
+
+.pm4d-topic-section {
+    margin: 2.5rem 0 2.75rem;
+}
+
+.pm4d-section-heading {
+    margin-bottom: 1.25rem;
+}
+
+.pm4d-section-heading h2 {
+    margin: 0;
+    font-size: 1.75rem;
+}
+
+.pm4d-topic-grid {
+    display: grid;
+    grid-template-columns: repeat(auto-fit, minmax(230px, 1fr));
+    gap: 1rem;
+}
+
+.pm4d-topic-card,
+.pm4d-topic-card:visited {
+    display: block;
+    min-height: 158px;
+    padding: 1.15rem;
+    color: var(--pm4d-ink);
+    text-decoration: none;
+    background: #ffffff;
+    border: 1px solid var(--pm4d-line);
+    border-left: 4px solid var(--nvidia-green);
+    border-radius: 8px;
+    box-shadow: 0 8px 20px rgba(17, 17, 17, 0.08);
+    transition: transform 160ms ease, box-shadow 160ms ease, border-color 160ms ease;
+}
+
+.pm4d-topic-card:hover {
+    color: var(--pm4d-ink);
+    text-decoration: none;
+    border-color: var(--nvidia-green);
+    box-shadow: 0 12px 26px rgba(17, 17, 17, 0.14);
+    transform: translateY(-3px);
+}
+
+.pm4d-topic-card:focus-visible {
+    color: var(--pm4d-ink);
+    text-decoration: none;
+    border-color: var(--nvidia-green);
+    outline: 3px solid #111111;
+    outline-offset: 3px;
+    box-shadow: 0 0 0 5px rgba(118, 185, 0, 0.35);
+}
+
+.pm4d-topic-card h3 {
+    margin: 0 0 0.6rem;
+    font-size: 1.05rem;
+}
+
+.pm4d-topic-card p {
+    margin: 0;
+    color: var(--pm4d-muted);
+    font-size: 0.92rem;
+    line-height: 1.5;
+}
+
+@media screen and (max-width: 768px) {
+    .pm4d-hero {
+        margin: -1.618em -1.618em 2rem;
+        padding: 3rem 1.618em 4.25rem;
+    }
+
+    .pm4d-hero h1 {
+        font-size: 2rem;
+    }
+
+    .pm4d-card-grid {
+        grid-template-columns: 1fr;
+    }
+}
+
+/* ==================== END ==================== */
diff --git a/docs/_static/nvidia-logo.svg b/docs/_static/nvidia-logo.svg
new file mode 100644
index 0000000..93a5848
--- /dev/null
+++ b/docs/_static/nvidia-logo.svg
@@ -0,0 +1,5 @@
+<svg xmlns="http://www.w3.org/2000/svg" width="328" height="60" viewBox="0 0 328 60" role="img" aria-labelledby="title">
+  <title id="title">NVIDIA</title>
+  <rect width="328" height="60" rx="0" fill="#76b900"/>
+  <text x="164" y="39" text-anchor="middle" font-family="Arial, Helvetica, sans-serif" font-size="30" font-weight="700" letter-spacing="4" fill="#111111">NVIDIA</text>
+</svg>
diff --git a/docs/conf.py b/docs/conf.py
index b1cb199..6b7e0fa 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -12,7 +12,7 @@
 # Add the source directory to the path
 sys.path.insert(0, os.path.abspath("../src"))
 
-# Suppress the ImportWarning emitted when cupy is absent (CPU-only docs build)
+# Suppress the ImportWarning emitted when CuPy is absent during docs builds.
 warnings.filterwarnings("ignore", category=UserWarning, message="CuPy is not installed")
 
 
@@ -71,6 +71,7 @@ def __getattr__(cls, name):
 # -- Options for HTML output -------------------------------------------------
 html_theme = "sphinx_rtd_theme"
 html_static_path = ["_static"]
+html_logo = "_static/nvidia-logo.svg"
 
 html_theme_options = {
     "logo_only": False,
@@ -78,7 +79,7 @@ def __getattr__(cls, name):
     "prev_next_buttons_location": "both",  # Show navigation at top and bottom
     "style_external_links": False,
     "vcs_pageview_mode": "",
-    "style_nav_header_background": "#2980B9",
+    "style_nav_header_background": "#111111",
     # Toc options - optimized for module browsing
     "collapse_navigation": False,  # Keep sidebar expanded
     "sticky_navigation": True,  # Sidebar follows scroll
diff --git a/docs/faq.rst b/docs/faq.rst
index d9266de..1bd0b6b 100644
--- a/docs/faq.rst
+++ b/docs/faq.rst
@@ -43,7 +43,6 @@ a ``UserWarning`` is emitted (visible by default in all standard Python runs):
    slow. Re-install with uv to get CuPy and CUDA-enabled PyTorch in one step
    (pip alone will not select the correct CUDA wheel):
      uv pip install 'physiomotion4d[cuda13]'  # CUDA 13
-     uv pip install 'physiomotion4d[cuda12]'  # CUDA 12
 
 CPU-only mode is suitable for evaluation and small datasets. For production
 workloads an NVIDIA GPU is strongly recommended.
@@ -51,25 +50,15 @@ workloads an NVIDIA GPU is strongly recommended.
 Which CUDA version is required?
 --------------------------------
 
-CUDA 13 and CUDA 12 are both supported. Install the extra that matches your
-system CUDA version:
+CUDA 13 is supported. Install the CUDA 13 extra for GPU acceleration:
 
 .. code-block:: bash
 
-   # CUDA 13 (recommended)
    uv pip install "physiomotion4d[cuda13]"
 
-   # CUDA 12
-   uv pip install "physiomotion4d[cuda12]"
-
-Each extra installs both CuPy and a CUDA-built PyTorch wheel in one step —
-there is no need to install PyTorch separately. The ``[cuda13]`` extra provides
-``cupy-cuda13x>=13.6.0`` and sources PyTorch, torchvision, and torchaudio from
-``https://download.pytorch.org/whl/cu130``. The ``[cuda12]`` extra provides
-``cupy-cuda12x>=12.0.0`` and sources them from
-``https://download.pytorch.org/whl/cu128``. PyTorch is listed in both extras
-so that uv's dependency resolver fetches the GPU wheel instead of the CPU wheel
-from PyPI.
+The extra installs CuPy. In uv-managed source environments, PyTorch,
+torchvision, and torchaudio are sourced from
+``https://download.pytorch.org/whl/cu130`` by default.
 
 What Python version is required?
 ---------------------------------
diff --git a/docs/index.rst b/docs/index.rst
index d919b99..f8840d9 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -1,79 +1,127 @@
 .. PhysioMotion4D documentation master file
 
-=====================================
-PhysioMotion4D Documentation
-=====================================
-
-**Generate anatomic models in Omniverse with physiological motion derived from 4D medical images.**
-
-PhysioMotion4D is a comprehensive medical imaging package that converts 3D and 4D medical scans (particularly heart and lung gated CT data) into dynamic 3D models for visualization in NVIDIA Omniverse. The package provides state-of-the-art deep learning-based image processing, segmentation, registration, and USD file generation capabilities.
-
-.. image:: https://img.shields.io/pypi/v/physiomotion4d.svg
-   :target: https://pypi.org/project/physiomotion4d/
-   :alt: PyPI Version
-
-.. image:: https://img.shields.io/pypi/pyversions/physiomotion4d.svg
-   :target: https://pypi.org/project/physiomotion4d/
-   :alt: Python Versions
-
-.. image:: https://img.shields.io/badge/license-Apache%202.0-blue.svg
-   :target: https://github.com/Project-MONAI/physiomotion4d/blob/main/LICENSE
-   :alt: License
-
-.. image:: https://img.shields.io/github/actions/workflow/status/Project-MONAI/physiomotion4d/nightly-health.yml?branch=main&label=Nightly%20CI%20Tests
-   :target: https://github.com/Project-MONAI/physiomotion4d/actions/workflows/nightly-health.yml
-   :alt: Nightly CI Tests
-
-.. image:: https://img.shields.io/badge/tests-Windows%20%7C%20Linux%20%7C%20Python%203.10--3.12-blue
-   :target: https://github.com/Project-MONAI/physiomotion4d/actions/workflows/ci.yml
-   :alt: Test Matrix: Windows, Linux, Python 3.10-3.12
-
-.. image:: https://codecov.io/gh/Project-MONAI/physiomotion4d/branch/main/graph/badge.svg
-   :target: https://codecov.io/gh/Project-MONAI/physiomotion4d
-   :alt: Test Coverage
-
-🚀 Key Features
-===============
-
-* **Complete 4D Medical Imaging Pipeline**: End-to-end processing from 4D CT/MR data to animated USD models
-* **Multiple AI Segmentation Methods**: TotalSegmentator and Simpleware cardiac segmentation
-* **Deep Learning Registration**: GPU-accelerated image registration using Icon algorithm
-* **NVIDIA Omniverse Integration**: Direct USD file export for medical visualization
-* **Physiological Motion Analysis**: Capture and visualize cardiac and respiratory motion
-* **Flexible Workflow Control**: Step-based processing with checkpoint management
-
-📋 Supported Applications
-==========================
-
-* **Cardiac Imaging**: Heart-gated CT/MR processing with cardiac motion analysis
-* **Pulmonary Imaging**: Lung 4D-CT/MR processing with respiratory motion tracking
-* **Medical Education**: Interactive 3D anatomical models with physiological motion
-* **Research Visualization**: Advanced medical imaging research in Omniverse
-* **Clinical Planning**: Dynamic anatomical models for treatment planning
-
-.. tip::
-
-   **Getting Started with Code Examples:**
-
-   This documentation uses examples from the CLI commands (``physiomotion4d-heart-gated-ct``,
-   ``physiomotion4d-create-statistical-model``, ``physiomotion4d-fit-statistical-model-to-patient``)
-   and their implementations in ``src/physiomotion4d/cli/``,
-   which contain production-ready workflows and proper library usage patterns. The repository also includes
-   an ``experiments/`` directory with research prototypes that can inspire adaptations to
-   new digital twin models and anatomical regions—see the experiments README for details on
-   how to adapt these conceptual patterns to your own applications.
+.. title:: PhysioMotion4D Documentation
+
+.. raw:: html
+
+   <section class="pm4d-hero">
+     <div class="pm4d-hero__brand">
+       <img src="_static/nvidia-logo.svg" alt="NVIDIA logo">
+     </div>
+     <p class="pm4d-kicker">PhysioMotion4D tutorials</p>
+     <h1>Build animated medical USD workflows for NVIDIA Omniverse</h1>
+     <p>
+       PhysioMotion4D converts 3D and 4D medical scans into dynamic OpenUSD
+       assets for NVIDIA Omniverse. Start with the tutorial cards, then use the
+       documentation sections below for installation, CLI workflows, API
+       references, developer notes, and contribution guidance.
+     </p>
+   </section>
+
+   <section class="pm4d-card-grid" aria-label="Tutorial cards">
+     <a class="pm4d-card" href="tutorials.html#tutorial-1-heart-gated-ct-to-animated-usd">
+       <span class="pm4d-card__number">01</span>
+       <h2>Heart-Gated CT to Animated USD</h2>
+       <p>Convert cardiac 4D CT frames into registered contours and an animated OpenUSD model.</p>
+       <span class="pm4d-card__meta">Slicer-Heart-CT</span>
+     </a>
+     <a class="pm4d-card" href="tutorials.html#tutorial-2-ct-segmentation-to-vtk-surfaces">
+       <span class="pm4d-card__number">02</span>
+       <h2>CT Segmentation to VTK Surfaces</h2>
+       <p>Segment one CT phase and export patient anatomy as VTK PolyData surfaces.</p>
+       <span class="pm4d-card__meta">Slicer-Heart-CT</span>
+     </a>
+     <a class="pm4d-card" href="tutorials.html#tutorial-3-fit-statistical-model-to-patient">
+       <span class="pm4d-card__number">03</span>
+       <h2>Fit Statistical Model to Patient</h2>
+       <p>Fit a PCA heart model to patient-specific anatomy for model-based reconstruction.</p>
+       <span class="pm4d-card__meta">KCL-Heart-Model</span>
+     </a>
+     <a class="pm4d-card" href="tutorials.html#tutorial-4-create-a-pca-shape-model">
+       <span class="pm4d-card__number">04</span>
+       <h2>Create a PCA Shape Model</h2>
+       <p>Build a statistical shape model from aligned cardiac meshes.</p>
+       <span class="pm4d-card__meta">KCL-Heart-Model</span>
+     </a>
+     <a class="pm4d-card" href="tutorials.html#tutorial-5-vtk-surface-series-to-animated-usd">
+       <span class="pm4d-card__number">05</span>
+       <h2>VTK Surface Series to Animated USD</h2>
+       <p>Convert VTK meshes into a time-sampled USD scene for Omniverse playback.</p>
+       <span class="pm4d-card__meta">Tutorial 2 output</span>
+     </a>
+     <a class="pm4d-card" href="tutorials.html#tutorial-6-reconstruct-high-resolution-4d-ct">
+       <span class="pm4d-card__number">06</span>
+       <h2>Reconstruct High-Resolution 4D CT</h2>
+       <p>Register respiratory CT phases and reconstruct a higher-resolution 4D volume series.</p>
+       <span class="pm4d-card__meta">DirLab-4DCT</span>
+     </a>
+   </section>
+
+   <section class="pm4d-topic-section" aria-label="Documentation topics">
+     <div class="pm4d-section-heading">
+       <p class="pm4d-kicker">Documentation</p>
+       <h2>Explore the rest of the docs</h2>
+     </div>
+     <div class="pm4d-topic-grid">
+       <a class="pm4d-topic-card" href="installation.html">
+         <h3>Installation</h3>
+         <p>Set up PhysioMotion4D with CUDA extras, CPU-only options, and required system tools.</p>
+       </a>
+       <a class="pm4d-topic-card" href="quickstart.html">
+         <h3>Getting Started</h3>
+         <p>Run your first workflow and understand the basic CT-to-USD processing path.</p>
+       </a>
+       <a class="pm4d-topic-card" href="examples.html">
+         <h3>Examples</h3>
+         <p>Review focused usage patterns for common cardiac, lung, segmentation, and USD tasks.</p>
+       </a>
+       <a class="pm4d-topic-card" href="cli_scripts/overview.html">
+         <h3>CLI Workflows</h3>
+         <p>Use production command-line workflows for conversion, reconstruction, modeling, and USD export.</p>
+       </a>
+       <a class="pm4d-topic-card" href="api/index.html">
+         <h3>API Reference</h3>
+         <p>Browse classes and modules for workflows, segmentation, registration, USD, and utilities.</p>
+       </a>
+       <a class="pm4d-topic-card" href="developer/architecture.html">
+         <h3>Developer Docs</h3>
+         <p>Understand architecture, extension points, coordinate transforms, and implementation boundaries.</p>
+       </a>
+       <a class="pm4d-topic-card" href="contributing.html">
+         <h3>Contributing</h3>
+         <p>Follow repository conventions for code style, testing, documentation, and pull requests.</p>
+       </a>
+       <a class="pm4d-topic-card" href="testing.html">
+         <h3>Testing</h3>
+         <p>Run the fast test suite, data-gated tutorial tests, and regression checks.</p>
+       </a>
+       <a class="pm4d-topic-card" href="troubleshooting.html">
+         <h3>Troubleshooting</h3>
+         <p>Diagnose environment, data, segmentation, registration, and USD playback issues.</p>
+       </a>
+     </div>
+   </section>
+
+Tutorial Details
+================
+
+See :doc:`tutorials` for the recommended run order, commands, datasets, and
+per-tutorial implementation details.
 
 .. toctree::
    :maxdepth: 2
    :caption: Getting Started
+   :hidden:
 
    installation
    quickstart
+   tutorials
    examples
 
 .. toctree::
    :maxdepth: 2
    :caption: CLI & Scripts Guide
+   :hidden:
 
    cli_scripts/overview
    cli_scripts/heart_gated_ct
@@ -88,6 +136,7 @@ PhysioMotion4D is a comprehensive medical imaging package that converts 3D and 4
 .. toctree::
    :maxdepth: 2
    :caption: API Reference
+   :hidden:
 
    api/index
    api/base
@@ -101,6 +150,7 @@ PhysioMotion4D is a comprehensive medical imaging package that converts 3D and 4
 .. toctree::
    :maxdepth: 2
    :caption: Developer Guides
+   :hidden:
 
    developer/architecture
    developer/extending
@@ -115,6 +165,7 @@ PhysioMotion4D is a comprehensive medical imaging package that converts 3D and 4
 .. toctree::
    :maxdepth: 1
    :caption: Contributing
+   :hidden:
 
    contributing
    testing
@@ -122,6 +173,7 @@ PhysioMotion4D is a comprehensive medical imaging package that converts 3D and 4
 .. toctree::
    :maxdepth: 1
    :caption: Additional Resources
+   :hidden:
 
    faq
    troubleshooting
diff --git a/docs/installation.rst b/docs/installation.rst
index 677759e..ad2e4cc 100644
--- a/docs/installation.rst
+++ b/docs/installation.rst
@@ -11,7 +11,7 @@ System Requirements
 -------------------
 
 * **Python**: 3.10, 3.11, or 3.12
-* **GPU**: NVIDIA GPU with CUDA 13 (default) or CUDA 12 — recommended for production use; CPU-only installation is supported but will be slow and will emit a runtime warning
+* **GPU**: NVIDIA GPU with CUDA 13 — recommended for production use; CPU-only PyPI installation is supported but will be slow and will emit a runtime warning
 * **RAM**: 16GB minimum (32GB+ recommended for large datasets)
 * **Storage**: 10GB+ for package and model weights
 * **Visualization**: NVIDIA Omniverse (optional, for USD visualization)
@@ -22,7 +22,7 @@ Software Dependencies
 PhysioMotion4D relies on several key packages:
 
 * **Medical Imaging**: ITK, TubeTK, MONAI, nibabel, PyVista
-* **AI/ML**: PyTorch, CuPy (CUDA 13 default; CUDA 12 via ``[cuda12]`` extra), transformers, MONAI
+* **AI/ML**: PyTorch, CuPy (CUDA 13), transformers, MONAI
 * **Registration**: icon-registration, unigradicon
 * **Visualization**: USD-core, PyVista
 * **Segmentation**: TotalSegmentator
@@ -35,7 +35,7 @@ Method 1: Install from PyPI (Recommended)
 
 The simplest way to install PhysioMotion4D is from PyPI.
 
-CPU-only install (evaluation / no GPU):
+CPU-only PyPI install (evaluation / no GPU):
 
 .. code-block:: bash
 
@@ -50,7 +50,6 @@ import time (visible by default in all standard Python runs):
    slow. Re-install with uv to get CuPy and CUDA-enabled PyTorch in one step
    (pip alone will not select the correct CUDA wheel):
      uv pip install 'physiomotion4d[cuda13]'  # CUDA 13
-     uv pip install 'physiomotion4d[cuda12]'  # CUDA 12
 
 CUDA 13 install (recommended for production):
 
@@ -58,17 +57,9 @@ CUDA 13 install (recommended for production):
 
    uv pip install "physiomotion4d[cuda13]"
 
-CUDA 12 install:
-
-.. code-block:: bash
-
-   uv pip install "physiomotion4d[cuda12]"
-
-The ``[cuda13]`` and ``[cuda12]`` extras install both CuPy and the correct
-CUDA-built PyTorch wheel in one step. PyTorch is listed inside the extras so
-that uv's dependency resolver fetches the GPU wheel from the PyTorch index
-rather than the CPU wheel from PyPI. There is no need to install PyTorch
-separately.
+The ``[cuda13]`` extra installs CuPy. In uv-managed source environments,
+PyTorch, torchvision, and torchaudio resolve from the CUDA 13.0 PyTorch wheel
+index. There is no need to install PyTorch separately.
 
 For development with NVIDIA NIM cloud services:
 
@@ -114,24 +105,19 @@ For development or to get the latest features:
 
 **Step 4: Install PhysioMotion4D**
 
-CPU-only (evaluation / no GPU):
+Default uv-managed source install:
 
 .. code-block:: bash
 
    uv pip install -e "."
 
-With uv (CUDA 13, recommended for production):
+This uses the CUDA 13.0 PyTorch wheel index by default. To add CuPy for CUDA 13
+GPU acceleration:
 
 .. code-block:: bash
 
    uv pip install -e ".[cuda13]"
 
-With uv (CUDA 12):
-
-.. code-block:: bash
-
-   uv pip install -e ".[cuda12]"
-
 Optional Dependencies
 =====================
 
@@ -209,11 +195,10 @@ GPU Setup
 CUDA Installation
 -----------------
 
-An NVIDIA GPU is strongly recommended. Two CUDA versions are supported via
-optional extras:
+An NVIDIA GPU is strongly recommended. CUDA 13 is supported via the optional
+extra:
 
 * **CUDA 13** — installed when you use the ``[cuda13]`` extra (recommended)
-* **CUDA 12** — installed when you use the ``[cuda12]`` extra
 
 A plain ``pip install physiomotion4d`` installs a CPU-only build. It runs
 without error but emits a ``UserWarning`` at import time and will be
@@ -230,10 +215,8 @@ If CUDA is not yet installed, download the CUDA Toolkit from
 PyTorch with CUDA
 -----------------
 
-A plain install pulls a CPU-only PyTorch wheel from PyPI. The ``[cuda13]``
-extra sources PyTorch, torchvision, and torchaudio from the
-``https://download.pytorch.org/whl/cu130`` index. The ``[cuda12]`` extra
-sources them from ``https://download.pytorch.org/whl/cu128``. To verify the
+uv-managed source environments source PyTorch, torchvision, and torchaudio from
+the ``https://download.pytorch.org/whl/cu130`` index by default. To verify the
 active version:
 
 .. code-block:: python
diff --git a/docs/quickstart.rst b/docs/quickstart.rst
index 0526594..9dee120 100644
--- a/docs/quickstart.rst
+++ b/docs/quickstart.rst
@@ -2,7 +2,41 @@
 Quick Start
 ===========
 
-This guide will help you get started with PhysioMotion4D quickly. We'll walk through a basic workflow for processing heart-gated CT data.
+This guide will help you get started with PhysioMotion4D quickly.
+
+.. _tutorials:
+
+Tutorials
+=========
+
+The ``tutorials/`` directory contains six end-to-end scripts, one per major
+workflow.  Each script is self-contained, includes its own ``argparse`` CLI, and
+can be imported as a module from the test suite.
+
+See :doc:`tutorials` for the NVIDIA-styled tutorial card index, dataset
+requirements, commands, and workflow details.
+
+After preparing the Slicer-Heart-CT data, run the first two tutorials:
+
+.. code-block:: bash
+
+   python tutorials/tutorial_01_heart_gated_ct_to_usd.py \
+       --data-dir ./data --output-dir ./output/tutorial_01 \
+       --registration-method ants
+
+   python tutorials/tutorial_02_ct_to_vtk.py \
+       --data-dir ./data --output-dir ./output/tutorial_02
+
+Each script prints the paths of outputs and screenshots it created.
+See ``tutorials/README.md`` for dataset download instructions and the
+recommended run order.
+
+Recommended run order:
+
+1. Tutorials 1 and 2 first, after preparing Slicer-Heart-CT data.
+2. Tutorial 5 after Tutorial 2 (consumes Tutorial 2 output).
+3. Tutorials 3 and 4 after downloading KCL-Heart-Model.
+4. Tutorial 6 after downloading DirLab-4DCT.
 
 Prerequisites
 =============
@@ -10,7 +44,7 @@ Prerequisites
 Before starting, ensure you have:
 
 * PhysioMotion4D installed (see :doc:`installation`)
-* NVIDIA GPU with CUDA 13 (default) or CUDA 12 — recommended for production performance; see :doc:`installation` for the ``[cuda13]`` and ``[cuda12]`` extras. A CPU-only install works for evaluation but is slow.
+* NVIDIA GPU with CUDA 13 - recommended for production performance; see :doc:`installation` for the ``[cuda13]`` extra. A CPU-only PyPI install works for evaluation but is slow.
 * 4D cardiac CT data or access to sample datasets
 
 Basic Workflow
@@ -216,7 +250,7 @@ In NVIDIA Omniverse
 
 1. Open NVIDIA Omniverse
 2. Launch USD Composer or USD Presenter
-3. File → Open → Select your generated `.usd` file
+3. File -> Open -> Select your generated `.usd` file
 4. Press Play to view the animation
 
 Using USD Viewer
diff --git a/docs/troubleshooting.rst b/docs/troubleshooting.rst
index b7bd81c..60f86eb 100644
--- a/docs/troubleshooting.rst
+++ b/docs/troubleshooting.rst
@@ -42,17 +42,14 @@ returning ``False``, or runtime messages indicating a CUDA library version confl
 **Cause**: The installed ``cupy`` or PyTorch wheel was built for a different CUDA
 version than the one present on the system.
 
-**Solution**: Install the extra that matches your system CUDA version:
+**Solution**: Install the CUDA 13 extra:
 
 .. code-block:: bash
 
-   # CUDA 13
    uv pip install "physiomotion4d[cuda13]"
 
-   # CUDA 12
-   uv pip install "physiomotion4d[cuda12]"
-
-Each extra installs both CuPy and the correct CUDA-built PyTorch wheel.
+The extra installs CuPy. In uv-managed source environments, PyTorch resolves
+from the CUDA 13.0 wheel index.
 
 Verify the active CUDA version before reinstalling:
 
diff --git a/docs/tutorials.rst b/docs/tutorials.rst
new file mode 100644
index 0000000..0da306c
--- /dev/null
+++ b/docs/tutorials.rst
@@ -0,0 +1,203 @@
+=========
+Tutorials
+=========
+
+.. raw:: html
+
+   <section class="pm4d-hero">
+     <div class="pm4d-hero__brand">
+       <img src="_static/nvidia-logo.svg" alt="NVIDIA logo">
+     </div>
+     <p class="pm4d-kicker">PhysioMotion4D tutorials</p>
+     <h1>Build animated medical USD workflows for NVIDIA Omniverse</h1>
+     <p>
+       Six focused tutorials walk through CT segmentation, registration,
+       statistical model fitting, high-resolution 4D reconstruction, and USD
+       export. Each card links to implementation details, datasets, and the
+       command used to run the workflow.
+     </p>
+   </section>
+
+   <section class="pm4d-card-grid" aria-label="Tutorial cards">
+     <a class="pm4d-card" href="#tutorial-1-heart-gated-ct-to-animated-usd">
+       <span class="pm4d-card__number">01</span>
+       <h2>Heart-Gated CT to Animated USD</h2>
+       <p>Convert cardiac 4D CT frames into registered contours and an animated OpenUSD model.</p>
+       <span class="pm4d-card__meta">Slicer-Heart-CT</span>
+     </a>
+     <a class="pm4d-card" href="#tutorial-2-ct-segmentation-to-vtk-surfaces">
+       <span class="pm4d-card__number">02</span>
+       <h2>CT Segmentation to VTK Surfaces</h2>
+       <p>Segment one CT phase and export patient anatomy as VTK PolyData surfaces.</p>
+       <span class="pm4d-card__meta">Slicer-Heart-CT</span>
+     </a>
+     <a class="pm4d-card" href="#tutorial-3-fit-statistical-model-to-patient">
+       <span class="pm4d-card__number">03</span>
+       <h2>Fit Statistical Model to Patient</h2>
+       <p>Fit a PCA heart model to patient-specific anatomy for model-based reconstruction.</p>
+       <span class="pm4d-card__meta">KCL-Heart-Model</span>
+     </a>
+     <a class="pm4d-card" href="#tutorial-4-create-a-pca-shape-model">
+       <span class="pm4d-card__number">04</span>
+       <h2>Create a PCA Shape Model</h2>
+       <p>Build a statistical shape model from aligned cardiac meshes.</p>
+       <span class="pm4d-card__meta">KCL-Heart-Model</span>
+     </a>
+     <a class="pm4d-card" href="#tutorial-5-vtk-surface-series-to-animated-usd">
+       <span class="pm4d-card__number">05</span>
+       <h2>VTK Surface Series to Animated USD</h2>
+       <p>Convert VTK meshes into a time-sampled USD scene for Omniverse playback.</p>
+       <span class="pm4d-card__meta">Tutorial 2 output</span>
+     </a>
+     <a class="pm4d-card" href="#tutorial-6-reconstruct-high-resolution-4d-ct">
+       <span class="pm4d-card__number">06</span>
+       <h2>Reconstruct High-Resolution 4D CT</h2>
+       <p>Register respiratory CT phases and reconstruct a higher-resolution 4D volume series.</p>
+       <span class="pm4d-card__meta">DirLab-4DCT</span>
+     </a>
+   </section>
+
+Recommended Run Order
+=====================
+
+1. Run Tutorials 1 and 2 after preparing Slicer-Heart-CT data.
+2. Run Tutorial 5 after Tutorial 2 because it consumes Tutorial 2 output.
+3. Run Tutorials 3 and 4 after downloading KCL-Heart-Model.
+4. Run Tutorial 6 after downloading DirLab-4DCT.
+
+Tutorial 1: Heart-Gated CT to Animated USD
+==========================================
+
+Script
+   ``tutorials/tutorial_01_heart_gated_ct_to_usd.py``
+
+Workflow
+   ``WorkflowConvertHeartGatedCTToUSD``
+
+Dataset
+   Slicer-Heart-CT, prepared before running the tutorial.
+
+Command
+   .. code-block:: bash
+
+      python tutorials/tutorial_01_heart_gated_ct_to_usd.py \
+          --data-dir ./data --output-dir ./output/tutorial_01 \
+          --registration-method ants
+
+Outputs
+   Registered phase images, transformed contours, preview screenshots, and an
+   animated USD model.
+
+Tutorial 2: CT Segmentation to VTK Surfaces
+===========================================
+
+Script
+   ``tutorials/tutorial_02_ct_to_vtk.py``
+
+Workflow
+   ``WorkflowConvertCTToVTK``
+
+Dataset
+   Slicer-Heart-CT, prepared before running the tutorial.
+
+Command
+   .. code-block:: bash
+
+      python tutorials/tutorial_02_ct_to_vtk.py \
+          --data-dir ./data --output-dir ./output/tutorial_02
+
+Outputs
+   Segmentation artifacts, VTK PolyData surfaces, and preview screenshots.
+
+Tutorial 3: Fit Statistical Model to Patient
+============================================
+
+Script
+   ``tutorials/tutorial_03_fit_statistical_model_to_patient.py``
+
+Workflow
+   ``WorkflowFitStatisticalModelToPatient``
+
+Dataset
+   KCL-Heart-Model, downloaded manually.
+
+Command
+   .. code-block:: bash
+
+      python tutorials/tutorial_03_fit_statistical_model_to_patient.py \
+          --data-dir ./data --output-dir ./output/tutorial_03
+
+Outputs
+   Patient-fitted statistical model surfaces and registration diagnostics.
+
+Tutorial 4: Create a PCA Shape Model
+====================================
+
+Script
+   ``tutorials/tutorial_04_create_statistical_model.py``
+
+Workflow
+   ``WorkflowCreateStatisticalModel``
+
+Dataset
+   KCL-Heart-Model, downloaded manually.
+
+Command
+   .. code-block:: bash
+
+      python tutorials/tutorial_04_create_statistical_model.py \
+          --data-dir ./data --output-dir ./output/tutorial_04
+
+Outputs
+   PCA model files, mean shape, and component diagnostics.
+
+Tutorial 5: VTK Surface Series to Animated USD
+==============================================
+
+Script
+   ``tutorials/tutorial_05_vtk_to_usd.py``
+
+Workflow
+   ``WorkflowConvertVTKToUSD``
+
+Dataset
+   Output from Tutorial 2.
+
+Command
+   .. code-block:: bash
+
+      python tutorials/tutorial_05_vtk_to_usd.py \
+          --data-dir ./data --output-dir ./output/tutorial_05 \
+          --vtk-file output/tutorial_02/patient_surfaces.vtp
+
+Outputs
+   Time-sampled USD scene and conversion logs for Omniverse inspection.
+
+Tutorial 6: Reconstruct High-Resolution 4D CT
+=============================================
+
+Script
+   ``tutorials/tutorial_06_reconstruct_highres_4d_ct.py``
+
+Workflow
+   ``WorkflowReconstructHighres4DCT``
+
+Dataset
+   DirLab-4DCT, downloaded manually.
+
+Command
+   .. code-block:: bash
+
+      python tutorials/tutorial_06_reconstruct_highres_4d_ct.py \
+          --data-dir ./data --output-dir ./output/tutorial_06
+
+Outputs
+   Registered respiratory phases, reconstructed high-resolution CT volumes,
+   and preview screenshots.
+
+Dataset Notes
+=============
+
+The repository-level ``tutorials/README.md`` has the most detailed dataset
+preparation notes. The tutorials are also exercised by ``tests/test_tutorials.py``
+behind the experiment marker.
diff --git a/pyproject.toml b/pyproject.toml
index 209e872..d7f4087 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -87,34 +87,19 @@ dependencies = [
 [tool.uv]
 link-mode = "copy"
 
-# torch/torchvision/torchaudio appear in both core deps (CPU fallback) and in
-# the cuda13/cuda12 extras. The extras listing is required by uv: it is what
-# triggers the CUDA-specific wheel index (cu130 / cu128) instead of PyPI.
-# Without it, uv would install CPU-only torch even when a CUDA extra is active.
+# PyTorch packages default to the CUDA 13.0 wheel index for uv-managed
+# environments. This avoids ambiguous resolution across multiple PyTorch CUDA
+# indexes.
 [tool.uv.sources]
-torch = [
-  { index = "pytorch-cu130", extra = "cuda13" },
-  { index = "pytorch-cu128", extra = "cuda12" },
-]
-torchvision = [
-  { index = "pytorch-cu130", extra = "cuda13" },
-  { index = "pytorch-cu128", extra = "cuda12" },
-]
-torchaudio = [
-  { index = "pytorch-cu130", extra = "cuda13" },
-  { index = "pytorch-cu128", extra = "cuda12" },
-]
+torch = { index = "pytorch-cu130" }
+torchvision = { index = "pytorch-cu130" }
+torchaudio = { index = "pytorch-cu130" }
 
 [[tool.uv.index]]
 name = "pytorch-cu130"
 url = "https://download.pytorch.org/whl/cu130"
 explicit = true
 
-[[tool.uv.index]]
-name = "pytorch-cu128"
-url = "https://download.pytorch.org/whl/cu128"
-explicit = true
-
 [project.optional-dependencies]
 cuda13 = [
     "cupy-cuda13x>=13.6.0",
@@ -122,12 +107,6 @@ cuda13 = [
     "torchvision",
     "torchaudio",
 ]
-cuda12 = [
-    "cupy-cuda12x>=12.0.0",
-    "torch>=2.0.0,<3.0.0",
-    "torchvision",
-    "torchaudio",
-]
 nim = [
     # Cloud and data handling
     "google-auth>=2.0.0",
@@ -259,6 +238,7 @@ addopts = [
     "--cov-report=xml"
 ]
 testpaths = ["tests"]
+pythonpath = ["."]
 markers = [
     "unit: marks tests as unit tests (fast, isolated)",
     "integration: marks tests as integration tests (slower, multiple components)",
diff --git a/src/physiomotion4d/test_tools.py b/src/physiomotion4d/test_tools.py
index d46dbe8..2961704 100644
--- a/src/physiomotion4d/test_tools.py
+++ b/src/physiomotion4d/test_tools.py
@@ -11,7 +11,7 @@
 import logging
 import shutil
 from pathlib import Path
-from typing import Any
+from typing import Any, Literal, Optional
 
 import itk
 import numpy as np
@@ -33,9 +33,9 @@ def set_create_baseline_if_missing(value: bool) -> None:
 
 class TestTools(PhysioMotion4DBase):
     """
-    Utilities for pytest image comparison: baseline directory, results directory,
+    Utilities for pytest image comparison: baseline directory, result directory,
     and comparison with configurable tolerances. Inherits from PhysioMotion4DBase
-    for logging. All image I/O uses ITK .mha with compression.
+    for logging. All image I/O uses ITK with compression where supported.
     """
 
     # Prevent pytest from collecting this as a test class
@@ -47,11 +47,28 @@ def __init__(
         baselines_dir: Path,
         class_name: str,
         *,
+        results_output_dir: Optional[Path] = None,
         log_level: int = logging.INFO,
     ) -> None:
+        """Initialize test helpers.
+
+        Args:
+            results_dir: Root directory for result artifacts.
+            baselines_dir: Root directory for baseline artifacts.
+            class_name: Subdirectory name for baselines and, by default, results.
+            results_output_dir: Optional exact directory for result artifacts.
+                When set, results are read and written there instead of
+                ``results_dir / class_name`` while baselines still use
+                ``baselines_dir / class_name``.
+            log_level: Logging level.
+        """
         super().__init__(class_name=class_name, log_level=log_level)
 
-        self._results_dir = results_dir / class_name
+        self._results_dir = (
+            results_output_dir
+            if results_output_dir is not None
+            else results_dir / class_name
+        )
         self._results_dir.mkdir(parents=True, exist_ok=True)
 
         self._baselines_dir = baselines_dir / class_name
@@ -149,11 +166,11 @@ def transform_difference(self) -> Any:
         return self._last_transform_difference_transform
 
     def write_result_image(self, image: Any, filename: str) -> None:
-        """Write the image to the results directory."""
+        """Write the image to the configured result artifact directory."""
         itk.imwrite(image, str(self._results_dir / filename), compression=True)
 
     def write_result_transform(self, transform: Any, filename: str) -> None:
-        """Write the transform to the results directory."""
+        """Write the transform to the configured result artifact directory."""
         itk.transformwrite(
             transform, str(self._results_dir / filename), compression=True
         )
@@ -342,3 +359,124 @@ def compare_result_to_baseline_image(
             )
 
         return passed
+
+    def save_screenshot_mesh(
+        self,
+        mesh: Any,  # pv.PolyData
+        filename: str,
+        *,
+        camera_position: Literal["xy", "xz", "yz", "yx", "zx", "zy", "iso"] = "iso",
+        window_size: tuple[int, int] = (800, 600),
+        color: str = "pink",
+        opacity: float = 0.9,
+    ) -> Path:
+        """Render a PyVista mesh off-screen and save a PNG.
+
+        Saves to the configured result artifact directory. On Linux headless
+        environments, calls pv.start_xvfb() before rendering (no-op when a
+        display is present).
+
+        Args:
+            mesh: PyVista PolyData or compatible mesh object.
+            filename: Output PNG filename, relative to the result artifact dir.
+            camera_position: PyVista camera preset, e.g. ``'iso'``, ``'xy'``, ``'xz'``.
+            window_size: Off-screen render size ``(width, height)`` in pixels.
+            color: Mesh color string accepted by PyVista.
+            opacity: Mesh opacity in [0, 1].
+
+        Returns:
+            Absolute Path to the saved PNG.
+        """
+        import pyvista as pv
+
+        xvfb_started = False
+        try:
+            pv.start_xvfb()
+            xvfb_started = True
+        except Exception:
+            pass
+
+        output_path = self._results_dir / filename
+        plotter = pv.Plotter(off_screen=True, window_size=list(window_size))
+        try:
+            plotter.add_mesh(mesh, color=color, opacity=opacity)
+            plotter.camera_position = camera_position
+            plotter.screenshot(str(output_path))
+        finally:
+            plotter.close()
+            if xvfb_started and hasattr(pv, "stop_xvfb"):
+                pv.stop_xvfb()
+        self.log_info("Screenshot saved: %s", output_path)
+        return output_path
+
+    def save_screenshot_image_slice(
+        self,
+        image: Any,  # itk.Image, axes X Y Z in RAS world space
+        filename: str,
+        *,
+        axis: int = 0,
+        slice_fraction: float = 0.5,
+        colormap: str = "gray",
+        vmin: Optional[float] = None,
+        vmax: Optional[float] = None,
+        overlay_mask: Optional[Any] = None,  # itk.Image same spatial extent
+        overlay_alpha: float = 0.4,
+    ) -> Path:
+        """Extract one slice from an ITK image and save a PNG via matplotlib.
+
+        The numpy array from ``itk.array_view_from_image`` has shape ``(Z, Y, X)``
+        (ITK stores X fastest; numpy reverses the axis order). Axis indices:
+        - axis=0: axial (constant-Z plane)
+        - axis=1: coronal (constant-Y plane)
+        - axis=2: sagittal (constant-X plane)
+
+        Saves to the configured result artifact directory.
+
+        Args:
+            image: 3-D ``itk.Image`` in RAS world space, axes X Y Z.
+            filename: Output PNG filename, relative to the result artifact dir.
+            axis: Numpy axis along which to slice (0=axial, 1=coronal, 2=sagittal).
+            slice_fraction: Fractional position along ``axis`` in [0, 1].
+            colormap: Matplotlib colormap name for the base image.
+            vmin: Lower clamp for display; None means data minimum.
+            vmax: Upper clamp for display; None means data maximum.
+            overlay_mask: Optional binary ITK mask rendered as a semi-transparent
+                overlay. Must have the same spatial extent as ``image``.
+            overlay_alpha: Opacity of the mask overlay in [0, 1].
+
+        Returns:
+            Absolute Path to the saved PNG.
+        """
+        import matplotlib.pyplot as plt
+
+        arr = np.asarray(itk.array_view_from_image(image), dtype=np.float64)
+        idx = int(arr.shape[axis] * slice_fraction)
+        idx = max(0, min(idx, arr.shape[axis] - 1))
+
+        slices: list[Any] = [slice(None)] * arr.ndim
+        slices[axis] = idx
+        slice_data = arr[tuple(slices)]
+
+        output_path = self._results_dir / filename
+        fig, ax = plt.subplots(figsize=(6, 6))
+        try:
+            ax.imshow(slice_data, cmap=colormap, vmin=vmin, vmax=vmax, origin="lower")
+
+            if overlay_mask is not None:
+                mask_arr = np.asarray(
+                    itk.array_view_from_image(overlay_mask), dtype=np.float64
+                )
+                mask_slice = mask_arr[tuple(slices)]
+                ax.imshow(
+                    np.ma.masked_where(mask_slice == 0, mask_slice),
+                    cmap="autumn",
+                    alpha=overlay_alpha,
+                    origin="lower",
+                )
+
+            ax.axis("off")
+            fig.savefig(str(output_path), bbox_inches="tight", dpi=100)
+        finally:
+            plt.close(fig)
+        self.log_info("Screenshot saved: %s", output_path)
+        return output_path
diff --git a/src/physiomotion4d/vtk_to_usd/usd_mesh_converter.py b/src/physiomotion4d/vtk_to_usd/usd_mesh_converter.py
index 31a5966..9ee0b8e 100644
--- a/src/physiomotion4d/vtk_to_usd/usd_mesh_converter.py
+++ b/src/physiomotion4d/vtk_to_usd/usd_mesh_converter.py
@@ -94,7 +94,8 @@ def create_mesh(
         mesh.CreateFaceVertexCountsAttr(face_counts_vt)
         mesh.CreateFaceVertexIndicesAttr(face_indices_vt)
 
-        # Set points (time-varying if time_code provided)
+        # Set points (time-varying if time_code provided). Also author a
+        # default value for readers that inspect the prim without a time code.
         points_attr = mesh.CreatePointsAttr()
         if time_code is not None:
             if points_attr.Get() is None:
diff --git a/src/physiomotion4d/workflow_convert_heart_gated_ct_to_usd.py b/src/physiomotion4d/workflow_convert_heart_gated_ct_to_usd.py
index 1610c2d..fda4b5b 100644
--- a/src/physiomotion4d/workflow_convert_heart_gated_ct_to_usd.py
+++ b/src/physiomotion4d/workflow_convert_heart_gated_ct_to_usd.py
@@ -169,18 +169,22 @@ def _load_time_series(self) -> None:
             self.converter.load_nrrd_3d(self.input_filenames)
 
         self._num_time_points = self.converter.get_number_of_3d_images()
+        if self._num_time_points <= 0:
+            raise ValueError("No time-series images were produced from input data")
 
         # Load all time series images into memory
         for i in range(self._num_time_points):
             self._time_series_images.append(self.converter.get_3d_image(i))
+        if not self._time_series_images:
+            raise ValueError("No time-series images were loaded from input data")
 
         # Load reference image
         if self.reference_image_filename:
             self._fixed_image = itk.imread(self.reference_image_filename)
         else:
-            # Use 70% frame as reference if none specified
-            mid_frame = int(self._num_time_points * 0.7)
-            self._fixed_image = self._time_series_images[mid_frame]
+            # Use 70% frame as reference if none specified.
+            reference_frame = int(self._num_time_points * 0.7)
+            self._fixed_image = self._time_series_images[reference_frame]
             itk.imwrite(
                 self._fixed_image,
                 os.path.join(self.output_directory, "fixed_image.mha"),
diff --git a/tests/test_tutorials.py b/tests/test_tutorials.py
new file mode 100644
index 0000000..3540cfd
--- /dev/null
+++ b/tests/test_tutorials.py
@@ -0,0 +1,423 @@
+"""Experiment tests that run each tutorial end-to-end and compare screenshots.
+
+Each test class maps to one tutorial script.  Tests are gated behind
+``--run-experiments`` (handled by conftest.py) and require the relevant dataset
+to be present (see data/README.md).
+
+Screenshot comparison uses the existing ITK-based baseline infrastructure:
+
+1. The tutorial's ``run_tutorial()`` saves PNGs directly to its ``output_dir``.
+2. ``TestTools.compare_result_to_baseline_image`` reads each PNG from that
+   directory and compares it against a stored baseline with loose tolerances.
+
+Run all tutorial tests::
+
+    pytest tests/test_tutorials.py --run-experiments -v
+
+Create baselines on first run::
+
+    pytest tests/test_tutorials.py --run-experiments --create-baselines -v
+"""
+
+from __future__ import annotations
+
+from pathlib import Path
+from types import SimpleNamespace
+from typing import Any, Optional
+
+import itk
+import numpy as np
+import pytest
+
+from physiomotion4d.test_tools import TestTools
+
+# Tolerances for screenshot comparison. Loose to survive minor rendering
+# differences across OS / GPU / driver versions.
+_PX_TOL = 10.0  # per-pixel absolute error (0-255 range)
+_MAX_PX = 2000  # maximum number of pixels allowed above _PX_TOL
+_TOT_TOL = float("inf")  # use the pixel-count criterion only
+
+
+class _FakeMesh:
+    """Minimal mesh double for tutorial screenshot selection tests."""
+
+    def __init__(self, n_points: int) -> None:
+        self.n_points = n_points
+
+
+class _FakeSurfaceExtractable:
+    """Minimal mesh double that records PyVista extraction options."""
+
+    def __init__(self) -> None:
+        self.algorithm: Optional[str] = None
+
+    def extract_surface(self, *, algorithm: str) -> "_FakeSurfaceExtractable":
+        self.algorithm = algorithm
+        return self
+
+
+def _compare_screenshots(
+    screenshots: list[Path],
+    tt: TestTools,
+) -> None:
+    """Read each PNG as itk.Image and compare against baseline."""
+    if not screenshots:
+        pytest.fail("No screenshots produced by run_tutorial")
+
+    for png_path in screenshots:
+        if not png_path.exists():
+            pytest.fail(f"Screenshot not created: {png_path}")
+        assert tt.compare_result_to_baseline_image(
+            png_path.name,
+            per_pixel_absolute_error_tol=_PX_TOL,
+            max_number_of_pixels_above_tol=_MAX_PX,
+            total_absolute_error_tol=_TOT_TOL,
+        ), f"Screenshot baseline mismatch: {png_path.name}"
+
+
+def test_testtools_results_output_dir_override(tmp_path: Path) -> None:
+    """Store result artifacts in an explicit directory when requested."""
+    results_root = tmp_path / "results"
+    output_dir = tmp_path / "tutorial_output"
+    baselines_root = tmp_path / "baselines"
+    tt = TestTools(
+        results_dir=results_root,
+        baselines_dir=baselines_root,
+        class_name="tutorial_example",
+        results_output_dir=output_dir,
+    )
+
+    image = itk.image_from_array(np.zeros((2, 2, 2), dtype=np.uint8))
+    tt.write_result_image(image, "current.mha")
+
+    assert (output_dir / "current.mha").exists()
+    assert not (results_root / "tutorial_example" / "current.mha").exists()
+    assert (baselines_root / "tutorial_example").exists()
+
+
+# -----------------------------------------------------------------------------
+# Tutorial 1 - Heart-Gated CT to Animated USD
+# -----------------------------------------------------------------------------
+
+
+def test_tutorial_01_contour_png_mesh_uses_current_run_results() -> None:
+    """Select current in-memory contours instead of disk VTP outputs."""
+    from tutorials.tutorial_01_heart_gated_ct_to_usd import (
+        _first_current_contour_mesh,
+    )
+
+    transformed_mesh = _FakeMesh(4)
+    reference_mesh = _FakeMesh(8)
+    workflow: Any = SimpleNamespace(
+        _transformed_contours={"all": [transformed_mesh]},
+        _reference_contours={"all": reference_mesh},
+    )
+
+    assert _first_current_contour_mesh(workflow) is transformed_mesh
+
+
+def test_tutorial_01_reference_png_uses_workflow_fixed_image(
+    tmp_path: Path,
+) -> None:
+    """Select the actual workflow reference image over cached slice images."""
+    from tutorials.tutorial_01_heart_gated_ct_to_usd import (
+        _current_reference_image,
+    )
+
+    fixed_image = object()
+    workflow: Any = SimpleNamespace(_fixed_image=fixed_image)
+    (tmp_path / "slice_000.mha").touch()
+
+    assert _current_reference_image(workflow, tmp_path) is fixed_image
+
+
+def test_tutorial_01_overlay_uses_workflow_fixed_segmentation(
+    tmp_path: Path,
+) -> None:
+    """Select the current fixed labelmap over cached slice labelmaps."""
+    from tutorials.tutorial_01_heart_gated_ct_to_usd import (
+        _current_reference_segmentation,
+    )
+
+    labelmap = object()
+    workflow: Any = SimpleNamespace(_fixed_segmentation={"labelmap": labelmap})
+    (tmp_path / "slice_000_labelmap.mha").touch()
+
+    assert _current_reference_segmentation(workflow, tmp_path) is labelmap
+
+
+def test_tutorial_01_overlay_falls_back_to_fixed_image_mask(
+    tmp_path: Path,
+) -> None:
+    """Read fixed_image_mask.mha before stale slice labelmap files."""
+    from tutorials.tutorial_01_heart_gated_ct_to_usd import (
+        _current_reference_segmentation,
+    )
+
+    arr = np.array(
+        [[[1, 0], [0, 1]], [[0, 1], [1, 0]]],
+        dtype=np.uint8,
+    )
+    mask = itk.image_from_array(arr)
+    fixed_mask_path = tmp_path / "fixed_image_mask.mha"
+    itk.imwrite(mask, str(fixed_mask_path), compression=True)
+    (tmp_path / "slice_000_labelmap.mha").touch()
+
+    workflow: Any = SimpleNamespace()
+    selected = _current_reference_segmentation(workflow, tmp_path)
+
+    assert selected is not None
+    assert tuple(selected.GetLargestPossibleRegion().GetSize()) == (2, 2, 2)
+
+
+@pytest.mark.experiment
+@pytest.mark.requires_data
+@pytest.mark.slow
+class TestTutorial01HeartGatedCTToUSD:
+    """End-to-end test for tutorial_01_heart_gated_ct_to_usd.py."""
+
+    _class_name = "tutorial_01_heart_gated_ct_to_usd"
+
+    def test_run(self, test_directories: dict[str, Path]) -> None:
+        from tutorials.tutorial_01_heart_gated_ct_to_usd import run_tutorial
+
+        out_dir = test_directories["output"] / self._class_name
+        results: dict[str, Any] = run_tutorial(
+            data_dir=test_directories["data"],
+            output_dir=out_dir,
+            registration_method="ants",
+        )
+        assert results["usd_file"], "USD file path should not be empty"
+        assert Path(results["usd_file"]).exists(), "USD file should exist"
+
+        tt = TestTools(
+            class_name=self._class_name,
+            results_dir=test_directories["output"],
+            baselines_dir=test_directories["baselines"],
+            results_output_dir=out_dir,
+        )
+        _compare_screenshots(results["screenshots"], tt)
+
+
+# -----------------------------------------------------------------------------
+# Tutorial 3 - Fit Statistical Model to Patient
+# -----------------------------------------------------------------------------
+
+
+def test_tutorial_03_extract_surface_uses_dataset_surface() -> None:
+    """Use the robust dataset_surface algorithm for VTK surface extraction."""
+    from tutorials.tutorial_03_fit_statistical_model_to_patient import (
+        _extract_surface,
+    )
+
+    mesh: Any = _FakeSurfaceExtractable()
+
+    assert _extract_surface(mesh) is mesh
+    assert mesh.algorithm == "dataset_surface"
+
+
+# -----------------------------------------------------------------------------
+# Tutorial 2 - CT Segmentation to VTK
+# -----------------------------------------------------------------------------
+
+
+@pytest.mark.experiment
+@pytest.mark.requires_data
+@pytest.mark.slow
+class TestTutorial02CTToVTK:
+    """End-to-end test for tutorial_02_ct_to_vtk.py."""
+
+    _class_name = "tutorial_02_ct_to_vtk"
+
+    def test_run(self, test_directories: dict[str, Path]) -> None:
+        from tutorials.tutorial_02_ct_to_vtk import run_tutorial
+
+        out_dir = test_directories["output"] / self._class_name
+        results: dict[str, Any] = run_tutorial(
+            data_dir=test_directories["data"],
+            output_dir=out_dir,
+        )
+        assert results["surface_file"].exists(), "Combined VTP surface should exist"
+        assert results["mesh_file"].exists(), "Combined VTU mesh should exist"
+
+        tt = TestTools(
+            class_name=self._class_name,
+            results_dir=test_directories["output"],
+            baselines_dir=test_directories["baselines"],
+            results_output_dir=out_dir,
+        )
+        _compare_screenshots(results["screenshots"], tt)
+
+
+# -----------------------------------------------------------------------------
+# Tutorial 3 - Fit Statistical Model to Patient
+# -----------------------------------------------------------------------------
+
+
+@pytest.mark.experiment
+@pytest.mark.requires_data
+@pytest.mark.slow
+class TestTutorial03FitStatisticalModelToPatient:
+    """End-to-end test for tutorial_03_fit_statistical_model_to_patient.py."""
+
+    _class_name = "tutorial_03_fit_statistical_model_to_patient"
+
+    def test_run(self, test_directories: dict[str, Path]) -> None:
+        kcl_dir = test_directories["data"] / "KCL-Heart-Model"
+        if not (kcl_dir / "pca_mean.vtu").exists():
+            pytest.skip(
+                "KCL-Heart-Model not downloaded. See data/README.md for instructions."
+            )
+
+        from tutorials.tutorial_03_fit_statistical_model_to_patient import run_tutorial
+
+        out_dir = test_directories["output"] / self._class_name
+        results: dict[str, Any] = run_tutorial(
+            data_dir=test_directories["data"],
+            output_dir=out_dir,
+        )
+        assert results["registered_file"].exists(), "Registered VTP should exist"
+
+        tt = TestTools(
+            class_name=self._class_name,
+            results_dir=test_directories["output"],
+            baselines_dir=test_directories["baselines"],
+            results_output_dir=out_dir,
+        )
+        _compare_screenshots(results["screenshots"], tt)
+
+
+# -----------------------------------------------------------------------------
+# Tutorial 4 - Create Statistical Shape Model
+# -----------------------------------------------------------------------------
+
+
+@pytest.mark.experiment
+@pytest.mark.requires_data
+@pytest.mark.slow
+class TestTutorial04CreateStatisticalModel:
+    """End-to-end test for tutorial_04_create_statistical_model.py."""
+
+    _class_name = "tutorial_04_create_statistical_model"
+
+    def test_run(self, test_directories: dict[str, Path]) -> None:
+        kcl_dir = test_directories["data"] / "KCL-Heart-Model"
+        if not (kcl_dir / "pca_mean.vtu").exists():
+            pytest.skip(
+                "KCL-Heart-Model not downloaded. See data/README.md for instructions."
+            )
+
+        from tutorials.tutorial_04_create_statistical_model import run_tutorial
+
+        out_dir = test_directories["output"] / self._class_name
+        results: dict[str, Any] = run_tutorial(
+            data_dir=test_directories["data"],
+            output_dir=out_dir,
+            pca_components=5,
+            max_samples=10,
+        )
+        assert results["model_file"].exists(), "pca_model.json should exist"
+        assert results["mean_surface_file"].exists(), "Mean surface VTP should exist"
+
+        tt = TestTools(
+            class_name=self._class_name,
+            results_dir=test_directories["output"],
+            baselines_dir=test_directories["baselines"],
+            results_output_dir=out_dir,
+        )
+        _compare_screenshots(results["screenshots"], tt)
+
+
+# -----------------------------------------------------------------------------
+# Tutorial 5 - VTK to USD
+# -----------------------------------------------------------------------------
+
+
+@pytest.mark.experiment
+@pytest.mark.requires_data
+@pytest.mark.slow
+class TestTutorial05VTKToUSD:
+    """End-to-end test for tutorial_05_vtk_to_usd.py."""
+
+    _class_name = "tutorial_05_vtk_to_usd"
+
+    def test_run(self, test_directories: dict[str, Path]) -> None:
+        # Prefer Tutorial 2 output; fall back to any .vtp in data
+        tutorial2_vtp = (
+            test_directories["output"]
+            / "tutorial_02_ct_to_vtk"
+            / "patient_surfaces.vtp"
+        )
+        vtk_file = tutorial2_vtp if tutorial2_vtp.exists() else None
+        if vtk_file is None:
+            found = list(test_directories["data"].rglob("*.vtp"))
+            if not found:
+                pytest.skip(
+                    "No VTK file available. Run Tutorial 2 first or place a .vtp "
+                    "file under data/."
+                )
+            vtk_file = found[0]
+
+        from tutorials.tutorial_05_vtk_to_usd import run_tutorial
+
+        out_dir = test_directories["output"] / self._class_name
+        results: dict[str, Any] = run_tutorial(
+            data_dir=test_directories["data"],
+            output_dir=out_dir,
+            vtk_file=vtk_file,
+        )
+        assert results["usd_file"], "USD file path should not be empty"
+        assert Path(results["usd_file"]).exists(), "USD file should exist"
+
+        tt = TestTools(
+            class_name=self._class_name,
+            results_dir=test_directories["output"],
+            baselines_dir=test_directories["baselines"],
+            results_output_dir=out_dir,
+        )
+        _compare_screenshots(results["screenshots"], tt)
+
+
+# -----------------------------------------------------------------------------
+# Tutorial 6 - Reconstruct High-Resolution 4D CT
+# -----------------------------------------------------------------------------
+
+
+@pytest.mark.experiment
+@pytest.mark.requires_data
+@pytest.mark.slow
+class TestTutorial06ReconstructHighres4DCT:
+    """End-to-end test for tutorial_06_reconstruct_highres_4d_ct.py."""
+
+    _class_name = "tutorial_06_reconstruct_highres_4d_ct"
+
+    def test_run(self, test_directories: dict[str, Path]) -> None:
+        dirlab_dir = test_directories["data"] / "DirLab-4DCT" / "Case1"
+        if not dirlab_dir.exists():
+            pytest.skip(
+                "DirLab-4DCT Case1 not downloaded. See data/README.md for instructions."
+            )
+
+        from tutorials.tutorial_06_reconstruct_highres_4d_ct import run_tutorial
+
+        out_dir = test_directories["output"] / self._class_name
+        results: dict[str, Any] = run_tutorial(
+            data_dir=test_directories["data"],
+            output_dir=out_dir,
+            case=1,
+            max_frames=3,
+            registration_method="ants",
+        )
+        assert results["reconstructed_files"], (
+            "At least one reconstructed frame expected"
+        )
+        for f in results["reconstructed_files"]:
+            assert f.exists(), f"Reconstructed frame missing: {f}"
+
+        tt = TestTools(
+            class_name=self._class_name,
+            results_dir=test_directories["output"],
+            baselines_dir=test_directories["baselines"],
+            results_output_dir=out_dir,
+        )
+        _compare_screenshots(results["screenshots"], tt)
diff --git a/tutorials/README.md b/tutorials/README.md
new file mode 100644
index 0000000..fab80df
--- /dev/null
+++ b/tutorials/README.md
@@ -0,0 +1,61 @@
+# PhysioMotion4D Tutorials
+
+End-to-end Python scripts covering each major workflow in the library.
+These are the recommended starting point for new users.
+
+## Before You Begin
+
+Each tutorial requires one or more public datasets.
+**See [../data/README.md](../data/README.md)** for download instructions,
+dataset licensing, and expected directory layout.
+
+## Tutorial Index
+
+| # | Script | Workflow Class | CLI Command | Dataset |
+|---|--------|---------------|-------------|---------|
+| 1 | [tutorial_01_heart_gated_ct_to_usd.py](tutorial_01_heart_gated_ct_to_usd.py) | `WorkflowConvertHeartGatedCTToUSD` | `physiomotion4d-heart-gated-ct` | Slicer-Heart-CT (prepare first) |
+| 2 | [tutorial_02_ct_to_vtk.py](tutorial_02_ct_to_vtk.py) | `WorkflowConvertCTToVTK` | `physiomotion4d-convert-ct-to-vtk` | Slicer-Heart-CT (prepare first) |
+| 3 | [tutorial_03_fit_statistical_model_to_patient.py](tutorial_03_fit_statistical_model_to_patient.py) | `WorkflowFitStatisticalModelToPatient` | `physiomotion4d-fit-statistical-model-to-patient` | KCL-Heart-Model (manual) |
+| 4 | [tutorial_04_create_statistical_model.py](tutorial_04_create_statistical_model.py) | `WorkflowCreateStatisticalModel` | `physiomotion4d-create-statistical-model` | KCL-Heart-Model (manual) |
+| 5 | [tutorial_05_vtk_to_usd.py](tutorial_05_vtk_to_usd.py) | `WorkflowConvertVTKToUSD` | `physiomotion4d-convert-vtk-to-usd` | Output of tutorial 2 |
+| 6 | [tutorial_06_reconstruct_highres_4d_ct.py](tutorial_06_reconstruct_highres_4d_ct.py) | `WorkflowReconstructHighres4DCT` | `physiomotion4d-reconstruct-highres-4d-ct` | DirLab-4DCT (manual) |
+
+## Running a Tutorial
+
+Each tutorial is a standalone Python script with a `run_tutorial(data_dir, output_dir)`
+function. Run from the repository root:
+
+```bash
+python tutorials/tutorial_01_heart_gated_ct_to_usd.py \
+    --data-dir ./data --output-dir ./output
+
+python tutorials/tutorial_02_ct_to_vtk.py \
+    --data-dir ./data --output-dir ./output
+```
+
+## Running as Pytest Experiment Tests
+
+All tutorials are wired into the test suite under the `experiment` marker.
+They run end-to-end and compare generated screenshots against baselines:
+
+```bash
+# Run all tutorial tests (requires data download first)
+pytest tests/test_tutorials.py --run-experiments -v
+
+# Create baselines on first run
+pytest tests/test_tutorials.py --run-experiments --create-baselines -v
+
+# Run a single tutorial test
+pytest tests/test_tutorials.py::TestTutorial01HeartGatedCTToUSD --run-experiments -v
+```
+
+## Recommended Order
+
+1. **Tutorial 1** and **Tutorial 2** use Slicer-Heart-CT - prepare it per `data/README.md`, then start here.
+2. **Tutorial 5** uses the VTK surfaces produced by Tutorial 2 - run Tutorial 2 first.
+3. **Tutorials 3 and 4** require the KCL-Heart-Model - download it per `data/README.md`.
+4. **Tutorial 6** requires DirLab-4DCT - download it per `data/README.md`.
+
+## For Contributors
+
+Class-level API reference: [../docs/API_MAP.md](../docs/API_MAP.md)
diff --git a/tutorials/__init__.py b/tutorials/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tutorials/tutorial_01_heart_gated_ct_to_usd.py b/tutorials/tutorial_01_heart_gated_ct_to_usd.py
new file mode 100644
index 0000000..143f262
--- /dev/null
+++ b/tutorials/tutorial_01_heart_gated_ct_to_usd.py
@@ -0,0 +1,297 @@
+"""
+Tutorial 1: Heart-Gated CT to Animated USD
+
+Purpose
+-------
+Convert a 4D cardiac CT scan (multiple gated time frames) into an animated USD
+model suitable for visualization in NVIDIA Omniverse. The workflow segments the
+heart and surrounding anatomy from a reference frame, registers all other frames
+to that reference using deep learning or classical registration, and assembles
+the resulting time-varying surface meshes into a single USD file with anatomical
+materials applied.
+
+Inputs
+------
+- A 4D NRRD sequence file (``*.seq.nrrd``) **or** a list of 3D CT volumes
+  (``*.mha`` / ``*.nrrd``) representing successive cardiac phases.
+  Expected location: ``data/Slicer-Heart-CT/TruncalValve_4DCT.seq.nrrd``
+- Optional: a reference frame image to fix the cardiac phase used as the
+  segmentation source.
+
+Outputs
+-------
+- ``output_dir/cardiac_model.dynamic_anatomy_painted.usd`` - animated USD with
+  anatomy materials
+- Screenshots (PNG) for documentation and regression testing:
+  - ``reference_frame_axial.png`` - axial slice of the reference CT frame
+  - ``segmentation_overlay.png`` - segmentation mask overlaid on reference
+  - ``contours_3d.png`` - 3-D isometric view of the current-run contours
+
+Strengths
+---------
+- Single call (``WorkflowConvertHeartGatedCTToUSD.process()``) runs the full pipeline.
+- Supports both GPU-accelerated ICON registration and CPU-capable ANTs registration.
+- Automatically detects contrast enhancement and adjusts segmentation thresholds.
+- Output is Omniverse-ready with anatomical materials (USDAnatomyTools).
+
+Weaknesses / Limitations
+------------------------
+- Requires a GPU for ICON registration (``registration_method='icon'``); use
+  ``registration_method='ants'`` for CPU-only environments (about 10x slower).
+- Segmentation quality depends on TotalSegmentator's training distribution;
+  unusual pathologies or pediatric anatomy may degrade results.
+- Large 4D datasets (>20 phases, high resolution) can require 32 GB+ RAM.
+
+Classes Used
+------------
+- WorkflowConvertHeartGatedCTToUSD (workflow_convert_heart_gated_ct_to_usd.py):
+    Orchestrates the full pipeline: 4D NRRD -> segmentation -> registration ->
+    contour extraction -> USD export.
+- SegmentChestTotalSegmentator (segment_chest_total_segmentator.py):
+    Deep-learning segmentation of 117 anatomical structures (used internally).
+- RegisterImagesICON / RegisterImagesANTs (register_images_icon.py / _ants.py):
+    Frame-to-frame image registration (used internally).
+- ContourTools (contour_tools.py):
+    Extracts and transforms surface meshes from segmentation masks (used internally).
+- USDAnatomyTools (usd_anatomy_tools.py):
+    Applies clinical material colours to USD prims (used internally).
+
+CLI Equivalent
+--------------
+The same main outputs (without screenshots) can be produced via the CLI::
+
+    physiomotion4d-heart-gated-ct \\
+        data/Slicer-Heart-CT/TruncalValve_4DCT.seq.nrrd \\
+        --contrast \\
+        --project-name cardiac_model \\
+        --registration-method ants \\
+        --registration-iterations 1 \\
+        --output-dir ./output/tutorial_01
+
+See ``src/physiomotion4d/cli/convert_heart_gated_ct_to_usd.py`` for full CLI
+documentation.
+
+Data Required
+-------------
+See data/README.md for download instructions and dataset licensing.
+Dataset: Slicer-Heart-CT - https://github.com/Slicer-Heart-CT/Slicer-Heart-CT
+This script expects the data to already exist at
+``data/Slicer-Heart-CT/TruncalValve_4DCT.seq.nrrd``. Run the repository data
+download notebook or download the file manually before running this tutorial.
+"""
+
+from __future__ import annotations
+
+import argparse
+import logging
+from pathlib import Path
+from typing import Any, Optional
+
+import itk
+
+from physiomotion4d.test_tools import TestTools
+from physiomotion4d.workflow_convert_heart_gated_ct_to_usd import (
+    WorkflowConvertHeartGatedCTToUSD,
+)
+
+
+def _first_current_contour_mesh(
+    workflow: WorkflowConvertHeartGatedCTToUSD,
+) -> Optional[Any]:
+    """Return a contour mesh produced by the current workflow run.
+
+    Prefer transformed all-anatomy contours, because those are the meshes passed
+    into USD conversion for the current run. Fall back to reference contours only
+    if transformed contours are unavailable.
+    """
+    transformed_contours = getattr(workflow, "_transformed_contours", {})
+    if isinstance(transformed_contours, dict):
+        for mesh in transformed_contours.get("all", []):
+            if getattr(mesh, "n_points", 0) > 0:
+                return mesh
+
+    reference_contours = getattr(workflow, "_reference_contours", {})
+    if isinstance(reference_contours, dict):
+        mesh = reference_contours.get("all")
+        if mesh is not None and getattr(mesh, "n_points", 0) > 0:
+            return mesh
+
+    return None
+
+
+def _current_reference_image(
+    workflow: WorkflowConvertHeartGatedCTToUSD,
+    output_dir: Path,
+) -> Optional[Any]:
+    """Return the reference image used by the current workflow run."""
+    fixed_image = getattr(workflow, "_fixed_image", None)
+    if fixed_image is not None:
+        return fixed_image
+
+    fixed_image_file = output_dir / "fixed_image.mha"
+    if fixed_image_file.exists():
+        return itk.imread(str(fixed_image_file))
+
+    ref_frames = sorted(output_dir.glob("slice_???.mha"))
+    if ref_frames:
+        return itk.imread(str(ref_frames[0]))
+
+    return None
+
+
+def _current_reference_segmentation(
+    workflow: WorkflowConvertHeartGatedCTToUSD,
+    output_dir: Path,
+) -> Optional[Any]:
+    """Return the labelmap for the current workflow reference image."""
+    fixed_segmentation = getattr(workflow, "_fixed_segmentation", None)
+    if isinstance(fixed_segmentation, dict):
+        labelmap = fixed_segmentation.get("labelmap")
+        if labelmap is not None:
+            return labelmap
+
+    fixed_mask_file = output_dir / "fixed_image_mask.mha"
+    if fixed_mask_file.exists():
+        return itk.imread(str(fixed_mask_file))
+
+    label_files = sorted(output_dir.glob("slice_???_labelmap*.mha"))
+    if label_files:
+        return itk.imread(str(label_files[0]))
+
+    return None
+
+
+def run_tutorial(
+    data_dir: Path,
+    output_dir: Path,
+    *,
+    registration_method: str = "ants",
+    log_level: int = logging.INFO,
+) -> dict[str, Any]:
+    """Run Tutorial 1: Heart-Gated CT to Animated USD.
+
+    Args:
+        data_dir: Root of the ``data/`` directory (see data/README.md).
+        output_dir: Directory to write outputs and screenshots.
+        registration_method: ``'ants'`` (CPU-capable, default) or ``'icon'`` (GPU).
+        log_level: Python logging level.
+
+    Returns:
+        dict with keys:
+
+        - ``'usd_file'`` (str): path to the final painted USD.
+        - ``'screenshots'`` (list[Path]): paths to saved PNG screenshots.
+          PNGs are rendered from data produced by this invocation, not from
+          previously saved VTK/VTP files in ``output_dir``. Reference-frame
+          screenshots use the workflow's selected fixed image.
+    """
+    output_dir.mkdir(parents=True, exist_ok=True)
+
+    nrrd_file = data_dir / "Slicer-Heart-CT" / "TruncalValve_4DCT.seq.nrrd"
+    if not nrrd_file.exists():
+        raise FileNotFoundError(
+            f"Slicer-Heart-CT data not found: {nrrd_file}\n"
+            "See data/README.md for download instructions."
+        )
+
+    workflow = WorkflowConvertHeartGatedCTToUSD(
+        input_filenames=[str(nrrd_file)],
+        contrast_enhanced=True,
+        output_directory=str(output_dir),
+        project_name="cardiac_model",
+        registration_method=registration_method,
+        number_of_registration_iterations=1,
+        log_level=log_level,
+    )
+
+    usd_file = output_dir / workflow.process()
+
+    # Screenshots
+    tt = TestTools(
+        results_dir=output_dir,
+        baselines_dir=output_dir / "baselines",
+        class_name="tutorial_01",
+        results_output_dir=output_dir,
+        log_level=log_level,
+    )
+
+    screenshots: list[Path] = []
+
+    # Reference frame: use the workflow's selected fixed image for this run.
+    ref_image = _current_reference_image(workflow, output_dir)
+    if ref_image is not None:
+        screenshots.append(
+            tt.save_screenshot_image_slice(
+                ref_image,
+                "reference_frame_axial.png",
+                axis=0,
+                slice_fraction=0.5,
+                colormap="gray",
+                vmin=-200,
+                vmax=600,
+            )
+        )
+
+        # Segmentation overlay: align with the selected fixed image.
+        overlay = _current_reference_segmentation(workflow, output_dir)
+        screenshots.append(
+            tt.save_screenshot_image_slice(
+                ref_image,
+                "segmentation_overlay.png",
+                axis=0,
+                slice_fraction=0.5,
+                colormap="gray",
+                vmin=-200,
+                vmax=600,
+                overlay_mask=overlay,
+            )
+        )
+
+    # 3-D contour view: render the current run's in-memory contours.
+    contour_mesh = _first_current_contour_mesh(workflow)
+    if contour_mesh is not None:
+        screenshots.append(
+            tt.save_screenshot_mesh(
+                contour_mesh,
+                "contours_3d.png",
+                camera_position="iso",
+                color="tomato",
+                opacity=0.85,
+            )
+        )
+
+    return {"usd_file": str(usd_file), "screenshots": screenshots}
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(
+        description=__doc__,
+        formatter_class=argparse.RawDescriptionHelpFormatter,
+    )
+    parser.add_argument(
+        "--data-dir",
+        type=Path,
+        default=Path("data"),
+        help="Root data directory (default: ./data)",
+    )
+    parser.add_argument(
+        "--output-dir",
+        type=Path,
+        default=Path("output") / "tutorial_01",
+        help="Output directory (default: ./output/tutorial_01)",
+    )
+    parser.add_argument(
+        "--registration-method",
+        default="ants",
+        choices=["ants", "icon"],
+        help="Registration method: ants (CPU) or icon (GPU). Default: ants",
+    )
+    args = parser.parse_args()
+
+    results = run_tutorial(
+        args.data_dir,
+        args.output_dir,
+        registration_method=args.registration_method,
+    )
+    print(f"USD file: {results['usd_file']}")
+    print(f"Screenshots: {[str(p) for p in results['screenshots']]}")
diff --git a/tutorials/tutorial_02_ct_to_vtk.py b/tutorials/tutorial_02_ct_to_vtk.py
new file mode 100644
index 0000000..bff9541
--- /dev/null
+++ b/tutorials/tutorial_02_ct_to_vtk.py
@@ -0,0 +1,216 @@
+"""
+Tutorial 2: CT Segmentation to VTK Surfaces
+
+Purpose
+-------
+Segment a 3D CT image into anatomical groups (heart, lungs, vessels, bone,
+soft tissue) and export per-group VTK surface and mesh files. Each output mesh
+is annotated with anatomy metadata and colour so it can be used directly in
+PyVista, ParaView, or the downstream USD pipeline (Tutorial 5).
+
+Inputs
+------
+- A single 3D CT image in any ITK-readable format (NIfTI, MHA, NRRD, etc.).
+  This tutorial uses one time frame from the Slicer-Heart-CT dataset.
+  Expected location: ``data/Slicer-Heart-CT/`` (any ``slice_???.mha`` frame).
+
+Outputs
+-------
+- ``output_dir/patient_surfaces.vtp`` - all anatomy surfaces in one file
+- ``output_dir/patient_meshes.vtu`` - all voxel meshes in one file
+- Screenshots (PNG):
+  - ``segmentation_overlay.png`` - segmentation mask overlaid on axial CT slice
+  - ``vtk_surfaces.png`` - 3-D isometric view of the combined surface
+
+Strengths
+---------
+- One call to ``WorkflowConvertCTToVTK.run_workflow()`` handles segmentation and
+  mesh extraction for all anatomy groups in a single pass.
+- Each output mesh carries field data (group name, label IDs, colour) for
+  downstream tools.
+- Combined-file output (default) produces a single VTP/VTU rather than one file
+  per group, simplifying downstream handling.
+
+Weaknesses / Limitations
+------------------------
+- TotalSegmentator requires ~8 GB GPU VRAM for full segmentation; CPU fallback
+  is available but much slower (~30 min per volume).
+- Small or unusual anatomical structures (e.g., pediatric heart) may be partially
+  missed by the default TotalSegmentator model.
+- Output mesh resolution is governed by the input CT voxel size; coarse scans
+  yield coarser meshes.
+
+Classes Used
+------------
+- WorkflowConvertCTToVTK (workflow_convert_ct_to_vtk.py):
+    Segments a CT image and extracts per-anatomy-group VTK surfaces and meshes.
+- SegmentChestTotalSegmentator (segment_chest_total_segmentator.py):
+    Deep-learning segmentation backend (used internally).
+- ContourTools (contour_tools.py):
+    Mesh extraction via marching cubes (used internally).
+- USDAnatomyTools (usd_anatomy_tools.py):
+    Provides anatomy group colours for mesh annotation (used internally).
+
+CLI Equivalent
+--------------
+The same main outputs (without screenshots) can be produced via the CLI::
+
+    physiomotion4d-convert-ct-to-vtk \\
+        --input-image data/Slicer-Heart-CT/slice_000.mha \\
+        --output-dir ./output/tutorial_02 \\
+        --output-prefix patient \\
+        --contrast
+
+See ``src/physiomotion4d/cli/convert_ct_to_vtk.py`` for full CLI documentation.
+
+Data Required
+-------------
+See data/README.md for download instructions and dataset licensing.
+Dataset: Slicer-Heart-CT - https://github.com/Slicer-Heart-CT/Slicer-Heart-CT
+Auto-download: the conftest fixture downloads
+``data/test/TruncalValve_4DCT.seq.nrrd`` and extracts frames as
+``data/test/slice_???.mha``. For this tutorial the full dataset is at
+``data/Slicer-Heart-CT/``.
+"""
+
+from __future__ import annotations
+
+import argparse
+import logging
+from pathlib import Path
+from typing import Any
+
+import itk
+import pyvista as pv
+
+from physiomotion4d.test_tools import TestTools
+from physiomotion4d.workflow_convert_ct_to_vtk import WorkflowConvertCTToVTK
+
+
+def run_tutorial(
+    data_dir: Path,
+    output_dir: Path,
+    *,
+    log_level: int = logging.INFO,
+) -> dict[str, Any]:
+    """Run Tutorial 2: CT Segmentation to VTK Surfaces.
+
+    Args:
+        data_dir: Root of the ``data/`` directory (see data/README.md).
+        output_dir: Directory to write outputs and screenshots.
+        log_level: Python logging level.
+
+    Returns:
+        dict with keys:
+
+        - ``'result'`` (dict): workflow result dict with ``'surfaces'`` and
+          ``'meshes'`` entries (pv.PolyData / pv.UnstructuredGrid per group).
+        - ``'surface_file'`` (Path): path to the combined ``.vtp`` file.
+        - ``'mesh_file'`` (Path): path to the combined ``.vtu`` file.
+        - ``'screenshots'`` (list[Path]): paths to saved PNG screenshots.
+    """
+    output_dir.mkdir(parents=True, exist_ok=True)
+
+    # Prefer full dataset; fall back to test cache
+    candidates = list((data_dir / "Slicer-Heart-CT").glob("slice_???.mha"))
+    if not candidates:
+        candidates = list((data_dir / "test").glob("slice_???_sml.mha"))
+    if not candidates:
+        raise FileNotFoundError(
+            "No CT frame found under data/Slicer-Heart-CT/ or data/test/.\n"
+            "See data/README.md for download instructions."
+        )
+    candidates.sort()
+    ct_file = candidates[0]
+
+    ct_image = itk.imread(str(ct_file))
+
+    workflow = WorkflowConvertCTToVTK(
+        segmentation_method="total_segmentator",
+        log_level=log_level,
+    )
+    result = workflow.run_workflow(
+        input_image=ct_image,
+        contrast_enhanced_study=True,
+    )
+
+    surface_file = output_dir / "patient_surfaces.vtp"
+    mesh_file = output_dir / "patient_meshes.vtu"
+    WorkflowConvertCTToVTK.save_combined_surface(
+        result["surfaces"], str(output_dir), prefix="patient"
+    )
+    WorkflowConvertCTToVTK.save_combined_mesh(
+        result["meshes"], str(output_dir), prefix="patient"
+    )
+
+    # Screenshots
+    tt = TestTools(
+        results_dir=output_dir,
+        baselines_dir=output_dir / "baselines",
+        class_name="tutorial_02",
+        results_output_dir=output_dir,
+        log_level=log_level,
+    )
+
+    screenshots: list[Path] = []
+
+    # Segmentation overlay on axial CT slice
+    labelmap = result.get("labelmap")
+    screenshots.append(
+        tt.save_screenshot_image_slice(
+            ct_image,
+            "segmentation_overlay.png",
+            axis=0,
+            slice_fraction=0.5,
+            colormap="gray",
+            vmin=-200,
+            vmax=600,
+            overlay_mask=labelmap,
+        )
+    )
+
+    # 3-D view of combined surface
+    surfaces = [s for s in result["surfaces"].values() if s is not None]
+    if surfaces:
+        combined = pv.merge(surfaces) if len(surfaces) > 1 else surfaces[0]
+        screenshots.append(
+            tt.save_screenshot_mesh(
+                combined,
+                "vtk_surfaces.png",
+                camera_position="iso",
+                color="lightblue",
+                opacity=0.85,
+            )
+        )
+
+    return {
+        "result": result,
+        "surface_file": surface_file,
+        "mesh_file": mesh_file,
+        "screenshots": screenshots,
+    }
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(
+        description=__doc__,
+        formatter_class=argparse.RawDescriptionHelpFormatter,
+    )
+    parser.add_argument(
+        "--data-dir",
+        type=Path,
+        default=Path("data"),
+        help="Root data directory (default: ./data)",
+    )
+    parser.add_argument(
+        "--output-dir",
+        type=Path,
+        default=Path("output") / "tutorial_02",
+        help="Output directory (default: ./output/tutorial_02)",
+    )
+    args = parser.parse_args()
+
+    results = run_tutorial(args.data_dir, args.output_dir)
+    print(f"Surface file: {results['surface_file']}")
+    print(f"Mesh file:    {results['mesh_file']}")
+    print(f"Screenshots:  {[str(p) for p in results['screenshots']]}")
diff --git a/tutorials/tutorial_03_fit_statistical_model_to_patient.py b/tutorials/tutorial_03_fit_statistical_model_to_patient.py
new file mode 100644
index 0000000..1a20c68
--- /dev/null
+++ b/tutorials/tutorial_03_fit_statistical_model_to_patient.py
@@ -0,0 +1,233 @@
+"""
+Tutorial 3: Fit Statistical Shape Model to Patient Data
+
+Purpose
+-------
+Register a generic anatomical template (e.g., a statistical heart model) to
+patient-specific surface meshes derived from medical imaging. The multi-stage
+pipeline performs ICP rough alignment, optional PCA-constrained shape fitting,
+mask-based deformable registration, and optional image-based refinement. The
+result is a patient-specific instance of the template mesh that matches the
+patient anatomy.
+
+Inputs
+------
+- Template model (``pv.UnstructuredGrid`` / ``.vtu``): generic anatomical mesh,
+  e.g., the KCL heart model.
+  Expected location: ``data/KCL-Heart-Model/pca_mean.vtu``
+- Patient surface models (list of ``pv.PolyData`` / ``.vtp``): anatomy surfaces
+  extracted from the patient CT (e.g., from Tutorial 2).
+  Expected location: ``data/KCL-Heart-Model/sample_meshes/*.vtu`` (used as
+  stand-in patient models for demonstration).
+- Optional patient CT image (``itk.Image``): used for image-based refinement.
+
+Outputs
+-------
+- ``output_dir/registered_template.vtp`` - template mesh fitted to patient
+- Screenshots (PNG):
+  - ``model_before_registration.png`` - template and patient overlaid (pre-ICP)
+  - ``model_after_registration.png`` - registered template on patient
+
+Strengths
+---------
+- Combines three complementary registration strategies in sequence, each
+  correcting different scales of misalignment.
+- PCA-constrained fitting (optional) prevents anatomically implausible
+  deformations by constraining shape variation to the training population.
+- Automatic mask generation means patient meshes are the only required input;
+  a CT image is optional.
+
+Weaknesses / Limitations
+------------------------
+- Requires the KCL-Heart-Model dataset (manual download; see data/README.md).
+- Deformable registration (ANTs) is the slowest stage (~5-15 min on CPU).
+- PCA mode is only beneficial when the template was trained on a population
+  that includes the patient's anatomical variant.
+- ICON-based image refinement requires a GPU.
+
+Classes Used
+------------
+- WorkflowFitStatisticalModelToPatient (workflow_fit_statistical_model_to_patient.py):
+    Orchestrates ICP -> (optional PCA) -> mask-to-mask -> (optional image) pipeline.
+- RegisterModelsICP (register_models_icp.py):
+    Centroid alignment followed by ICP affine registration (used internally).
+- RegisterModelsDistanceMaps (register_models_distance_maps.py):
+    ANTs deformable registration via signed distance maps (used internally).
+- ContourTools (contour_tools.py):
+    Creates reference images and masks from meshes (used internally).
+
+CLI Equivalent
+--------------
+The same main outputs (without screenshots) can be produced via the CLI::
+
+    physiomotion4d-fit-statistical-model-to-patient \\
+        --template-model data/KCL-Heart-Model/pca_mean.vtu \\
+        --patient-models data/KCL-Heart-Model/sample_meshes/sample_000.vtu \\
+        --output-dir ./output/tutorial_03
+
+See ``src/physiomotion4d/cli/fit_statistical_model_to_patient.py`` for full
+CLI documentation.
+
+Data Required
+-------------
+See data/README.md for download instructions and dataset licensing.
+Dataset: KCL-Heart-Model - manual download required.
+Place files under ``data/KCL-Heart-Model/`` as described in data/README.md.
+"""
+
+from __future__ import annotations
+
+import argparse
+import logging
+from pathlib import Path
+from typing import Any, cast
+
+import pyvista as pv
+
+from physiomotion4d.workflow_fit_statistical_model_to_patient import (
+    WorkflowFitStatisticalModelToPatient,
+)
+
+
+def _extract_surface(mesh: pv.DataSet) -> pv.PolyData:
+    """Extract a PolyData surface using the library's standard VTK algorithm."""
+    if isinstance(mesh, pv.PolyData):
+        return mesh
+    return cast(pv.PolyData, mesh.extract_surface(algorithm="dataset_surface"))
+
+
+def run_tutorial(
+    data_dir: Path,
+    output_dir: Path,
+    *,
+    log_level: int = logging.INFO,
+) -> dict[str, Any]:
+    """Run Tutorial 3: Fit Statistical Shape Model to Patient Data.
+
+    Args:
+        data_dir: Root of the ``data/`` directory (see data/README.md).
+        output_dir: Directory to write outputs and screenshots.
+        log_level: Python logging level.
+
+    Returns:
+        dict with keys:
+
+        - ``'registered_model'`` (pv.PolyData): fitted template surface.
+        - ``'registered_file'`` (Path): path to saved ``.vtp``.
+        - ``'screenshots'`` (list[Path]): paths to saved PNG screenshots.
+    """
+    output_dir.mkdir(parents=True, exist_ok=True)
+
+    kcl_dir = data_dir / "KCL-Heart-Model"
+    template_file = kcl_dir / "pca_mean.vtu"
+    if not template_file.exists():
+        raise FileNotFoundError(
+            f"KCL-Heart-Model template not found: {template_file}\n"
+            "See data/README.md for manual download instructions."
+        )
+
+    template_model = pv.read(str(template_file))
+    template_model = _extract_surface(template_model)
+
+    # Use a subset of sample meshes as stand-in patient models
+    sample_files = sorted((kcl_dir / "sample_meshes").glob("*.vtu"))[:3]
+    if not sample_files:
+        sample_files = sorted(kcl_dir.glob("*.vtu"))[:3]
+    if not sample_files:
+        raise FileNotFoundError(
+            f"No sample meshes found under {kcl_dir}.\n"
+            "See data/README.md for manual download instructions."
+        )
+    patient_models = []
+    for sample_file in sample_files:
+        sample_model = pv.read(str(sample_file))
+        patient_models.append(_extract_surface(sample_model))
+
+    workflow = WorkflowFitStatisticalModelToPatient(
+        template_model=template_model,
+        patient_models=patient_models,
+        log_level=log_level,
+    )
+    result = workflow.run_workflow()
+
+    registered_surface: pv.PolyData = result["registered_template_model_surface"]
+    registered_file = output_dir / "registered_template.vtp"
+    registered_surface.save(str(registered_file))
+
+    screenshots: list[Path] = []
+
+    patient_combined = (
+        pv.merge(patient_models) if len(patient_models) > 1 else patient_models[0]
+    )
+
+    # Before: template (blue) + patient (red)
+    try:
+        pv.start_xvfb()
+    except Exception:
+        pass
+    plotter = pv.Plotter(off_screen=True, window_size=[800, 600])
+    plotter.add_mesh(
+        template_model,
+        color="dodgerblue",
+        opacity=0.6,
+        label="Template",
+    )
+    plotter.add_mesh(
+        _extract_surface(patient_combined),
+        color="tomato",
+        opacity=0.6,
+        label="Patient",
+    )
+    plotter.camera_position = "iso"
+    before_path = output_dir / "model_before_registration.png"
+    before_path.parent.mkdir(parents=True, exist_ok=True)
+    plotter.screenshot(str(before_path))
+    plotter.close()
+    screenshots.append(before_path)
+
+    # After: registered template (green) + patient (red)
+    plotter2 = pv.Plotter(off_screen=True, window_size=[800, 600])
+    plotter2.add_mesh(
+        registered_surface, color="limegreen", opacity=0.7, label="Registered"
+    )
+    plotter2.add_mesh(
+        _extract_surface(patient_combined),
+        color="tomato",
+        opacity=0.4,
+        label="Patient",
+    )
+    plotter2.camera_position = "iso"
+    after_path = output_dir / "model_after_registration.png"
+    plotter2.screenshot(str(after_path))
+    plotter2.close()
+    screenshots.append(after_path)
+
+    return {
+        "registered_model": registered_surface,
+        "registered_file": registered_file,
+        "screenshots": screenshots,
+    }
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(
+        description=__doc__,
+        formatter_class=argparse.RawDescriptionHelpFormatter,
+    )
+    parser.add_argument(
+        "--data-dir",
+        type=Path,
+        default=Path("data"),
+        help="Root data directory (default: ./data)",
+    )
+    parser.add_argument(
+        "--output-dir",
+        type=Path,
+        default=Path("output") / "tutorial_03",
+        help="Output directory (default: ./output/tutorial_03)",
+    )
+    args = parser.parse_args()
+
+    results = run_tutorial(args.data_dir, args.output_dir)
+    print(f"Registered model: {results['registered_file']}")
+    print(f"Screenshots:      {[str(p) for p in results['screenshots']]}")
diff --git a/tutorials/tutorial_04_create_statistical_model.py b/tutorials/tutorial_04_create_statistical_model.py
new file mode 100644
index 0000000..9002b93
--- /dev/null
+++ b/tutorials/tutorial_04_create_statistical_model.py
@@ -0,0 +1,275 @@
+"""
+Tutorial 4: Create a PCA Statistical Shape Model
+
+Purpose
+-------
+Build a PCA (Principal Component Analysis) statistical shape model from a
+population of anatomical meshes aligned to a reference. The model captures
+the mean shape and the principal modes of geometric variation across the
+population. The resulting model can be used in Tutorial 3 to constrain
+patient-specific fitting to anatomically plausible shapes.
+
+Inputs
+------
+- A reference mesh (``pv.DataSet`` / ``.vtu``): defines the template topology.
+  Expected location: ``data/KCL-Heart-Model/pca_mean.vtu``
+- A collection of sample meshes (list of ``pv.DataSet`` / ``.vtu``):
+  population shapes to learn from.
+  Expected location: ``data/KCL-Heart-Model/sample_meshes/*.vtu``
+
+Outputs
+-------
+- ``output_dir/pca_model.json`` - PCA model (eigenvectors, eigenvalues, mean)
+- ``output_dir/pca_mean_surface.vtp`` - mean shape as a surface
+- Screenshots (PNG):
+  - ``pca_mean_model.png`` - 3-D view of the PCA mean surface
+  - ``pca_mode_01.png`` - mean +/- 2 sigma for the first PCA mode (side-by-side)
+  - ``pca_mode_02.png`` - mean +/- 2 sigma for the second PCA mode
+
+Strengths
+---------
+- Single call to ``WorkflowCreateStatisticalModel.run_workflow()`` covers the
+  full pipeline: ICP alignment, deformable correspondence, and PCA.
+- Returns a pure-Python dict (eigenvectors as numpy arrays) compatible with
+  ``WorkflowFitStatisticalModelToPatient.set_use_pca_registration()``.
+- Surface-mode PCA (``solve_for_surface_pca=True``, default) is faster and
+  sufficient for most cardiac applications.
+
+Weaknesses / Limitations
+------------------------
+- Requires the KCL-Heart-Model dataset (manual download; see data/README.md).
+- Population size directly affects model quality; small populations (<20 meshes)
+  produce unreliable high-order modes.
+- ICP alignment (step 2) assumes all sample meshes share a common approximate
+  orientation; large pose variation may degrade correspondence.
+- Deformable correspondence (step 3) uses ANTs, which is slow on CPU.
+
+Classes Used
+------------
+- WorkflowCreateStatisticalModel (workflow_create_statistical_model.py):
+    Runs the full pipeline: ICP -> deformable correspondence -> PCA.
+- RegisterModelsICP (register_models_icp.py):
+    Aligns each sample to the reference (used internally).
+- RegisterModelsDistanceMaps (register_models_distance_maps.py):
+    Dense deformable correspondence via signed distance maps (used internally).
+
+CLI Equivalent
+--------------
+The same main outputs (without screenshots) can be produced via the CLI::
+
+    physiomotion4d-create-statistical-model \\
+        --sample-meshes-dir data/KCL-Heart-Model/sample_meshes \\
+        --reference-mesh data/KCL-Heart-Model/pca_mean.vtu \\
+        --pca-components 10 \\
+        --output-dir ./output/tutorial_04
+
+See ``src/physiomotion4d/cli/create_statistical_model.py`` for full CLI
+documentation.
+
+Data Required
+-------------
+See data/README.md for download instructions and dataset licensing.
+Dataset: KCL-Heart-Model - manual download required.
+Place files under ``data/KCL-Heart-Model/`` as described in data/README.md.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import logging
+from pathlib import Path
+from typing import Any, cast
+
+import numpy as np
+import pyvista as pv
+
+from physiomotion4d.test_tools import TestTools
+from physiomotion4d.workflow_create_statistical_model import (
+    WorkflowCreateStatisticalModel,
+)
+
+
+def run_tutorial(
+    data_dir: Path,
+    output_dir: Path,
+    *,
+    pca_components: int = 10,
+    max_samples: int = 20,
+    log_level: int = logging.INFO,
+) -> dict[str, Any]:
+    """Run Tutorial 4: Create a PCA Statistical Shape Model.
+
+    Args:
+        data_dir: Root of the ``data/`` directory (see data/README.md).
+        output_dir: Directory to write outputs and screenshots.
+        pca_components: Number of PCA modes to retain.
+        max_samples: Maximum number of sample meshes to use (cap for speed).
+        log_level: Python logging level.
+
+    Returns:
+        dict with keys:
+
+        - ``'pca_model'`` (dict): PCA model dict (eigenvectors, eigenvalues, mean).
+        - ``'mean_surface'`` (pv.PolyData): mean shape surface.
+        - ``'model_file'`` (Path): path to saved ``pca_model.json``.
+        - ``'mean_surface_file'`` (Path): path to saved ``pca_mean_surface.vtp``.
+        - ``'screenshots'`` (list[Path]): paths to saved PNG screenshots.
+    """
+    output_dir.mkdir(parents=True, exist_ok=True)
+
+    kcl_dir = data_dir / "KCL-Heart-Model"
+    reference_file = kcl_dir / "pca_mean.vtu"
+    if not reference_file.exists():
+        raise FileNotFoundError(
+            f"KCL-Heart-Model reference mesh not found: {reference_file}\n"
+            "See data/README.md for manual download instructions."
+        )
+
+    sample_dir = kcl_dir / "sample_meshes"
+    sample_files = sorted(sample_dir.glob("*.vtu"))
+    if not sample_files:
+        sample_files = sorted(kcl_dir.glob("*.vtu"))
+    sample_files = [f for f in sample_files if f.name != "pca_mean.vtu"]
+    sample_files = sample_files[:max_samples]
+    if len(sample_files) < 3:
+        raise FileNotFoundError(
+            f"Need at least 3 non-reference sample meshes under {sample_dir} "
+            f"or {kcl_dir}.\n"
+            "See data/README.md for manual download instructions."
+        )
+
+    reference_mesh = cast(pv.DataSet, pv.read(str(reference_file)))
+    sample_meshes = [cast(pv.DataSet, pv.read(str(f))) for f in sample_files]
+
+    workflow = WorkflowCreateStatisticalModel(
+        sample_meshes=sample_meshes,
+        reference_mesh=reference_mesh,
+        pca_number_of_components=pca_components,
+        log_level=log_level,
+    )
+    result = workflow.run_workflow()
+
+    mean_surface: pv.PolyData = result["mean_surface"]
+    mean_surface_file = output_dir / "pca_mean_surface.vtp"
+    mean_surface.save(str(mean_surface_file))
+
+    # Serialise the JSON-safe parts of the PCA model
+    pca_model: dict[str, Any] = result["pca_model"]
+    model_file = output_dir / "pca_model.json"
+    json_safe: dict[str, Any] = {}
+    for k, v in pca_model.items():
+        if isinstance(v, np.ndarray):
+            json_safe[k] = v.tolist()
+        elif isinstance(v, (int, float, str, bool, list)):
+            json_safe[k] = v
+    with open(model_file, "w") as fh:
+        json.dump(json_safe, fh, indent=2)
+
+    # Screenshots
+    tt = TestTools(
+        results_dir=output_dir,
+        baselines_dir=output_dir / "baselines",
+        class_name="tutorial_04",
+        results_output_dir=output_dir,
+        log_level=log_level,
+    )
+
+    screenshots: list[Path] = []
+
+    # Mean model
+    screenshots.append(
+        tt.save_screenshot_mesh(
+            mean_surface,
+            "pca_mean_model.png",
+            camera_position="iso",
+            color="steelblue",
+            opacity=0.9,
+        )
+    )
+
+    # First two PCA modes: show mean +/- 2 sigma side-by-side
+    eigenvectors: Any = pca_model.get("eigenvectors")
+    eigenvalues: Any = pca_model.get("eigenvalues")
+    mean_points = np.asarray(mean_surface.points)
+
+    try:
+        pv.start_xvfb()
+    except Exception:
+        pass
+
+    for mode_idx in range(min(2, pca_components)):
+        if eigenvectors is None or eigenvalues is None:
+            break
+
+        sigma = float(np.sqrt(eigenvalues[mode_idx]))
+        ev = np.asarray(eigenvectors[:, mode_idx]).reshape(-1, 3)
+
+        minus_mesh = mean_surface.copy()
+        minus_mesh.points = mean_points - 2 * sigma * ev
+        plus_mesh = mean_surface.copy()
+        plus_mesh.points = mean_points + 2 * sigma * ev
+
+        plotter = pv.Plotter(off_screen=True, window_size=[1200, 500], shape=(1, 3))
+        plotter.subplot(0, 0)
+        plotter.add_mesh(minus_mesh, color="royalblue", opacity=0.9)
+        plotter.add_text("mean - 2 sigma", font_size=10)
+        plotter.camera_position = "iso"
+        plotter.subplot(0, 1)
+        plotter.add_mesh(mean_surface, color="steelblue", opacity=0.9)
+        plotter.add_text("mean", font_size=10)
+        plotter.camera_position = "iso"
+        plotter.subplot(0, 2)
+        plotter.add_mesh(plus_mesh, color="coral", opacity=0.9)
+        plotter.add_text("mean + 2 sigma", font_size=10)
+        plotter.camera_position = "iso"
+
+        png_name = f"pca_mode_{mode_idx + 1:02d}.png"
+        png_path = output_dir / png_name
+        png_path.parent.mkdir(parents=True, exist_ok=True)
+        plotter.screenshot(str(png_path))
+        plotter.close()
+        screenshots.append(png_path)
+
+    return {
+        "pca_model": pca_model,
+        "mean_surface": mean_surface,
+        "model_file": model_file,
+        "mean_surface_file": mean_surface_file,
+        "screenshots": screenshots,
+    }
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(
+        description=__doc__,
+        formatter_class=argparse.RawDescriptionHelpFormatter,
+    )
+    parser.add_argument(
+        "--data-dir",
+        type=Path,
+        default=Path("data"),
+        help="Root data directory (default: ./data)",
+    )
+    parser.add_argument(
+        "--output-dir",
+        type=Path,
+        default=Path("output") / "tutorial_04",
+        help="Output directory (default: ./output/tutorial_04)",
+    )
+    parser.add_argument(
+        "--pca-components",
+        type=int,
+        default=10,
+        help="Number of PCA modes to retain (default: 10)",
+    )
+    args = parser.parse_args()
+
+    results = run_tutorial(
+        args.data_dir,
+        args.output_dir,
+        pca_components=args.pca_components,
+    )
+    print(f"PCA model:    {results['model_file']}")
+    print(f"Mean surface: {results['mean_surface_file']}")
+    print(f"Screenshots:  {[str(p) for p in results['screenshots']]}")
diff --git a/tutorials/tutorial_05_vtk_to_usd.py b/tutorials/tutorial_05_vtk_to_usd.py
new file mode 100644
index 0000000..0d5cf45
--- /dev/null
+++ b/tutorials/tutorial_05_vtk_to_usd.py
@@ -0,0 +1,183 @@
+"""
+Tutorial 5: VTK Surface Series to Animated USD
+
+Purpose
+-------
+Convert one or more VTK surface files into a USD file suitable for NVIDIA
+Omniverse. Supports a single static mesh, a time-series (animated) set of
+meshes, or a mesh with scalar data visualised via a colormap. This tutorial
+uses the surface files produced by Tutorial 2 (CT Segmentation to VTK) as
+input, but any VTK/VTP/VTU files will work.
+
+Inputs
+------
+- One or more VTK-compatible surface files (``.vtp`` / ``.vtk`` / ``.vtu``).
+  This tutorial looks for ``output/tutorial_02/patient_surfaces.vtp`` (output
+  of Tutorial 2). If that file does not exist, it falls back to any ``.vtp``
+  under ``data/``.
+
+Outputs
+-------
+- ``output_dir/surfaces.usd`` - USD file with anatomy materials applied
+- Screenshots (PNG):
+  - ``usd_mesh_rendering.png`` - PyVista off-screen render of the mesh
+
+Strengths
+---------
+- Supports time-varying USD for animated sequences (one VTK file per frame).
+- Three appearance modes: ``solid`` (flat colour), ``anatomy`` (clinical
+  material by anatomy type), and ``colormap`` (scalar field visualisation).
+- Coordinate system is automatically converted from RAS to Omniverse Y-up.
+
+Weaknesses / Limitations
+------------------------
+- Requires Tutorial 2 output (or any VTK file) as input; not standalone.
+- Time-series ordering relies on a filename regex pattern (``\.t\d+\.vtp$``);
+  non-conforming filenames are treated as static single-frame input.
+- USD materials use UsdPreviewSurface; advanced Omniverse MDL materials require
+  additional post-processing.
+
+Classes Used
+------------
+- WorkflowConvertVTKToUSD (workflow_convert_vtk_to_usd.py):
+    Loads VTK files, splits meshes, converts to USD, and applies appearance.
+- ConvertVTKToUSD (convert_vtk_to_usd.py):
+    High-level PyVista-to-USD converter with colormap support (used internally).
+- USDAnatomyTools (usd_anatomy_tools.py):
+    Applies clinical material colours to USD prims (used internally).
+
+CLI Equivalent
+--------------
+The same main outputs (without screenshots) can be produced via the CLI::
+
+    physiomotion4d-convert-vtk-to-usd \\
+        --input output/tutorial_02/patient_surfaces.vtp \\
+        --output-usd ./output/tutorial_05/surfaces.usd \\
+        --appearance anatomy \\
+        --anatomy-type heart
+
+See ``src/physiomotion4d/cli/convert_vtk_to_usd.py`` for full CLI documentation.
+
+Data Required
+-------------
+This tutorial uses the output of Tutorial 2 (``output/tutorial_02/patient_surfaces.vtp``).
+Run Tutorial 2 first, or provide any VTK surface file via the ``--vtk-file`` flag.
+"""
+
+from __future__ import annotations
+
+import argparse
+import logging
+from pathlib import Path
+from typing import Any, Optional
+
+import pyvista as pv
+
+from physiomotion4d.test_tools import TestTools
+from physiomotion4d.workflow_convert_vtk_to_usd import WorkflowConvertVTKToUSD
+
+
+def run_tutorial(
+    data_dir: Path,
+    output_dir: Path,
+    *,
+    vtk_file: Optional[Path] = None,
+    log_level: int = logging.INFO,
+) -> dict[str, Any]:
+    """Run Tutorial 5: VTK Surface Series to Animated USD.
+
+    Args:
+        data_dir: Root of the ``data/`` directory (see data/README.md).
+        output_dir: Directory to write outputs and screenshots.
+        vtk_file: Explicit path to a VTK file; overrides auto-discovery.
+        log_level: Python logging level.
+
+    Returns:
+        dict with keys:
+
+        - ``'usd_file'`` (str): path to the output USD file.
+        - ``'screenshots'`` (list[Path]): paths to saved PNG screenshots.
+    """
+    output_dir.mkdir(parents=True, exist_ok=True)
+
+    if vtk_file is None:
+        # Prefer Tutorial 2 output
+        project_root = data_dir.parent
+        candidate = project_root / "output" / "tutorial_02" / "patient_surfaces.vtp"
+        if candidate.exists():
+            vtk_file = candidate
+        else:
+            # Fall back to any .vtp under data/
+            found = list(data_dir.rglob("*.vtp"))
+            if not found:
+                raise FileNotFoundError(
+                    "No VTK file found. Run Tutorial 2 first, or specify "
+                    "--vtk-file <path>."
+                )
+            vtk_file = found[0]
+
+    output_usd = output_dir / "surfaces.usd"
+
+    workflow = WorkflowConvertVTKToUSD(
+        vtk_files=[vtk_file],
+        output_usd=output_usd,
+        appearance="anatomy",
+        anatomy_type="heart",
+        separate_by_connectivity=True,
+        log_level=log_level,
+    )
+    usd_path = workflow.run()
+
+    # Screenshots
+    tt = TestTools(
+        results_dir=output_dir,
+        baselines_dir=output_dir / "baselines",
+        class_name="tutorial_05",
+        results_output_dir=output_dir,
+        log_level=log_level,
+    )
+
+    screenshots: list[Path] = []
+
+    mesh = pv.read(str(vtk_file))
+    screenshots.append(
+        tt.save_screenshot_mesh(
+            mesh,
+            "usd_mesh_rendering.png",
+            camera_position="iso",
+            color="lightcoral",
+            opacity=0.9,
+        )
+    )
+
+    return {"usd_file": usd_path, "screenshots": screenshots}
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(
+        description=__doc__,
+        formatter_class=argparse.RawDescriptionHelpFormatter,
+    )
+    parser.add_argument(
+        "--data-dir",
+        type=Path,
+        default=Path("data"),
+        help="Root data directory (default: ./data)",
+    )
+    parser.add_argument(
+        "--output-dir",
+        type=Path,
+        default=Path("output") / "tutorial_05",
+        help="Output directory (default: ./output/tutorial_05)",
+    )
+    parser.add_argument(
+        "--vtk-file",
+        type=Path,
+        default=None,
+        help="Explicit VTK input file (default: auto-discover from Tutorial 2 output)",
+    )
+    args = parser.parse_args()
+
+    results = run_tutorial(args.data_dir, args.output_dir, vtk_file=args.vtk_file)
+    print(f"USD file:    {results['usd_file']}")
+    print(f"Screenshots: {[str(p) for p in results['screenshots']]}")
diff --git a/tutorials/tutorial_06_reconstruct_highres_4d_ct.py b/tutorials/tutorial_06_reconstruct_highres_4d_ct.py
new file mode 100644
index 0000000..586c882
--- /dev/null
+++ b/tutorials/tutorial_06_reconstruct_highres_4d_ct.py
@@ -0,0 +1,242 @@
+"""
+Tutorial 6: Reconstruct High-Resolution 4D CT
+
+Purpose
+-------
+Reconstruct a high-resolution dynamic 4D CT volume from a time series of
+lower-resolution or sparse CT frames. The workflow registers each time frame
+to a high-resolution reference image, producing a sequence of reconstructed
+volumes that share the spatial resolution of the reference. This is useful for
+respiratory-gated lung CT (DirLab-4DCT) where breath-hold reference scans are
+available alongside lower-quality respiratory-phase images.
+
+Inputs
+------
+- A list of 3D CT images (``itk.Image``): the time series to reconstruct.
+  Expected location: ``data/DirLab-4DCT/Case1/`` (T00-T90 phases).
+- A high-resolution fixed reference image (``itk.Image``):
+  the target space for reconstruction.
+  Expected location: ``data/DirLab-4DCT/Case1/`` (any phase used as reference).
+
+Note
+----
+The DirLab-4DCT sample data used by this tutorial does not include a separate
+high-resolution breath-hold reference image. For demonstration and regression
+testing, the tutorial uses one available DirLab respiratory phase as the fixed
+reference, so the reconstructed outputs inherit that phase image's resolution
+rather than true higher-resolution reference spacing.
+
+Outputs
+-------
+- ``output_dir/reconstructed_frame_<N>.mha`` - one reconstructed 3D image per frame
+- Screenshots (PNG):
+  - ``reference_frame.png`` - axial slice of the high-resolution reference image
+  - ``reconstructed_frame.png`` - axial slice of the first reconstructed frame
+
+Strengths
+---------
+- Bidirectional propagation of registration from the reference frame reduces
+  accumulated error for frames far from the reference.
+- Temporal smoothing via ``prior_weight`` parameter reduces frame-to-frame jitter.
+- Supports ``'ants'``, ``'icon'``, and ``'ants_icon'`` (two-stage) registration.
+
+Weaknesses / Limitations
+------------------------
+- Requires the DirLab-4DCT dataset (manual download; see data/README.md).
+- ICON registration (default part of ``'ants_icon'``) requires a GPU.
+- Reconstruction quality is bounded by the accuracy of the registration; large
+  respiratory excursion between phases can cause residual artefacts.
+- Runtime is proportional to the number of frames times registration cost.
+
+Classes Used
+------------
+- WorkflowReconstructHighres4DCT (workflow_reconstruct_highres_4d_ct.py):
+    Registers each time frame to the fixed reference and reconstructs volumes.
+- RegisterTimeSeriesImages (register_time_series_images.py):
+    Chains frame-to-frame registration with optional temporal smoothing (used
+    internally).
+- RegisterImagesANTs / RegisterImagesICON (register_images_ants.py / _icon.py):
+    Individual frame registration backends (used internally).
+
+CLI Equivalent
+--------------
+The same main outputs (without screenshots) can be produced via the CLI::
+
+    physiomotion4d-reconstruct-highres-4d-ct \\
+        --time-series-dir data/DirLab-4DCT/Case1 \\
+        --fixed-image data/DirLab-4DCT/Case1/case1_T00.mhd \\
+        --registration-method ants \\
+        --output-dir ./output/tutorial_06
+
+See ``src/physiomotion4d/cli/reconstruct_highres_4d_ct.py`` for full CLI
+documentation.
+
+Data Required
+-------------
+See data/README.md for download instructions and dataset licensing.
+Dataset: DirLab 4D-CT - https://www.dir-lab.com/ReferenceData.html
+Manual download required. Place files under ``data/DirLab-4DCT/`` as described
+in data/README.md.
+"""
+
+from __future__ import annotations
+
+import argparse
+import logging
+from pathlib import Path
+from typing import Any
+
+import itk
+
+from physiomotion4d.test_tools import TestTools
+from physiomotion4d.workflow_reconstruct_highres_4d_ct import (
+    WorkflowReconstructHighres4DCT,
+)
+
+
+def run_tutorial(
+    data_dir: Path,
+    output_dir: Path,
+    *,
+    case: int = 1,
+    max_frames: int = 4,
+    registration_method: str = "ants",
+    log_level: int = logging.INFO,
+) -> dict[str, Any]:
+    """Run Tutorial 6: Reconstruct High-Resolution 4D CT.
+
+    Args:
+        data_dir: Root of the ``data/`` directory (see data/README.md).
+        output_dir: Directory to write outputs and screenshots.
+        case: DirLab case number (1-10). Default: 1.
+        max_frames: Maximum number of time frames to reconstruct (for speed).
+        registration_method: ``'ants'`` (CPU-capable) or ``'icon'`` (GPU) or
+            ``'ants_icon'`` (two-stage). Default: ``'ants'``.
+        log_level: Python logging level.
+
+    Returns:
+        dict with keys:
+
+        - ``'reconstructed_images'`` (list[itk.Image]): reconstructed volumes.
+        - ``'reconstructed_files'`` (list[Path]): saved ``.mha`` paths.
+        - ``'screenshots'`` (list[Path]): paths to saved PNG screenshots.
+    """
+    output_dir.mkdir(parents=True, exist_ok=True)
+
+    dirlab_dir = data_dir / "DirLab-4DCT"
+    case_dir = dirlab_dir / f"Case{case}"
+
+    # Discover phase images (MetaImage .mhd or .mha)
+    phase_files = sorted(list(case_dir.glob("*.mhd")) + list(case_dir.glob("*.mha")))
+    if not phase_files:
+        phase_files = sorted(
+            list(dirlab_dir.glob(f"Case{case}Pack_T*.mhd"))
+            + list(dirlab_dir.glob(f"Case{case}Pack_T*.mha"))
+        )
+    if not phase_files:
+        raise FileNotFoundError(
+            f"No .mhd / .mha files found for DirLab-4DCT case {case} under "
+            f"{case_dir} or {dirlab_dir}.\n"
+            "See data/README.md for manual download instructions."
+        )
+
+    phase_files = phase_files[:max_frames]
+    time_series = [itk.imread(str(f)) for f in phase_files]
+    fixed_image = time_series[0]  # use first phase as high-res reference
+
+    workflow = WorkflowReconstructHighres4DCT(
+        time_series_images=time_series,
+        fixed_image=fixed_image,
+        reference_frame=0,
+        registration_method=registration_method,
+        log_level=log_level,
+    )
+    workflow.set_modality("ct")
+    result = workflow.run_workflow()
+
+    reconstructed: list[itk.Image] = result["reconstructed_images"]
+    reconstructed_files: list[Path] = []
+    for i, vol in enumerate(reconstructed):
+        out_path = output_dir / f"reconstructed_frame_{i:03d}.mha"
+        itk.imwrite(vol, str(out_path), compression=True)
+        reconstructed_files.append(out_path)
+
+    # Screenshots
+    tt = TestTools(
+        results_dir=output_dir,
+        baselines_dir=output_dir / "baselines",
+        class_name="tutorial_06",
+        results_output_dir=output_dir,
+        log_level=log_level,
+    )
+
+    screenshots: list[Path] = []
+
+    screenshots.append(
+        tt.save_screenshot_image_slice(
+            fixed_image,
+            "reference_frame.png",
+            axis=0,
+            slice_fraction=0.5,
+            colormap="gray",
+        )
+    )
+
+    if reconstructed:
+        screenshots.append(
+            tt.save_screenshot_image_slice(
+                reconstructed[0],
+                "reconstructed_frame.png",
+                axis=0,
+                slice_fraction=0.5,
+                colormap="gray",
+            )
+        )
+
+    return {
+        "reconstructed_images": reconstructed,
+        "reconstructed_files": reconstructed_files,
+        "screenshots": screenshots,
+    }
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(
+        description=__doc__,
+        formatter_class=argparse.RawDescriptionHelpFormatter,
+    )
+    parser.add_argument(
+        "--data-dir",
+        type=Path,
+        default=Path("data"),
+        help="Root data directory (default: ./data)",
+    )
+    parser.add_argument(
+        "--output-dir",
+        type=Path,
+        default=Path("output") / "tutorial_06",
+        help="Output directory (default: ./output/tutorial_06)",
+    )
+    parser.add_argument(
+        "--case",
+        type=int,
+        default=1,
+        choices=list(range(1, 11)),
+        help="DirLab case number 1-10 (default: 1)",
+    )
+    parser.add_argument(
+        "--registration-method",
+        default="ants",
+        choices=["ants", "icon", "ants_icon"],
+        help="Registration method (default: ants)",
+    )
+    args = parser.parse_args()
+
+    results = run_tutorial(
+        args.data_dir,
+        args.output_dir,
+        case=args.case,
+        registration_method=args.registration_method,
+    )
+    print(f"Reconstructed frames: {[str(p) for p in results['reconstructed_files']]}")
+    print(f"Screenshots:          {[str(p) for p in results['screenshots']]}")