GhostTracts

Physics-Aware Diffusion MRI Reconstruction & Tractography Stress Testing

“When MRI acquisition and reconstruction go wrong, how do those errors propagate into diffusion biomarkers and apparent white-matter tracts?”

GhostTracts is compact research software for controlled diffusion-MRI error-propagation experiments. It begins with analytic fiber anatomy, simulates diffusion contrast and slice-wise multi-coil k-space, applies interpretable MRI/EPI perturbations, reconstructs with conventional algorithms, and measures the consequences for DTI, CSD, fiber orientation, and tractography. The synthetic phantom supplies the known ground truth needed to distinguish an image that looks acceptable from one that preserves the intended diffusion geometry.

The core pipeline is deterministic, CPU compatible, and designed to run within 8 GB RAM. Volumes are synthesized and acquired sequentially; k-space is stored one volume at a time as complex64 instead of materializing all coils and directions together.

Motivation

MRI reconstruction is often ranked with image-space quantities such as RMSE, NRMSE, or PSNR. Diffusion MRI, however, is used to infer tensor-derived tissue properties, local fiber orientations, and long-range pathways. GhostTracts asks whether image-space reconstruction quality predicts the validity of those downstream scientific conclusions.

The framework uses controlled synthetic anatomy so every acquisition starts from the same fiber geometry. Errors can then be traced through signal formation, encoding, reconstruction, diffusion modeling, and tracking without treating a noisy estimate as truth. The aim is not a vendor scanner replica or clinical validation; it is isolation of mechanisms in a reproducible methods experiment.

Pipeline

flowchart LR
    A["Analytic fiber phantom"] --> B["Diffusion signal"]
    B --> C["Complex coil images"]
    C --> D["Sampled k-space"]
    D --> E["MRI / EPI artifacts"]
    E --> F["Conventional reconstruction"]
    F --> G["DTI / CSD"]
    G --> H["Fiber orientation"]
    H --> I["Tractography"]
    I --> J["Image, diffusion, and tract metrics"]

Loading

The default phantom is 96 × 96 × 32 at 2 mm isotropic resolution, with one b=0 volume and 32 approximately uniform b=1000 s/mm² directions. A 64-direction shell is also supported. The quick study uses 24 × 24 × 16, retains 32 directions for DTI/CSD, and exercises 13 conditions.

Physics and signal model

Each single-fiber compartment uses an axisymmetric diffusion tensor with configurable axial and radial diffusivities. The defaults are 1.7 × 10⁻³ and 0.3 × 10⁻³ mm²/s. For a unit gradient direction g, diffusion weighting b, and tensor D, the Stejskal–Tanner signal is

$$ S(b,\mathbf g)=S_0\exp!\left(-b,\mathbf g^T\mathbf D\mathbf g\right). $$

A crossing is not represented by a fictitious single tensor. Its signal is a normalized, positive multi-compartment mixture:

$$ S(b,\mathbf g)=S_0\sum_i w_i\exp!\left(-b,\mathbf g^T\mathbf D_i\mathbf g\right), \qquad w_i>0,\quad \sum_i w_i=1. $$

The analytic phantom contains a straight left–right bundle, an oblique bundle, a curved bundle, and a known straight/oblique crossing. Ground truth includes the tissue mask, bundle-label bitmask, up to two fiber directions and weights per voxel, compartment tensor coefficients, expected single-fiber FA/MD, and analytic reference streamlines. Expected tensor metrics are undefined in the background and crossing voxels rather than being assigned a misleading value.

For receive coil c, the coil image and sampled data are

$$ I_c(\mathbf r)=C_c(\mathbf r)S(\mathbf r), $$

$$ y_c=P,\mathcal F{I_c}+\varepsilon, $$

where $C_c$ is a smooth complex sensitivity, $\mathcal F$ is a centered orthonormal 2D FFT, $P$ is the Cartesian sampling operator, and $\varepsilon$ is complex receiver noise. Spatial phase/distortion operators are inserted before encoding; line-dependent ghost and relaxation operators act during k-space formation.

NIfTI images use a centered RAS+ affine with millimetre units. FSL-style b-vectors are saved as 3 × N rows in image x/y/z axes; b=0 vectors are zero and diffusion vectors are unit length.

Simulated acquisition and artifacts

GhostTracts describes its acquisition stage as a physics-informed EPI-style forward model for controlled reconstruction experiments. It uses independent axial 2D encoding rather than a complete pulse-sequence, gradient-system, or vendor raw-data model.

