Spherical Polar Fourier Imaging

Theory

Spherical Polar Fourier Imaging (SPFI) is a reconstruction method to estimate continuous diffusion signal, Ensemble Average Propagator (EAP), diffusion Orientation Distribution Function (dODF), some meaningful quantities (RTO, MSD, PFA) from samples of DWI signals.

SPFI represents the signal using 3D SPF basis [Assemlal2009]:

E(q\uu|\zeta) = \sum_{n=0}^N\sum_{l=0}^L\sum_{m=-l}^l a_{nlm} G_n(q|\zeta) Y_l^m(\uu)

G_n(q|\zeta)=\left[ \frac{2n!}{\zeta ^{3/2}\Gamma(n+3/2)}\right ]^{1/2}\exp\left(-\frac{q^2}{2\zeta}\right) L_n^{1/2}(\frac{q^2}{\zeta})

where \q=q\uu, \uu\in \mathbb{S}^2, \zeta is the scale parameter and Y_l^m(\uu) is the real spherical harmonic basis.

The reconstruction in SPFI has two steps:

  1. Estimate the coefficients of SPF basis from measurements of diffusion signals using compressed sensing, dictionary learning, etc.
  2. Analytically reconstruct EAP, dODF, scalar maps from the SPF coefficients [Cheng2010a] [Cheng2010b].

SPFI can be seen as a generalization of Q-Ball Imaging [Tuch2004]. Compared to QBI which works only for single shell data, SPFI works for arbitrarily sampled data. Compared to Diffusion spectrum imaging (DSI), SPFI requires relatively less number of samples and lower b values. SPFI uses analytical transforms, avoiding numerical Fourier transform and numerical integral in DSI.

A typical sampling scheme for SPFI has 2 or 3 b-values, maximal b value more than 3000 s/mm^2, and each shell has more than 30 samples. For example, b values are 1500 and 3000, 40 samples per shell. The samples in different shell are better to be staggered. See the tutorial on uniform sampling scheme.

Synthetic data Experiment

DWI data generation

Generate DWI data in 3 shells. See the tutorial on DWI data simulation.

export DMRITOOL_EXAMPLE_DIR=${DMRITOOL_SOURCE_DIR}/Examples

b=1000,2000,3000
DWISimulator ${DMRITOOL_EXAMPLE_DIR}/dwi_circle_crossing.txt --outdwi dwi.nii.gz --outodf odfTrue.nii.gz --outeap eapTrue_r0.015.nii.gz --qorientations ${DMRITOOL_EXAMPLE_DIR}/Elec060.txt --bvalues ${b} --rorientations ${DMRITOOL_EXAMPLE_DIR}/directions_t4.txt --rvalues 0.015 --noisesigma 0.0 --outb0 dwi_diagonal_b0.nii.gz --outputdwitype EACHSHELL

DL-SPFI reconstruction

Reconstruction of SPF coefficients using DL-SPFI (with the default scale).

SphericalPolarFourierImaging dwi.txt --sh 8 --ra 4 --signal signalSPF.nii.gz --radius 0.015 --estimation L1_DL --lambdaL1 1e-7
  • The above commend is to perform DL-SPFI using SH rank 8 and radial rank 4, regularization parameter lambda 1e-7.
  • It uses default scale \zeta = 1.0 / (8 \pi^2 \tau D_0), default mean diffusivity D_0=0.7\times 10^{-3} mm^2/s for all voxels. You can set default mean diffusivity in --md0
  • You can try different regularization lamdaL1 around 1e-6 in range [1e-8, 1e-5]
  • For real data or general synthetic data with unknown mean diffusivity, It is better to estimate mean diffusivity and adaptively set scale for each voxel.

We use a generalized high order tensor (GHOT) model to estimate the mean diffusivity [ChengThesis2012]. Then use DL-SPFI with adaptive scale for adaptive dictionary [Cheng2013].

