Implicit Neural Representation for Spherical Harmonics Reconstruction of Motion-Corrupted Fetal Diffusion MRI

Conference: ECCV 2026
Paper: CVF Open Access
Code: https://github.com/baby-MedIA/INR-dMRI-Recon
Area: Medical Imaging
Keywords: Diffusion MRI, Fetal Brain Imaging, Implicit Neural Representation, Spherical Harmonics, Slice-to-Volume Reconstruction

TL;DR

Addressing severe motion corruption and single-orientation stripe artifacts in fetal diffusion MRI, this paper proposes the first self-supervised implicit neural representation framework parameterized by continuous spherical harmonics, achieving high-fidelity 3D reconstruction and sharp tractography via physics-based acquisition modeling and through-plane continuity regularization.

Background & Motivation

In-vivo fetal diffusion Magnetic Resonance Imaging (dMRI) probes the restricted microscopic displacement of water molecules within cellular compartments, providing unprecedented potential for non-invasively tracking early axonal growth and white matter pathway maturation in the womb. However, acquiring diagnostic-quality 3D fetal dMRI volumes is notoriously plagued by unpredictable and rapid fetal head motion. While ultrafast single-shot Echo Planar Imaging (EPI) effectively freezes motion during each 2D slice readout, inter-slice motion induces severe spatial dislocations and geometric distortions, corrupting whole 3D stacks and rendering them fundamentally unusable for quantitative microstructural modeling.

Slice-to-volume reconstruction (SVR) frameworks have been widely explored to reconstruct motion-corrected continuous 3D volumes from scattered 2D slices. Conventional model-based iterative approaches (such as IRTK) alternate between multi-stage slice-to-volume registration and scattered discrete data interpolation. These classical algorithms not only suffer from crippling computational bottlenecks, but their fixed voxel grids also struggle to represent continuous directional diffusion signals without resolution degradation. Supervised deep learning methods offer accelerated inference, but their reliance on motion-free ground truth volumes makes them impractical in in-utero fetal regimes where clean reference volumes cannot be acquired. While self-supervised Implicit Neural Representations (INR) have recently overcome data dependencies in structural MRI (e.g., T2-weighted SVR), existing architectures model only scalar intensity fields and fail on diffusion MRI.

Diffusion MRI presents unique physical complexities: each spatial coordinate represents a spherical continuous diffusion function parameterized by the orientation of the applied diffusion-encoding gradient vector, and fetal head rotations dynamically disrupt the physical alignment between the tissue's microscopic architecture and the scanner gradient axes. Moreover, scan time constraints in pregnant subjects typically allow only a single slice orientation per gradient direction, unlike the orthogonal multi-stack geometries common in structural MRI. Lacking cross-sectional geometric intersections, implicit networks readily overfit inter-slice discontinuities as genuine high-frequency anatomy, manifesting as severe stripe artifacts aligned with the slice-acquisition plane. The core idea is to model the continuous diffusion signal as a continuous field of even-degree spherical harmonics coefficients parameterized by a coordinate-based hash grid, optimize motion, gradient reorientation, and bias fields via a forward physical acquisition model, and penalize anisotropic through-plane gradient surges to eliminate striping artifacts.

Method

Overall Architecture

The framework reconstructs a spatially and angularly continuous, motion-free 3D diffusion volume \(V(\mathbf{x}, \mathbf{g})\) from a set of motion-corrupted, slice-wise 2D stacks acquired across multiple gradient directions \(\mathbf{g}_k\). The inputs are 2D slice coordinates and their nominal gradient directions, where each slice has undergone an independent rigid transformation \(T_{k,j}\). Continuous spatial coordinates are encoded via a multi-resolution hash grid and passed into the primary network \(\text{MLP}_V\) to predict direction-independent spherical harmonics (SH) coefficients and latent spatial features. Conditioned on the motion-corrected local diffusion gradient direction, the continuous diffusion signal is synthesized via an analytic linear combination of real-valued orthonormal SH basis functions. To match acquired scanner slices, the synthetic volume undergoes a physics-based forward acquisition pipeline comprising low-frequency multiplicative bias field modulation, 3D Gaussian point spread function (PSF) slice-profile blurring, and uncertainty weighting within an end-to-end negative log-likelihood objective. Finally, through-plane continuity regularization is enforced on the canonical anatomical space to explicitly suppress striping artifacts arising from single-orientation acquisitions.