• Receive coils: analytic circular arrays with smooth magnitude and phase; sensitivities are normalized pointwise so $\sum_c|C_c|^2=1$.
• Sampling: full Cartesian, uniform R=2/R=3 acceleration with optional central calibration lines, and 6/8 or 7/8 partial Fourier.
• Motion: known deterministic volume-wise rigid rotations/translations before encoding. Rotation matrices are recorded so spatial correction and b-vector rotation can be evaluated.
• Eddy-like distortion: a transparent gradient-dependent affine PE scale, shear, and translation. It is a parametric distortion, not a physical eddy-current field solver.
• B0/off-resonance: a smooth synthetic field map in Hz produces PE displacement proportional to field offset and total readout time, $\Delta p=\Delta f,T_\mathrm{readout}$ pixels.
• Nyquist ghost: an odd/even ky phase mismatch produces an N/2 phase-encode ghost.
• T2* blurring: exponential line attenuation across the simplified echo train.
• Receiver noise: independent zero-mean Gaussian noise in real and imaginary acquired k-space channels, configured by sigma or target SNR.

Every artifact has an independent enable/severity setting and deterministic seed. Provenance records coil count, sampling mask, acceleration, partial-Fourier fraction, phase-encoding direction, field-map parameters, motion/eddy transforms, artifact settings, and source phantom.

Reconstruction

All reconstruction methods operate slice by slice and volume by volume.

• Reference FFT: inverse centered FFT per coil followed by root-sum-of-squares, or complex sensitivity combination when the simulated sensitivities are available.

• Zero-filled: missing Cartesian samples remain zero before inverse FFT. This deliberately naive baseline exposes aliasing and partial-Fourier truncation.

• Iterative SENSE: a matrix-free forward operator $A(x)_c=P\mathcal F(C_cx)$ and matching adjoint are used in conjugate gradients to solve

  $$ \underset{x}{\operatorname{argmin}};|Ax-y|_2^2+\lambda|x|_2^2. $$

  No encoding matrix is constructed; residual histories, tolerance, iterations, and achieved residual are saved. The known analytic sensitivity maps make this an idealized SENSE test.

• Partial Fourier: POCS estimates phase from symmetric central k-space, alternates an image-domain phase constraint with Fourier transformation, and restores acquired samples exactly after every iteration.

Reconstruction metrics against the source DWI include RMSE, NRMSE, MAE, and PSNR per volume and in aggregate. Image-space quality is kept distinct from downstream diffusion validity.

Diffusion analysis

DIPY is the reproducible Python reference implementation.

• Weighted least-squares DTI produces FA, MD, axial diffusivity, radial diffusivity, eigenvalues, eigenvectors, tensor coefficients, and a mask as NIfTI images.
• Tensor errors are evaluated only where the phantom has a meaningful single-tensor ground truth. Principal-direction angular error uses $|\hat v\cdot v|$ so eigenvector sign is irrelevant.
• Constrained spherical deconvolution uses DIPY's official CSD implementation and a known synthetic single-shell response. Spherical-harmonic coefficients, ODF peak directions, and amplitudes are saved.
• Crossing peaks are matched to true directions by the permutation minimizing antipodally symmetric angular error. Missed-peak, spurious-peak, and crossing-specific statistics are reported. A tensor principal eigenvector is not treated as a valid crossing-fiber model.
• QC includes b0 identity/SNR, per-volume signal, outlier scores, gradient coverage, and reconstruction error when the source exists.

Tractography

Reference bundles are generated analytically from the exact phantom vector fields; they are never inferred from corrupted DWI. Estimated tracts use deterministic local tracking of the DIPY CSD field, known seed regions, a tissue mask, configurable 0.75 mm default steps, and curvature and length limits. Runs remain on the order of thousands—not millions—of streamlines.

Tractography outputs include .trk, compact .npz, and per-bundle visitation NIfTI maps. Metrics include visitation Dice/precision/recall, endpoint displacement in millimetres, bundle recovery, crossing-path recovery, and reversal-invariant symmetric nearest-neighbour MDF distance after fixed-point resampling.

A seeded streamline is accepted only when both orientation-independent endpoints enter the expected endpoint spheres and at least the configured fraction (0.85 by default) of its subsegmented points remains in the intended bundle-label corridor. It is marked spurious if it is too short, misses that endpoint pair, leaves the corridor, or connects another bundle's endpoint pair. A bundle is recovered only if it meets both the minimum accepted-streamline count and connection-rate rules. These definitions are saved with every metrics file.