MeanDiffusivityEstimator dwi.txt D_sh4_ra1.nii.gz --sh 4 --ra 1
SphericalPolarFourierImaging dwi.txt --sh 8 --ra 4 --signal signalSPF.nii.gz --radius 0.015 --estimation L1_DL --lambdaL1 1e-7  --mdImage D_sh4_ra1.nii.gz

Note

You cannot use different ranks from SH rank 8 and radial rank 4 for DL-SPFI, because the used dictionary was learned using rank (8,4).

L1-SPFI reconstruction

Without using learned dictionary, you can try L1-SPFI which uses least squares with L1 norm regularization based on compressed sensing [Cheng2011].

SphericalPolarFourierImaging dwi.txt --sh 8 --ra 4 --lambdaSH 1e-9 --lambdaRA 1e-9 --signal signalSPF.nii.gz --radius 0.015 --estimation L1_2 --mdImage D_sh4_ra1.nii.gz
  • Without --mdImage D_sh4_ra1.nii.gz, it uses default scale for all voxels.
  • You can try different ranks, and different regularizations lambdaSH and lambdaRA around 1e-8.

L2-SPFI reconstruction

Besides DL-SPFI and L1-SPFI, you can also try L2-SPFI which uses least squares with L2 norm regularization [Cheng2010a].

SphericalPolarFourierImaging dwi.txt --sh 6 --ra 2 --lambdaSH 1e-9 --lambdaRA 1e-9 --signal signalSPF.nii.gz --radius 0.015 --estimation LS --mdImage D_sh4_ra1.nii.gz
  • Without --mdImage D_sh4_ra1.nii.gz, it uses the default scale for all voxels based on the default mean diffusivity.
  • You can try different ranks, and different regularizations lambdaSH and lambdaRA around 1e-8.
  • You may need to use lower ranks in L2-SPFI than L1-SPFI and DL-SPFI.

Note

We encourage you to use DL-SPFI with adaptive scales other than L2-SPFI or L1-SPFI.

Analytical reconstruction of EAP and ODF

SPFI has analytical relationship between SPF coefficients and EAP profiles, ODFs. See [Cheng2010a] and [Cheng2010b]. You can use different methods (L2-SPFI, L1-SPFI or DL-SPFI) to reconstruct SPF coefficients, then you always can efficiently obtain the EAP profiles and ODFs from the SPF coefficients. But please set the scale in SPF basis correctly based the mean diffusivity used in SPF coefficients reconstruction, also set the ranks correctly.

  • Analytically obtain EAP profiles and ODFs (using default scale based on default mean diffusivity).
SPFToProfile signalSPF.nii.gz eap_r0.015.nii.gz --sh 8 --ra 4 --radius 0.015 --fourier
SPFToODF signalSPF.nii.gz odf.nii.gz --sh 8 --ra 4
  • Analytically obtain EAP profiles and ODFs (using adaptive scale based on adaptive mean diffusivity).
SPFToProfile signalSPF.nii.gz eap_r0.015.nii.gz --sh 8 --ra 4 --radius 0.015 --fourier --mdImage D_sh4_ra1.nii.gz
SPFToODF signalSPF.nii.gz odf.nii.gz --sh 8 --ra 4 --mdImage D_sh4_ra1.nii.gz
  • Visualization of EAP profiles and ODFs.
MeshFromSHCoefficients eap_r0.015.nii.gz eap_r0.015_vis.vtk --tessorder 4 --scale 8e-6
VTKPolyData.py --vtk eap_r0.015_vis.vtk --png synthetic_eap_r0.015.png --zoom 1.3
MeshFromSHCoefficients odf.nii.gz odf_vis.vtk --tessorder 4 --scale 1.5
VTKPolyData.py --vtk odf_vis.vtk --png synthetic_odf.png --zoom 1.3
synthetic_eap_r0.015.png

eap_r0.015

synthetic_odf.png

odf

Analytical reconstruction of scalar maps

