Spherical Polar Fourier Imaging¶

Table of Contents

Theory
Synthetic data Experiment
Real data Experiment for Human Connectome Project

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:

Estimate the coefficients of SPF basis from measurements of diffusion signals using compressed sensing, dictionary learning, etc.
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 -o 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 -o odf_vis.vtk --tessorder 4 --scale 1.5
VTKPolyData.py --vtk odf_vis.vtk --png synthetic_odf.png --zoom 1.3

eap_r0.015

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

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 -o 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 -o 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.


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.


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