Neuroimaging tools

Core I/O uses NiBabel and modeling uses DIPY. FSL and MRtrix3 are optional validation targets and are never imported as Python dependencies.

ghosttracts tools check
ghosttracts tools validate runs/recon_clean --dipy-analysis runs/analysis_clean

tools check reports paths and versions for dtifit, mrconvert, dwi2tensor, tensor2metric, dwi2response, dwi2fod, and tckgen. tools validate runs available tools, captures commands/stdout/stderr/version metadata, compares FA/MD (or ADC) and eigenvectors where available, writes JSON/CSV and difference statistics, and records an explicit skipped status when software is absent.

Neuroimaging tool validation

FSL validation prepares NIfTI/bval/bvec/mask inputs and invokes the equivalent of:

dtifit -k dwi.nii.gz -o PREFIX -m mask.nii.gz -r dwi.bvec -b dwi.bval

MRtrix3 validation uses mrconvert -fslgrad, dwi2tensor, and tensor2metric; response/FOD and tracking commands are optional flags. The wrappers never hard-code an installation path.

Reproducibility

Python 3.11 or newer is supported. CPU execution is the reference path; no GPU or learned model is required.

Reviewer path

git clone https://github.com/jiff3/ghosttracts.git
cd ghosttracts
python -m pip install -e .
ghosttracts experiment run configs/study_quick.yaml --output runs/demo
ghosttracts report runs/demo --output reports/demo

The same workflow with the default-resolution study is:

ghosttracts experiment run configs/study.yaml --output runs/study
ghosttracts report runs/study --output reports/study

The experiment runner hashes canonical stage configuration plus upstream hashes. Identical phantom/acquisition/reconstruction products are cached; a config change invalidates the affected stage and its descendants. A failed condition remains an explicit failed row without corrupting the rest of the table.

The report contains:

reports/demo/
├── report.md
├── metrics.csv
├── summary.json
├── reproducibility_manifest.json
└── figures/
    ├── end_to_end_pipeline.png
    ├── artifact_comparison_grid.png
    ├── nrmse_vs_angular_error.png
    ├── nrmse_vs_tract_error.png
    ├── motion_bvec_ablation.png
    ├── reconstruction_method_comparison.png
    └── clean_vs_corrupted_tractography.png

The manifest records the Git commit when available, dirty state, UTC timestamp, Python and package versions, platform, exact YAML, all seeds, input size, and detected FSL/MRtrix versions. Generated runs, downloaded data, NIfTI, k-space NPZ files, and reports are ignored by Git.

Individual stages

ghosttracts phantom generate --config configs/baseline.yaml --output runs/baseline
ghosttracts acquire runs/baseline --config configs/acquisition_clean.yaml --output runs/acq_clean
ghosttracts reconstruct runs/acq_clean --method reference --output runs/recon_clean
ghosttracts analyze runs/recon_clean --output runs/analysis_clean
ghosttracts tractography runs/analysis_clean --output runs/tract_clean

Useful development commands are make install, make test, make study-quick, make report-quick, and make check. make check runs Ruff, the full test suite, the quick study, and report generation.

Real-data demonstration

The optional real-data path uses DIPY's official Stanford HARDI single-subject fetch/read API. Stanford HARDI starts from reconstructed magnitude DWI, not vendor raw k-space. GhostTracts selects a compact gradient/spatial subset and then performs:

reconstructed real DWI
→ forward-simulated coil/k-space acquisition
→ controlled undersampling/artifacts
→ conventional reconstruction
→ reference-relative diffusion analysis

It therefore tests the pipeline on realistic anatomy and diffusion contrast but does not reproduce the scanner's original acquisition. There is no true fiber geometry, so metrics are labelled reference-relative: image NRMSE, FA/MD difference, high-FA principal-direction disagreement, ODF peak disagreement where available, and tract visitation similarity—not absolute tract accuracy.

The compact preprocessing keeps one b0 and 32 main-shell directions. It starts from the direction farthest from the current antipodal set and greedily maximizes the minimum antipodal angle $\arccos(|g_i\cdot g_j|)$, rather than taking the first volumes blindly. Data are cropped to at most 96 × 96 × 32, stored as float32, masked by robust thresholding/morphology, and processed per volume/slice.