%%{init: {'flowchart': {'rankSpacing': 24, 'nodeSpacing': 28, 'padding': 6, 'wrappingWidth': 400}}}%%
flowchart TD
    A["Corrupted Slice Stacks Input<br/>Gradient directions and slice coordinates"] --> B["Slice Motion and Gradient Dynamic Reorientation<br/>6-DoF pose estimation and inverse gradient rotation"]
    B --> C["Continuous Spherical Harmonics Field Representation<br/>Hash grid encoding and SH coefficient prediction"]
    C --> D["Physics-Based Forward Slice Acquisition<br/>Low-frequency bias modulation and 3D Gaussian PSF blur"]
    D --> E["Through-Plane Continuity Regularization<br/>Ratio penalization of through-plane vs 3D spatial gradient"]
    E --> F["High-Fidelity Continuous Diffusion Output<br/>Direct tractography and fiber orientation distribution"]

Key Designs

1. Slice Motion and Gradient Dynamic Reorientation: Preserving Physical Microstructure Alignment
Unpredictable fetal movement introduces an independent rigid transformation \(T_{k,j}(\mathbf{x}') = \mathbf{R}_{k,j}\mathbf{x}' + \mathbf{t}_{k,j}\) for each slice \(S_{k,j}\), mapping the canonical anatomical space \(\mathbf{x}'\) to the scanner coordinate space \(\mathbf{x}\). Because white matter fiber bundles rotate rigidly alongside the fetal head, the nominal diffusion gradient direction \(\mathbf{g}_k\) applied by the scanner becomes misaligned with the underlying local tissue orientation. To maintain physical consistency, each slice is assigned six learnable Lie algebra parameters (axis-angle rotation and translation), and the diffusion gradient vector is dynamically reoriented into the anatomical frame: \(\mathbf{g}'_{k,j} = \mathbf{R}_{k,j}^{-1}\mathbf{g}_k\). This allows the continuous representation to accurately fit intrinsic microstructural anisotropy in a canonical coordinate system across all acquired slices.

2. Continuous Spherical Harmonics Field Representation: Continuous Coordinate-to-Sphere Mapping
High Angular Resolution Diffusion Imaging (HARDI) typically incurs discretization errors and dimensional explosion when stored on discrete volumetric grids. The proposed architecture maps any continuous spatial coordinate \(\mathbf{x}'\) to direction-independent real-valued orthonormal Spherical Harmonics (SH) basis coefficients \(\mathbf{c}(\mathbf{x}')\). Exploiting the antipodal symmetry of water diffusion in fibrous tissue, only even degrees (\(l = 0, 2, \dots, L\)) are preserved, resulting in coefficient dimension \(M = (L+1)(L+2)/2\). An \(\text{MLP}_V\) powered by a multi-resolution hash grid maps coordinate \(\mathbf{x}'\) to coefficients \(c_l^m(\mathbf{x}')\) and latent code \(\mathbf{z}(\mathbf{x}')\): $\(V(\mathbf{x}', \mathbf{g}') = \sum_{l=0,2,\dots}^{L} \sum_{m=-l}^{l} c_l^m(\mathbf{x}') Y_l^m(\mathbf{g}')\)$ This formulation decouples spatial representation from directional sampling: the network evaluates continuous diffusion signals across arbitrary continuous gradient angles via closed-form SH synthesis, providing native compatibility with Constrained Spherical Deconvolution (CSD) for tractography without resampling artifacts.

3. Physics-Based Forward Slice Acquisition: Disentangling Low-Frequency Bias and Finite Resolution Blur
Directly comparing continuous volume predictions with raw slices introduces severe confounding from hardware and motion artifacts. Magnetic field inhomogeneities and magnetization interactions with fetal movement (spin-history effects) induce slowly varying spatial intensity modulations across slices. To isolate these effects, a dedicated bias branch \(\text{MLPB}\) takes only the lowest-resolution hash grid components \(\gamma_{\text{low}}\) and slice embedding \(\mathbf{e}_{k,j}\), creating an architectural bottleneck that restricts high frequencies and models exclusively smooth multiplicative bias fields \(B_{k,j}(\mathbf{x}')\). The intermediate clean signal \(I'_{k,j}(\mathbf{x}', \mathbf{g}'_{k,j}) = B_{k,j}(\mathbf{x}') V(\mathbf{x}', \mathbf{g}'_{k,j})\) is subsequently convolved with a 3D Gaussian point spread function (PSF) \(P_{k,j}\) in scanner space to account for through-plane slice profile blurring and finite in-plane resolution.

4. Through-Plane Continuity Regularization: Anisotropic Gradient Ratio Penalty Against Stripe Artifacts
Because fetal diffusion MRI acquires single-orientation slice stacks per gradient, the reconstruction problem is ill-posed along the through-plane direction \(\mathbf{n}_{k,j} = \mathbf{R}_{k,j}^{-1}(\mathbf{n}_k)\), causing implicit networks to overfit slice-to-slice intensity steps into pronounced striping. The Through-Plane Continuity (TPC) regularizer detects and suppresses these non-anatomical discontinuities. For each voxel, 1D finite-difference gradient magnitudes are computed over continuous random offsets \(\boldsymbol{\delta}_z = u \mathbf{n}_{k,j}\) bounded by slice thickness \(s\), yielding through-plane variation \(G_z\); concurrently, an all-direction 3D neighborhood gradient \(G_{\text{all}}\) is computed across the local PSF kernel: $\(G_z = \frac{1}{|\mathcal{N}_z|} \sum_{\boldsymbol{\delta}_z \in \mathcal{N}_z} \frac{|V(\mathbf{x}' + \boldsymbol{\delta}_z, \mathbf{g}'_{k,j}) - V(\mathbf{x}', \mathbf{g}'_{k,j})|}{\|\boldsymbol{\delta}_z\|}, \quad G_{\text{all}} = \frac{1}{|\mathcal{N}_{\text{PSF}}|} \sum_{\boldsymbol{\delta} \in \mathcal{N}_{\text{PSF}}} \frac{|V(\mathbf{x}' + \boldsymbol{\delta}, \mathbf{g}'_{k,j}) - V(\mathbf{x}', \mathbf{g}'_{k,j})|}{\|\boldsymbol{\delta}\|}\)$ The ratio \(R_{k,j,i} = \frac{G_z}{G_{\text{all}} + \epsilon}\) flags artifacts when directional through-plane variation drastically exceeds local 3D spatial variation, penalizing disproportionate gradients beyond threshold \(\tau\) while preserving sharp, multidirectional anatomical boundaries.

Loss & Training

The entire pipeline is trained end-to-end in a self-supervised regime. The total loss combines Gaussian data fidelity \(\mathcal{L}_D\), total variation regularization \(\mathcal{L}_R\), and through-plane continuity regularization \(\mathcal{L}_T\): $\(\mathcal{L} = \mathcal{L}_D + \mathcal{L}_R + 0.5 \mathcal{L}_T\)$ Data fidelity \(\mathcal{L}_D\) incorporates voxel-wise standard deviation \(\sigma_{k,j}(\mathbf{x}')\) predicted by \(\text{MLP}_U\) and learnable slice-level outlier parameter \(\nu_{k,j}\) to down-weight corrupted samples: $\(\mathcal{L}_D = \frac{1}{N} \sum_{k,j,i} \left( \frac{\left( \hat{S}_{k,j,i} - S_{k,j}(\mathbf{x}_{k,j,i}) \right)^2}{2 \sigma_{k,j}^2(\mathbf{x}'_{k,j,i})} + \log \left( \sigma_{k,j}(\mathbf{x}'_{k,j,i}) + \nu_{k,j} \right) \right)\)$ The through-plane regularization is formulated as: $\(\mathcal{L}_T = \frac{1}{N} \sum_{k,j,i} \text{ReLU}\left( \log(R_{k,j,i}) - \tau \right)\)$ with margin threshold \(\tau = 80\). Implemented in PyTorch with custom CUDA acceleration via tiny-cuda-nn, optimization converges in 6000 iterations on a single NVIDIA Tesla V100 GPU (32GB).

Key Experimental Results

Main Results

The method was evaluated against IRTK (state-of-the-art classical SVR baseline) on three datasets: 9 real in-vivo fetal subjects (3 clinical protocols, 6 dHCP protocols), 20 simulated neonatal subjects with realistic fetal motion and noise, and an isotropic 1.0 mm neonatal template simulation. Metrics include PSNR, SSIM, NRMSE, NCC, and execution time (hours).

Dataset	Metric	Ours	Prev. SOTA (IRTK)	Gain / Change
In-vivo Fetal Data	PSNR ↑ (dB)	23.54 (±1.23)	23.34 (±0.98)	+0.86%
In-vivo Fetal Data	SSIM ↑	0.9014 (±0.0246)	0.7516 (±0.1263)	+19.93% (+0.1498)
In-vivo Fetal Data	NRMSE ↓	0.1788 (±0.0284)	0.1828 (±0.0274)	-2.19%
In-vivo Fetal Data	NCC ↑	0.8194 (±0.0374)	0.8085 (±0.0350)	+1.35%
In-vivo Fetal Data	Run Time ↓ (h)	5.52 (±0.16)	8.77 (±4.93)	-37.06% speedup
Simulated Neonatal Data	PSNR ↑ (dB)	25.20 (±0.72)	23.99 (±0.78)	+5.04% (+1.21 dB)
Simulated Neonatal Data	SSIM ↑	0.9828 (±0.0055)	0.9639 (±0.1000)	+1.96%
Simulated Neonatal Data	NRMSE ↓	0.1325 (±0.0105)	0.1520 (±0.0133)	-12.83%
Simulated Neonatal Data	NCC ↑	0.8412 (±0.0312)	0.7654 (±0.0540)	+9.90%
Simulated Neonatal Data	Run Time ↓ (h)	6.75 (±0.11)	47.80 (±3.74)	-85.88% speedup
Simulated Neonatal Template	PSNR ↑ (dB)	24.88 (±0.90)	24.74 (±1.24)	+0.57%
Simulated Neonatal Template	SSIM ↑	0.9688 (±0.0040)	0.9554 (±0.0048)	+1.40%
Simulated Neonatal Template	NRMSE ↓	0.1379 (±0.0115)	0.1393 (±0.0055)	-1.01%
Simulated Neonatal Template	NCC ↑	0.9504 (±0.0090)	0.9469 (±0.0039)	+0.37%
Simulated Neonatal Template	Run Time ↓ (h)	6.03 (±0.09)	102.15 (±2.86)	94.10% time reduction (17× speedup)

Ablation Study

1. Spherical Harmonics Maximum Order (\(L\)) on Real Fetal Data
The impact of truncating SH orders \(L \in \{2, 4, 6, 8\}\) on real fetal reconstructions:

Config (SH Order)	PSNR ↑ (dB)	SSIM ↑	NRMSE ↓	NCC ↑	Note
\(L = 2\)	23.59 (±1.03)	0.9013 (±0.0291)	0.1779 (±0.0287)	0.8205 (±0.0319)	Smooth approximation, minimal computation
\(L = 4\) (Selected)	23.54 (±1.23)	0.9014 (±0.0246)	0.1788 (±0.0284)	0.8194 (±0.0374)	Preserves fiber crossings with identical fidelity
\(L = 6\)	23.14 (±1.27)	0.8926 (±0.0277)	0.1871 (±0.0290)	0.7989 (±0.0454)	Moderate overfitting to noise and motion artifacts
\(L = 8\)	22.41 (±1.94)	0.8776 (±0.0359)	0.2045 (±0.0421)	0.7572 (±0.0900)	Severe degradation from overfitting residual noise

2. Artifact Mitigation Modules (Bias Field & TPC Regularization) on Fetal Data
Evaluating individual and combined contributions of the bias-field network and Through-Plane Continuity:

Config	Bias Field	TPC Reg.	PSNR ↑ (dB)	SSIM ↑	NRMSE ↓	NCC ↑	Note
Baseline w/o both	✗	✗	22.77 (±1.31)	0.8898 (±0.0293)	0.1954 (±0.0297)	0.7841 (±0.0466)	Severe striping and heavy intensity drifts
Bias field only	✓	✗	23.25 (±1.24)	0.8964 (±0.0273)	0.1848 (±0.0286)	0.8071 (±0.0410)	Mitigates field bias, residual striping remains
TPC reg. only	✗	✓	23.07 (±1.29)	0.8949 (±0.0263)	0.1889 (±0.0301)	0.7972 (±0.0434)	Suppresses stripes, uncorrected slow bias
Full model	✓	✓	23.54 (±1.23)	0.9014 (±0.0246)	0.1788 (±0.0284)	0.8194 (±0.0374)	Optimal quality, stripe-free with sharp anatomy

Key Findings

• On in-vivo fetal data, the proposed method demonstrates a massive +19.93% jump in SSIM over IRTK (0.9014 vs. 0.7516), showing remarkable preservation of fine cortical folding and deep white matter structures along with tightly bounded variance across subjects.
• Truncating spherical harmonics at \(L=4\) represents the optimal operational point: expanding to \(L=8\) degrades PSNR by 1.13 dB due to fitting motion-induced noise, while \(L=4\) provides sufficient angular degrees of freedom for estimating complex crossing fibers via Constrained Spherical Deconvolution (CSD).
• In terms of computational efficiency, neural hash-grid optimization on GPU reduces runtime by up to 94.1% compared to CPU multi-threading (102.15 h down to 6.03 h on neonatal templates), demonstrating dramatic scalability for volumetric diffusion reconstruction.

Highlights & Insights

• Dual Spatial-Angular Continuous Formulation: Coupling coordinate-based multi-resolution hash encoding with even-degree real spherical harmonics maps continuous 3D space to continuous spherical orientations, eliminating discrete voxel grid sampling bottlenecks.
• Physical Dynamic Gradient Reorientation: Rotating nominal scanner gradients into the moving fetal frame simultaneously with 6-DoF pose optimization guarantees strict adherence to diffusion anisotropy physics.
• Selective Anisotropic Ratio Regularization: TPC penalizes directional gradient surges relative to 3D local contexts, successfully eliminating single-orientation slice striping while preserving true cortical and anatomical boundaries.

Limitations & Future Work

• Rigid Motion Assumption: The model assumes rigid transformations per slice; non-rigid fetal deformations (such as neck bending or maternal contraction pressure) violate this assumption and would benefit from local non-rigid deformation modeling.
• Offline Reconstruction Latency: Although 6-hour training represents a dramatic improvement over multi-day classical methods, it remains an offline procedure; integrating 3D Gaussian Splatting primitives or learned flow matching priors could enable near-real-time reconstruction.

Related Work & Insights

• vs IRTK (Deprez et al., TMI 2020): IRTK performs iterative two-stage registration and discrete SH voxel fitting on CPU, taking up to days and suffering from unconstrained slice striping. This work provides continuous implicit representations with 17× speedup and eliminates stripes via TPC regularization.
• vs NeSVoR (Xu et al., TMI 2023): NeSVoR pioneered self-supervised INR for fetal structural MRI but is limited to scalar intensity fields with orthogonal multi-view constraints. This work extends INR to high-dimensional spherical harmonics diffusion fields under ill-posed single-stack acquisitions.

Rating

• Novelty: ⭐⭐⭐⭐⭐ First learning-based framework formulating continuous spherical harmonics fields for motion-corrupted fetal diffusion SVR with tailored TPC regularization.
• Experimental Thoroughness: ⭐⭐⭐⭐⭐ Rigorously tested across in-vivo fetal scans, simulated subjects, and template benchmarks with detailed ablations.
• Writing Quality: ⭐⭐⭐⭐⭐ Elegant mathematical derivations, sound physical grounding, and coherent structural presentation.
• Value: ⭐⭐⭐⭐⭐ Resolves a long-standing roadblock in prenatal neuroimaging, unlocking reliable fetal white matter connectome analysis with open-source code.