We can analytically generate some scalar maps using the estimated SPF coefficients. See [Wu2007], [Cheng2010a], [ChengThesis2012].

  • Return-To-Origin probability (RTO) is the EAP with r=0, i.e. P(0).
SPFToScalarMap signalSPF.nii.gz rto.nii.gz --mapType RTO  --sh 8 --ra 4 --mdImage D_sh4_ra1.nii.gz
  • Mean Squared Displacement (MSD) is the variance of the probability, i.e. \int_{\mathbb{R}^3} P(R) R^TR d R.
SPFToScalarMap signalSPF.nii.gz msd.nii.gz --mapType MSD  --sh 8 --ra 4 --mdImage D_sh4_ra1.nii.gz
  • Propagator Fractional Anisotropy (PFA) is a generalization of FA for non-Gaussian EAP. It is defined as the normalized L2 distance between the EAP to its nearest EAP. See [ChengThesis2012].
SPFToScalarMap signalSPF.nii.gz pfa.nii.gz --mapType PFA  --sh 8 --ra 4 --mdImage D_sh4_ra1.nii.gz
  • From the above estimated EAP profiles and ODFs, you can also generate generalized FA (GFA). See [Tuch2004].
SHCoefficientsToGFA eap_r0.015.nii.gz eap_r0.015_gfa.nii.gz
SHCoefficientsToGFA odf.nii.gz odf_gfa.nii.gz
  • We can visualize eap profile using scalar maps as background. --image-range is to control the contrast in visualization. If it is not given, it maps the minimal value to black and the maximal value to white.
VTKPolyData.py --vtk eap_r0.015_vis.vtk --image eap_r0.015_gfa.nii.gz --image-range 0,1 --png synthetic_eap_r0.015_withgfa.png --zoom 1.3
synthetic_eap_r0.015_withgfa.png

eap_r0.015

Real data Experiment for Human Connectome Project

HCP data

We downloaded a preprocessed subject data (ID: 100307) from HCP data. The data has three shells with b values of 1000, 2000, 3000 s/mm^2, 90 samples per shell. It also contains 18 volumes with b=0, thus totally it has 288 volumes. Considering the data size is big (145, 174, 145, 288), we cropped only one slice for this tutorial. You can download it from this link.

Preprocess for DMRITool

The data has already preprocessed. Here we just use DWIPreprocess to normalize the DWI data using the baseline image.

echo b_raw.txt grad_raw.txt dwi_c88.nii.gz > data_c88_raw.txt
DWIPreprocess data_c88_raw.txt  data_c88_normalize.txt --oEachShell --bThreshold 15 --mask mask_c88.nii.gz  --odwi dwi_c88_normalize.nii.gz  --ograd grad_normalize.txt --ob0Image dwi_c88_b0.nii.gz
  • The b values in the HCP data all larger than 0. There is no exact b=0. Thus we consider b values smaller than the bThreshold as b=0.
  • The b values are not exactly equal 1000, 2000, or 3000. DWIPreprocess can group b values and replace b values using the mean of the group.
  • It outputs data_c88_normalize.txt which contains the b values, gradients and the normalized DWI data.

DL-SPFI Reconstruction of SPF coefficients

Estimate the mean diffusivity first.

MeanDiffusivityEstimator data_c88_normalize.txt D_sh4_ra1.nii.gz --sh 4 --ra 1

Perform DL-SPFI to estimate the SPF coefficients.

SphericalPolarFourierImaging data_c88_normalize.txt --estimation L1_DL --sh 8 --ra 4 --lambdaL1 1e-6  --signal signalSPF.nii.gz --solver SPAMS --mask mask_c88.nii.gz --mdImage D_sh4_ra1.nii.gz

Reconstruction of EAPs

Obtain EAP profiles from the SPF coefficients, and obtain GFA of EAP profiles.

SPFToProfile signalSPF.nii.gz eap_r0.015.nii.gz  --radius 0.015 --ra 4 --sh 8 --mdImage D_sh4_ra1.nii.gz --fourier
SHCoefficientsToGFA eap_r0.015.nii.gz eap_r0.015_gfa.nii.gz