ghosttracts realdata path
ghosttracts realdata fetch --data-home data/dipy
ghosttracts realdata stanford --data-home data/dipy --output runs/stanford_hardi

Downloads remain in the local DIPY cache and are not committed. Network or fetch failures give a direct, non-destructive error message and do not affect the synthetic pipeline.

Results

The following values were generated from configs/study_quick.yaml with seed 2025 on 2026-08-07; they are not hand-selected targets. All 13 conditions completed. This compact study reported strong positive rank correlations between reconstruction NRMSE and median principal-direction error (Spearman ρ=0.916, p=1.09×10⁻⁵) and between NRMSE and spurious tract rate (ρ=0.959, p=2.41×10⁻⁷). In this run, image error generally tracked downstream error; the results do not support a broad claim that NRMSE and tract validity are unrelated.

The report's predeclared illustrative-pair rule—NRMSE separation no more than 10% of the observed range and spurious-rate separation at least 0.15—identified clean versus SNR 20: ΔNRMSE=0.0664 and Δspurious rate=0.210. This is reported as a rule-qualified comparison, not evidence of clinical divergence.

Condition	NRMSE	Median direction error	Bundle Dice	Spurious rate
clean	2.11×10⁻⁷	0.000001°	0.786	0.373
SNR 20	0.0664	1.26°	0.750	0.583
R2 zero-fill	0.216	0.163°	0.705	0.701
R2 SENSE	6.30×10⁻⁶	0.000056°	0.786	0.373
R3 zero-fill	0.273	0.113°	0.708	0.689
R3 SENSE	0.000618	0.00244°	0.786	0.375
6/8 zero-fill	0.0921	0.0255°	0.740	0.508
6/8 POCS	0.0182	0.00464°	0.781	0.381
high motion	0.683	50.8°	0.618	0.996
severe combined artifacts	0.743	44.8°	0.530	1.000

Known-motion image correction without b-vector rotation left 4.10° median direction error; rotating the gradients reduced it to 1.79° at the same image NRMSE (0.258). The b-vector-aware case did not improve every tract summary—mean bundle Dice was 0.743 versus 0.760 without rotation in this small tracking run—so the report preserves both diffusion and tract results rather than claiming uniform improvement. Even the clean deterministic tracking baseline has a nonzero spurious rate (0.373), which reflects the encoded tracking/endpoint rules and underscores that tractography itself is an imperfect inference stage.

Complete per-condition values and automatically generated observations are written by ghosttracts report; they are intentionally not committed as large run products.

Limitations

• Slice-wise Cartesian encoding is an EPI-style abstraction, not a complete scanner or pulse-sequence simulation.
• Coil sensitivities are analytical and known exactly; no vendor calibration, coil compression, noise covariance estimation, gradient nonlinearity, or scanner raw-data format is modeled.
• Field maps and eddy behavior are synthetic parameterizations. Susceptibility, eddy, motion, relaxation, and readout interactions are not a full coupled physical system.
• SENSE receives ideal sensitivity maps, and POCS phase is estimated under controlled partial Fourier sampling.
• The phantom has limited bundle complexity and tissue compartments. DTI truth is deliberately excluded from crossings.
• CSD response and tractography parameters are controlled. Tractography is an indirect model of white matter, can fail on clean data, and is not anatomical proof.
• Stanford HARDI has no known tract ground truth and no original scanner k-space; every real-data comparison is relative to its selected reconstructed DWI.
• Results and correlations apply to the saved experiment matrix and seeds, not patients or clinical decision-making.

Repository layout

ghosttracts/
├── phantom/          # analytic geometry, tensors, signal formation
├── acquisition/      # coils, centered FFT, masks, forward encoding
├── artifacts/        # motion, off-resonance, eddy, ghost, T2*, noise
├── reconstruction/   # RSS/combine, zero-fill, SENSE, POCS
├── neuroimaging/     # gradients, DIPY DTI/CSD, QC, external tools
├── tractography/     # seeds, analytic truth, tracking, explicit metrics
├── experiments/      # matrix expansion, caching, correction, aggregation
├── realdata/         # Stanford HARDI fetching and compact preprocessing
└── reporting/        # grounded summaries, manifests, publication figures
configs/              # baseline, acquisition, reconstruction, study YAML
scripts/              # plotting and optional validation entry points
tests/                # physics, numerical, neuroimaging, and integration tests

GhostTracts is released under the MIT License. Citation metadata are provided in CITATION.cff.