diff --git a/README.md b/README.md index bb8286e..441ba96 100644 --- a/README.md +++ b/README.md @@ -7,14 +7,17 @@ # PhantomKit -**PhantomKit** is a Python toolkit for automated quality assurance (QA) of medical imaging scanners using physical phantoms. It provides pydra-based workflows that register phantom scans to a reference template, extract per-vial signal statistics across multiple contrast types, and generate publication-quality plots — supporting both MRI and PET phantom protocols. +**PhantomKit** is a Python toolkit for automated quality assurance (QA) of MRI scanners using physical phantoms. It provides a high-level processing engine and pydra-based workflows that register phantom scans to a reference template, extract per-vial signal statistics across multiple contrast types, and generate publication-quality plots — with full support for DWI preprocessing and diffusion metric (ADC, FA) extraction. ## Features -- **Template-based registration** — iterative ANTs SyN registration with automatic orientation search across a rotation library -- **Vial metric extraction** — per-vial mean, median, std, min and max across all contrast images, written to CSV -- **Plotting** — scatter plots of vial intensity and parametric map plots (T1/IR, T2/TE) with mrview ROI overlays -- **Protocol support** — extensible `protocols` sub-package for phantom- and project-specific workflow configurations +- **End-to-end pipeline** — single command processes a raw DICOM session directory through DWI preprocessing, phantom QC in DWI space, and native contrast QC +- **Automatic series classification** — DWI, reverse phase-encode, T1, IR, and TE series are detected and paired automatically from folder names and DICOM sidecar metadata; no manual configuration required +- **DWI preprocessing** — FSL `dwifslpreproc` with automatic phase-encoding correction mode selection (`rpe_none`, `rpe_pair`, `rpe_all`), optional denoising/Gibbs correction, tensor fitting, and T1-to-DWI co-registration via FLIRT +- **Template-based registration** — iterative ANTs rigid registration with automatic orientation search across a rotation library; vial masks propagated to subject space via inverse transform +- **Vial metric extraction** — per-vial mean, median, std, min and max across all contrast images (T1, IR, TE, ADC, FA), written to CSV +- **Plotting** — ADC/FA scatter plots with SPIRIT reference values, T1/IR and T2/TE parametric map plots with mrview ROI overlays and Monte Carlo 95% CI bands +- **Checkpoint-based resumption** — re-running the pipeline skips stages whose outputs already exist - **Parallel batch processing** — pydra-native splitting and combining for multi-session datasets ## Installation @@ -23,37 +26,100 @@ python -m pip install phantomkit ``` +### External dependencies + +The pipeline requires FSL, MRtrix3, ANTs, and dcm2niix to be available on `PATH`: + +- [FSL](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/) — `dwifslpreproc`, `flirt`, `convert_xfm` +- [MRtrix3](https://www.mrtrix.org/) — `mrconvert`, `dwi2tensor`, `tensor2metric`, `dwidenoise`, `mrdegibbs`, `mrstats`, `mrview` +- [ANTs](http://stnava.github.io/ANTs/) — `antsRegistrationSyN.sh`, `antsApplyTransforms` +- [dcm2niix](https://github.com/rordenlab/dcm2niix) — DICOM to NIfTI conversion + ## Basic usage +### End-to-end pipeline + +Point the pipeline at a session directory containing DICOM subdirectories: + +```bash +phantomkit pipeline \ + --input-dir /data/session01 \ + --output-dir /results/session01 \ + --phantom SPIRIT +``` + +Optional flags: + +``` +--denoise-degibbs Apply dwidenoise + mrdegibbs before preprocessing +--gradcheck Run dwigradcheck to verify gradient orientations +--nocleanup Keep intermediate tmp/ directories after completion +--readout-time Override TotalReadoutTime (seconds) for dwifslpreproc +--eddy-options Override FSL eddy options string +--dry-run Print the planned workflow without executing +``` + +Output structure: + +``` +/results/session01/ + / + DWI_preproc_biascorr.mif.gz + ADC.nii.gz + FA.nii.gz + T1_in_DWI_space.nii.gz + metrics/ ← per-vial CSVs and QA plots + vial_segmentations/ + native_contrasts_staging/ + metrics/ ← per-vial CSVs and parametric map plots + vial_segmentations/ + images_template_space/ +``` + +### Python API + +```python +from phantomkit.phantom_processor import PhantomProcessor + +processor = PhantomProcessor( + template_dir="/templates/SPIRIT", + output_base_dir="/results", + rotation_library_file="/templates/rotations.txt", +) +results = processor.process_session("/data/session01/T1.nii.gz") +``` + +Or via the pydra workflow API: + ```python -from phantomkit.protocols.gsp_spirit import GspSpiritAnalysis +from phantomkit.analyses.vial_signal import VialSignalAnalysis -wf = GspSpiritAnalysis( - input_image="/data/session01/t1_mprage.nii.gz", - template_dir="/templates/gsp_spirit", - rotation_library_file="/templates/gsp_spirit/rotations.txt", +wf = VialSignalAnalysis( + input_image="/data/session01/T1.nii.gz", + template_dir="/templates/SPIRIT", + rotation_library_file="/templates/rotations.txt", ) -outputs = wf(cache_root="/data/cache-root") +outputs = wf(cache_root="/data/pydra-cache") ``` -Or via the command line: +### CLI — run a protocol directly ```bash # Single session -phantom-process run gsp-spirit /data/session01/t1_mprage.nii.gz \ - --template-dir /templates/gsp_spirit \ - --rotation-library-file /templates/gsp_spirit/rotations.txt \ - --output-base-dir /results - -# Batch — process every matching image found under /data/ -phantom-process run gsp-spirit /data/ \ - --template-dir /templates/gsp_spirit \ - --rotation-library-file /templates/gsp_spirit/rotations.txt \ - --output-base-dir /results \ - --pattern "*t1*mprage*.nii.gz" +phantomkit run vial-signal /data/session01/T1.nii.gz \ + --template-dir /templates/SPIRIT \ + --rotation-library-file /templates/rotations.txt \ + --output-base-dir /results + +# Batch mode — process every matching image under /data/ +phantomkit run vial-signal /data/ \ + --template-dir /templates/SPIRIT \ + --rotation-library-file /templates/rotations.txt \ + --output-base-dir /results \ + --pattern "*T1*.nii.gz" # List available protocols -phantom-process list +phantomkit list ``` ### Plotting @@ -61,22 +127,30 @@ phantom-process list Generate QA plots from existing CSV metric files: ```bash -# Vial intensity scatter plot for one contrast -phantom-process plot vial-intensity \ - /results/session01/metrics/session01_t1_mprage_mean_matrix.csv scatter \ - --std_csv /results/session01/metrics/session01_t1_mprage_std_matrix.csv \ - --output /results/session01/metrics/session01_t1_PLOTmeanstd.png +# Generic scatter plot +phantomkit plot vial-intensity \ + /results/session01/metrics/session01_T1_mean_matrix.csv scatter \ + --std-csv /results/session01/metrics/session01_T1_std_matrix.csv \ + --output /results/session01/metrics/session01_T1_PLOTmeanstd.png + +# ADC scatter plot with SPIRIT reference values +phantomkit plot vial-intensity \ + /results/session01/metrics/session01_ADC_mean_matrix.csv scatter \ + --std-csv /results/session01/metrics/session01_ADC_std_matrix.csv \ + --phantom SPIRIT \ + --template-dir /templates \ + --output /results/session01/metrics/session01_ADC_PLOTmeanstd.png # T1 inversion-recovery parametric map plot -phantom-process plot maps-ir \ - /results/session01/images_template_space/ir_*.nii.gz \ - --metric_dir /results/session01/metrics \ +phantomkit plot maps-ir \ + /results/session01/images_template_space/se_ir_*.nii.gz \ + --metric-dir /results/session01/metrics \ --output /results/session01/metrics/session01_T1map_plot.png # T2 spin-echo parametric map plot -phantom-process plot maps-te \ - /results/session01/images_template_space/te_*.nii.gz \ - --metric_dir /results/session01/metrics \ +phantomkit plot maps-te \ + /results/session01/images_template_space/t2_se_*.nii.gz \ + --metric-dir /results/session01/metrics \ --output /results/session01/metrics/session01_T2map_plot.png ``` @@ -84,4 +158,4 @@ See the [CLI documentation](https://australian-imaging-service.github.io/phantom ## License -Copyright 2026 Australian Imaging Service. Released under the [Apache License 2.0](LICENSE). +Copyright 2026 Australian Imaging Service. Released under the [Apache License 2.0](LICENSE). \ No newline at end of file diff --git a/docs/source/api.rst b/docs/source/api.rst index ec8ee58..44ea4b4 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -1,19 +1,83 @@ API === +This page documents the public Python API of PhantomKit. The package is +organised into five sub-packages: + +- :mod:`phantomkit.analyses` — top-level pydra workflows (user-facing entry + points) +- :mod:`phantomkit.phantom_processor` — imperative processing engine + (:class:`~phantomkit.phantom_processor.PhantomProcessor`) +- :mod:`phantomkit.registration` — ANTs registration tasks and workflows +- :mod:`phantomkit.metrics` — vial mask transform and metric extraction + workflows +- :mod:`phantomkit.plotting` — QA plot generation functions and pydra tasks + + Analyses --------- +--------- + +The ``analyses`` sub-package contains the pydra workflow classes that serve as +the primary programmatic entry points. + +Vial signal analysis +~~~~~~~~~~~~~~~~~~~~ .. autoclass:: phantomkit.analyses.vial_signal.VialSignalAnalysis -.. autoclass:: phantomkit.protocols.diffusion_metrics.DiffusionMetricsAnalysis +.. autoclass:: phantomkit.analyses.vial_signal.VialSignalAnalysisBatch + +Diffusion metrics analysis +~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. autoclass:: phantomkit.analyses.diffusion_metrics.DiffusionMetricsAnalysis + +.. autoclass:: phantomkit.analyses.diffusion_metrics.DiffusionMetricsAnalysisBatch + + +Phantom processor +------------------ + +:class:`~phantomkit.phantom_processor.PhantomProcessor` is the imperative +counterpart to :class:`~phantomkit.analyses.vial_signal.VialSignalAnalysis`. +It is used internally by the pipeline orchestrator (Stages 2 and 3) and can +also be used directly for scripted single-session processing. + +.. autoclass:: phantomkit.phantom_processor.PhantomProcessor + :members: process_session + + +Registration +------------ + +.. autoclass:: phantomkit.registration.IterativeRegistration + +.. autoclass:: phantomkit.registration.CheckRegistration + +.. autoclass:: phantomkit.registration.SaveTemplateInScannerSpace + +.. autoclass:: phantomkit.registration.RegistrationSynN + + +Metrics +------- + +.. autoclass:: phantomkit.metrics.TransformVialsToSubjectSpace + +.. autoclass:: phantomkit.metrics.ExtractMetricsFromContrasts + +.. autoclass:: phantomkit.metrics.TransformContrastsToTemplateSpace Plotting -------- +.. autofunction:: phantomkit.plotting.vial_intensity.plot_vial_intensity + .. autofunction:: phantomkit.plotting.maps_ir.plot_vial_ir_means_std .. autofunction:: phantomkit.plotting.maps_te.plot_vial_te_means_std -.. autofunction:: phantomkit.plotting.vial_intensity.plot_vial_intensity +.. autoclass:: phantomkit.plotting.visualization.GeneratePlots + +.. autoclass:: phantomkit.plotting.visualization.BuildRoiOverlay \ No newline at end of file diff --git a/docs/source/cli.rst b/docs/source/cli.rst index 3608f63..fb5d103 100644 --- a/docs/source/cli.rst +++ b/docs/source/cli.rst @@ -1,34 +1,98 @@ Command-line interface ====================== -PhantomKit installs the ``phantom-process`` command, which groups two -subcommand families: +PhantomKit installs the ``phantomkit`` command, which groups four subcommand +families: -- ``phantom-process run`` — execute a full QA protocol workflow -- ``phantom-process plot`` — generate individual QA plots from existing CSVs +- ``phantomkit pipeline`` — end-to-end phantom QC + DWI orchestrator +- ``phantomkit run`` — execute a single pydra protocol workflow +- ``phantomkit plot`` — generate individual QA plots from existing CSVs +- ``phantomkit list`` — list all available protocol workflows Use ``--help`` at any level for a full option listing:: - phantom-process --help - phantom-process run --help - phantom-process run gsp-spirit --help - phantom-process plot --help - phantom-process plot vial-intensity --help + phantomkit --help + phantomkit pipeline --help + phantomkit run --help + phantomkit run vial-signal --help + phantomkit plot --help + phantomkit plot vial-intensity --help + + +End-to-end pipeline +-------------------- + +The ``pipeline`` command is the primary entry point for processing a complete +scanner session. It scans ``--input-dir`` for acquisition subdirectories, +classifies them automatically, and runs up to three stages: + +- **Stage 1** — DWI preprocessing (runs if DWI acquisitions are found) +- **Stage 2** — Phantom QC in DWI space (runs per DWI series after Stage 1) +- **Stage 3** — Phantom QC on native T1/IR/TE contrasts (runs if those series + are present) + +Stages 1 and 3 execute in parallel where possible. + +.. code-block:: console + + $ phantomkit pipeline \ + --input-dir /data/session01 \ + --output-dir /results/session01 \ + --phantom SPIRIT + +Options: + +.. code-block:: text + + --input-dir PATH Root directory containing acquisition subdirectories + (DICOM folders). [required] + --output-dir PATH Top-level output directory. All results are written + here. [required] + --phantom TEXT Phantom name, e.g. SPIRIT. Used to locate + template_data//. [required] + --denoise-degibbs Apply dwidenoise + mrdegibbs before preprocessing. + --gradcheck Run dwigradcheck to verify gradient orientations. + --nocleanup Keep DWI tmp/ intermediate directories. + --readout-time FLOAT Override TotalReadoutTime (seconds) for dwifslpreproc. + --eddy-options TEXT Override FSL eddy options string. + --dry-run Plan and print commands; do not execute any processing. + +Output structure: + +.. code-block:: text + + /results/session01/ + / + DWI_preproc_biascorr.mif.gz + ADC.nii.gz + FA.nii.gz + T1_in_DWI_space.nii.gz + metrics/ + vial_segmentations/ + native_contrasts_staging/ + metrics/ + vial_segmentations/ + images_template_space/ Discovering available protocols --------------------------------- +--------------------------------- .. code-block:: console - $ phantom-process list + $ phantomkit list Available protocols: - gsp-spirit (batch supported) + vial-signal (batch supported) Pydra workflow for processing a single GSP SPIRIT phantom MRI session. + diffusion-metrics (batch supported) + Full SPIRIT phantom DWI analysis workflow. + +Running a protocol workflow +---------------------------- -Running a protocol ------------------- +The ``run`` subgroup executes pydra protocol workflows directly, bypassing the +pipeline orchestrator. Single session ~~~~~~~~~~~~~~ @@ -37,10 +101,10 @@ Pass a path to a NIfTI file as ``INPUT`` to process one session: .. code-block:: console - $ phantom-process run gsp-spirit /data/session01/t1_mprage.nii.gz \ - --template-dir /templates/gsp_spirit \ - --rotation-library-file /templates/gsp_spirit/rotations.txt \ - --output-base-dir /results + $ phantomkit run vial-signal /data/session01/T1.nii.gz \ + --template-dir /templates/SPIRIT \ + --rotation-library-file /templates/rotations.txt \ + --output-base-dir /results The session name is taken from the parent directory of the image (``session01`` in the example above). All ``*.nii.gz`` files in the same @@ -54,23 +118,23 @@ recursively inside it: .. code-block:: console - $ phantom-process run gsp-spirit /data \ - --template-dir /templates/gsp_spirit \ - --rotation-library-file /templates/gsp_spirit/rotations.txt \ - --output-base-dir /results \ - --pattern "*t1*mprage*.nii.gz" \ - --plugin cf + $ phantomkit run vial-signal /data \ + --template-dir /templates/SPIRIT \ + --rotation-library-file /templates/rotations.txt \ + --output-base-dir /results \ + --pattern "*T1*.nii.gz" \ + --worker cf ``--pattern`` controls which files are selected (default ``*.nii.gz``). -``--plugin`` selects the pydra execution backend: ``cf`` (concurrent -futures, default) for parallel execution or ``serial`` for sequential. +``--worker`` selects the pydra execution backend: ``cf`` (concurrent futures, +default) for parallel execution or ``serial`` for sequential. Protocol options ~~~~~~~~~~~~~~~~ -All parameters of the underlying workflow are exposed as options. Run -``phantom-process run --help`` to see the full list. For -``gsp-spirit``: +All parameters of the underlying workflow are exposed as CLI options. Run +``phantomkit run --help`` to see the full list. For +``vial-signal``: .. code-block:: text @@ -80,55 +144,68 @@ All parameters of the underlying workflow are exposed as options. Run --rotation-library-file PATH Path to the rotation library text file. [required] --output-base-dir PATH Root output directory (default: cwd). - --plugin [cf|serial] Pydra execution plugin. [default: cf] + --worker [cf|serial] Pydra execution backend. [default: cf] --pattern TEXT Glob pattern for batch mode. [default: *.nii.gz] Generating plots ----------------- +----------------- -The ``plot`` subcommands generate standalone QA figures from CSV metric -files already written to disk (e.g. after a completed ``run``). +The ``plot`` subcommands generate standalone QA figures from CSV metric files +already written to disk (e.g. after a completed ``run`` or ``pipeline``). Vial intensity scatter plot ~~~~~~~~~~~~~~~~~~~~~~~~~~~ Plots per-vial mean ± std as a scatter, line, or bar chart for a single -contrast: +contrast. In ADC mode (auto-detected from filename), measured values are +overlaid against SPIRIT reference values: .. code-block:: console - $ phantom-process plot vial-intensity \ - /results/session01/metrics/session01_t1_mprage_mean_matrix.csv \ + $ phantomkit plot vial-intensity \ + /results/session01/metrics/session01_T1_mean_matrix.csv \ + scatter \ + --std-csv /results/session01/metrics/session01_T1_std_matrix.csv \ + --output /results/session01/metrics/session01_T1_PLOTmeanstd.png + + # ADC mode — reference values overlaid automatically + $ phantomkit plot vial-intensity \ + /results/session01/metrics/session01_ADC_mean_matrix.csv \ scatter \ - --std_csv /results/session01/metrics/session01_t1_mprage_std_matrix.csv \ - --output /results/session01/metrics/session01_t1_PLOTmeanstd.png + --std-csv /results/session01/metrics/session01_ADC_std_matrix.csv \ + --phantom SPIRIT \ + --template-dir /templates \ + --output /results/session01/metrics/session01_ADC_PLOTmeanstd.png Options: .. code-block:: text Arguments: - CSV_FILE Mean-value CSV produced by ExtractMetricsFromContrasts. + CSV_FILE Mean-value CSV produced by metric extraction. {line|bar|scatter} Plot style. Options: - --std_csv TEXT CSV file containing standard deviations. - --roi_image TEXT PNG screenshot of the contrast with vial ROI overlay. - --annotate Annotate each point with mean ± std. - --output TEXT Output filename. [default: vial_subplot.png] + --std-csv PATH CSV file containing standard deviations. + --roi-image PATH PNG screenshot of the contrast with vial ROI overlay. + --annotate Annotate each point with mean ± std. + --phantom TEXT Phantom name (required in ADC mode). + --template-dir PATH TemplateData directory (required in ADC mode). + --output PATH Output filename. [default: vial_subplot.png] Inversion-recovery (T1) map plot ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -Plots vial mean ± std across multiple IR contrasts with T₁ curve fitting: +Plots vial mean ± std across multiple IR contrasts with T₁ inversion-recovery +curve fitting and 95% Monte Carlo confidence intervals: .. code-block:: console - $ phantom-process plot maps-ir \ - /results/session01/images_template_space/ir_*.nii.gz \ - --metric_dir /results/session01/metrics \ + $ phantomkit plot maps-ir \ + /results/session01/images_template_space/se_ir_*.nii.gz \ + --metric-dir /results/session01/metrics \ --output /results/session01/metrics/session01_T1map_plot.png Options: @@ -136,25 +213,26 @@ Options: .. code-block:: text Arguments: - CONTRAST_FILES... IR NIfTI images in template space. + CONTRAST_FILES... IR NIfTI images in template space (one per TI). Options: - -m, --metric_dir TEXT Directory containing mean/std CSV files. [required] - -o, --output TEXT Output filename. [default: vial_summary_T1.png] + -m, --metric-dir PATH Directory containing mean/std CSV files. [required] + -o, --output PATH Output filename. [default: vial_summary_T1.png] --annotate Annotate each point with mean ± std. - --roi_image TEXT PNG ROI overlay for the extra subplot. + --roi-image PATH PNG ROI overlay for the extra subplot. Spin-echo (T2/TE) map plot ~~~~~~~~~~~~~~~~~~~~~~~~~~~ -Plots vial mean ± std across multiple TE contrasts with mono-exponential -T₂ fitting: +Plots vial mean ± std across multiple TE contrasts with mono-exponential T₂ +decay fitting and 95% Monte Carlo confidence intervals: .. code-block:: console - $ phantom-process plot maps-te \ - /results/session01/images_template_space/te_*.nii.gz \ - --metric_dir /results/session01/metrics \ + $ phantomkit plot maps-te \ + /results/session01/images_template_space/t2_se_*.nii.gz \ + --metric-dir /results/session01/metrics \ --output /results/session01/metrics/session01_T2map_plot.png -Options are identical to ``maps-ir``. +Options are identical to ``maps-ir``, except the fitted parameter is T₂ and +the output CSV records T₂ values. \ No newline at end of file diff --git a/docs/source/index.rst b/docs/source/index.rst index 94df8c7..1d0c075 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -14,24 +14,36 @@ PhantomKit :target: https://pypi.python.org/pypi/phantomkit/ :alt: Latest Version -*PhantomKit* is a Python toolkit for automated quality assurance (QA) of medical imaging -scanners using physical phantoms. It provides pydra-based workflows that register phantom -scans to a reference template, extract per-vial signal statistics across multiple contrast -types, and generate publication-quality plots — supporting both MRI and PET phantom -protocols. +*PhantomKit* is a Python toolkit for automated quality assurance (QA) of MRI +scanners using physical phantoms. It provides pydra-based workflows and a +high-level processing engine that register phantom scans to a reference +template, extract per-vial signal statistics across multiple contrast types, +and generate publication-quality plots — with full support for DWI +preprocessing and diffusion metric (ADC, FA) extraction. Key features: -- **Template-based registration** — iterative ANTs SyN registration with automatic - orientation search across a rotation library -- **Vial metric extraction** — per-vial mean, median, std, min and max across all - contrast images, written to CSV -- **Plotting** — scatter plots of vial intensity and parametric map plots (T1/IR, T2/TE) - with mrview ROI overlays -- **Protocol support** — extensible ``protocols`` sub-package for phantom- and - project-specific workflow configurations -- **Parallel batch processing** — pydra-native splitting and combining for multi-session - datasets +- **End-to-end pipeline** — single ``phantomkit pipeline`` command processes a + raw DICOM session directory through DWI preprocessing (Stage 1), phantom QC + in DWI space (Stage 2), and native contrast QC (Stage 3) +- **Automatic series classification** — DWI, reverse phase-encode, T1, IR, and + TE series are detected from folder names and DICOM sidecar metadata; no manual + configuration required +- **DWI preprocessing** — FSL ``dwifslpreproc`` with automatic phase-encoding + correction mode selection (``rpe_none``, ``rpe_pair``, ``rpe_all``), optional + denoising/Gibbs correction, tensor fitting, and T1-to-DWI co-registration via + FLIRT +- **Template-based registration** — iterative ANTs rigid registration with + automatic orientation search across a rotation library; vial masks propagated + to subject space via inverse transform +- **Vial metric extraction** — per-vial mean, median, std, min and max across + all contrast images, written to CSV +- **Plotting** — ADC/FA scatter plots with SPIRIT reference values, T1/IR and + T2/TE parametric map plots with mrview ROI overlays +- **Parallel batch processing** — pydra-native splitting and combining for + multi-session datasets +- **Checkpoint-based resumption** — re-running the pipeline skips stages whose + outputs already exist Installation @@ -43,6 +55,17 @@ Installation $ python3 -m pip install phantomkit +External dependencies (must be on ``PATH``): + +- `FSL `_ — ``dwifslpreproc``, + ``flirt``, ``convert_xfm`` +- `MRtrix3 `_ — ``mrconvert``, ``dwi2tensor``, + ``tensor2metric``, ``dwidenoise``, ``mrdegibbs``, ``mrstats``, ``mrview`` +- `ANTs `_ — ``antsRegistrationSyN.sh``, + ``antsApplyTransforms`` +- `dcm2niix `_ — DICOM to NIfTI + conversion + License ------- @@ -59,4 +82,4 @@ Copyright 2026 Australian Imaging Service. quick_start cli - api + api \ No newline at end of file diff --git a/docs/source/quick_start.rst b/docs/source/quick_start.rst index 5130ef4..55ae3ab 100644 --- a/docs/source/quick_start.rst +++ b/docs/source/quick_start.rst @@ -1,18 +1,119 @@ Quick start =========== -PhantomKit analyses are implemented as `Pydra `__ workflows. To -execute them, you first parameterise the workflow by instantiate the Pydra workflow class, -then call the workflow object with execution args:: +PhantomKit can be used through the ``phantomkit`` command-line interface or +directly via its Python API. + + +End-to-end pipeline (CLI) +-------------------------- + +The fastest way to process a session is the ``pipeline`` command. Point it at +a directory of DICOM subdirectories, specify an output location, and name the +phantom: + +.. code-block:: console + + $ phantomkit pipeline \ + --input-dir /data/session01 \ + --output-dir /results/session01 \ + --phantom SPIRIT + +PhantomKit will: + +1. Scan ``--input-dir`` and classify subdirectories into T1, IR, TE, and DWI + series automatically. +2. Run DWI preprocessing (Stage 1) if DWI series are found. +3. Run phantom QC in DWI space (Stage 2) on each preprocessed DWI output. +4. Run phantom QC on native T1/IR/TE contrasts (Stage 3) if those series are + present. + +See :doc:`cli` for the full option reference. + + +Python API — single session +---------------------------- + +The :class:`~phantomkit.phantom_processor.PhantomProcessor` class orchestrates +phantom QC for a single session end-to-end without the pipeline scaffolding: .. code-block:: python - >>> from pydra.utils import get_fields - >>> from phantomkit.analyses.vial_signal import VialSignalAnalysis - >>> vial_signal_wf = VialSignalAnalysis( - input_image="/path/to/input/image.nii.gz", - template_dir="/path/to/template/dir", - rotation_library_file="/path/to/rotation/library/file.txt" - ) - >>> outputs = vial_signal_wf(cache_root="/path/to/cache/root") - >>> print("\n".join(f.name for f in outputs.metrics_dir.contents)) + from phantomkit.phantom_processor import PhantomProcessor + + processor = PhantomProcessor( + template_dir="/templates/SPIRIT", + output_base_dir="/results", + rotation_library_file="/templates/rotations.txt", + ) + results = processor.process_session("/data/session01/T1.nii.gz") + + print(results["metrics_dir"]) + print(results["vial_dir"]) + print(results["images_template_space_dir"]) + +All ``*.nii.gz`` files in the same directory as the input image are treated as +additional contrast images. + + +Python API — pydra workflow +---------------------------- + +For finer-grained control, :class:`~phantomkit.analyses.vial_signal.VialSignalAnalysis` +exposes the same processing as a pydra workflow, enabling lazy execution and +caching: + +.. code-block:: python + + from phantomkit.analyses.vial_signal import VialSignalAnalysis + + wf = VialSignalAnalysis( + input_image="/data/session01/T1.nii.gz", + template_dir="/templates/SPIRIT", + rotation_library_file="/templates/rotations.txt", + output_base_dir="/results", + ) + outputs = wf(cache_root="/data/pydra-cache") + + print(outputs.metrics_dir) + +Run ``phantomkit list`` to discover all available protocol workflows. Batch +processing is available via :class:`~phantomkit.analyses.vial_signal.VialSignalAnalysisBatch`. + + +Output structure +----------------- + +A completed session produces: + +.. code-block:: text + + /results/session01/ + metrics/ + session01__mean_matrix.csv + session01__std_matrix.csv + session01__PLOTmeanstd.png + session01_ir_map_PLOTmeanstd_T1mapping.png (if IR contrasts present) + session01_TE_map_PLOTmeanstd_TEmapping.png (if TE contrasts present) + vial_segmentations/ + A.nii.gz … T.nii.gz + images_template_space/ + .nii.gz … + TemplatePhantom_ScannerSpace.nii.gz + +When called through the pipeline, DWI series each get their own subdirectory: + +.. code-block:: text + + /results/session01/ + / + DWI_preproc_biascorr.mif.gz + ADC.nii.gz + FA.nii.gz + T1_in_DWI_space.nii.gz + metrics/ + vial_segmentations/ + native_contrasts_staging/ + metrics/ + vial_segmentations/ + images_template_space/ \ No newline at end of file diff --git a/phantomkit/cli.py b/phantomkit/cli.py index 126a236..436df62 100644 --- a/phantomkit/cli.py +++ b/phantomkit/cli.py @@ -124,7 +124,7 @@ def _build_command(slug: str, single_cls, batch_cls) -> click.Command: def _callback(**kwargs): input_val: str = kwargs.pop("input") - plugin: str = kwargs.pop("plugin") + worker: str = kwargs.pop("worker") pattern: str = kwargs.pop("pattern") # Convert click tuples (from multiple=True) to lists wf_kwargs = { @@ -161,7 +161,7 @@ def _callback(**kwargs): from pydra.engine import Submitter - with Submitter(plugin=plugin) as sub: + with Submitter(worker=worker) as sub: sub(wf) logger.info("Done.") @@ -189,11 +189,11 @@ def _callback(**kwargs): # Common options click_params += [ click.Option( - ["--plugin"], + ["--worker"], type=click.Choice(["cf", "serial"]), default="cf", show_default=True, - help="Pydra execution plugin.", + help="Pydra worker type.", ), click.Option( ["--pattern"], @@ -301,5 +301,176 @@ def _register_plot_commands() -> None: _register_plot_commands() -if __name__ == "__main__": - main() +# --------------------------------------------------------------------------- +# Pipeline command — end-to-end phantom QC + DWI orchestrator +# --------------------------------------------------------------------------- + + +@main.command("pipeline") +@click.option( + "--input-dir", + required=True, + help="Root directory containing acquisition subdirectories (DICOM folders).", +) +@click.option( + "--output-dir", + required=True, + help="Top-level output directory. All results are written here.", +) +@click.option( + "--phantom", + required=True, + help="Phantom name, e.g. SPIRIT. Used to locate template_data//.", +) +@click.option( + "--denoise-degibbs", + is_flag=True, + default=False, + help="Apply dwidenoise + mrdegibbs before preprocessing.", +) +@click.option( + "--gradcheck", + is_flag=True, + default=False, + help="Run dwigradcheck to verify gradient orientations.", +) +@click.option( + "--nocleanup", + is_flag=True, + default=False, + help="Keep DWI tmp/ intermediate directories.", +) +@click.option( + "--readout-time", + type=float, + default=None, + help="Override TotalReadoutTime (seconds) for dwifslpreproc.", +) +@click.option( + "--eddy-options", + type=str, + default=None, + help="Override FSL eddy options string.", +) +@click.option( + "--dry-run", + is_flag=True, + default=False, + help="Plan and print commands; do not execute any processing.", +) +def run_pipeline( + input_dir, + output_dir, + phantom, + denoise_degibbs, + gradcheck, + nocleanup, + readout_time, + eddy_options, + dry_run, +): + """Run the end-to-end phantom QC + DWI processing pipeline. + + Orchestrates DWI preprocessing (Stage 1), phantom QC in DWI space + (Stage 2), and phantom QC on native T1/IR/TE contrasts (Stage 3). + Stages 1 and 3 run in parallel; Stage 2 follows Stage 1. + + INPUT_DIR should contain acquisition subdirectories (DICOM folders). + """ + import sys + from pathlib import Path + from phantomkit.pipeline import ( + scan_input_dir, + validate_inputs, + run_stage1, + run_stage2, + run_stage3, + print_header, + TEMPLATE_DATA_ROOT, + ROTATIONS_FILE, + ) + import concurrent.futures + import threading + + # Build a minimal args-like namespace so validate_inputs() can be reused + class _Args: + pass + + _args = _Args() + _args.input_dir = input_dir + _args.output_dir = output_dir + _args.phantom = phantom + + validate_inputs(_args) + + input_path = Path(input_dir).resolve() + output_path = Path(output_dir).resolve() + template_dir = TEMPLATE_DATA_ROOT / phantom + + output_path.mkdir(parents=True, exist_ok=True) + + dwi_cfg = { + "denoise_degibbs": denoise_degibbs, + "gradcheck": gradcheck, + "nocleanup": nocleanup, + "readout_time": readout_time, + "eddy_options": eddy_options, + } + + print_header("Input Directory Scan") + scan_info = scan_input_dir(input_path) + + run_stage1_flag = scan_info["has_dwi"] + run_stage3_flag = bool(scan_info["t1_dirs"]) and ( + scan_info["has_native_contrasts"] or not run_stage1_flag + ) + + _print_lock = threading.Lock() + + def _locked_header(title): + with _print_lock: + print_header(title) + + dwi_output_dirs = [] + stage1_error = stage3_error = None + + def _s1(): + nonlocal dwi_output_dirs, stage1_error + try: + if run_stage1_flag: + dwi_output_dirs[:] = run_stage1( + input_path, output_path, dwi_cfg, dry_run + ) + else: + _locked_header("STAGE 1 — DWI Processing") + print(" Skipped: no DWI acquisitions found.\n") + except Exception as exc: + stage1_error = exc + + def _s3(): + nonlocal stage3_error + try: + if run_stage3_flag: + run_stage3(input_path, output_path, template_dir, scan_info, dry_run) + else: + _locked_header("STAGE 3 — Phantom QC on Native Contrasts") + print(" Skipped.\n") + except Exception as exc: + stage3_error = exc + + with concurrent.futures.ThreadPoolExecutor(max_workers=2) as executor: + concurrent.futures.wait([executor.submit(_s1), executor.submit(_s3)]) + + if stage1_error: + raise click.ClickException(f"Stage 1 failed: {stage1_error}") + if stage3_error: + raise click.ClickException(f"Stage 3 failed: {stage3_error}") + + if run_stage1_flag: + run_stage2(dwi_output_dirs, output_path, template_dir, dry_run) + else: + print_header("STAGE 2 — Phantom QC in DWI Space") + print(" Skipped: Stage 1 did not run.\n") + + print_header("Pipeline Complete") + print(f" All outputs written to: {output_path}\n") diff --git a/phantomkit/dwi_processing.py b/phantomkit/dwi_processing.py new file mode 100644 index 0000000..b95da37 --- /dev/null +++ b/phantomkit/dwi_processing.py @@ -0,0 +1,2067 @@ +#!/usr/bin/env python3 +""" +DWI processing and tensor metrics pipeline (Pydra 0.23+) + +Usage: + python Functions/pydra_dwi_processing.py --config pydra_dwi_processing.yaml + python Functions/pydra_dwi_processing.py --scans-dir /path/to/scans --output-dir /path/to/out [options] + +Options: + --nocleanup Keep the tmp/ intermediate file directory after processing + (default: tmp/ is removed once final outputs are copied) + +Output structure: + / + / + ADC.nii.gz + FA.nii.gz + DWI_preproc_biascorr.mif.gz (or DWI_denoise_gibbs_preproc_biascorr.mif.gz) + T1_in_DWI_space.nii.gz + tmp/ (intermediate files — removed unless --nocleanup) +""" + +import argparse +import json +import os +import re +import shutil +import subprocess +from pathlib import Path + +from pydra.compose import python, workflow +from pydra.engine import Submitter +import yaml + + +# ============================================================================= +# Utility helpers +# ============================================================================= + + +def run_cmd(cmd: list, cwd: str = None): + """Run a shell command, raising on failure.""" + print(f" >> {' '.join(str(c) for c in cmd)}") + subprocess.run([str(c) for c in cmd], check=True, cwd=cwd) + + +def sanitise_name(name: str) -> str: + """Convert a folder name to a valid Python identifier for Pydra task names.""" + return re.sub(r"[^a-zA-Z0-9_]", "_", name) + + +def strip_series_number(name: str) -> str: + """Strip leading series number, e.g. '87-ep2d_diff_...' -> 'ep2d_diff_...'""" + return re.sub(r"^\d+-", "", name) + + +def get_series_number(name: str) -> int: + """Extract leading series number from folder name, e.g. '87-ep2d...' -> 87.""" + m = re.match(r"^(\d+)-", name) + return int(m.group(1)) if m else 0 + + +def common_prefix_len(a: str, b: str) -> int: + i = 0 + while i < len(a) and i < len(b) and a[i] == b[i]: + i += 1 + return i + + +def get_nvols(path: str) -> int: + """Return number of volumes using mrinfo.""" + result = subprocess.run( + ["mrinfo", path, "-ndim"], capture_output=True, text=True, check=True + ) + ndim = int(result.stdout.strip()) + if ndim == 3: + return 1 + result2 = subprocess.run( + ["mrinfo", path, "-size"], capture_output=True, text=True, check=True + ) + return int(result2.stdout.strip().split()[-1]) + + +def read_bvals(bval_path: str) -> list: + """Read bval file and return list of float values. Returns [] if missing/empty.""" + try: + with open(bval_path) as f: + content = f.read().strip() + if not content: + return [] + return [float(v) for v in content.split()] + except Exception: + return [] + + +def has_nonzero_bvals(bval_path: str) -> bool: + """Return True if bval file exists, is non-empty, and contains non-zero values.""" + vals = read_bvals(bval_path) + return bool(vals) and any(v > 10 for v in vals) # threshold >10 to ignore rounding + + +def all_zero_bvals(bval_path: str) -> bool: + """Return True if bval file exists and all values are effectively zero.""" + vals = read_bvals(bval_path) + return bool(vals) and all(v <= 10 for v in vals) + + +def get_readout_time(json_path: str, fallback: float = 0.0342002) -> float: + """Extract TotalReadoutTime from dcm2niix JSON sidecar.""" + try: + with open(json_path) as f: + data = json.load(f) + val = data.get("TotalReadoutTime") + if val is not None: + return float(val) + except Exception: + pass + print(f" Warning: TotalReadoutTime not found in {json_path} -- using {fallback}") + return fallback + + +def get_pe_from_json(json_path: str) -> tuple: + """ + (#9) Fallback: read PhaseEncodingDirection from JSON sidecar and map to + MRtrix3-compatible direction string. + Returns (pe_dir, rpe_dir) or raises ValueError if not found/unrecognised. + """ + mapping = { + "j-": ("AP", "PA"), + "j": ("PA", "AP"), + "i": ("LR", "RL"), + "i-": ("RL", "LR"), + "k": ("SI", "IS"), + "k-": ("IS", "SI"), + } + try: + with open(json_path) as f: + data = json.load(f) + ped = data.get("PhaseEncodingDirection", "").strip() + if ped in mapping: + return mapping[ped] + except Exception: + pass + raise ValueError(f"PhaseEncodingDirection not found or unrecognised in {json_path}") + + +def detect_pe_direction(folder_name: str) -> tuple: + """ + Return (pe_dir, rpe_dir) from a folder name. + Raises ValueError if not found — caller should fall back to get_pe_from_json (#9). + """ + pairs = [ + (r"_AP(_|$)", "AP", "PA"), + (r"_A_P(_|$)", "AP", "PA"), + (r"_PA(_|$)", "PA", "AP"), + (r"_P_A(_|$)", "PA", "AP"), + (r"_LR(_|$)", "LR", "RL"), + (r"_L_R(_|$)", "LR", "RL"), + (r"_RL(_|$)", "RL", "LR"), + (r"_R_L(_|$)", "RL", "LR"), + (r"_SI(_|$)", "SI", "IS"), + (r"_S_I(_|$)", "SI", "IS"), + (r"_IS(_|$)", "IS", "SI"), + (r"_I_S(_|$)", "IS", "SI"), + ] + for pattern, pe, rpe in pairs: + if re.search(pattern, folder_name, re.IGNORECASE): + return pe, rpe + raise ValueError( + f"Could not determine phase encoding direction from folder name: {folder_name}" + ) + + +# ============================================================================= +# Directory scanning constants +# ============================================================================= + +FWD_DIRS = ["AP", "LR", "SI", "A_P", "L_R", "S_I"] +RPE_DIRS = ["PA", "RL", "IS", "P_A", "R_L", "I_S"] +ALL_PE_DIRS = FWD_DIRS + RPE_DIRS + + +def has_terminal_pe(name: str, directions: list) -> bool: + """ + Returns True if a PE direction tag appears as the final token in the folder + name (nothing follows, or only a single short alphanumeric token with no + underscores). This identifies dedicated b0 PE images vs full DWI series. + """ + for d_ in directions: + for pat in [rf"_{d_}$", rf"_{d_}_[A-Za-z0-9]{{1,10}}$"]: + if re.search(pat, name, re.IGNORECASE): + return True + # Underscore-separated variant (e.g. A_P, P_A) + if len(d_) == 3 and d_[1] == "_": + for pat in [ + rf"_{d_[0]}_{d_[2]}$", + rf"_{d_[0]}_{d_[2]}_[A-Za-z0-9]{{1,10}}$", + ]: + if re.search(pat, name, re.IGNORECASE): + return True + return False + + +def get_acq_stem_and_suffix(folder_name: str) -> tuple: + """ + Split acquisition name into (prefix, suffix) around the PE direction tag. + Both must match exactly for two series to be considered valid AP/PA partners + (#14 — prevents pairing of MONO/BIPOLAR variants or different parameters). + + e.g. '12-ep2d_diff__FREE30DIR_b1000_x7b0_Ghost_R_AP_ESpt54_BW2264_BIPOLAR' + -> prefix: 'ep2d_diff__FREE30DIR_b1000_x7b0_Ghost_R' + suffix: 'ESpt54_BW2264_BIPOLAR' + + Returns (full_name_no_series, "", "") if no PE direction tag found. + """ + name = strip_series_number(folder_name) + # Try longer variants first (A_P before AP) to avoid partial matches + all_dirs_ordered = [ + "A_P", + "P_A", + "L_R", + "R_L", + "S_I", + "I_S", + "AP", + "PA", + "LR", + "RL", + "SI", + "IS", + ] + for d in all_dirs_ordered: + m = re.search(rf"_({re.escape(d)})(_|$)", name, re.IGNORECASE) + if m: + prefix = name[: m.start()] + # suffix is everything after the direction tag (and the following _) + suffix_start = m.end() + suffix = name[suffix_start:] if suffix_start < len(name) else "" + return prefix, suffix + return name, "" + + +# ============================================================================= +# DICOM conversion helper (used in planning stage) +# ============================================================================= + + +def convert_dicom_to_nii(dicom_dir: str, out_dir: str) -> dict: + """ + Run dcm2niix on a DICOM directory, returning a dict with paths to: + nii, json, bvec, bval + Returns empty strings for missing files. + """ + os.makedirs(out_dir, exist_ok=True) + subprocess.run( + ["dcm2niix", "-o", out_dir, "-f", "%p", "-z", "y", dicom_dir], + check=True, + capture_output=True, + ) + niis = sorted(Path(out_dir).glob("*.nii.gz")) + if not niis: + raise FileNotFoundError(f"No NIfTI produced from DICOM: {dicom_dir}") + nii = str(niis[0]) + base = nii.replace(".nii.gz", "") + return { + "nii": nii, + "json": base + ".json" if Path(base + ".json").exists() else "", + "bvec": base + ".bvec" if Path(base + ".bvec").exists() else "", + "bval": base + ".bval" if Path(base + ".bval").exists() else "", + } + + +# ============================================================================= +# Directory scanning and candidate classification +# ============================================================================= + + +def scan_directory(scans_dir: str) -> dict: + """ + First pass: classify subdirectories into t1_dirs, candidate_dwi (all series + containing _diff_ or _DWI_), and ignored. PE classification is PROVISIONAL + at this stage — final assignment happens after DICOM conversion in + classify_candidates(). + """ + scans_path = Path(scans_dir) + t1_dirs = [] + candidate_dwi = [] # all _diff_ / _DWI_ series, including PE images + ignored = [] + + for d in sorted(scans_path.iterdir()): + if not d.is_dir(): + continue + name = d.name + + # T1 + if re.search(r"t1", name, re.IGNORECASE): + t1_dirs.append(str(d)) + continue + + # Skip scanner-derived maps + if re.search(r"_(ADC|FA)$", name, re.IGNORECASE): + ignored.append(str(d)) + continue + + # Only consider DWI-like acquisitions + if not re.search(r"(_diff_|_DWI_)", name, re.IGNORECASE): + ignored.append(str(d)) + continue + + candidate_dwi.append(str(d)) + + if not t1_dirs: + raise ValueError(f"Could not identify any T1 directory in {scans_dir}") + if not candidate_dwi: + raise ValueError(f"Could not identify any DWI directories in {scans_dir}") + + return { + "t1_dirs": t1_dirs, + "candidate_dwi": candidate_dwi, + "ignored": ignored, + } + + +def convert_all_candidates(candidate_dwi: list, output_dir: str) -> dict: + """ + (#18 ordering step 2) Convert all candidate DWI DICOMs upfront. + Returns a dict: {dicom_dir: {nii, json, bvec, bval}} + """ + conversions = {} + for dicom_dir in candidate_dwi: + name = Path(dicom_dir).name + conv_dir = str(Path(output_dir) / name / "tmp" / "dwi_nii") + print(f" Converting: {name}") + try: + result = convert_dicom_to_nii(dicom_dir, conv_dir) + conversions[dicom_dir] = result + except Exception as e: + print(f" WARNING: DICOM conversion failed for {name}: {e}") + conversions[dicom_dir] = None + return conversions + + +def classify_candidates(candidate_dwi: list, conversions: dict) -> dict: + """ + (#18) Final classification of candidate DWI series into: + - dwi_dirs: full DWI series eligible for tensor processing (non-zero bvals) + - fwd_pe_dirs: forward PE correction images (all b-values zero) + - rpe_dirs: reverse PE correction images (all b-values zero) + - pending_fwd: series with non-zero bvals and a fwd PE direction tag anywhere + in name — reclassified to dwi_dirs or fwd_pe_dirs after pairing + - pending_rpe: same but reverse PE direction tag + - skipped: list of (path, reason) tuples + + Classification rules (applied in order after DICOM conversion): + 1. No bvec/bval files → skip (#16) + 2. Empty bval file → skip (#16) + 3. All b-values zero → PE correction image candidate, classified by PE + direction tag found ANYWHERE in folder name (direction tag position + does not affect classification — only bval content matters) + 4. Has non-zero b-values AND a PE direction tag → pending (reclassified + after rpe_all pairing check) + 5. Has non-zero b-values, no PE direction tag → dwi_dirs + """ + dwi_dirs = [] + fwd_pe_dirs = [] + rpe_dirs = [] + skipped = [] + pending_fwd = [] # non-zero bvals + fwd PE tag — needs pairing check + pending_rpe = [] # non-zero bvals + rpe PE tag — needs pairing check + + for dicom_dir in candidate_dwi: + name = Path(dicom_dir).name + conv = conversions.get(dicom_dir) + + if conv is None: + skipped.append((dicom_dir, "DICOM conversion failed")) + continue + + bvec = conv["bvec"] + bval = conv["bval"] + + # #16 — check bvec/bval exist and are non-empty + if not bvec or not bval or not Path(bvec).exists() or not Path(bval).exists(): + skipped.append( + (dicom_dir, "no bvec/bval files — not a raw DWI acquisition") + ) + continue + if not read_bvals(bval): + skipped.append((dicom_dir, "empty bval file — not a raw DWI acquisition")) + continue + + is_b0_only = all_zero_bvals(bval) + + # Detect PE direction tag anywhere in the folder name (not just terminal) + has_fwd_tag = any( + re.search(rf"_{re.escape(d)}(_|$)", name, re.IGNORECASE) for d in FWD_DIRS + ) + has_rpe_tag = any( + re.search(rf"_{re.escape(d)}(_|$)", name, re.IGNORECASE) for d in RPE_DIRS + ) + + if is_b0_only: + # All b-values zero: classify as PE correction image by direction tag + if has_fwd_tag: + fwd_pe_dirs.append(dicom_dir) + elif has_rpe_tag: + rpe_dirs.append(dicom_dir) + else: + skipped.append( + ( + dicom_dir, + "all b-values zero and no phase encoding direction in folder name — " + "cannot use for correction", + ) + ) + else: + # Has non-zero b-values: pending if PE tag present, dwi_dirs otherwise + if has_fwd_tag: + pending_fwd.append(dicom_dir) + elif has_rpe_tag: + pending_rpe.append(dicom_dir) + else: + dwi_dirs.append(dicom_dir) + + return { + "dwi_dirs": dwi_dirs, + "fwd_pe_dirs": fwd_pe_dirs, + "rpe_dirs": rpe_dirs, + "pending_fwd": pending_fwd, + "pending_rpe": pending_rpe, + "skipped": skipped, + } + + +# ============================================================================= +# AP/PA pairing +# ============================================================================= + + +def match_ap_pa_pairs( + dwi_dirs: list, pending_fwd: list, pending_rpe: list, conversions: dict +) -> tuple: + """ + (#14, #18) Detect AP/PA (rpe_all/rpe_split) pairs across dwi_dirs, + pending_fwd, and pending_rpe. + + Pairing requires: + - Same acquisition prefix (everything before the PE direction tag) + - Same acquisition suffix (everything after the PE direction tag) — #14 + - Opposing PE directions + + Pending series with non-zero bvals that find a partner are moved to dwi_dirs. + Pending series that find no partner are moved to fwd_pe_dirs / rpe_dirs. + + Returns: + (dwi_dirs, fwd_pe_dirs, rpe_dirs, rpe_all_map, tie_warnings) + rpe_all_map: {fwd_path: rpe_path} + tie_warnings: list of (pe_image_path, [matched_dwi_paths]) for #20 + """ + fwd_set = {"AP", "LR", "SI"} + rpe_set = {"PA", "RL", "IS"} + + # Build stem map over all candidates (dwi + pending) + all_candidates = dwi_dirs + pending_fwd + pending_rpe + stem_map = {} # (prefix, suffix) -> list of (path, pe_dir) + for d in all_candidates: + name = Path(d).name + try: + pe, _ = detect_pe_direction(name) + except ValueError: + pe = None + prefix, suffix = get_acq_stem_and_suffix(name) + key = (prefix, suffix) + stem_map.setdefault(key, []).append((d, pe)) + + rpe_all_map = {} + tie_warnings = [] + paired_paths = set() + + for (prefix, suffix), entries in stem_map.items(): + fwds = [(p, pe) for p, pe in entries if pe in fwd_set] + rpes = [(p, pe) for p, pe in entries if pe in rpe_set] + + if not (fwds and rpes): + continue + + for fwd_path, _ in fwds: + fwd_num = get_series_number(Path(fwd_path).name) + # Find closest RPE by series number + best_rpe_path = min( + rpes, key=lambda x: abs(get_series_number(Path(x[0]).name) - fwd_num) + )[0] + rpe_all_map[fwd_path] = best_rpe_path + paired_paths.add(fwd_path) + paired_paths.add(best_rpe_path) + + # #20 — check for tie (multiple equal-prefix matches) + if len(fwds) > 1 or len(rpes) > 1: + tie_warnings.append((fwd_path, [p for p, _ in rpes])) + + # Move paired pending series into dwi_dirs; unpaired pending into PE dirs + final_dwi = list(dwi_dirs) + final_fwd = list(conversions.get("fwd_pe_dirs", [])) + final_rpe = list(conversions.get("rpe_dirs", [])) + + for p in pending_fwd: + if p in paired_paths: + final_dwi.append(p) + else: + final_fwd.append(p) + + for p in pending_rpe: + if p in paired_paths: + # PA partner of a paired series — recorded in rpe_all_map, not in dwi_dirs + pass + else: + final_rpe.append(p) + + return sorted(final_dwi), final_fwd, final_rpe, rpe_all_map, tie_warnings + + +# ============================================================================= +# T1 and PE assignment helpers +# ============================================================================= + + +def assign_t1(dwi_name: str, t1_dirs: list) -> str: + """ + Assign the nearest preceding T1 by series number. + Falls back to the first T1 if none precedes the DWI. + """ + dwi_num = get_series_number(dwi_name) + best = None + best_num = -1 + for t1 in t1_dirs: + t1_num = get_series_number(Path(t1).name) + if t1_num < dwi_num and t1_num > best_num: + best_num = t1_num + best = t1 + return best if best is not None else t1_dirs[0] + + +def build_pe_assignment_map( + dwi_dirs: list, fwd_pe_dirs: list, rpe_dirs: list, rpe_all_map: dict +) -> dict: + """ + Build a map of {dwi_dir: {"rpe": rpe_dir_or_None, "fwd": fwd_dir_or_None}} + for all DWI series not already in rpe_all_map. + + A single b0 PE correction image is assigned to ALL DWI series sharing the + longest common prefix with it — not just the nearest one. This correctly + handles the case where one b0 image is intended to correct an entire block + of acquisitions that vary only in a trailing parameter (e.g. bandwidth sweep). + + If multiple PE images share the same best prefix length against a DWI block, + the one with the nearest series number to the block is chosen, with a warning. + + Priority: rpe_all (already in rpe_all_map) > rpe_pair > rpe_none. + """ + unpaired = [d for d in dwi_dirs if d not in rpe_all_map] + assignment = {d: {"rpe": None, "fwd": None, "tie_warning": False} for d in unpaired} + + def assign_pe_images(pe_dirs, key): + """ + For each PE image, find all DWI series sharing the longest common prefix. + Assign this PE image to that entire group. + """ + # Build prefix lengths: pe_dir -> {dwi_dir: prefix_len} + pe_scores = {} + for pe_dir in pe_dirs: + pe_stem = strip_series_number(Path(pe_dir).name) + scores = {} + for dwi_dir in unpaired: + dwi_stem = strip_series_number(Path(dwi_dir).name) + scores[dwi_dir] = common_prefix_len(dwi_stem, pe_stem) + pe_scores[pe_dir] = scores + + # For each PE image, identify the group of DWI series with max prefix length + # Only assign if prefix length is meaningful (>10 chars) + pe_groups = {} # pe_dir -> set of dwi_dirs it should serve + for pe_dir, scores in pe_scores.items(): + max_len = max(scores.values()) if scores else 0 + if max_len <= 0: + continue # no common prefix at all — skip + group = {d for d, l in scores.items() if l == max_len} + pe_groups[pe_dir] = group + + # For each DWI series, find which PE image has the best (longest) prefix + # If multiple PE images tie for a DWI series, pick nearest by series number + dwi_to_best_pe = {} + for dwi_dir in unpaired: + candidates = [] + best_len = 0 + for pe_dir, group in pe_groups.items(): + if dwi_dir in group: + pe_stem = strip_series_number(Path(pe_dir).name) + dwi_stem = strip_series_number(Path(dwi_dir).name) + length = common_prefix_len(dwi_stem, pe_stem) + if length > best_len: + best_len = length + candidates = [pe_dir] + elif length == best_len: + candidates.append(pe_dir) + if not candidates: + continue + tie = len(candidates) > 1 + dwi_num = get_series_number(Path(dwi_dir).name) + best_pe = min( + candidates, key=lambda p: abs(get_series_number(Path(p).name) - dwi_num) + ) + dwi_to_best_pe[dwi_dir] = (best_pe, tie) + + # Now assign: use the winning PE image for each DWI, extend to whole group + # Group DWI series that share the same winning PE image + pe_to_dwi_group = {} + for dwi_dir, (pe_dir, tie) in dwi_to_best_pe.items(): + pe_to_dwi_group.setdefault(pe_dir, []).append((dwi_dir, tie)) + + for pe_dir, dwi_list in pe_to_dwi_group.items(): + any_tie = any(tie for _, tie in dwi_list) + for dwi_dir, _ in dwi_list: + if assignment[dwi_dir][key] is None: + assignment[dwi_dir][key] = pe_dir + assignment[dwi_dir]["tie_warning"] = any_tie + + assign_pe_images(rpe_dirs, "rpe") + assign_pe_images(fwd_pe_dirs, "fwd") + + return assignment + + +# ============================================================================= +# PE table header-based mode inference +# ============================================================================= + +# Mapping from (i, j, k) unit vector (rounded) to direction string. +# MRtrix convention: j- = AP, j = PA, i = LR, i- = RL, k = SI, k- = IS +_PETABLE_VEC_TO_DIR = { + (0, 1, 0): "PA", + (0, -1, 0): "AP", + (1, 0, 0): "LR", + (-1, 0, 0): "RL", + (0, 0, 1): "SI", + (0, 0, -1): "IS", +} + + +def _petable_vec_to_dir(vec: tuple) -> str: + """ + Map a (i, j, k) unit vector from mrinfo -petable to a direction string. + Returns 'UNKNOWN' if the vector doesn't match any known direction. + """ + rounded = tuple(round(float(v)) for v in vec) + return _PETABLE_VEC_TO_DIR.get(rounded, f"UNKNOWN{rounded}") + + +def get_petable_mode( + nii: str, bvec: str, bval: str, json_path: str, tmp_dir: str, has_rpe: bool +) -> tuple: + """ + Infer the preproc mode for a DWI series by reading the phase encoding table + embedded in the MIF header produced by mrconvert. + + Steps: + 1. mrconvert NIfTI+bvec+bval+JSON → temporary MIF + 2. mrinfo -petable on the MIF + 3. Parse rows to extract PE directions per volume + 4. Apply mode inference rules: + - All same direction, no RPE available → rpe_none + - All same direction, RPE available → rpe_pair + - Mixed directions, equal counts → rpe_all + - Mixed directions, unequal counts → rpe_split + + Returns: + (mode, directions_found, petable_rows) + mode: one of 'rpe_none', 'rpe_pair', 'rpe_all', 'rpe_split' + directions_found: list of direction strings per volume (e.g. ['AP','AP','PA']) + petable_rows: raw petable as list of lists (for diagnostics) + + Raises RuntimeError if mrconvert or mrinfo fails. + """ + os.makedirs(tmp_dir, exist_ok=True) + tmp_mif = str(Path(tmp_dir) / "petable_check.mif") + + # Build mrconvert command — include JSON sidecar if available (embeds PE info) + cmd = ["mrconvert", nii, tmp_mif, "-fslgrad", bvec, bval, "-force", "-quiet"] + if json_path and Path(json_path).exists(): + cmd += ["-json_import", json_path] + + try: + subprocess.run( + [str(c) for c in cmd], check=True, capture_output=True, text=True + ) + except subprocess.CalledProcessError as e: + raise RuntimeError(f"mrconvert failed during petable check: {e.stderr.strip()}") + + # Read petable — columns: i j k readout_time (one row per volume) + try: + result = subprocess.run( + ["mrinfo", tmp_mif, "-petable"], capture_output=True, text=True, check=True + ) + except subprocess.CalledProcessError as e: + raise RuntimeError(f"mrinfo -petable failed: {e.stderr.strip()}") + finally: + # Always clean up the temporary MIF + try: + Path(tmp_mif).unlink(missing_ok=True) + except Exception: + pass + + raw = result.stdout.strip() + if not raw: + raise RuntimeError( + "mrinfo -petable returned no output — PE information " + "may not be embedded in this image." + ) + + petable_rows = [] + directions = [] + for line in raw.splitlines(): + parts = line.split() + if len(parts) < 3: + continue + vec = (parts[0], parts[1], parts[2]) + petable_rows.append(parts) + directions.append(_petable_vec_to_dir(vec)) + + if not directions: + raise RuntimeError("Could not parse any rows from mrinfo -petable output.") + + unique_dirs = set(directions) + + if len(unique_dirs) == 1: + # All volumes have the same PE direction + mode = "rpe_pair" if has_rpe else "rpe_none" + else: + # Mixed PE directions in a single series + counts = {d: directions.count(d) for d in unique_dirs} + if len(set(counts.values())) == 1: + mode = "rpe_all" # equal counts + else: + mode = "rpe_split" # unequal counts + + return mode, directions, petable_rows + + +# ============================================================================= +# Plan construction +# ============================================================================= + + +def resolve_pe_direction(dwi_name: str, dwi_json: str) -> tuple: + """ + (#9) Try to determine PE direction from folder name first. + Fall back to JSON sidecar if not found. Returns (pe_dir, rpe_dir, source) + where source is 'folder_name' or 'json_sidecar'. + Raises ValueError if neither source yields a direction. + """ + try: + pe_dir, rpe_dir = detect_pe_direction(dwi_name) + return pe_dir, rpe_dir, "folder_name" + except ValueError: + pass + + if dwi_json and Path(dwi_json).exists(): + try: + pe_dir, rpe_dir = get_pe_from_json(dwi_json) + return pe_dir, rpe_dir, "json_sidecar" + except ValueError: + pass + + raise ValueError( + f"Could not determine phase encoding direction for {dwi_name} " + f"from folder name or JSON sidecar." + ) + + +def plan_workflow( + dwi_dir: str, + t1_dirs: list, + fwd_pe_dirs: list, + rpe_dirs: list, + rpe_all_map: dict, + pe_assignment_map: dict, + conversions: dict, + cfg: dict, +) -> dict: + """ + Build a complete plan dict for a single DWI series. + All DICOM conversion has already been done — we use results from conversions. + PE correction image assignments come from pe_assignment_map (pre-built so that + one b0 image can be shared across a whole block of DWI series). + """ + dwi_name = Path(dwi_dir).name + conv = conversions[dwi_dir] + dwi_nii = conv["nii"] + dwi_json = conv["json"] + dwi_bvec = conv["bvec"] + dwi_bval = conv["bval"] + dwi_nvols = get_nvols(dwi_nii) + + out_dir = str(Path(cfg["output_dir"]) / dwi_name) + tmp_dir = str(Path(out_dir) / "tmp") + os.makedirs(tmp_dir, exist_ok=True) + + # T1 assignment + t1_dir = assign_t1(dwi_name, t1_dirs) + + # PE direction (#9) + warnings = [] + try: + pe_dir, rpe_dir, pe_source = resolve_pe_direction(dwi_name, dwi_json) + except ValueError as e: + raise ValueError(str(e)) + + if pe_source == "json_sidecar": + warnings.append( + f"Phase encoding direction ({pe_dir}) inferred from JSON sidecar " + f"— not found in folder name. Please verify." + ) + + # RPE assignment — priority: rpe_all > rpe_pair > rpe_none + rpe_all_partner = rpe_all_map.get(dwi_dir) + + rpe_nii = None + rpe_nvols = 0 + fwd_pe_nii = None + fwd_pe_dir_used = None + rpe_dir_path = None + + if rpe_all_partner: + # rpe_all: full-volume PA partner — most preferred + rpe_dir_path = rpe_all_partner + rpe_conv = conversions.get(rpe_all_partner) + if rpe_conv: + rpe_nii = rpe_conv["nii"] + rpe_nvols = get_nvols(rpe_nii) + + if rpe_nvols == dwi_nvols: + warnings.append( + f"{dwi_name} and {Path(rpe_all_partner).name} have equal volumes " + f"({dwi_nvols}). Assuming rpe_all (full repeats). If directions are " + f"split across PE directions, set mode: rpe_split in config YAML." + ) + fwd_pe_dir_used = None # b0 PE images not needed for rpe_all + + else: + # Use pe_assignment_map for b0 PE correction (rpe_pair) + assignment = pe_assignment_map.get(dwi_dir, {}) + rpe_dir_path = assignment.get("rpe") + fwd_pe_dir_used = assignment.get("fwd") + tie_warning = assignment.get("tie_warning", False) + + # #20 — tie warning + if tie_warning: + warnings.append( + f"PE correction image matched equally to multiple DWI series. " + f"Assigned by nearest series number. If incorrect, consider adding " + f"a distinguishing suffix (e.g. _repeat) to the second acquisition " + f"block and its correction pair." + ) + + # Ensure RPE conversion is available + if rpe_dir_path: + rpe_conv = conversions.get(rpe_dir_path) + if rpe_conv is None: + rpe_conv_dir = str(Path(tmp_dir) / "rpe_nii") + rpe_conv = convert_dicom_to_nii(rpe_dir_path, rpe_conv_dir) + conversions[rpe_dir_path] = rpe_conv + rpe_nii = rpe_conv["nii"] + rpe_nvols = get_nvols(rpe_nii) + + # Ensure FWD PE conversion is available + if fwd_pe_dir_used: + fwd_conv = conversions.get(fwd_pe_dir_used) + if fwd_conv is None: + fwd_conv_dir = str(Path(tmp_dir) / "fwd_pe_nii") + fwd_conv = convert_dicom_to_nii(fwd_pe_dir_used, fwd_conv_dir) + conversions[fwd_pe_dir_used] = fwd_conv + fwd_pe_nii = fwd_conv["nii"] + + # #17 — orphaned RPE handling (no forward PE partner) + if rpe_dir_path and not fwd_pe_dir_used: + if rpe_nvols == 1: + warnings.append( + f"{Path(rpe_dir_path).name} has no matching forward PE image. " + f"Will use mean b0 from main DWI as forward PE image (rpe_pair)." + ) + elif rpe_nvols > 1: + warnings.append( + f"{Path(rpe_dir_path).name} has no matching forward PE image " + f"and RPE has multiple volumes. Cannot use for correction. " + f"Falling back to rpe_none." + ) + rpe_dir_path = None + rpe_nii = None + rpe_nvols = 0 + + # Determine preproc mode + if rpe_nii is None: + preproc_mode = "rpe_none" + elif rpe_nvols == dwi_nvols: + preproc_mode = "rpe_all" + elif rpe_nvols > 1 and fwd_pe_nii is not None: + preproc_mode = "rpe_split" + else: + preproc_mode = "rpe_pair" + + # ── Header-based mode sanity check ────────────────────────────────────── + # Run mrconvert + mrinfo -petable on the main DWI (and RPE if present) to + # derive the mode independently from the MIF header PE table. If the result + # disagrees with the filename-derived mode, warn and override. + # + # For rpe_all: the check requires a concatenated AP+PA MIF which does not + # exist at planning time. We skip the header check for rpe_all and note it. + petable_tmp = str(Path(tmp_dir) / "petable_check") + if preproc_mode == "rpe_all": + notes_petable = ( + "Header PE check: skipped for rpe_all " + "(requires concatenated AP+PA MIF — verified during processing)" + ) + else: + try: + header_mode, header_dirs, _ = get_petable_mode( + nii=dwi_nii, + bvec=dwi_bvec, + bval=dwi_bval, + json_path=dwi_json, + tmp_dir=petable_tmp, + has_rpe=rpe_nii is not None, + ) + + # Also check RPE series if present + rpe_header_mode = None + if rpe_nii: + rpe_conv_for_check = conversions.get(rpe_dir_path, {}) + try: + rpe_header_mode, _, _ = get_petable_mode( + nii=rpe_nii, + bvec=rpe_conv_for_check.get("bvec", ""), + bval=rpe_conv_for_check.get("bval", ""), + json_path=rpe_conv_for_check.get("json", ""), + tmp_dir=petable_tmp + "_rpe", + has_rpe=False, # checking the RPE series itself in isolation + ) + except RuntimeError as e: + warnings.append( + f"Header check on RPE series could not be completed: {e} " + f"— RPE series assignment retained." + ) + + if header_mode != preproc_mode: + warnings.append( + f"Preproc mode mismatch — filename-derived: {preproc_mode}, " + f"header-derived: {header_mode}. " + f"Overriding with header-derived mode. " + f"PE directions found in header: {sorted(set(header_dirs))}. " + f"Please verify series naming and acquisition protocol." + ) + preproc_mode = header_mode + + notes_petable = ( + f"Header PE check: {header_mode} " + f"(directions: {', '.join(sorted(set(header_dirs)))})" + ) + if rpe_header_mode: + notes_petable += f"; RPE header: {rpe_header_mode}" + + except RuntimeError as e: + warnings.append( + f"Header-based PE mode check could not be completed: {e} " + f"— falling back to filename-derived mode ({preproc_mode})." + ) + notes_petable = "Header PE check: unavailable" + + # #21 — note unused PE images when rpe_all is used + notes = [notes_petable] + if rpe_all_partner: + unused = [d for d in fwd_pe_dirs + rpe_dirs if d != rpe_all_partner] + for d in unused: + notes.append(f"{Path(d).name} — not used (DWI series has rpe_all partner)") + + # Resolve readout time + readout_time_override = cfg.get("readout_time", None) + if readout_time_override is not None: + readout_time = float(readout_time_override) + readout_time_source = "config override" + elif dwi_json and Path(dwi_json).exists(): + readout_time = get_readout_time(dwi_json) + readout_time_source = f"JSON ({Path(dwi_json).name})" + else: + readout_time = 0.0342002 + readout_time_source = "fallback (no JSON found)" + + do_denoise = cfg.get("denoise_degibbs", False) + do_gradcheck = cfg.get("gradcheck", False) + eddy_options = cfg.get("eddy_options", " --slm=linear") + keep_tmp = cfg.get("keep_tmp", False) + + dwi_preproc_name = ( + "DWI_denoise_gibbs_preproc_biascorr.mif.gz" + if do_denoise + else "DWI_preproc_biascorr.mif.gz" + ) + + # T1 output filename depends on whether registration was performed. + t1_output_name = "T1_in_DWI_space.nii.gz" + + rpe_conv = conversions.get(rpe_dir_path) if rpe_dir_path else None + fwd_conv = conversions.get(fwd_pe_dir_used) if fwd_pe_dir_used else None + + return { + "dwi_name": dwi_name, + "dwi_dir": dwi_dir, + "t1_dir": t1_dir, + "fwd_pe_dir": fwd_pe_dir_used, + "rpe_dir_path": rpe_dir_path, + "pe_dir": pe_dir, + "rpe_dir": rpe_dir, + "pe_source": pe_source, + "dwi_nii": dwi_nii, + "dwi_json": dwi_json, + "dwi_bvec": dwi_bvec, + "dwi_bval": dwi_bval, + "dwi_nvols": dwi_nvols, + "rpe_nii": rpe_nii, + "rpe_json": rpe_conv["json"] if rpe_conv else "", + "rpe_bvec": rpe_conv["bvec"] if rpe_conv else "", + "rpe_bval": rpe_conv["bval"] if rpe_conv else "", + "rpe_nvols": rpe_nvols, + "fwd_pe_nii": fwd_pe_nii, + "fwd_pe_json": fwd_conv["json"] if fwd_conv else "", + "fwd_pe_bvec": fwd_conv["bvec"] if fwd_conv else "", + "fwd_pe_bval": fwd_conv["bval"] if fwd_conv else "", + "preproc_mode": preproc_mode, + "readout_time": readout_time, + "readout_time_src": readout_time_source, + "do_denoise": do_denoise, + "do_gradcheck": do_gradcheck, + "eddy_options": eddy_options, + "keep_tmp": keep_tmp, + "dwi_preproc_name": dwi_preproc_name, + "t1_output_name": t1_output_name, + "out_dir": out_dir, + "tmp_dir": tmp_dir, + "warnings": warnings, + "notes": notes, + } + + +def print_plan(plans: list, skipped: list): + """Print workflow summary including all warnings and notes.""" + print("\n" + "=" * 70) + print("WORKFLOW PLAN") + print("=" * 70) + + if skipped: + print("\nSKIPPED SERIES:") + for path, reason in skipped: + print(f" SKIPPED: {Path(path).name}") + print(f" ({reason})") + + for p in plans: + rpe_name = Path(p["rpe_dir_path"]).name if p["rpe_dir_path"] else "not found" + fwd_name = Path(p["fwd_pe_dir"]).name if p["fwd_pe_dir"] else "not used" + + print(f"\nDWI series: {p['dwi_name']}") + print(f" T1: {Path(p['t1_dir']).name}") + + if p["preproc_mode"] == "rpe_all": + print(f" AP/PA pair: {p['dwi_name']} ({p['dwi_nvols']} vols)") + print(f" {rpe_name} ({p['rpe_nvols']} vols)") + print(f" Forward PE b0: not used (rpe_all)") + print(f" Reverse PE b0: not used (rpe_all)") + else: + print(f" Forward PE: {fwd_name}") + print(f" Reverse PE: {rpe_name}") + print(f" DWI volumes: {p['dwi_nvols']}") + print(f" RPE volumes: {p['rpe_nvols']}") + + pe_str = p["pe_dir"] + if p["pe_source"] == "json_sidecar": + pe_str += ( + " (WARNING: inferred from JSON sidecar — not found in folder name)" + ) + print(f" PE direction: {pe_str}") + print(f" Preproc mode: {p['preproc_mode']}") + print(f" Readout time: {p['readout_time']} ({p['readout_time_src']})") + print(f" Denoise/Gibbs: {p['do_denoise']}") + print(f" Gradcheck: {p['do_gradcheck']}") + + if p["do_gradcheck"]: + targets = [f"DWI ({p['pe_dir']})"] + if p["preproc_mode"] == "rpe_all": + targets.append(f"DWI ({p['rpe_dir']})") + elif p["fwd_pe_nii"]: + targets.append("Forward PE b0") + if p["rpe_nii"] and p["preproc_mode"] != "rpe_all": + targets.append("Reverse PE b0") + print(f" Gradcheck on: {', '.join(targets)}") + + if p["do_denoise"] and p["preproc_mode"] == "rpe_all": + print( + f" Denoise on: DWI ({p['pe_dir']}) and DWI ({p['rpe_dir']}) separately" + ) + + print(f" Output: {p['out_dir']}") + + if p["preproc_mode"] == "rpe_none": + print( + f" NOTE: preproc_mode is rpe_none — no distortion correction applied.\n" + f" B0-to-T1 co-registration will still be performed, but EPI\n" + f" distortion may reduce registration accuracy.\n" + f" Consider acquiring a reverse PE image for optimal results." + ) + + for w in p["warnings"]: + print(f" WARNING: {w}") + for n in p["notes"]: + print(f" NOTE: {n}") + + print("\n" + "=" * 70) + print("Proceed? (Ctrl+C to abort)") + print("=" * 70 + "\n") + + +# ============================================================================= +# Pydra tasks +# ============================================================================= + + +@python.define(outputs=["nii", "json_file", "bvec", "bval", "sentinel"]) +def convert_dicoms(dicom_dir: str, out_dir: str) -> tuple[str, str, str, str, str]: + """Convert a DICOM directory to NIfTI. Returns NIfTI + sidecar paths.""" + os.makedirs(out_dir, exist_ok=True) + run_cmd(["dcm2niix", "-o", out_dir, "-f", "%p", "-z", "y", dicom_dir]) + niis = sorted(Path(out_dir).glob("*.nii.gz")) + if not niis: + raise FileNotFoundError(f"No NIfTI files produced in {out_dir}") + nii = str(niis[0]) + base = nii.replace(".nii.gz", "") + return ( + nii, + base + ".json" if Path(base + ".json").exists() else "", + base + ".bvec" if Path(base + ".bvec").exists() else "", + base + ".bval" if Path(base + ".bval").exists() else "", + nii, # sentinel: forces true lazy dependency in downstream tasks + ) + + +@python.define +def convert_to_mif_initial( + nii: str, json_file: str, bvec: str, bval: str, out_path: str +) -> str: + """Step 1 of gradcheck: mrconvert NIfTI -> MIF with original gradients embedded.""" + run_cmd( + ["mrconvert", nii, out_path, "-json_import", json_file, "-fslgrad", bvec, bval] + ) + return out_path + + +@python.define(outputs=["bvec", "bval"]) +def run_gradcheck(mif_path: str, out_dir: str, prefix: str) -> tuple[str, str]: + """Run dwigradcheck on a MIF, export corrected bvec/bval.""" + out_bvec = str(Path(out_dir) / f"{prefix}_corrected.bvec") + out_bval = str(Path(out_dir) / f"{prefix}_corrected.bval") + run_cmd(["dwigradcheck", mif_path, "-export_grad_fsl", out_bvec, out_bval]) + return (out_bvec, out_bval) + + +@python.define +def convert_to_mif_final( + nii: str, + json_file: str, + bvec: str, + bval: str, + out_path: str, + sentinel: str = "", +) -> str: + """mrconvert NIfTI -> MIF with (optionally corrected) gradients.""" + run_cmd( + ["mrconvert", nii, out_path, "-json_import", json_file, "-fslgrad", bvec, bval] + ) + result = subprocess.run( + ["mrinfo", out_path, "-petable"], capture_output=True, text=True + ) + if not result.stdout.strip(): + print(f" Warning: PE table not found in MIF header: {out_path}") + return out_path + + +@python.define +def run_dwidenoise(in_mif: str, out_mif: str) -> str: + run_cmd(["dwidenoise", in_mif, out_mif]) + return out_mif + + +@python.define +def run_mrdegibbs(in_mif: str, out_mif: str) -> str: + run_cmd(["mrdegibbs", in_mif, out_mif]) + return out_mif + + +@python.define +def concatenate_ap_pa(ap_mif: str, pa_mif: str, out_mif: str) -> str: + """Concatenate AP and PA DWI series (AP first) for rpe_all.""" + run_cmd(["dwicat", ap_mif, pa_mif, out_mif]) + return out_mif + + +@python.define +def build_se_epi( + dwi_mif: str, + rpe_nii: str, + rpe_json: str, + rpe_bvec: str, + rpe_bval: str, + fwd_pe_nii: str, + fwd_pe_json: str, + fwd_pe_bvec: str, + fwd_pe_bval: str, + pe_dir: str, + rpe_dir: str, + preproc_mode: str, + tmp_dir: str, +) -> str: + """Build se_epi pair for rpe_pair / rpe_split. Returns bzero_pair path.""" + if preproc_mode not in ("rpe_pair", "rpe_split"): + return "" + + rpe_mif = str(Path(tmp_dir) / f"DWI_ref_{rpe_dir}.mif.gz") + run_cmd( + [ + "mrconvert", + rpe_nii, + rpe_mif, + "-json_import", + rpe_json, + "-fslgrad", + rpe_bvec, + rpe_bval, + ] + ) + + if fwd_pe_nii: + fwd_mif = str(Path(tmp_dir) / f"DWI_ref_{pe_dir}.mif.gz") + run_cmd( + [ + "mrconvert", + fwd_pe_nii, + fwd_mif, + "-json_import", + fwd_pe_json, + "-fslgrad", + fwd_pe_bvec, + fwd_pe_bval, + ] + ) + else: + # #17 — use mean b0 from main DWI as forward PE + print(" No forward PE image — computing mean b0 from DWI...") + mean_b0 = str(Path(tmp_dir) / f"mean_bzero_{pe_dir}.mif.gz") + extract = subprocess.Popen( + ["dwiextract", dwi_mif, "-", "-bzero"], stdout=subprocess.PIPE + ) + subprocess.run( + ["mrmath", "-", "mean", mean_b0, "-axis", "3"], + stdin=extract.stdout, + check=True, + ) + extract.stdout.close() + extract.wait() + fwd_mif = mean_b0 + + bzero_pair = str(Path(tmp_dir) / "bzero_pair.mif.gz") + run_cmd(["dwicat", fwd_mif, rpe_mif, bzero_pair]) + return bzero_pair + + +@python.define +def run_dwifslpreproc( + dwi_mif: str, + out_mif: str, + pe_dir: str, + preproc_mode: str, + se_epi: str, + readout_time: float, + eddy_options: str, + scratch_dir: str = "", +) -> str: + scratch_dir = scratch_dir or str(Path(out_mif).parent / "dwifslpreproc_scratch") + cmd = [ + "dwifslpreproc", + dwi_mif, + out_mif, + "-pe_dir", + pe_dir, + f"-{preproc_mode}", + "-eddy_options", + eddy_options, + "-scratch", + scratch_dir, + "-nocleanup", + ] + if preproc_mode in ("rpe_pair", "rpe_split"): + cmd += ["-se_epi", se_epi, "-readout_time", str(readout_time), "-align_seepi"] + elif preproc_mode == "rpe_all": + cmd += ["-readout_time", str(readout_time)] + run_cmd(cmd) + # Clean up scratch after successful completion (dwifslpreproc -nocleanup keeps it) + import shutil as _shutil + + if Path(scratch_dir).exists(): + _shutil.rmtree(scratch_dir) + return out_mif + + +@python.define +def run_dwi2mask(in_mif: str, out_mif: str) -> str: + run_cmd(["dwi2mask", in_mif, out_mif]) + return out_mif + + +@python.define +def run_dwibiascorrect(in_mif: str, out_mif: str, mask_mif: str, bias_mif: str) -> str: + run_cmd( + [ + "dwibiascorrect", + "ants", + in_mif, + out_mif, + "-mask", + mask_mif, + "-bias", + bias_mif, + ] + ) + return out_mif + + +@python.define +def extract_mean_b0(dwi_biascorr_mif: str, out_nii: str) -> str: + extract = subprocess.Popen( + ["dwiextract", dwi_biascorr_mif, "-", "-bzero"], stdout=subprocess.PIPE + ) + subprocess.run( + ["mrmath", "-", "mean", out_nii, "-axis", "3"], stdin=extract.stdout, check=True + ) + extract.stdout.close() + extract.wait() + return out_nii + + +@python.define +def register_b0_to_t1(b0_nii: str, t1_nii: str, out_b0_in_t1: str, out_mat: str) -> str: + run_cmd( + [ + "flirt", + "-in", + b0_nii, + "-ref", + t1_nii, + "-out", + out_b0_in_t1, + "-omat", + out_mat, + "-dof", + "6", + ] + ) + return out_mat + + +@python.define(outputs=["t1_in_dwi", "mrtrix_xfm"]) +def invert_and_apply_transform( + b02t1_mat: str, t1_nii: str, b0_nii: str, tmp_dir: str +) -> tuple[str, str]: + t12b0_mat = str(Path(tmp_dir) / "T12b0.mat") + t1_in_dwi = str(Path(tmp_dir) / "T1_in_DWI_space.nii.gz") + mrtrix_txt = str(Path(tmp_dir) / "struct2diff_mrtrix.txt") + run_cmd(["convert_xfm", "-omat", t12b0_mat, "-inverse", b02t1_mat]) + run_cmd( + [ + "flirt", + "-in", + t1_nii, + "-ref", + b0_nii, + "-out", + t1_in_dwi, + "-init", + t12b0_mat, + "-applyxfm", + ] + ) + run_cmd(["transformconvert", t12b0_mat, t1_nii, b0_nii, "flirt_import", mrtrix_txt]) + return (t1_in_dwi, mrtrix_txt) + + +@python.define +def copy_t1_as_dwi_space(t1_nii: str, tmp_dir: str) -> str: + """ + rpe_none fallback: copy the T1 directly to + T1_in_DWI_space.nii.gz without performing any registration. + + B0-to-T1 co-registration is unreliable when no distortion correction + is applied (rpe_none), so the native T1 is used as-is. Stage 2 of the + orchestrator (phantom QC in DWI space) will therefore operate on the + native T1 geometry rather than a registered approximation. + """ + + t1_out = str(Path(tmp_dir) / "T1.nii.gz") + shutil.copy2(t1_nii, t1_out) + print(f" Copied T1 (no registration): {Path(t1_out).name}") + return t1_out + + +@python.define(outputs=["adc", "fa"]) +def compute_tensor_metrics(dwi_biascorr_mif: str, tmp_dir: str) -> tuple[str, str]: + tensor_mif = str(Path(tmp_dir) / "tensor.mif.gz") + adc_mif = str(Path(tmp_dir) / "ADC.mif.gz") + fa_mif = str(Path(tmp_dir) / "FA.mif.gz") + adc_nii = str(Path(tmp_dir) / "ADC.nii.gz") + fa_nii = str(Path(tmp_dir) / "FA.nii.gz") + run_cmd(["dwi2tensor", dwi_biascorr_mif, tensor_mif]) + run_cmd(["tensor2metric", "-adc", adc_mif, "-fa", fa_mif, tensor_mif]) + run_cmd(["mrconvert", adc_mif, adc_nii]) + run_cmd(["mrconvert", fa_mif, fa_nii]) + return (adc_nii, fa_nii) + + +@python.define(outputs=["out"]) +def copy_final_outputs_task( + dwi_biascorr_mif: str, + t1_in_dwi: str, + adc_nii: str, + fa_nii: str, + out_dir: str, + dwi_preproc_name: str, + t1_output_name: str, +): + os.makedirs(out_dir, exist_ok=True) + dst_dwi = str(Path(out_dir) / dwi_preproc_name) + dst_t1 = str(Path(out_dir) / t1_output_name) + dst_adc = str(Path(out_dir) / "ADC.nii.gz") + dst_fa = str(Path(out_dir) / "FA.nii.gz") + for src, dst in [ + (dwi_biascorr_mif, dst_dwi), + (t1_in_dwi, dst_t1), + (adc_nii, dst_adc), + (fa_nii, dst_fa), + ]: + shutil.copy2(src, dst) + print(f" Copied: {Path(dst).name}") + return dst_dwi # sentinel: used by cleanup ordering + + +def cleanup_tmp(sentinel: str, tmp_dir: str, keep_tmp: bool) -> str: + """ + Remove the tmp directory after all outputs have been copied. + Skipped if keep_tmp is True (--nocleanup flag). + sentinel: any upstream output — used only to enforce task ordering. + """ + if keep_tmp: + print(f" --nocleanup set: retaining {tmp_dir}") + return tmp_dir + if Path(tmp_dir).exists(): + shutil.rmtree(tmp_dir) + print(f" Removed tmp directory: {tmp_dir}") + return tmp_dir + + +# ============================================================================= +# Per-DWI workflow builder +# ============================================================================= + + +def _exists(*paths) -> bool: + """Return True only if every path exists on disk.""" + return all(Path(p).exists() for p in paths) + + +def build_dwi_workflow(plan: dict): + """ + Build and execute a pydra 1.0 workflow for a single DWI series. + + Design decisions: + - cache_root is set per-series to /.pydra_cache — isolated, + avoids cross-series contamination. + - rerun=True ensures pydra re-executes rather than replaying stale cache. + - Conditional branching (preproc_mode, do_denoise, do_gradcheck) is + resolved at workflow-class-definition time using concrete plan values + captured in the closure. Each unique combination gets a distinctly-named + class so pydra's cache correctly distinguishes them. + - convert_dicoms returns a sentinel output (=nii) so downstream tasks that + need T1-before-DWI ordering have a true lazy dependency edge, not a + resolved string literal that pydra treats as a known constant. + - CP1 (all final outputs exist) is checked before building the workflow — + fast early-return with no pydra overhead. + - cleanup_tmp runs as plain Python after sub() returns, guaranteeing it + only executes after all outputs are confirmed copied. + """ + tmp_dir = plan["tmp_dir"] + out_dir = plan["out_dir"] + pe_dir = plan["pe_dir"] + rpe_dir = plan["rpe_dir"] + preproc_mode = plan["preproc_mode"] + do_denoise = plan["do_denoise"] + do_gradcheck = plan["do_gradcheck"] + readout_time = plan["readout_time"] + eddy_options = plan["eddy_options"] + keep_tmp = plan["keep_tmp"] + dwi_preproc_name = plan["dwi_preproc_name"] + t1_output_name = plan["t1_output_name"] + t1_dir = plan["t1_dir"] + dwi_nii = plan["dwi_nii"] + dwi_json = plan["dwi_json"] + dwi_bvec = plan["dwi_bvec"] + dwi_bval = plan["dwi_bval"] + rpe_nii = plan["rpe_nii"] + rpe_json = plan["rpe_json"] + rpe_bvec = plan["rpe_bvec"] + rpe_bval = plan["rpe_bval"] + fwd_pe_nii = plan["fwd_pe_nii"] + fwd_pe_json = plan["fwd_pe_json"] + fwd_pe_bvec = plan["fwd_pe_bvec"] + fwd_pe_bval = plan["fwd_pe_bval"] + + t1_conv_dir = str(Path(tmp_dir) / "t1_nii") + dwi_raw_mif = str(Path(tmp_dir) / f"DWI_raw_{pe_dir}.mif.gz") + preproc_mif = str(Path(tmp_dir) / "DWI_preproc.mif.gz") + biascorr_mif = str(Path(tmp_dir) / dwi_preproc_name) + t1_in_dwi = str(Path(tmp_dir) / t1_output_name) + + out_dwi = str(Path(out_dir) / dwi_preproc_name) + out_adc = str(Path(out_dir) / "ADC.nii.gz") + out_fa = str(Path(out_dir) / "FA.nii.gz") + out_t1 = str(Path(out_dir) / t1_output_name) + + # ── CP1: all final outputs exist — nothing to do ───────────────────────── + if _exists(out_dwi, out_adc, out_fa, out_t1): + print(f" [skip] All outputs already exist in {Path(out_dir).name}") + return + + # ── CP2: bias-corrected MIF exists in out_dir — skip DWI preprocessing ─── + # dwifslpreproc + dwibiascorrect are by far the most expensive steps. + # If the bias-corrected MIF already exists in out_dir (e.g. from a prior + # successful run where only downstream steps failed), copy it into tmp/ so + # the pydra workflow can use it directly and skip straight to coregistration + # and tensor metrics. + skip_preproc = _exists(out_dwi) + if skip_preproc: + print(f" [skip] DWI preprocessing already done: {Path(out_dwi).name}") + os.makedirs(tmp_dir, exist_ok=True) + if not _exists(biascorr_mif): + shutil.copy2(out_dwi, biascorr_mif) + + # ── Workflow class name encodes the branch combination ──────────────────── + # Each unique combination gets a distinct class name so pydra caches them + # independently. skip_preproc adds an extra tag to distinguish the two paths. + dn_tag = "dn" if do_denoise else "nodn" + gc_tag = "gc" if do_gradcheck else "nogc" + sp_tag = "skippreproc" if skip_preproc else "fullpreproc" + wf_class_name = f"DwiWf_{preproc_mode}_{dn_tag}_{gc_tag}_{sp_tag}" + + if skip_preproc: + # ── Reduced workflow: T1 conversion + registration + tensor + copy ─── + # biascorr_mif is already in tmp/ (copied from out_dir above). + @workflow.define(outputs=["result"]) + def DwiWf(x: int) -> str: + + # T1 DICOM conversion (always needed for registration) + conv_t1 = workflow.add( + convert_dicoms(dicom_dir=t1_dir, out_dir=t1_conv_dir), + name="convert_t1", + ) + + # ── T1-to-DWI registration ──────────────────────────────────────── + b0 = workflow.add( + extract_mean_b0( + dwi_biascorr_mif=biascorr_mif, + out_nii=str(Path(tmp_dir) / "bzero_f.nii.gz"), + ), + name="mean_b0", + ) + flirt = workflow.add( + register_b0_to_t1( + b0_nii=b0.out, + t1_nii=conv_t1.nii, + out_b0_in_t1=str(Path(tmp_dir) / "b0_to_T1.nii.gz"), + out_mat=str(Path(tmp_dir) / "b02T1.mat"), + ), + name="register", + ) + xfm = workflow.add( + invert_and_apply_transform( + b02t1_mat=flirt.out, + t1_nii=conv_t1.nii, + b0_nii=b0.out, + tmp_dir=tmp_dir, + ), + name="invert_xfm", + ) + t1_result = xfm.t1_in_dwi + + # ── Tensor metrics ──────────────────────────────────────────────── + tensor = workflow.add( + compute_tensor_metrics( + dwi_biascorr_mif=biascorr_mif, + tmp_dir=tmp_dir, + ), + name="tensor", + ) + + # ── Copy final outputs ──────────────────────────────────────────── + copy = workflow.add( + copy_final_outputs_task( + dwi_biascorr_mif=biascorr_mif, + t1_in_dwi=t1_result, + adc_nii=tensor.adc, + fa_nii=tensor.fa, + out_dir=out_dir, + dwi_preproc_name=dwi_preproc_name, + t1_output_name=t1_output_name, + ), + name="copy_outputs", + ) + return copy.out + + else: + # ── Full workflow: all stages ───────────────────────────────────────── + @workflow.define(outputs=["result"]) + def DwiWf(x: int) -> str: + + # ── T1 DICOM conversion ─────────────────────────────────────────── + # Always run — needed for registration regardless of skip_preproc. + conv_t1 = workflow.add( + convert_dicoms(dicom_dir=t1_dir, out_dir=t1_conv_dir), + name="convert_t1", + ) + # conv_t1.sentinel is a true lazy output — used to gate DWI chain + + # ── DWI → MIF (with optional gradcheck) ────────────────────────── + if do_gradcheck: + mif_init = workflow.add( + convert_to_mif_initial( + nii=dwi_nii, + json_file=dwi_json, + bvec=dwi_bvec, + bval=dwi_bval, + out_path=str(Path(tmp_dir) / f"DWI_raw_{pe_dir}_init.mif.gz"), + ), + name="dwi_to_mif_init", + ) + gc_dwi = workflow.add( + run_gradcheck(mif_path=mif_init.out, out_dir=tmp_dir, prefix="dwi"), + name="gradcheck_dwi", + ) + dwi_bvec_lazy = gc_dwi.bvec + dwi_bval_lazy = gc_dwi.bval + + if preproc_mode == "rpe_all" and rpe_nii: + rpe_init = workflow.add( + convert_to_mif_initial( + nii=rpe_nii, + json_file=rpe_json, + bvec=rpe_bvec, + bval=rpe_bval, + out_path=str( + Path(tmp_dir) / f"DWI_raw_{rpe_dir}_init.mif.gz" + ), + ), + name="rpe_to_mif_init", + ) + gc_rpe = workflow.add( + run_gradcheck( + mif_path=rpe_init.out, out_dir=tmp_dir, prefix="rpe" + ), + name="gradcheck_rpe", + ) + rpe_bvec_lazy = gc_rpe.bvec + rpe_bval_lazy = gc_rpe.bval + else: + rpe_bvec_lazy = rpe_bvec + rpe_bval_lazy = rpe_bval + + if fwd_pe_nii and preproc_mode in ("rpe_pair", "rpe_split"): + fwd_init = workflow.add( + convert_to_mif_initial( + nii=fwd_pe_nii, + json_file=fwd_pe_json, + bvec=fwd_pe_bvec, + bval=fwd_pe_bval, + out_path=str( + Path(tmp_dir) / f"DWI_ref_{pe_dir}_init.mif.gz" + ), + ), + name="fwd_to_mif_init", + ) + gc_fwd = workflow.add( + run_gradcheck( + mif_path=fwd_init.out, out_dir=tmp_dir, prefix="fwd" + ), + name="gradcheck_fwd", + ) + fwd_bvec_lazy = gc_fwd.bvec + fwd_bval_lazy = gc_fwd.bval + else: + fwd_bvec_lazy = fwd_pe_bvec + fwd_bval_lazy = fwd_pe_bval + else: + dwi_bvec_lazy = dwi_bvec + dwi_bval_lazy = dwi_bval + rpe_bvec_lazy = rpe_bvec + rpe_bval_lazy = rpe_bval + fwd_bvec_lazy = fwd_pe_bvec + fwd_bval_lazy = fwd_pe_bval + + dwi_mif = workflow.add( + convert_to_mif_final( + nii=dwi_nii, + json_file=dwi_json, + bvec=dwi_bvec_lazy, + bval=dwi_bval_lazy, + out_path=dwi_raw_mif, + sentinel=conv_t1.sentinel, + ), + name="dwi_to_mif", + ) + + if preproc_mode == "rpe_all" and rpe_nii: + rpe_mif = workflow.add( + convert_to_mif_final( + nii=rpe_nii, + json_file=rpe_json, + bvec=rpe_bvec_lazy, + bval=rpe_bval_lazy, + out_path=str(Path(tmp_dir) / f"DWI_raw_{rpe_dir}.mif.gz"), + sentinel=conv_t1.sentinel, + ), + name="rpe_to_mif", + ) + rpe_mif_out = rpe_mif.out + else: + rpe_mif_out = None + + # ── Denoise + Gibbs ─────────────────────────────────────────────── + if do_denoise: + dn_dwi = workflow.add( + run_dwidenoise( + in_mif=dwi_mif.out, + out_mif=str(Path(tmp_dir) / f"DWI_raw_{pe_dir}_denoise.mif.gz"), + ), + name="denoise_dwi", + ) + dg_dwi = workflow.add( + run_mrdegibbs( + in_mif=dn_dwi.out, + out_mif=str( + Path(tmp_dir) / f"DWI_raw_{pe_dir}_denoise_gibbs.mif.gz" + ), + ), + name="degibbs_dwi", + ) + dwi_for_preproc = dg_dwi.out + + if preproc_mode == "rpe_all" and rpe_nii: + dn_rpe = workflow.add( + run_dwidenoise( + in_mif=rpe_mif_out, + out_mif=str( + Path(tmp_dir) / f"DWI_raw_{rpe_dir}_denoise.mif.gz" + ), + ), + name="denoise_rpe", + ) + dg_rpe = workflow.add( + run_mrdegibbs( + in_mif=dn_rpe.out, + out_mif=str( + Path(tmp_dir) + / f"DWI_raw_{rpe_dir}_denoise_gibbs.mif.gz" + ), + ), + name="degibbs_rpe", + ) + rpe_for_preproc = dg_rpe.out + else: + rpe_for_preproc = rpe_mif_out + else: + dwi_for_preproc = dwi_mif.out + rpe_for_preproc = rpe_mif_out + + # ── Concatenate AP+PA for rpe_all ───────────────────────────────── + if preproc_mode == "rpe_all" and rpe_nii: + cat = workflow.add( + concatenate_ap_pa( + ap_mif=dwi_for_preproc, + pa_mif=rpe_for_preproc, + out_mif=str(Path(tmp_dir) / "DWI_AP_PA_concat.mif.gz"), + ), + name="concat_ap_pa", + ) + dwi_input = cat.out + else: + dwi_input = dwi_for_preproc + + # ── se_epi pair (rpe_pair / rpe_split) ─────────────────────────── + se = workflow.add( + build_se_epi( + dwi_mif=dwi_for_preproc, + rpe_nii=rpe_nii or "", + rpe_json=rpe_json, + rpe_bvec=rpe_bvec_lazy, + rpe_bval=rpe_bval_lazy, + fwd_pe_nii=fwd_pe_nii or "", + fwd_pe_json=fwd_pe_json, + fwd_pe_bvec=fwd_bvec_lazy, + fwd_pe_bval=fwd_bval_lazy, + pe_dir=pe_dir, + rpe_dir=rpe_dir, + preproc_mode=preproc_mode, + tmp_dir=tmp_dir, + ), + name="se_epi", + ) + + # ── dwifslpreproc ───────────────────────────────────────────────── + fslpreproc = workflow.add( + run_dwifslpreproc( + dwi_mif=dwi_input, + out_mif=preproc_mif, + pe_dir=pe_dir, + preproc_mode=preproc_mode, + se_epi=se.out, + readout_time=readout_time, + eddy_options=eddy_options, + scratch_dir=str(Path(tmp_dir) / "dwifslpreproc_scratch"), + ), + name="fslpreproc", + ) + + # ── Mask + DWI bias correction ──────────────────────────────────── + dwi_mask = workflow.add( + run_dwi2mask( + in_mif=fslpreproc.out, + out_mif=str(Path(tmp_dir) / "DWI_preproc_mask.mif.gz"), + ), + name="dwi_mask", + ) + biascorr = workflow.add( + run_dwibiascorrect( + in_mif=fslpreproc.out, + out_mif=biascorr_mif, + mask_mif=dwi_mask.out, + bias_mif=str(Path(tmp_dir) / "DWI_preproc_bias.mif.gz"), + ), + name="dwi_biascorr", + ) + biascorr_out = biascorr.out + + # ── T1-to-DWI registration ──────────────────────────────────────── + b0 = workflow.add( + extract_mean_b0( + dwi_biascorr_mif=biascorr_out, + out_nii=str(Path(tmp_dir) / "bzero_f.nii.gz"), + ), + name="mean_b0", + ) + flirt = workflow.add( + register_b0_to_t1( + b0_nii=b0.out, + t1_nii=conv_t1.nii, + out_b0_in_t1=str(Path(tmp_dir) / "b0_to_T1.nii.gz"), + out_mat=str(Path(tmp_dir) / "b02T1.mat"), + ), + name="register", + ) + xfm = workflow.add( + invert_and_apply_transform( + b02t1_mat=flirt.out, + t1_nii=conv_t1.nii, + b0_nii=b0.out, + tmp_dir=tmp_dir, + ), + name="invert_xfm", + ) + t1_result = xfm.t1_in_dwi + + # ── Tensor metrics ──────────────────────────────────────────────── + tensor = workflow.add( + compute_tensor_metrics( + dwi_biascorr_mif=biascorr_out, + tmp_dir=tmp_dir, + ), + name="tensor", + ) + + # ── Copy final outputs ──────────────────────────────────────────── + copy = workflow.add( + copy_final_outputs_task( + dwi_biascorr_mif=biascorr_out, + t1_in_dwi=t1_result, + adc_nii=tensor.adc, + fa_nii=tensor.fa, + out_dir=out_dir, + dwi_preproc_name=dwi_preproc_name, + t1_output_name=t1_output_name, + ), + name="copy_outputs", + ) + return copy.out + + # Rename the class so pydra treats each branch combination as a distinct type + DwiWf.__name__ = wf_class_name + DwiWf.__qualname__ = wf_class_name + + # ── Submit ──────────────────────────────────────────────────────────────── + cache_dir = str(Path(out_dir) / ".pydra_cache") + wf = DwiWf(x=1) + with Submitter(worker="cf", cache_root=cache_dir) as sub: + sub(wf, rerun=True) + + # Cleanup runs after submission completes — guaranteed post-copy + cleanup_tmp(sentinel="", tmp_dir=tmp_dir, keep_tmp=keep_tmp) + + +# ============================================================================= +# Top-level pipeline +# ============================================================================= + + +def run_pipeline(cfg: dict): + scans_dir = cfg["scans_dir"] + output_dir = cfg.get("output_dir", str(Path(scans_dir).name)) + cfg["output_dir"] = output_dir + os.makedirs(output_dir, exist_ok=True) + + print(f"Scans directory: {scans_dir}") + print(f"Denoise/Degibbs: {cfg.get('denoise_degibbs', False)}") + print(f"Gradcheck: {cfg.get('gradcheck', False)}") + print(f"Keep tmp: {cfg.get('keep_tmp', False)}") + print(f"Output: {Path(output_dir).resolve()}") + + # Step 1: classify directories (provisional) + dirs = scan_directory(scans_dir) + print(f"\nFound:") + print(f" {len(dirs['t1_dirs'])} T1 series") + print(f" {len(dirs['candidate_dwi'])} candidate DWI series (before filtering)") + if dirs["ignored"]: + print(f" {len(dirs['ignored'])} ignored (non-DWI)") + + # Step 2: convert all candidate DICOMs + print("\nConverting DICOMs for all candidate series...") + conversions = convert_all_candidates(dirs["candidate_dwi"], output_dir) + + # Steps 3-6: validate bvec/bval, check b-values, finalise classification + classified = classify_candidates(dirs["candidate_dwi"], conversions) + + # Steps 7: detect AP/PA pairs (across dwi_dirs and pending series) + print("\nDetecting AP/PA pairs...") + # Pass existing fwd/rpe dirs into match so they are preserved + conversions["fwd_pe_dirs"] = classified["fwd_pe_dirs"] + conversions["rpe_dirs"] = classified["rpe_dirs"] + dwi_dirs, fwd_pe_dirs, rpe_dirs, rpe_all_map, tie_warnings = match_ap_pa_pairs( + classified["dwi_dirs"], + classified["pending_fwd"], + classified["pending_rpe"], + conversions, + ) + + # Step 8: build PE assignment map — one b0 PE image → all DWI series in same block + print("Assigning PE correction images to DWI blocks...") + pe_assignment_map = build_pe_assignment_map( + dwi_dirs, fwd_pe_dirs, rpe_dirs, rpe_all_map + ) + + print(f"\nAfter classification:") + print(f" {len(dwi_dirs)} DWI series to process") + print(f" {len(fwd_pe_dirs)} forward PE images") + print(f" {len(rpe_dirs)} reverse PE images") + print(f" {len(rpe_all_map)} rpe_all pairs detected") + if classified["skipped"]: + print(f" {len(classified['skipped'])} series skipped") + + if not dwi_dirs: + raise ValueError("No processable DWI series found after filtering.") + + # Steps 9-12: plan each workflow + print("\nPlanning workflows...") + plans = [] + for dwi_dir in dwi_dirs: + plan = plan_workflow( + dwi_dir=dwi_dir, + t1_dirs=dirs["t1_dirs"], + fwd_pe_dirs=fwd_pe_dirs, + rpe_dirs=rpe_dirs, + rpe_all_map=rpe_all_map, + pe_assignment_map=pe_assignment_map, + conversions=conversions, + cfg=cfg, + ) + plans.append(plan) + + # Step 13: print plan and pause + print_plan(plans, classified["skipped"]) + + for plan in plans: + build_dwi_workflow(plan) + + print("\n=== All done ===") + print(f"Outputs in: {Path(output_dir).resolve()}") + + +# ============================================================================= +# CLI +# ============================================================================= + + +def load_config(args) -> dict: + cfg = {} + if args.config: + with open(args.config) as f: + cfg = yaml.safe_load(f) or {} + if args.scans_dir: + cfg["scans_dir"] = args.scans_dir + if args.output_dir: + cfg["output_dir"] = args.output_dir + if args.denoise_degibbs: + cfg["denoise_degibbs"] = True + if args.gradcheck: + cfg["gradcheck"] = True + if args.nocleanup: + cfg["keep_tmp"] = True + if args.readout_time is not None: + cfg["readout_time"] = args.readout_time + if args.eddy_options is not None: + cfg["eddy_options"] = args.eddy_options + if "scans_dir" not in cfg: + raise ValueError("scans_dir must be provided via --scans-dir or config YAML") + return cfg + + +def main(): + parser = argparse.ArgumentParser( + description="DWI processing and tensor metrics pipeline" + ) + parser.add_argument("--config", type=str, help="Path to YAML config file") + parser.add_argument("--scans-dir", type=str, help="Path to scans directory") + parser.add_argument("--output-dir", type=str, help="Path to output directory") + parser.add_argument("--denoise-degibbs", action="store_true", default=None) + parser.add_argument("--gradcheck", action="store_true", default=None) + parser.add_argument( + "--nocleanup", + action="store_true", + default=False, + help="Keep the tmp/ directory after processing (default: remove)", + ) + parser.add_argument("--readout-time", type=float, default=None) + parser.add_argument("--eddy-options", type=str, default=None) + args = parser.parse_args() + cfg = load_config(args) + run_pipeline(cfg) + + +if __name__ == "__main__": + main() diff --git a/phantomkit/phantom_processor.py b/phantomkit/phantom_processor.py new file mode 100644 index 0000000..ffc1c4b --- /dev/null +++ b/phantomkit/phantom_processor.py @@ -0,0 +1,847 @@ +#!/usr/bin/env python3 +""" +phantom_processor.py +==================== +Core phantom processing engine for the phantomkit package. + +Contains ``PhantomProcessor``, a class that orchestrates: + + 1. ANTs registration + 2. Vial mask inverse-transform to subject space + 3. Per-vial metric extraction from all contrast images + 4. Plot generation (per-contrast scatter plots, T1/T2 parametric maps) + 5. Forward transform of all contrasts to template space + +This module is the functional equivalent of ``pydra_phantom_iterative.py`` +from the personal repo, re-implemented to integrate with the ``phantomkit`` +package structure. Plotting calls use the ``phantomkit.plotting`` API +directly rather than shelling out to external scripts. + +Path conventions (shared repo): + template_data//ImageTemplate.nii.gz + template_data//vials_labelled/*.nii.gz + template_data//adc_reference.json +""" + +import matplotlib + +matplotlib.use("Agg") # non-interactive backend; required when plotting runs +# in a background thread (e.g. ThreadPoolExecutor on macOS) + +import json +import re +import shutil +import subprocess +from pathlib import Path +from typing import Dict, List, Optional, Tuple + +import pandas as pd + + +# ============================================================================= +# PhantomProcessor +# ============================================================================= + + +class PhantomProcessor: + """ + Orchestrate phantom QC processing for a single session. + + Parameters + ---------- + template_dir: + Directory for this phantom type, e.g. + ``/template_data/SPIRIT``. + Must contain ``ImageTemplate.nii.gz`` and ``vials_labelled/``. + output_base_dir: + Top-level output directory. Results are written to a + ``/`` subdirectory within this path. + """ + + def __init__( + self, + template_dir: str, + output_base_dir: str, + ): + self.template_dir = Path(template_dir) + self.output_base_dir = Path(output_base_dir) + + # Phantom name is the last component of template_dir (e.g. "SPIRIT") + self.phantom_name = self.template_dir.name + + self.template_phantom = self.template_dir / "ImageTemplate.nii.gz" + self.vial_dir = self.template_dir / "vials_labelled" + self.vial_masks = sorted(self.vial_dir.glob("*.nii.gz")) + + # Load ADC vials from adc_reference.json + adc_ref_path = self.template_dir / "adc_reference.json" + if not adc_ref_path.exists(): + raise FileNotFoundError(f"adc_reference.json not found: {adc_ref_path}") + with open(adc_ref_path) as fh: + adc_ref = json.load(fh) + if "vials" not in adc_ref: + raise KeyError(f"'vials' key missing from: {adc_ref_path}") + self.adc_vials = {v.upper() for v in adc_ref["vials"]} + + if not self.template_phantom.exists(): + raise FileNotFoundError(f"Template not found: {self.template_phantom}") + if len(self.vial_masks) == 0: + raise FileNotFoundError(f"No vial masks found in: {self.vial_dir}") + + # ------------------------------------------------------------------------- + # Private helpers + # ------------------------------------------------------------------------- + + def _run_ants_registration( + self, input_image: str, output_prefix: str + ) -> Tuple[str, str, str]: + """Run ANTs rigid registration. Returns (warped, transform, inverse_warped).""" + cmd = [ + "antsRegistrationSyN.sh", + "-d", + "3", + "-f", + str(self.template_phantom), + "-m", + input_image, + "-o", + output_prefix, + "-t", + "r", + "-n", + "8", + "-j", + "1", + ] + print("Running ANTs registration...") + result = subprocess.run(cmd, capture_output=True, text=True) + if result.returncode != 0: + raise RuntimeError(f"ANTs failed: {result.stderr}") + + warped = f"{output_prefix}Warped.nii.gz" + transform = f"{output_prefix}0GenericAffine.mat" + inverse_warped = f"{output_prefix}InverseWarped.nii.gz" + return warped, transform, inverse_warped + + def _register( + self, input_image: str, session_name: str, tmp_dir: Path + ) -> Tuple[str, str, str]: + """ + Register the input image to the template phantom. + + Returns + ------- + (warped, transform, inverse_warped) + """ + output_prefix = str(tmp_dir / f"{session_name}_Transformed_") + warped, transform, inverse_warped = self._run_ants_registration( + input_image, output_prefix + ) + print(f" ✓ Registration complete") + return warped, transform, inverse_warped + + def _transform_vials_to_subject_space( + self, + reference_image: str, + transform_matrix: str, + output_vial_dir: Path, + ) -> List[str]: + """Transform vial masks from template space into subject (scanner) space.""" + output_vial_dir.mkdir(parents=True, exist_ok=True) + tmp_vial_dir = output_vial_dir.parent / "tmp_vials" + tmp_vial_dir.mkdir(parents=True, exist_ok=True) + + transformed_vials = [] + + for vial_mask in self.vial_masks: + vial_name = ( + Path(vial_mask) + .name.replace(".nii.gz", "") + .replace(".nii", "") + .split(".")[0] + ) + + tmp_vial = str(tmp_vial_dir / f"{vial_name}.nii") + output_vial = str(output_vial_dir / f"{vial_name}.nii.gz") + + # Inverse ANTs transform + cmd = [ + "antsApplyTransforms", + "-d", + "3", + "-i", + str(vial_mask), + "-r", + reference_image, + "-o", + tmp_vial, + "-t", + f"[{transform_matrix}, 1]", + "-n", + "NearestNeighbor", + ] + result = subprocess.run(cmd, capture_output=True, text=True) + if result.returncode != 0: + raise RuntimeError(f"Transform failed for {vial_name}: {result.stderr}") + + cmd = ["mrconvert", "-quiet", tmp_vial, output_vial, "-force"] + result = subprocess.run(cmd, capture_output=True, text=True) + if result.returncode != 0: + raise RuntimeError( + f"Final conversion failed for {vial_name}: {result.stderr}" + ) + + transformed_vials.append(output_vial) + + return transformed_vials + + def _transform_contrast_to_template_space( + self, + contrast_file: Path, + transform_matrix: str, + tmp_dir: Path, + ) -> str: + """ + Forward-transform a contrast image from subject space to template space. + + Handles 4D series and single-slice (degenerate z) images, which cause + ``antsApplyTransforms -d 3`` to fail silently without padding. + + Returns + ------- + str + Path to the contrast image in template space. + """ + tmp_dir.mkdir(parents=True, exist_ok=True) + contrast_name = contrast_file.stem.replace(".nii", "") + + source_image = str(contrast_file) + + # Detect 4D and single-slice + detect = subprocess.run( + ["mrinfo", "-size", source_image], capture_output=True, text=True + ) + size_parts = detect.stdout.strip().split() + is_4d = len(size_parts) >= 4 and int(size_parts[3]) > 1 + is_single_slice = len(size_parts) >= 3 and int(size_parts[2]) == 1 + + if is_single_slice: + padded = str(tmp_dir / f"{contrast_name}_padded.nii.gz") + cmd = [ + "mrgrid", + source_image, + "pad", + "-axis", + "2", + "1,1", + padded, + "-force", + ] + result = subprocess.run(cmd, capture_output=True, text=True) + if result.returncode != 0: + raise RuntimeError( + f"Padding of single-slice contrast {contrast_name} failed: " + f"{result.stderr}" + ) + transform_input = padded + else: + transform_input = source_image + + warped_tmp = str(tmp_dir / f"{contrast_name}_template_space_tmp.nii.gz") + warped_contrast = str(tmp_dir / f"{contrast_name}_template_space.nii.gz") + + cmd = [ + "antsApplyTransforms", + "-d", + "3", + "-e", + "3" if is_4d else "0", + "-i", + transform_input, + "-r", + str(self.template_phantom), + "-o", + warped_tmp, + "-t", + transform_matrix, + "-n", + "Linear", + ] + result = subprocess.run(cmd, capture_output=True, text=True) + if result.returncode != 0: + raise RuntimeError( + f"Forward transform of contrast {contrast_name} failed: {result.stderr}" + ) + + # Validate output was actually produced + verify = subprocess.run( + ["mrinfo", "-size", warped_tmp], capture_output=True, text=True + ) + if verify.returncode != 0 or not verify.stdout.strip(): + raise RuntimeError( + f"antsApplyTransforms produced no valid output for {contrast_name}. " + f"stderr: {result.stderr}" + ) + + shutil.move(warped_tmp, warped_contrast) + return warped_contrast + + @staticmethod + def _classify_contrast(contrast_file: Path) -> Optional[str]: + """ + Classify a contrast image by its filename stem. + + Returns + ------- + "adc" – filename contains 'ADC' (case-insensitive) + "fa" – filename contains standalone 'FA' (whole-word, case-insensitive) + None – no special classification + """ + stem = contrast_file.stem + if re.search(r"ADC", stem, re.IGNORECASE): + return "adc" + if re.search(r"(? Dict: + """Extract per-vial statistics from one contrast image across all vials.""" + # Strip both .nii.gz and .nii from contrast name + contrast_name = contrast_file.name + for ext in (".nii.gz", ".nii"): + if contrast_name.endswith(ext): + contrast_name = contrast_name[: -len(ext)] + break + + clean_contrast_name = contrast_name + + contrast_type = self._classify_contrast(contrast_file) + if contrast_type == "adc": + original_count = len(vial_masks) + vial_masks = [ + m + for m in vial_masks + if Path(m) + .name.replace(".nii.gz", "") + .replace(".nii", "") + .split(".")[0] + .upper() + in self.adc_vials + ] + print( + f" [ADC mode] Restricting to {len(self.adc_vials)} ADC vials " + f"({len(vial_masks)} of {original_count})" + ) + + # Get number of volumes + cmd = ["mrinfo", "-size", str(contrast_file)] + result = subprocess.run(cmd, capture_output=True, text=True) + size_info = result.stdout.strip().split() + nvols = ( + int(size_info[3]) if len(size_info) >= 4 and int(size_info[3]) > 0 else 1 + ) + + print( + f" Processing {clean_contrast_name} " + f"({nvols} volume{'s' if nvols > 1 else ''})" + ) + + metrics_data: Dict[str, Dict[str, List[float]]] = { + "mean": {}, + "median": {}, + "std": {}, + "min": {}, + "max": {}, + } + + tmp_vol_dir = output_metrics_dir.parent / "tmp_vols" + tmp_vol_dir.mkdir(parents=True, exist_ok=True) + + for vial_mask in vial_masks: + vial_name = ( + Path(vial_mask) + .name.replace(".nii.gz", "") + .replace(".nii", "") + .split(".")[0] + ) + for metric in metrics_data: + metrics_data[metric][vial_name] = [] + + regridded_mask = str(tmp_vol_dir / f"{contrast_name}_{vial_name}.nii") + cmd = [ + "mrgrid", + "-template", + str(contrast_file), + vial_mask, + "regrid", + regridded_mask, + "-force", + ] + result = subprocess.run(cmd, capture_output=True, text=True) + if result.returncode != 0: + raise RuntimeError( + f"mrgrid regrid failed for vial {vial_name}: {result.stderr}" + ) + + for vol_idx in range(nvols): + if nvols == 1: + vol_file = str(contrast_file) + else: + vol_file = str(tmp_vol_dir / f"{contrast_name}_vol{vol_idx}.nii.gz") + cmd = [ + "mrconvert", + str(contrast_file), + "-coord", + "3", + str(vol_idx), + vol_file, + "-quiet", + "-force", + ] + subprocess.run(cmd, check=True, capture_output=True) + + cmd = [ + "mrstats", + "-quiet", + vol_file, + "-output", + "mean", + "-output", + "median", + "-output", + "std", + "-output", + "min", + "-output", + "max", + "-mask", + regridded_mask, + ] + result = subprocess.run(cmd, capture_output=True, text=True) + if result.returncode != 0 or not result.stdout.strip(): + raise RuntimeError( + f"mrstats failed for vial {vial_name}: {result.stderr}" + ) + values = result.stdout.strip().split() + metrics_data["mean"][vial_name].append(float(values[0])) + metrics_data["median"][vial_name].append(float(values[1])) + metrics_data["std"][vial_name].append(float(values[2])) + metrics_data["min"][vial_name].append(float(values[3])) + metrics_data["max"][vial_name].append(float(values[4])) + + # Write CSV files + for metric_name, vial_data in metrics_data.items(): + csv_file = ( + output_metrics_dir + / f"{session_name}_{contrast_name}_{metric_name}_matrix.csv" + ) + rows = [ + { + "vial": vn, + **{f"{clean_contrast_name}_vol{i}": v for i, v in enumerate(vals)}, + } + for vn, vals in vial_data.items() + ] + pd.DataFrame(rows).to_csv(csv_file, index=False) + print(f" Saved: {csv_file.name}") + + return metrics_data + + def _generate_mrview_screenshot( + self, + contrast_file: Path, + roi_overlay: str, + output_image: str, + intensity_range: Optional[Tuple[float, float]] = None, + ) -> Optional[str]: + """ + Capture an mrview screenshot with a vial ROI overlay. + + Parameters + ---------- + intensity_range: + If provided, passes ``-intensity_range min,max`` to mrview. + Use (0, 1) for FA maps and (0, 0.005) for ADC maps. + """ + cmd = [ + "mrview", + str(contrast_file), + "-mode", + "1", + "-plane", + "2", + "-interpolation", + "0", + "-roi.load", + roi_overlay, + "-roi.colour", + "1,0,0", + "-roi.opacity", + "1", + "-comments", + "0", + "-noannotations", + "-fullscreen", + ] + + if intensity_range is not None: + cmd.extend( + ["-intensity_range", f"{intensity_range[0]},{intensity_range[1]}"] + ) + + cmd.extend( + [ + "-capture.folder", + str(Path(output_image).parent), + "-capture.prefix", + Path(output_image).stem, + "-capture.grab", + "-exit", + ] + ) + + result = subprocess.run(cmd, capture_output=True, text=True) + if result.returncode != 0: + print(f" ⚠ mrview screenshot failed: {result.stderr}") + return None + + # mrview appends "0000" to the prefix + actual_file = str( + Path(output_image).parent / f"{Path(output_image).stem}0000.png" + ) + return actual_file if Path(actual_file).exists() else None + + def _build_roi_overlay( + self, + contrast_file: Path, + vial_masks_list: List[Path], + prefix: str, + tmp_vial_dir: Path, + ) -> Optional[str]: + """Regrid each vial mask to contrast space and combine into one overlay NIfTI.""" + roi_overlay = str(tmp_vial_dir / f"{prefix}_VialsCombined.nii.gz") + regridded_vials = [] + + for vial_mask in vial_masks_list: + vial_name = vial_mask.name.replace(".nii.gz", "").replace(".nii", "") + regridded = str(tmp_vial_dir / f"{prefix}_{vial_name}.nii") + cmd = [ + "mrgrid", + "-template", + str(contrast_file), + str(vial_mask), + "regrid", + regridded, + "-force", + ] + subprocess.run(cmd, check=True, capture_output=True) + regridded_vials.append(regridded) + + if not regridded_vials: + return None + + cmd_cat = ["mrcat"] + regridded_vials + ["-", "-axis", "3"] + cmd_math = ["mrmath", "-", "max", roi_overlay, "-axis", "3", "-force"] + + proc_cat = subprocess.Popen( + cmd_cat, stdout=subprocess.PIPE, stderr=subprocess.PIPE + ) + proc_math = subprocess.Popen( + cmd_math, + stdin=proc_cat.stdout, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + ) + proc_cat.stdout.close() + proc_math.communicate() + + return roi_overlay if Path(roi_overlay).exists() else None + + def _generate_plots( + self, + contrast_files: List[Path], + metrics_dir: Path, + vial_dir: Path, + session_name: str, + ): + """ + Generate per-contrast scatter plots and parametric map plots (IR / TE). + + Plotting functions are imported directly from ``phantomkit.plotting`` + rather than shelled out to external scripts. + """ + from phantomkit.plotting.vial_intensity import plot_vial_intensity + from phantomkit.plotting.maps_ir import plot_vial_ir_means_std + from phantomkit.plotting.maps_te import plot_vial_te_means_std + + tmp_vial_dir = vial_dir / "tmp" + tmp_vial_dir.mkdir(exist_ok=True) + vial_masks_list = list(vial_dir.glob("*.nii.gz")) + + # ── Per-contrast scatter plots ──────────────────────────────────────── + for contrast_file in contrast_files: + contrast_name = contrast_file.name + for ext in (".nii.gz", ".nii"): + if contrast_name.endswith(ext): + contrast_name = contrast_name[: -len(ext)] + break + + mean_csv = metrics_dir / f"{session_name}_{contrast_name}_mean_matrix.csv" + std_csv = metrics_dir / f"{session_name}_{contrast_name}_std_matrix.csv" + + if not mean_csv.exists(): + print(f" ⚠ Mean CSV not found, skipping plot for {contrast_name}") + continue + + output_plot = str( + metrics_dir / f"{session_name}_{contrast_name}_PLOTmeanstd.png" + ) + + # Build ROI overlay and screenshot + roi_image = None + if vial_masks_list: + roi_overlay = self._build_roi_overlay( + contrast_file, vial_masks_list, contrast_name, tmp_vial_dir + ) + if roi_overlay: + contrast_type = self._classify_contrast(contrast_file) + intensity_range = ( + (0, 1) + if contrast_type == "fa" + else (0, 0.005) if contrast_type == "adc" else None + ) + screenshot_base = str( + tmp_vial_dir / f"{contrast_name}_roi_overlay.png" + ) + actual_screenshot = self._generate_mrview_screenshot( + contrast_file, + roi_overlay, + screenshot_base, + intensity_range=intensity_range, + ) + if actual_screenshot and Path(actual_screenshot).exists(): + roi_image = actual_screenshot + + try: + plot_vial_intensity( + csv_file=str(mean_csv), + plot_type="scatter", + std_csv=str(std_csv) if std_csv.exists() else None, + roi_image=roi_image, + output=output_plot, + phantom=self.phantom_name, + # template_dir is the parent of the phantom-specific dir + template_dir=str(self.template_dir.parent), + ) + print(f" ✓ Generated plot: {Path(output_plot).name}") + except Exception as e: + print(f" ✗ Plot generation failed for {contrast_name}: {e}") + + # ── Parametric map plots (IR and TE) ────────────────────────────────── + def _matches(stem: str, token: str) -> bool: + return bool(re.search(rf"(? Dict: + """ + Process a single phantom session end-to-end. + + Runs both the subject-space pipeline (vial mask transform, metric + extraction, plot generation) and the template-space pipeline (forward + transform of all contrasts). + + Parameters + ---------- + input_image: + Path to the primary input image (T1 MPRAGE or T1 in DWI space). + + Returns + ------- + dict + Processing results and output paths. + """ + input_path = Path(input_image) + session_name = input_path.parent.name + + output_dir = self.output_base_dir / session_name + tmp_dir = output_dir / "tmp" + vial_dir = output_dir / "vial_segmentations" + metrics_dir = output_dir / "metrics" + images_template_space_dir = output_dir / "images_template_space" + + for d in [tmp_dir, vial_dir, metrics_dir, images_template_space_dir]: + d.mkdir(parents=True, exist_ok=True) + + print(f"\n{'=' * 60}") + print(f"Processing Session: {session_name}") + print(f"Input: {input_image}") + print(f"Output: {output_dir}") + print(f"{'=' * 60}\n") + + # Step 1: Registration + print("Step 1: Registration") + warped, transform, inverse_warped = self._register( + str(input_image), session_name, tmp_dir + ) + + # Save template in scanner space + template_scanner_space = str(output_dir / "TemplatePhantom_ScannerSpace.nii.gz") + cmd = [ + "mrconvert", + "-quiet", + inverse_warped, + template_scanner_space, + "-force", + ] + subprocess.run(cmd, check=True, capture_output=True) + print(f" ✓ Saved template in scanner space") + + # Gather all contrast images — exclude derived template scanner space file + contrast_files = [ + f + for f in input_path.parent.glob("*.nii.gz") + if f.name != "TemplatePhantom_ScannerSpace.nii.gz" + ] + + # Step 2: Transform vials to subject space + print("\nStep 2: Transforming vials to subject space") + transformed_vials = self._transform_vials_to_subject_space( + reference_image=str(input_image), + transform_matrix=transform, + output_vial_dir=vial_dir, + ) + print(f" ✓ Transformed {len(transformed_vials)} vial masks") + + # Step 3: Extract metrics from all contrasts + print("\nStep 3: Extracting metrics from all contrasts") + print(f" Found {len(contrast_files)} contrast image(s)") + all_metrics = {} + for contrast_file in contrast_files: + metrics_data = self._extract_metrics_from_contrast( + contrast_file=contrast_file, + vial_masks=transformed_vials, + output_metrics_dir=metrics_dir, + session_name=session_name, + ) + all_metrics[contrast_file.name] = metrics_data + + # Step 4: Generate plots + print("\nStep 4: Generating plots") + self._generate_plots( + contrast_files=contrast_files, + metrics_dir=metrics_dir, + vial_dir=vial_dir, + session_name=session_name, + ) + + # Step 5: Forward-transform all contrasts to template space + print("\nStep 5: Transforming all contrasts to template space") + tmp_template_space_dir = tmp_dir / "template_space_contrasts" + tmp_template_space_dir.mkdir(parents=True, exist_ok=True) + + for contrast_file in contrast_files: + print(f" Transforming: {contrast_file.name}") + warped_path = self._transform_contrast_to_template_space( + contrast_file=contrast_file, + transform_matrix=transform, + tmp_dir=tmp_template_space_dir, + ) + shutil.copy2( + warped_path, + str(images_template_space_dir / contrast_file.name), + ) + print(f" ✓ {contrast_file.name} → template space") + print(f" ✓ All contrasts saved to: {images_template_space_dir}") + + # Step 6: Cleanup + print("\nStep 6: Cleaning up temporary directories") + for temp_dir in [ + tmp_dir, + output_dir / "tmp_vials", + output_dir / "tmp_vols", + vial_dir / "tmp", + ]: + if temp_dir.exists(): + try: + shutil.rmtree(temp_dir) + print(f" ✓ Removed: {temp_dir.name}") + except Exception as e: + print(f" ⚠ Warning: Could not remove {temp_dir.name}: {e}") + + print(f"\n{'=' * 60}") + print(f"✓ Session {session_name} complete!") + print(f" Metrics: {metrics_dir}") + print(f" Vial masks: {vial_dir}") + print(f" Template-space images: {images_template_space_dir}") + print(f"{'=' * 60}\n") + + return { + "session": session_name, + "output_dir": str(output_dir), + "metrics_dir": str(metrics_dir), + "vial_dir": str(vial_dir), + "images_template_space_dir": str(images_template_space_dir), + "space_image": str(output_dir / "TemplatePhantom_ScannerSpace.nii.gz"), + "metrics": all_metrics, + } diff --git a/phantomkit/pipeline.py b/phantomkit/pipeline.py new file mode 100644 index 0000000..432b21e --- /dev/null +++ b/phantomkit/pipeline.py @@ -0,0 +1,622 @@ +#!/usr/bin/env python3 +""" +pipeline.py +=========== +End-to-end orchestrator for the phantomkit phantom QC pipeline. + +Stages +------ + Stage 1 DWI processing (phantomkit.dwi_processing) + Runs if DWI acquisitions are found in --input-dir. + Outputs per DWI series: + T1_in_DWI_space.nii.gz, ADC.nii.gz, FA.nii.gz, + DWI_preproc_biascorr.mif.gz + + Stage 2 Phantom QC in DWI space (phantomkit.phantom_processor) + Runs once per DWI series produced by Stage 1. + Input image: T1_in_DWI_space.nii.gz from Stage 1. + ADC and FA files in the same folder are picked up automatically. + + Stage 3 Phantom QC on native contrasts (phantomkit.phantom_processor) + Runs if the input directory contains IR and/or TE series. + All contrast DICOMs (T1, IR, TE) are converted to NIfTI into a + temporary staging folder; the phantom processor then registers + the T1 and extracts vial metrics from every NIfTI present. + +Usage +----- + phantomkit pipeline \\ + --input-dir /path/to/patient/scans \\ + --output-dir /path/to/outputs \\ + --phantom SPIRIT + + Optional flags: + --denoise-degibbs pass to DWI pipeline + --gradcheck pass to DWI pipeline + --nocleanup keep DWI tmp/ directories + --readout-time override TotalReadoutTime (seconds) + --eddy-options override FSL eddy options string + --dry-run plan and print; do not execute + +Path conventions (shared repo) +------------------------------- + template_data//ImageTemplate.nii.gz + template_data//vials_labelled/*.nii.gz + template_data//adc_reference.json + template_data/rotations.txt + +The template_data/ directory is resolved relative to the installed +package's location (``phantomkit/`` directory → ``../template_data/``). +""" + +import argparse +import re +import shutil +import subprocess +import sys +from pathlib import Path + + +# --------------------------------------------------------------------------- +# Path resolution +# +# template_data/ sits at /template_data/ — one level above the +# phantomkit/ package directory. +# --------------------------------------------------------------------------- + + +def _find_template_data_root() -> Path: + """ + Locate the template_data/ directory relative to this file. + Walks up from phantomkit/ to the repo root. + """ + candidate = Path(__file__).resolve().parent # phantomkit/ + for _ in range(3): + td = candidate / "template_data" + if td.is_dir(): + return td + candidate = candidate.parent + raise RuntimeError( + f"Could not locate template_data/ within 3 levels of {Path(__file__).resolve()}" + ) + + +TEMPLATE_DATA_ROOT = _find_template_data_root() + + +def _find_dwi_processing_script() -> Path: + """Locate dwi_processing.py within the package.""" + return Path(__file__).resolve().parent / "dwi_processing.py" + + +DWI_SCRIPT = _find_dwi_processing_script() + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + + +def print_header(title: str): + bar = "=" * 70 + print(f"\n{bar}") + print(f" {title}") + print(f"{bar}\n") + + +def run_cmd(cmd: list, label: str): + """Run a subprocess command, streaming output and raising on failure.""" + print(f" >> {' '.join(str(c) for c in cmd)}\n") + result = subprocess.run([str(c) for c in cmd]) + if result.returncode != 0: + raise RuntimeError(f"{label} failed (exit {result.returncode}).") + + +def convert_dicom_dir(dicom_dir: Path, out_dir: Path) -> list: + """ + Run dcm2niix on dicom_dir, placing NIfTIs in out_dir. + Returns list of produced .nii.gz paths. + """ + out_dir.mkdir(parents=True, exist_ok=True) + subprocess.run( + ["dcm2niix", "-o", str(out_dir), "-f", "%p", "-z", "y", str(dicom_dir)], + check=True, + capture_output=True, + ) + niis = sorted(out_dir.glob("*.nii.gz")) + if not niis: + raise FileNotFoundError(f"dcm2niix produced no NIfTI from {dicom_dir}") + return [str(p) for p in niis] + + +def scan_input_dir(input_dir: Path) -> dict: + """ + Classify subdirectories in input_dir into: + t1_dirs – folders matching 't1' (case-insensitive) + ir_dirs – folders matching 'se_ir' or standalone 'ir' + te_dirs – folders matching 't2_se' or standalone 'te' + dwi_dirs – folders matching '_diff_' or '_DWI_' (case-insensitive) + AND not ending in _ADC or _FA + other_dirs – everything else + """ + t1_dirs, ir_dirs, te_dirs, dwi_dirs, other_dirs = [], [], [], [], [] + + for d in sorted(input_dir.iterdir()): + if not d.is_dir(): + continue + name = d.name + + if re.search(r"_(ADC|FA)$", name, re.IGNORECASE): + other_dirs.append(d) + continue + + if re.search(r"t1", name, re.IGNORECASE): + t1_dirs.append(d) + elif re.search(r"se_ir|(? str: + """ + Derive a clean session name from the input directory. + Strips any leading series-number prefix (e.g. '87-PatientID' → 'PatientID'). + """ + name = input_dir.name + stripped = re.sub(r"^\d+-", "", name) + return stripped if stripped else name + + +# --------------------------------------------------------------------------- +# Stage 1: DWI processing +# --------------------------------------------------------------------------- + + +def run_stage1(input_dir: Path, output_dir: Path, cfg: dict, dry_run: bool) -> list: + """ + Run dwi_processing.py on input_dir. + Returns list of DWI output subdirectory Paths (one per processed series). + """ + print_header("STAGE 1 — DWI Processing") + + cmd = [ + sys.executable, + str(DWI_SCRIPT), + "--scans-dir", + str(input_dir), + "--output-dir", + str(output_dir), + ] + + if cfg.get("denoise_degibbs"): + cmd.append("--denoise-degibbs") + if cfg.get("gradcheck"): + cmd.append("--gradcheck") + if cfg.get("nocleanup"): + cmd.append("--nocleanup") + if cfg.get("readout_time") is not None: + cmd += ["--readout-time", str(cfg["readout_time"])] + if cfg.get("eddy_options") is not None: + cmd += ["--eddy-options", cfg["eddy_options"]] + + if dry_run: + print(" [DRY RUN] Would execute:") + print(f" {' '.join(str(c) for c in cmd)}\n") + return [] + + run_cmd(cmd, "Stage 1 (DWI processing)") + + _t1_names = {"T1_in_DWI_space.nii.gz", "T1.nii.gz"} + dwi_output_dirs = [ + d + for d in sorted(output_dir.iterdir()) + if d.is_dir() and any((d / n).exists() for n in _t1_names) + ] + + print(f"\n Stage 1 complete. Found {len(dwi_output_dirs)} processed DWI series.") + for d in dwi_output_dirs: + print(f" {d.name}") + + return dwi_output_dirs + + +# --------------------------------------------------------------------------- +# Stage 2: Phantom QC in DWI space +# --------------------------------------------------------------------------- + + +def run_stage2( + dwi_output_dirs: list, + output_dir: Path, + template_dir: Path, + dry_run: bool, +): + """ + For each DWI output directory, run PhantomProcessor on + T1_in_DWI_space.nii.gz (or T1.nii.gz for rpe_none). + ADC.nii.gz and FA.nii.gz in the same folder are picked up automatically. + """ + print_header("STAGE 2 — Phantom QC in DWI Space") + + if not dwi_output_dirs: + print(" No DWI output directories to process — skipping Stage 2.\n") + return + + for dwi_dir in dwi_output_dirs: + t1_in_dwi = next( + ( + dwi_dir / n + for n in ("T1_in_DWI_space.nii.gz", "T1.nii.gz") + if (dwi_dir / n).exists() + ), + None, + ) + + if t1_in_dwi is None: + print(f" WARNING: No T1 image found in {dwi_dir.name} — skipping.") + continue + + print(f" Processing: {dwi_dir.name}") + print(f" Input image: {t1_in_dwi}") + print(f" Contrast images found alongside T1:") + for nii in sorted(dwi_dir.glob("*.nii.gz")): + print(f" {nii.name}") + + if dry_run: + print(f" [DRY RUN] Would run PhantomProcessor on {t1_in_dwi.name}\n") + continue + + from phantomkit.phantom_processor import PhantomProcessor + + processor = PhantomProcessor( + template_dir=str(template_dir), + output_base_dir=str(output_dir), + ) + processor.process_session(str(t1_in_dwi)) + print() + + print(" Stage 2 complete.\n") + + +# --------------------------------------------------------------------------- +# Stage 3: Phantom QC on native contrasts +# --------------------------------------------------------------------------- + + +def run_stage3( + input_dir: Path, + output_dir: Path, + template_dir: Path, + scan_info: dict, + dry_run: bool, +): + """ + Convert T1, IR, and TE DICOMs to NIfTI into a staging folder, then + run PhantomProcessor on the T1. The processor picks up all NIfTIs + in the staging folder automatically. + """ + print_header("STAGE 3 — Phantom QC on Native Contrasts") + + t1_dirs = scan_info["t1_dirs"] + ir_dirs = scan_info["ir_dirs"] + te_dirs = scan_info["te_dirs"] + + if not t1_dirs: + print(" No T1 directory found — cannot run Stage 3.\n") + return + + session_name = derive_session_name(input_dir) + staging_dir = output_dir / "native_contrasts_staging" + + print(f" Session name: {session_name}") + print(f" Staging folder: {staging_dir}") + print(f" T1 dirs: {[d.name for d in t1_dirs]}") + print(f" IR dirs: {[d.name for d in ir_dirs]}") + print(f" TE dirs: {[d.name for d in te_dirs]}") + print() + + # Pick the T1 with the lowest series number + def _series_num(d: Path) -> int: + m = re.match(r"^(\d+)-", d.name) + return int(m.group(1)) if m else 0 + + primary_t1_dir = sorted(t1_dirs, key=_series_num)[0] + contrast_dirs_to_convert = [primary_t1_dir] + ir_dirs + te_dirs + + if dry_run: + print(" [DRY RUN] Would convert DICOMs to NIfTI and place in staging folder:") + for d in contrast_dirs_to_convert: + print(f" {d.name} → {staging_dir}/") + print("\n [DRY RUN] Would run PhantomProcessor on converted T1 NIfTI.\n") + return + + # Convert DICOMs to NIfTI + staging_dir.mkdir(parents=True, exist_ok=True) + t1_nii_path = None + + for dicom_dir in contrast_dirs_to_convert: + print(f" Converting DICOM: {dicom_dir.name}") + try: + produced = convert_dicom_dir(dicom_dir, staging_dir) + print(f" Produced: {[Path(p).name for p in produced]}") + if dicom_dir == primary_t1_dir: + t1_nii_path = Path(produced[0]) + except Exception as e: + print(f" WARNING: DICOM conversion failed for {dicom_dir.name}: {e}") + if dicom_dir == primary_t1_dir: + print(" Cannot proceed with Stage 3 without a T1 image.") + shutil.rmtree(staging_dir, ignore_errors=True) + return + + if t1_nii_path is None or not t1_nii_path.exists(): + print(" ERROR: T1 NIfTI was not produced — aborting Stage 3.") + shutil.rmtree(staging_dir, ignore_errors=True) + return + + print(f"\n All NIfTIs in staging folder:") + for nii in sorted(staging_dir.glob("*.nii.gz")): + print(f" {nii.name}") + + # Run phantom processor + print(f"\n Running PhantomProcessor:") + print(f" Input image: {t1_nii_path}") + print(f" Output base: {output_dir}") + + from phantomkit.phantom_processor import PhantomProcessor + + processor = PhantomProcessor( + template_dir=str(template_dir), + output_base_dir=str(output_dir), + ) + processor.process_session(str(t1_nii_path)) + + # Clean up staging NIfTIs only — leave processed outputs in place + print(f"\n Removing temporary NIfTIs from staging folder: {staging_dir}") + for nii in staging_dir.glob("*.nii.gz"): + nii.unlink(missing_ok=True) + for jsn in staging_dir.glob("*.json"): + jsn.unlink(missing_ok=True) + try: + staging_dir.rmdir() + print(f" Removed empty staging directory: {staging_dir.name}") + except OSError: + print( + f" Staging directory retained (contains phantom outputs): " + f"{staging_dir.name}" + ) + + print("\n Stage 3 complete.\n") + + +# --------------------------------------------------------------------------- +# Validation +# --------------------------------------------------------------------------- + + +def validate_inputs(args): + """Validate paths and phantom name before any processing.""" + errors = [] + + input_dir = Path(args.input_dir) + if not input_dir.is_dir(): + errors.append(f"--input-dir does not exist or is not a directory: {input_dir}") + + phantom_dir = TEMPLATE_DATA_ROOT / args.phantom + if not phantom_dir.is_dir(): + errors.append( + f"Phantom template directory not found: {phantom_dir}\n" + f" Expected: template_data/{args.phantom}/" + ) + else: + template_img = phantom_dir / "ImageTemplate.nii.gz" + vials_dir = phantom_dir / "vials_labelled" + if not template_img.exists(): + errors.append(f"ImageTemplate.nii.gz not found in: {phantom_dir}") + if not vials_dir.is_dir() or not list(vials_dir.glob("*.nii.gz")): + errors.append( + f"vials_labelled/ with .nii.gz masks not found in: {phantom_dir}" + ) + + if not DWI_SCRIPT.exists(): + errors.append(f"DWI processing script not found: {DWI_SCRIPT}") + + if errors: + print("\nValidation errors:") + for e in errors: + print(f" ✗ {e}") + sys.exit(1) + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + + +def main(): + parser = argparse.ArgumentParser( + description="End-to-end phantomkit phantom QC + DWI processing pipeline.", + formatter_class=argparse.RawDescriptionHelpFormatter, + epilog=__doc__, + ) + + parser.add_argument( + "--input-dir", + required=True, + help="Root directory containing acquisition subdirectories (DICOM folders).", + ) + parser.add_argument( + "--output-dir", + required=True, + help="Top-level output directory. All results are written here.", + ) + parser.add_argument( + "--phantom", + required=True, + help="Phantom name, e.g. SPIRIT. Used to locate template_data//.", + ) + parser.add_argument( + "--denoise-degibbs", + action="store_true", + default=False, + help="Apply dwidenoise + mrdegibbs before preprocessing.", + ) + parser.add_argument( + "--gradcheck", + action="store_true", + default=False, + help="Run dwigradcheck to verify gradient orientations.", + ) + parser.add_argument( + "--nocleanup", + action="store_true", + default=False, + help="Keep DWI tmp/ intermediate directories.", + ) + parser.add_argument( + "--readout-time", + type=float, + default=None, + help="Override TotalReadoutTime (seconds) for dwifslpreproc.", + ) + parser.add_argument( + "--eddy-options", + type=str, + default=None, + help="Override FSL eddy options string.", + ) + parser.add_argument( + "--dry-run", + action="store_true", + default=False, + help="Plan and print commands; do not execute any processing.", + ) + + args = parser.parse_args() + + input_dir = Path(args.input_dir).resolve() + output_dir = Path(args.output_dir).resolve() + phantom = args.phantom + template_dir = TEMPLATE_DATA_ROOT / phantom + + validate_inputs(args) + output_dir.mkdir(parents=True, exist_ok=True) + + dwi_cfg = { + "denoise_degibbs": args.denoise_degibbs, + "gradcheck": args.gradcheck, + "nocleanup": args.nocleanup, + "readout_time": args.readout_time, + "eddy_options": args.eddy_options, + } + + # ── Discover what's in the input directory ─────────────────────────────── + print_header("Input Directory Scan") + scan_info = scan_input_dir(input_dir) + + print(f" Input directory: {input_dir}") + print(f" Phantom: {phantom}") + print(f" Template dir: {template_dir}") + print(f" Output dir: {output_dir}") + print() + print(f" T1 directories: {len(scan_info['t1_dirs'])}") + for d in scan_info["t1_dirs"]: + print(f" {d.name}") + print(f" IR directories: {len(scan_info['ir_dirs'])}") + for d in scan_info["ir_dirs"]: + print(f" {d.name}") + print(f" TE directories: {len(scan_info['te_dirs'])}") + for d in scan_info["te_dirs"]: + print(f" {d.name}") + print(f" DWI candidate dirs: {len(scan_info['dwi_dirs'])}") + for d in scan_info["dwi_dirs"]: + print(f" {d.name}") + print(f" Other (ignored): {len(scan_info['other_dirs'])}") + print() + + run_stage1_flag = scan_info["has_dwi"] + run_stage3_flag = bool(scan_info["t1_dirs"]) and ( + scan_info["has_native_contrasts"] or not run_stage1_flag + ) + + print( + f" Stage 1 (DWI): " + f"{'YES' if run_stage1_flag else 'NO (no DWI found)'}" + ) + print( + f" Stage 2 (phantom QC, DWI space): " + f"{'YES (runs per DWI series)' if run_stage1_flag else 'NO'}" + ) + print( + f" Stage 3 (phantom QC, native T1): " + f"{'YES' if run_stage3_flag else 'NO (no T1/IR/TE found)'}" + ) + print() + + if args.dry_run: + print(" NOTE: --dry-run is active. No processing will be performed.\n") + + # ── Execution ───────────────────────────────────────────────────────────── + # Stages run sequentially: 1 → 2 → 3. + # (Stage 1 shells out to dwifslpreproc/eddy which are already multi-threaded; + # Stage 3 runs ANTs which is also multi-threaded. Parallel execution of + # stages gives no wall-time benefit and produces interleaved output.) + + # Stage 1 — DWI processing + dwi_output_dirs = [] + if run_stage1_flag: + dwi_output_dirs = run_stage1(input_dir, output_dir, dwi_cfg, args.dry_run) + else: + print_header("STAGE 1 — DWI Processing") + print(" Skipped: no DWI acquisitions found.\n") + + # Stage 2 — Phantom QC in DWI space (follows Stage 1) + if run_stage1_flag: + run_stage2(dwi_output_dirs, output_dir, template_dir, args.dry_run) + else: + print_header("STAGE 2 — Phantom QC in DWI Space") + print(" Skipped: Stage 1 did not run.\n") + + # Stage 3 — Phantom QC on native contrasts + if run_stage3_flag: + run_stage3(input_dir, output_dir, template_dir, scan_info, args.dry_run) + else: + print_header("STAGE 3 — Phantom QC on Native Contrasts") + if not scan_info["t1_dirs"]: + print(" Skipped: no T1 directory found.\n") + else: + print(" Skipped: no IR or TE series found, and DWI pipeline was run.\n") + + # ── Summary ─────────────────────────────────────────────────────────────── + print_header("Pipeline Complete") + print(f" All outputs written to: {output_dir}\n") + + if not args.dry_run: + print(" Output structure:") + for item in sorted(output_dir.iterdir()): + if item.is_dir(): + print(f" {item.name}/") + for sub in sorted(item.iterdir()): + if sub.is_dir(): + print(f" {sub.name}/") + else: + print(f" {sub.name}") + print() + + +if __name__ == "__main__": + main() diff --git a/phantomkit/plotting/maps_ir.py b/phantomkit/plotting/maps_ir.py index eca601b..e692728 100644 --- a/phantomkit/plotting/maps_ir.py +++ b/phantomkit/plotting/maps_ir.py @@ -13,6 +13,8 @@ ================================================================================ """ +import matplotlib + import logging import os import re @@ -23,6 +25,9 @@ import pandas as pd from scipy.optimize import curve_fit +matplotlib.use("Agg") # non-interactive backend; required when plotting runs +# in a background thread (e.g. ThreadPoolExecutor on macOS) + logger = logging.getLogger(__name__) @@ -30,6 +35,11 @@ def find_csv_file(metric_dir, contrast_name, suffix): """ Find a CSV file in a directory that matches the contrast name and suffix. + Uses an exact token match: the contrast_name must appear in the filename + immediately followed by the suffix (no extra characters between them). + This prevents ambiguous substring hits such as 'se_ir_100' matching + 'se_ir_1000_mean_matrix.csv' or 'se_ir_50' matching 'se_ir_500_mean_matrix.csv'. + Args: metric_dir: Directory containing CSV files contrast_name: Base name of the contrast to search for @@ -41,8 +51,9 @@ def find_csv_file(metric_dir, contrast_name, suffix): Raises: FileNotFoundError: If no matching file is found """ + exact_tail = contrast_name + suffix for f in os.listdir(metric_dir): - if contrast_name in f and f.endswith(suffix): + if f.endswith(exact_tail): return os.path.join(metric_dir, f) raise FileNotFoundError( f"No CSV file found for contrast '{contrast_name}' with suffix '{suffix}' in {metric_dir}" @@ -509,29 +520,14 @@ def plot_vial_ir_means_std( return os.path.abspath(output_file) -@click.command() +@click.command("maps-ir") @click.argument("contrast_files", nargs=-1, required=True) -@click.option( - "-m", - "--metric_dir", - required=True, - help="Directory containing the mean/std CSV files.", -) -@click.option( - "-o", - "--output", - default="vial_summary_T1.png", - show_default=True, - help="Output filename for the plot.", -) -@click.option("--annotate", is_flag=True, help="Annotate each point with mean ± std.") -@click.option( - "--roi_image", - default=None, - help="Path to ROI overlay PNG image for the extra subplot.", -) +@click.option("-m", "--metric-dir", required=True, help="Directory with mean/std CSVs.") +@click.option("-o", "--output", default="vial_summary_T1.png", show_default=True) +@click.option("--annotate", is_flag=True, default=False) +@click.option("--roi-image", default=None) def main(contrast_files, metric_dir, output, annotate, roi_image): - """Plot vial mean ± std for inversion recovery with T₁ fitting and save fit metrics.""" + """Plot vial mean ± std for inversion recovery with T1 fitting.""" logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s") plot_vial_ir_means_std( list(contrast_files), diff --git a/phantomkit/plotting/maps_te.py b/phantomkit/plotting/maps_te.py index c74a58f..2e9632f 100644 --- a/phantomkit/plotting/maps_te.py +++ b/phantomkit/plotting/maps_te.py @@ -13,23 +13,30 @@ ================================================================================ """ -import logging -import os -import re +import matplotlib -import click +matplotlib.use("Agg") # non-interactive backend; required when plotting runs +# in a background thread (e.g. ThreadPoolExecutor on macOS) + +import pandas as pd import matplotlib.pyplot as plt +import os import numpy as np -import pandas as pd +import argparse +import click +import re from scipy.optimize import curve_fit -logger = logging.getLogger(__name__) - def find_csv_file(metric_dir, contrast_name, suffix): """ Find a CSV file in a directory that matches the contrast name and suffix. + Uses an exact token match: the contrast_name must appear in the filename + immediately followed by the suffix (no extra characters between them). + This prevents ambiguous substring hits such as 't2_se_TE_14' matching + 't2_se_TE_140_mean_matrix.csv'. + Args: metric_dir: Directory containing CSV files contrast_name: Base name of the contrast to search for @@ -41,8 +48,9 @@ def find_csv_file(metric_dir, contrast_name, suffix): Raises: FileNotFoundError: If no matching file is found """ + exact_tail = contrast_name + suffix for f in os.listdir(metric_dir): - if contrast_name in f and f.endswith(suffix): + if f.endswith(exact_tail): return os.path.join(metric_dir, f) raise FileNotFoundError( f"No CSV file found for contrast '{contrast_name}' with suffix '{suffix}' in {metric_dir}" @@ -103,12 +111,12 @@ def calc_r2(y_true, y_pred): def plot_vial_te_means_std( - contrast_files: list[str], - metric_dir: str, - output_file: str = "vial_summary_pub.png", - annotate: bool = False, - roi_image: str | None = None, -) -> str: + contrast_files, + metric_dir, + output_file="vial_summary_pub.png", + annotate=False, + roi_image=None, +): """ Create publication-quality plots of vial intensity data with T2 curve fitting. @@ -121,9 +129,6 @@ def plot_vial_te_means_std( output_file: Output filename for the plot (default: 'vial_summary_pub.png') annotate: Whether to annotate points with mean ± std (not currently used) roi_image: Optional path to ROI overlay image for extra subplot - - Returns: - Absolute path to the saved plot file. """ # ======================================================================== @@ -296,9 +301,7 @@ def plot_vial_te_means_std( except (np.linalg.LinAlgError, ValueError) as e: # Covariance matrix might be singular or ill-conditioned - logger.warning( - "Could not calculate 95%% CI for vial %s: %s", vial, e - ) + print(f"[WARN] Could not calculate 95% CI for vial {vial}: {e}") # ============================================================ # *** FITTED CURVE AND CI BAND PLOTTING *** @@ -327,7 +330,7 @@ def plot_vial_te_means_std( ) except RuntimeError: # Curve fitting failed for this vial - logger.warning("Could not fit T\u2082 for vial %s", vial) + print(f"[WARN] Could not fit T₂ for vial {vial}") # -------------------------------------------------------------------- # CASE 2: Multiple vials in one subplot @@ -495,41 +498,25 @@ def plot_vial_te_means_std( # ======================================================================== plt.tight_layout(rect=[0, 0, 1, 1]) plt.savefig(output_file, dpi=300, bbox_inches="tight") - logger.info("Publication-ready plot saved to %s", output_file) + print(f"[INFO] Publication-ready plot saved to {output_file}") # Save fitted T2 values to CSV csv_output = os.path.splitext(output_file)[0] + "_T2_fits.csv" pd.DataFrame(fit_results).to_csv(csv_output, index=False) - logger.info("Fitted parameters saved to %s", csv_output) + print(f"[INFO] Fitted parameters saved to {csv_output}") plt.close(fig) - return os.path.abspath(output_file) + return output_file -@click.command() +@click.command("maps-te") @click.argument("contrast_files", nargs=-1, required=True) -@click.option( - "-m", - "--metric_dir", - required=True, - help="Directory containing the mean/std CSV files.", -) -@click.option( - "-o", - "--output", - default="vial_summary_pub.png", - show_default=True, - help="Output filename for the plot.", -) -@click.option("--annotate", is_flag=True, help="Annotate each point with mean ± std.") -@click.option( - "--roi_image", - default=None, - help="Path to ROI overlay PNG image for the extra subplot.", -) +@click.option("-m", "--metric-dir", required=True, help="Directory with mean/std CSVs.") +@click.option("-o", "--output", default="vial_summary_pub.png", show_default=True) +@click.option("--annotate", is_flag=True, default=False) +@click.option("--roi-image", default=None) def main(contrast_files, metric_dir, output, annotate, roi_image): - """Plot grouped vial mean ± std with mono-exponential T₂ fitting and save fit metrics.""" - logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s") + """Plot vial mean ± std for T2 spin-echo with mono-exponential fitting.""" plot_vial_te_means_std( list(contrast_files), metric_dir=metric_dir, diff --git a/phantomkit/plotting/vial_intensity.py b/phantomkit/plotting/vial_intensity.py index ee269f6..f9ccfb3 100644 --- a/phantomkit/plotting/vial_intensity.py +++ b/phantomkit/plotting/vial_intensity.py @@ -1,16 +1,44 @@ #!/usr/bin/env python3 -import csv -import logging -import os -import sys +""" +plot_vial_intensity.py +Plot vial vs intensity (mean ± std) for 3D or 4D contrasts, with optional ROI overlay. + +Plot mode is detected automatically from the csv_file filename: + + ADC mode – filename contains 'ADC' (case-insensitive) + · Y-axis: ADC ×10⁻³ mm²/s + · Reference values + ±5%/±10% tolerance bands drawn (blue crosshairs / shading) + · Measured values plotted as filled red circles + · Requires --phantom and --template_dir so the reference JSON can be loaded + + FA mode – filename contains standalone 'FA' (case-insensitive, whole-word) + · Y-axis: Fractional Anisotropy, fixed range 0–1 + + Generic – all other filenames → standard intensity plot +""" import click -import matplotlib.image as mpimg +import matplotlib + +matplotlib.use("Agg") # non-interactive backend; required when plotting runs +# in a background thread (e.g. ThreadPoolExecutor on macOS) + +import pandas as pd import matplotlib.pyplot as plt +import matplotlib.image as mpimg +import argparse +import sys +import csv +import os +import json +import re import numpy as np -import pandas as pd +from pathlib import Path -logger = logging.getLogger(__name__) + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- def detect_separator(file_path): @@ -25,27 +53,201 @@ def detect_separator(file_path): sys.exit(f"Error detecting separator: {e}") +def load_adc_reference(template_dir: str, phantom: str) -> dict: + """ + Load ADC reference values from TemplateData//adc_reference.json. + + Returns a dict with keys: + 'vials' : list of vial labels in order + 'adc_mm2_per_s': dict mapping vial label -> reference ADC (mm²/s) + """ + ref_path = os.path.join(template_dir, phantom, "adc_reference.json") + if not os.path.isfile(ref_path): + sys.exit( + f"Error: ADC reference file not found at '{ref_path}'.\n" + f"Expected structure: //adc_reference.json" + ) + with open(ref_path) as fh: + data = json.load(fh) + required_keys = {"vials", "adc_mm2_per_s"} + if not required_keys.issubset(data.keys()): + sys.exit(f"Error: adc_reference.json must contain keys: {required_keys}") + return data + + +def overlay_adc_reference(ax: plt.Axes, ref_data: dict): + """ + Draw per-vial reference ADC values as open circles onto *ax*. + Values are converted to ×10⁻³ for display. + + Open circles (no fill) are drawn after the measured scatter so both + are visible when they coincide. + """ + vials = ref_data["vials"] + ref_vals = np.array([ref_data["adc_mm2_per_s"][v] for v in vials]) + x = np.arange(len(vials)) + ref_display = ref_vals * 1e3 # convert to ×10⁻³ for display + + # Open circle — no fill, steelblue edge, 1.5× the measured markersize (7→10.5) + ax.scatter( + x, + ref_display, + marker="o", + s=110, # approx 10.5² ≈ 110 + facecolors="none", + edgecolors="steelblue", + linewidths=1.5, + zorder=6, + ) + + +def build_adc_legend(has_measured: bool) -> list: + """Return legend handles for the ADC overlay.""" + ref_handle = plt.Line2D( + [0], + [0], + marker="o", + color="none", + markeredgecolor="steelblue", + markeredgewidth=1.5, + linestyle="None", + markersize=10, + label="Reference ADC", + ) + handles = [ref_handle] + if has_measured: + meas_handle = plt.Line2D( + [0], + [0], + marker="o", + color="crimson", + linestyle="None", + markersize=7, + label="Mean (SD) ADC", + ) + handles.append(meas_handle) + return handles + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + + +def main(): + parser = argparse.ArgumentParser( + description=( + "Plot vial vs intensity (mean ± std) for 3D or 4D contrasts, " + "with optional ROI overlay.\n\n" + "Plot mode (ADC / FA / generic) is detected automatically from " + "the csv_file filename — no flags required." + ) + ) + # ---- positional -------------------------------------------------------- + parser.add_argument( + "csv_file", + help="CSV file containing mean values (vial, mean[, vol1, vol2...]).", + ) + parser.add_argument( + "plot_type", + choices=["line", "bar", "scatter"], + help="Type of plot to generate.", + ) + # ---- standard options -------------------------------------------------- + parser.add_argument( + "--std_csv", + help="Optional CSV file containing standard deviations (same shape as mean CSV).", + ) + parser.add_argument( + "--roi_image", + help="Optional PNG image (e.g. mrview screenshot with ROI overlay).", + ) + parser.add_argument( + "--annotate", + action="store_true", + help="Annotate points with mean ± std values.", + ) + parser.add_argument( + "--output", + default="vial_subplot.png", + help="Filename for saving the plot (default: vial_subplot.png).", + ) + # ---- phantom options (used when ADC mode is auto-detected) ------------- + parser.add_argument( + "--phantom", + default=None, + help=( + "Phantom name (e.g. SPIRIT). Used to locate ADC reference data " + "and to label plot titles. Required when the csv_file name " + "contains 'ADC'." + ), + ) + parser.add_argument( + "--template_dir", + default=None, + help=( + "Path to the TemplateData directory. " + "ADC reference is read from //adc_reference.json. " + "Required when the csv_file name contains 'ADC'." + ), + ) + + args = parser.parse_args() + plot_vial_intensity( + csv_file=args.csv_file, + plot_type=args.plot_type, + std_csv=args.std_csv, + roi_image=args.roi_image, + annotate=args.annotate, + output=args.output, + phantom=args.phantom, + template_dir=args.template_dir, + ) + + +# --------------------------------------------------------------------------- +# Public API +# --------------------------------------------------------------------------- + + def plot_vial_intensity( csv_file: str, - plot_type: str, + plot_type: str = "scatter", std_csv: str | None = None, roi_image: str | None = None, annotate: bool = False, output: str = "vial_subplot.png", -) -> str: - """Plot vial vs intensity (mean ± std) for 3D or 4D contrasts, with optional ROI overlay. - - Args: - csv_file: CSV file containing mean values (vial, mean[, vol1, vol2...]). - plot_type: Type of plot to generate: 'line', 'bar', or 'scatter'. - std_csv: Optional CSV file containing standard deviations (same shape as mean). - roi_image: Optional PNG image (e.g. mrview screenshot with ROI overlay). - annotate: Annotate points with mean ± std values. - output: Filename for saving the plot. - - Returns: - Absolute path to the saved plot file. + phantom: str | None = None, + template_dir: str | None = None, +): + """Plot vial vs intensity. Mode (ADC/FA/generic) auto-detected from csv_file name. + + Can be called programmatically or via the CLI entry point (main()). """ + + # ---- Auto-detect contrast mode from csv_file filename ------------------ + csv_stem = Path(csv_file).stem + if re.search(r"ADC", csv_stem, re.IGNORECASE): + contrast_mode = "adc" + elif re.search(r"(? 1 else "Mean (SD) ADC" + else: + color = cmap(vol_idx % 10) + fmt_line = "-o" + fmt_scatter = "o" + marker_kw = {} + vol_label = f"Vol {vol_idx}" if plot_type == "line": ax.errorbar( vials, means, yerr=stds, - fmt="-o", + fmt=fmt_line, capsize=5, color=color, - label=f"Vol {vol_idx}", + label=vol_label, + **marker_kw, ) elif plot_type == "bar": x = np.arange(len(vials)) + (vol_idx - n_vols / 2) * 0.1 ax.bar( - x, - means, - yerr=stds, - capsize=5, - color=color, - width=0.1, - label=f"Vol {vol_idx}", + x, means, yerr=stds, capsize=5, color=color, width=0.1, label=vol_label ) ax.set_xticks(np.arange(len(vials))) ax.set_xticklabels(vials) @@ -106,10 +344,11 @@ def plot_vial_intensity( vials, means, yerr=stds, - fmt="o", + fmt=fmt_scatter, capsize=5, color=color, - label=f"Vol {vol_idx}", + label=vol_label, + **marker_kw, ) if annotate: @@ -119,18 +358,67 @@ def plot_vial_intensity( ax.text( vial, mean + (max(means) * 0.02), - f"{mean:.1f}±{std:.1f}" if stds is not None else f"{mean:.1f}", + f"{mean:.3f}±{std:.3f}" if stds is not None else f"{mean:.3f}", ha="center", fontsize=8, color=color, ) - ax.set_xlabel("Vial") - ax.set_ylabel("Intensity") - ax.set_title("Vial vs Intensity (Mean ± Std)") + # ---- ADC reference overlay (drawn after measured values so crosshairs + # are visible even when they coincide with measured data) ----------- + if ref_data is not None: + overlay_adc_reference(ax, ref_data) + + # ---- Axis labels ------------------------------------------------------- + ax.set_xlabel("Vial", fontsize=12) + + phantom_label = f"{phantom} Phantom – " if phantom else "" + + if contrast_mode == "adc": + ax.set_ylabel("ADC ×10⁻³ mm²/s", fontsize=12) + ax.set_title( + f"{phantom_label}Measured vs Reference ADC", fontsize=13, fontweight="bold" + ) + # Plain numeric ticks (values are already in ×10⁻³) + ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda val, _: f"{val:.2f}")) + elif contrast_mode == "fa": + ax.set_ylabel("Fractional Anisotropy", fontsize=12) + ax.set_title( + f"{phantom_label}Fractional Anisotropy (Mean ± Std)", + fontsize=13, + fontweight="bold", + ) + ax.set_ylim(0, 1) + else: + ax.set_ylabel("Intensity", fontsize=12) + ax.set_title(f"{phantom_label}Vial vs Intensity (Mean ± Std)", fontsize=13) + ax.grid(True, linestyle="--", alpha=0.6) - ax.legend(title="Volumes") + # ---- Legend ------------------------------------------------------------ + if contrast_mode == "adc": + extra_handles = build_adc_legend(has_measured=True) + if n_vols > 1: + vol_handles = [ + plt.Line2D( + [0], + [0], + marker="o", + color="crimson", + linestyle="-" if plot_type == "line" else "None", + markersize=7, + label=f"Vol {i}", + alpha=0.4 + 0.6 * i / max(n_vols - 1, 1), + ) + for i in range(n_vols) + ] + ax.legend(handles=extra_handles + vol_handles, fontsize=10) + else: + ax.legend(handles=extra_handles, fontsize=10) + else: + ax.legend(title="Volumes", fontsize=10) + + # ---- ROI overlay subplot ----------------------------------------------- if roi_image: ax_img = axes[1] img = mpimg.imread(roi_image) @@ -141,30 +429,24 @@ def plot_vial_intensity( plt.tight_layout() output_file = os.path.abspath(output) plt.savefig(output_file, bbox_inches="tight", dpi=300) - logger.info("Plot saved to: %s", output_file) + print(f"[INFO] Plot saved to: {output_file}") plt.close(fig) return output_file -@click.command() +@click.command("vial-intensity") @click.argument("csv_file") -@click.argument("plot_type", type=click.Choice(["line", "bar", "scatter"])) -@click.option( - "--std_csv", default=None, help="Optional CSV file containing standard deviations." -) -@click.option("--roi_image", default=None, help="Optional PNG image with ROI overlay.") -@click.option( - "--annotate", is_flag=True, help="Annotate points with mean ± std values." -) -@click.option( - "--output", - default="vial_subplot.png", - show_default=True, - help="Filename for saving the plot.", -) -def main(csv_file, plot_type, std_csv, roi_image, annotate, output): - """Plot vial vs intensity (mean ± std) for 3D or 4D contrasts, with optional ROI overlay.""" - logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s") +@click.argument("plot_type", default="scatter") +@click.option("--std-csv", default=None) +@click.option("--roi-image", default=None) +@click.option("--annotate", is_flag=True, default=False) +@click.option("--output", default="vial_subplot.png") +@click.option("--phantom", default=None) +@click.option("--template-dir", default=None) +def main( + csv_file, plot_type, std_csv, roi_image, annotate, output, phantom, template_dir +): + """Plot vial vs intensity (mean ± std) for 3D or 4D contrasts.""" plot_vial_intensity( csv_file=csv_file, plot_type=plot_type, @@ -172,6 +454,8 @@ def main(csv_file, plot_type, std_csv, roi_image, annotate, output): roi_image=roi_image, annotate=annotate, output=output, + phantom=phantom, + template_dir=template_dir, ) diff --git a/phantomkit/plotting/visualization.py b/phantomkit/plotting/visualization.py index 2dba83f..4c9cab5 100644 --- a/phantomkit/plotting/visualization.py +++ b/phantomkit/plotting/visualization.py @@ -120,6 +120,14 @@ class MrView(shell.Task["MrView.Outputs"]): argstr="-capture.prefix", help="Filename prefix for the screenshot.", ) + intensity_range: str | None = shell.arg( + default=None, + argstr="-intensity_range", + help=( + "Display intensity range as 'min,max'. " + "Use '0,1' for FA maps and '0,0.005' for ADC maps." + ), + ) capture_grab: bool = shell.arg( default=True, argstr="-capture.grab", @@ -211,6 +219,35 @@ def GeneratePlots( vial_masks_list = sorted(vial_path.glob("*.nii.gz")) # triggers runtime fallback + def _classify_contrast(stem: str) -> str | None: + """Detect ADC / FA contrast types from filename stem.""" + import re as _re2 + if _re2.search(r"ADC", stem, _re2.IGNORECASE): + return "adc" + if _re2.search(r"(? str | None: + if contrast_type == "fa": + return "0,1" + if contrast_type == "adc": + return "0,0.005" + return None + + # Resolve template_dir (parent of phantom-specific dir) from metrics_dir. + # metrics_dir is //metrics/ — template_data/ is + # resolved from the installed package location. + import importlib.util as _ilu + _pkg_spec = _ilu.find_spec("phantomkit") + _pkg_dir = _Path(_pkg_spec.origin).parent if _pkg_spec else _Path(metrics_dir) + _template_data_root = None + for _candidate in [_pkg_dir.parent / "template_data", + _pkg_dir / "template_data"]: + if _candidate.is_dir(): + _template_data_root = str(_candidate) + break + # Per-contrast scatter plots for contrast_file in contrast_files: contrast_path = _Path(contrast_file) @@ -224,6 +261,9 @@ def GeneratePlots( metrics_path / f"{session_name}_{contrast_name}_PLOTmeanstd.png" ) + contrast_type = _classify_contrast(contrast_path.stem) + intensity_range_str = _intensity_range_for(contrast_type) + roi_image: File | None = None if vial_masks_list: overlay = workflow.add( @@ -244,6 +284,7 @@ def GeneratePlots( comments=0, noannotations=True, fullscreen=True, + intensity_range=intensity_range_str, capture_folder=tmp_vial_dir, capture_prefix=f"{contrast_name}_roi_overlay", ), @@ -251,6 +292,10 @@ def GeneratePlots( ) roi_image = screenshot.out_png + # Determine phantom name from metrics_dir path structure: + # //metrics/ — we don't know phantom here, + # so pass None and let plot_vial_intensity skip ADC reference if absent. + # For ADC mode to work fully, callers should use PhantomProcessor directly. workflow.add( python.define(plot_vial_intensity)( csv_file=mean_csv, @@ -258,6 +303,8 @@ def GeneratePlots( std_csv=std_csv, roi_image=roi_image, output=output_plot, + phantom=None, + template_dir=_template_data_root, ), name=f"scatter_{contrast_name}", ) @@ -267,7 +314,7 @@ def _matches(stem: str, token: str) -> bool: return bool(_re.search(rf"(? bool: """ Given per-vial [mean, std] pairs, check intensity-ranking criteria: - - No vial std > 50 + - No vial std > 60 - High-intensity vials A, O, Q are in top-5 by mean - Low-intensity vials S, D, P are in bottom-5 by mean """ vial_means = {name: vals[0] for name, vals in zip(vial_names, means_stds)} - high_std = [(n, vals[1]) for n, vals in zip(vial_names, means_stds) if vals[1] > 50] + high_std = [(n, vals[1]) for n, vals in zip(vial_names, means_stds) if vals[1] > 60] failures: list[str] = [] if high_std: diff --git a/phantomkit/tests/test_cli.py b/phantomkit/tests/test_cli.py index ec48569..1f4ff9b 100644 --- a/phantomkit/tests/test_cli.py +++ b/phantomkit/tests/test_cli.py @@ -60,7 +60,7 @@ def test_run_vial_signal_help(runner) -> None: assert "--template-dir" in result.output assert "--rotation-library-file" in result.output assert "--output-base-dir" in result.output - assert "--plugin" in result.output + assert "--worker" in result.output assert "--pattern" in result.output @@ -100,13 +100,13 @@ def test_run_vial_signal_single(runner, tmp_path) -> None: str(rot_lib), "--output-base-dir", str(tmp_path / "output"), - "--plugin", + "--worker", "serial", ], ) assert result.exit_code == 0, result.output - mock_sub.assert_called_once_with(plugin="serial") + mock_sub.assert_called_once_with(worker="serial") mock_submitter.assert_called_once() # sub(wf) @@ -142,7 +142,7 @@ def test_run_vial_signal_batch_from_dir(runner, tmp_path) -> None: str(rot_lib), "--output-base-dir", str(tmp_path / "output"), - "--plugin", + "--worker", "serial", "--pattern", "*.nii.gz", @@ -150,7 +150,7 @@ def test_run_vial_signal_batch_from_dir(runner, tmp_path) -> None: ) assert result.exit_code == 0, result.output - mock_sub.assert_called_once_with(plugin="serial") + mock_sub.assert_called_once_with(worker="serial") mock_submitter.assert_called_once() # sub(wf) @@ -202,11 +202,11 @@ def test_plot_maps_ir_help(runner) -> None: result = runner.invoke(main, ["plot", "maps-ir", "--help"]) assert result.exit_code == 0 assert "CONTRAST_FILES" in result.output - assert "--metric_dir" in result.output + assert "--metric-dir" in result.output def test_plot_maps_te_help(runner) -> None: result = runner.invoke(main, ["plot", "maps-te", "--help"]) assert result.exit_code == 0 assert "CONTRAST_FILES" in result.output - assert "--metric_dir" in result.output + assert "--metric-dir" in result.output diff --git a/template_data/120E/ImageTemplate.nii.gz b/template_data/120E/ImageTemplate.nii.gz new file mode 100644 index 0000000..9eab1c8 Binary files /dev/null and b/template_data/120E/ImageTemplate.nii.gz differ diff --git a/template_data/120E/VialSegmentationsCombined.nii.gz b/template_data/120E/VialSegmentationsCombined.nii.gz new file mode 100644 index 0000000..cd97947 Binary files /dev/null and b/template_data/120E/VialSegmentationsCombined.nii.gz differ diff --git a/template_data/120E/adc_reference.json b/template_data/120E/adc_reference.json new file mode 100644 index 0000000..8e4f0ac --- /dev/null +++ b/template_data/120E/adc_reference.json @@ -0,0 +1,29 @@ +{ + "vials": ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M", "N", "O", "P", "Q", "R", "S", "T", "U", "V", "W", "X"], +"adc_mm2_per_s": { + "A": 1.0e-3, + "B": 0.84e-3, + "C": 0.54e-3, + "D": 0.32e-3, + "E": 0.84e-3, + "F": 0.54e-3, + "G": 1.82e-3, + "H": 1.59e-3, + "I": 1.37e-3, + "J": 1.18e-3, + "K": 2.02e-3, + "L": 1.91e-3, + "M": 1.82e-3, + "N": 1.59e-3, + "O": 1.37e-3, + "P": 1.18e-3, + "Q": 1.0e-3, + "R": 0.84e-3, + "S": 0.54e-3, + "T": 0.32e-3, + "U": 2.02e-3, + "V": 1.91e-3, + "W": 1.18e-3, + "X": 1.0e-3 + } +} diff --git a/template_data/120E/vials_labelled/A.nii.gz b/template_data/120E/vials_labelled/A.nii.gz new file mode 100644 index 0000000..4bb824f Binary files /dev/null and b/template_data/120E/vials_labelled/A.nii.gz differ diff --git a/template_data/120E/vials_labelled/B.nii.gz b/template_data/120E/vials_labelled/B.nii.gz new file mode 100644 index 0000000..b528efc Binary files /dev/null and b/template_data/120E/vials_labelled/B.nii.gz differ diff --git a/template_data/120E/vials_labelled/C.nii.gz b/template_data/120E/vials_labelled/C.nii.gz new file mode 100644 index 0000000..3a61671 Binary files /dev/null and b/template_data/120E/vials_labelled/C.nii.gz differ diff --git a/template_data/120E/vials_labelled/D.nii.gz b/template_data/120E/vials_labelled/D.nii.gz new file mode 100644 index 0000000..1ab992d Binary files /dev/null and b/template_data/120E/vials_labelled/D.nii.gz differ diff --git a/template_data/120E/vials_labelled/E.nii.gz b/template_data/120E/vials_labelled/E.nii.gz new file mode 100644 index 0000000..98f446d Binary files /dev/null and b/template_data/120E/vials_labelled/E.nii.gz differ diff --git a/template_data/120E/vials_labelled/F.nii.gz b/template_data/120E/vials_labelled/F.nii.gz new file mode 100644 index 0000000..6af536f Binary files /dev/null and b/template_data/120E/vials_labelled/F.nii.gz differ diff --git a/template_data/120E/vials_labelled/G.nii.gz b/template_data/120E/vials_labelled/G.nii.gz new file mode 100644 index 0000000..1220ba7 Binary files /dev/null and b/template_data/120E/vials_labelled/G.nii.gz differ diff --git a/template_data/120E/vials_labelled/H.nii.gz b/template_data/120E/vials_labelled/H.nii.gz new file mode 100644 index 0000000..3b36bfb Binary files /dev/null and b/template_data/120E/vials_labelled/H.nii.gz differ diff --git a/template_data/120E/vials_labelled/I.nii.gz b/template_data/120E/vials_labelled/I.nii.gz new file mode 100644 index 0000000..a636517 Binary files /dev/null and b/template_data/120E/vials_labelled/I.nii.gz differ diff --git a/template_data/120E/vials_labelled/J.nii.gz b/template_data/120E/vials_labelled/J.nii.gz new file mode 100644 index 0000000..8f191a6 Binary files /dev/null and b/template_data/120E/vials_labelled/J.nii.gz differ diff --git a/template_data/120E/vials_labelled/K.nii.gz b/template_data/120E/vials_labelled/K.nii.gz new file mode 100644 index 0000000..f0894f6 Binary files /dev/null and b/template_data/120E/vials_labelled/K.nii.gz differ diff --git a/template_data/120E/vials_labelled/L.nii.gz b/template_data/120E/vials_labelled/L.nii.gz new file mode 100644 index 0000000..321ceb8 Binary files /dev/null and b/template_data/120E/vials_labelled/L.nii.gz differ diff --git a/template_data/120E/vials_labelled/M.nii.gz b/template_data/120E/vials_labelled/M.nii.gz new file mode 100644 index 0000000..a044865 Binary files /dev/null and b/template_data/120E/vials_labelled/M.nii.gz differ diff --git a/template_data/120E/vials_labelled/N.nii.gz b/template_data/120E/vials_labelled/N.nii.gz new file mode 100644 index 0000000..bfe641d Binary files /dev/null and b/template_data/120E/vials_labelled/N.nii.gz differ diff --git a/template_data/120E/vials_labelled/O.nii.gz b/template_data/120E/vials_labelled/O.nii.gz new file mode 100644 index 0000000..d79be31 Binary files /dev/null and b/template_data/120E/vials_labelled/O.nii.gz differ diff --git a/template_data/120E/vials_labelled/P.nii.gz b/template_data/120E/vials_labelled/P.nii.gz new file mode 100644 index 0000000..b01d4ce Binary files /dev/null and b/template_data/120E/vials_labelled/P.nii.gz differ diff --git a/template_data/120E/vials_labelled/Q.nii.gz b/template_data/120E/vials_labelled/Q.nii.gz new file mode 100644 index 0000000..0e0ccb4 Binary files /dev/null and b/template_data/120E/vials_labelled/Q.nii.gz differ diff --git a/template_data/120E/vials_labelled/R.nii.gz b/template_data/120E/vials_labelled/R.nii.gz new file mode 100644 index 0000000..eff1c61 Binary files /dev/null and b/template_data/120E/vials_labelled/R.nii.gz differ diff --git a/template_data/120E/vials_labelled/S.nii.gz b/template_data/120E/vials_labelled/S.nii.gz new file mode 100644 index 0000000..1bfd133 Binary files /dev/null and b/template_data/120E/vials_labelled/S.nii.gz differ diff --git a/template_data/120E/vials_labelled/T.nii.gz b/template_data/120E/vials_labelled/T.nii.gz new file mode 100644 index 0000000..d2274a0 Binary files /dev/null and b/template_data/120E/vials_labelled/T.nii.gz differ diff --git a/template_data/120E/vials_labelled/U.nii.gz b/template_data/120E/vials_labelled/U.nii.gz new file mode 100644 index 0000000..744ea92 Binary files /dev/null and b/template_data/120E/vials_labelled/U.nii.gz differ diff --git a/template_data/120E/vials_labelled/V.nii.gz b/template_data/120E/vials_labelled/V.nii.gz new file mode 100644 index 0000000..87abe39 Binary files /dev/null and b/template_data/120E/vials_labelled/V.nii.gz differ diff --git a/template_data/120E/vials_labelled/W.nii.gz b/template_data/120E/vials_labelled/W.nii.gz new file mode 100644 index 0000000..ed1a95b Binary files /dev/null and b/template_data/120E/vials_labelled/W.nii.gz differ diff --git a/template_data/120E/vials_labelled/X.nii.gz b/template_data/120E/vials_labelled/X.nii.gz new file mode 100644 index 0000000..1f64b49 Binary files /dev/null and b/template_data/120E/vials_labelled/X.nii.gz differ diff --git a/template_data/SPIRIT/GSP_ADCvalues.xlsx b/template_data/SPIRIT/GSP_ADCvalues.xlsx new file mode 100644 index 0000000..7c25fc1 Binary files /dev/null and b/template_data/SPIRIT/GSP_ADCvalues.xlsx differ diff --git a/template_data/ImageTemplate.nii.gz b/template_data/SPIRIT/ImageTemplate.nii.gz similarity index 100% rename from template_data/ImageTemplate.nii.gz rename to template_data/SPIRIT/ImageTemplate.nii.gz diff --git a/template_data/VialSegmentationsCombined.nii.gz b/template_data/SPIRIT/VialSegmentationsCombined.nii.gz similarity index 100% rename from template_data/VialSegmentationsCombined.nii.gz rename to template_data/SPIRIT/VialSegmentationsCombined.nii.gz diff --git a/template_data/SPIRIT/adc_reference.json b/template_data/SPIRIT/adc_reference.json new file mode 100644 index 0000000..5a66fed --- /dev/null +++ b/template_data/SPIRIT/adc_reference.json @@ -0,0 +1,13 @@ +{ + "vials": ["E", "F", "G", "H", "I", "J", "K", "L"], +"adc_mm2_per_s": { + "E": 1.82e-3, + "F": 1.59e-3, + "G": 1.37e-3, + "H": 1.18e-3, + "I": 1.00e-3, + "J": 0.84e-3, + "K": 0.54e-3, + "L": 0.32e-3 + } +} \ No newline at end of file diff --git a/template_data/vials_labelled/A.nii.gz b/template_data/SPIRIT/vials_labelled/A.nii.gz similarity index 100% rename from template_data/vials_labelled/A.nii.gz rename to template_data/SPIRIT/vials_labelled/A.nii.gz diff --git a/template_data/vials_labelled/B.nii.gz b/template_data/SPIRIT/vials_labelled/B.nii.gz similarity index 100% rename from template_data/vials_labelled/B.nii.gz rename to template_data/SPIRIT/vials_labelled/B.nii.gz diff --git a/template_data/vials_labelled/C.nii.gz b/template_data/SPIRIT/vials_labelled/C.nii.gz similarity index 100% rename from template_data/vials_labelled/C.nii.gz rename to template_data/SPIRIT/vials_labelled/C.nii.gz diff --git a/template_data/vials_labelled/D.nii.gz b/template_data/SPIRIT/vials_labelled/D.nii.gz similarity index 100% rename from template_data/vials_labelled/D.nii.gz rename to template_data/SPIRIT/vials_labelled/D.nii.gz diff --git a/template_data/vials_labelled/E.nii.gz b/template_data/SPIRIT/vials_labelled/E.nii.gz similarity index 100% rename from template_data/vials_labelled/E.nii.gz rename to template_data/SPIRIT/vials_labelled/E.nii.gz diff --git a/template_data/vials_labelled/F.nii.gz b/template_data/SPIRIT/vials_labelled/F.nii.gz similarity index 100% rename from template_data/vials_labelled/F.nii.gz rename to template_data/SPIRIT/vials_labelled/F.nii.gz diff --git a/template_data/vials_labelled/G.nii.gz b/template_data/SPIRIT/vials_labelled/G.nii.gz similarity index 100% rename from template_data/vials_labelled/G.nii.gz rename to template_data/SPIRIT/vials_labelled/G.nii.gz diff --git a/template_data/vials_labelled/H.nii.gz b/template_data/SPIRIT/vials_labelled/H.nii.gz similarity index 100% rename from template_data/vials_labelled/H.nii.gz rename to template_data/SPIRIT/vials_labelled/H.nii.gz diff --git a/template_data/vials_labelled/I.nii.gz b/template_data/SPIRIT/vials_labelled/I.nii.gz similarity index 100% rename from template_data/vials_labelled/I.nii.gz rename to template_data/SPIRIT/vials_labelled/I.nii.gz diff --git a/template_data/vials_labelled/J.nii.gz b/template_data/SPIRIT/vials_labelled/J.nii.gz similarity index 100% rename from template_data/vials_labelled/J.nii.gz rename to template_data/SPIRIT/vials_labelled/J.nii.gz diff --git a/template_data/vials_labelled/K.nii.gz b/template_data/SPIRIT/vials_labelled/K.nii.gz similarity index 100% rename from template_data/vials_labelled/K.nii.gz rename to template_data/SPIRIT/vials_labelled/K.nii.gz diff --git a/template_data/vials_labelled/L.nii.gz b/template_data/SPIRIT/vials_labelled/L.nii.gz similarity index 100% rename from template_data/vials_labelled/L.nii.gz rename to template_data/SPIRIT/vials_labelled/L.nii.gz diff --git a/template_data/vials_labelled/M.nii.gz b/template_data/SPIRIT/vials_labelled/M.nii.gz similarity index 100% rename from template_data/vials_labelled/M.nii.gz rename to template_data/SPIRIT/vials_labelled/M.nii.gz diff --git a/template_data/vials_labelled/N.nii.gz b/template_data/SPIRIT/vials_labelled/N.nii.gz similarity index 100% rename from template_data/vials_labelled/N.nii.gz rename to template_data/SPIRIT/vials_labelled/N.nii.gz diff --git a/template_data/vials_labelled/O.nii.gz b/template_data/SPIRIT/vials_labelled/O.nii.gz similarity index 100% rename from template_data/vials_labelled/O.nii.gz rename to template_data/SPIRIT/vials_labelled/O.nii.gz diff --git a/template_data/vials_labelled/P.nii.gz b/template_data/SPIRIT/vials_labelled/P.nii.gz similarity index 100% rename from template_data/vials_labelled/P.nii.gz rename to template_data/SPIRIT/vials_labelled/P.nii.gz diff --git a/template_data/vials_labelled/Q.nii.gz b/template_data/SPIRIT/vials_labelled/Q.nii.gz similarity index 100% rename from template_data/vials_labelled/Q.nii.gz rename to template_data/SPIRIT/vials_labelled/Q.nii.gz diff --git a/template_data/vials_labelled/R.nii.gz b/template_data/SPIRIT/vials_labelled/R.nii.gz similarity index 100% rename from template_data/vials_labelled/R.nii.gz rename to template_data/SPIRIT/vials_labelled/R.nii.gz diff --git a/template_data/vials_labelled/S.nii.gz b/template_data/SPIRIT/vials_labelled/S.nii.gz similarity index 100% rename from template_data/vials_labelled/S.nii.gz rename to template_data/SPIRIT/vials_labelled/S.nii.gz diff --git a/template_data/vials_labelled/T.nii.gz b/template_data/SPIRIT/vials_labelled/T.nii.gz similarity index 100% rename from template_data/vials_labelled/T.nii.gz rename to template_data/SPIRIT/vials_labelled/T.nii.gz diff --git a/vial_subplot.png b/vial_subplot.png deleted file mode 100644 index ac4e996..0000000 Binary files a/vial_subplot.png and /dev/null differ