Generate a coarse mesh (--tessorder 3) for visualization the EAP profiles in a whole slice.

MeshFromSHCoefficients eap_r0.015.nii.gz eap_r0.015_visall.vtk --tessorder 3 --scale 1e-5

Generate a fine mesh (--tessorder 4) for visualization the EAP profiles in a ROI.

MeshFromSHCoefficients eap_r0.015.nii.gz eap_r0.015_vis.vtk --tessorder 4 --scale 1e-5 --box 80,100,0,0,70,90

You can use VTKPolyData.py or paraview to visualize the vtk files.

VTKPolyData.py --vtk eap_r0.015_visall.vtk
VTKPolyData.py --vtk eap_r0.015_vis.vtk

You can also put scalar map like GFA map in the background. See the following results obtained by paraview.

HCP_eap_r0.015_whole HCP_eap_r0.015_ROI
eap_r0.015_visall.vtk eap_r0.015_vis.vtk

Reconstruction of scalar maps

Analytical reconstruction of RTO, MSD and PFA maps.

SPFToScalarMap signalSPF.nii.gz rto.nii.gz --mapType RTO  --sh 8 --ra 4 --mdImage D_sh4_ra1.nii.gz
SPFToScalarMap signalSPF.nii.gz msd.nii.gz --mapType MSD  --sh 8 --ra 4 --mdImage D_sh4_ra1.nii.gz
SPFToScalarMap signalSPF.nii.gz pfa.nii.gz --mapType PFA  --sh 8 --ra 4 --mdImage D_sh4_ra1.nii.gz

You can use many tools to visualize the 3D volume, e.g. VTKPolyData.py. Here are pictures for the scalar maps I got by fslview in FSL.

HCP_RTO HCP_MSD HCP_PFA
RTO map MSD map PFA map
[Assemlal2009]Haz-Edine Assemlal, David Tschumperle, Luc Brun, Efficient and robust computation of PDF features from diffusion MR signal, Medical Image Analysis, vol 13, p. 715-729, 2009
[Tuch2004](1, 2) David S. Tuch, Q-Ball Imaging, Magnetic Resonance in Medicine, vol 52, p. 1358-1372, 2004.
[Wu2007]Yu-Chien Wu, Andrew L. Alexanderb, Hybrid diffusion imaging, NeuroImage, vol 36, p. 617-629, 2007.
[Cheng2010a](1, 2, 3, 4) Jian Cheng, Aurobrata Ghosh, Rachid Deriche, Tianzi Jiang, Model-free and Analytical EAP Reconstruction via Spherical Polar Fourier Diffusion MRI, 13th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI‘10), Beijing, September 20-24, 2010
[Cheng2010b](1, 2) Jian Cheng, Aurobrata Ghosh, Rachid Deriche, Tianzi Jiang, Model-free, regularized, fast, and robust analytical orientation distribution function estimation, 13th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI‘10), Beijing, September 20-24, 2010
[Cheng2011]Jian Cheng, Sylvain Merlet, Emmanuel Caruyer, Aurobrata Ghosh, Rachid Deriche, Tianzi Jiang, Compressive Sensing Ensemble Average Propagator Estimation via L1 Spherical Polar Fourier Imaging, MICCAI Workshop on Computational Diffusion MRI (CDMRI‘11), Toronto, Canada, September, 2011
[ChengThesis2012](1, 2, 3) Jian Cheng, Estimation and Processing of Ensemble Average Propagator and Its Features in Diffusion MRI, Universite Nice Sophia Antipolis, 2012
[Cheng2013]Jian Cheng, Tianzi Jiang, Rachid Deriche, Dinggang Shen, Pew-Thian Yap, Regularized Spherical Polar Fourier Diffusion MRI with Optimal Dictionary Learning, 16th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI‘13), Nagoya, September 22-26, 2013