Review of eddy’s slice-to-volume correction (introduced v. 5.0.11)
date written: 06/02/2019
Overview
This analysis compares (i) our existing eddy correction, which corrects for inter-volume subject movement [Andersson, 2016], to (ii) eddy correction with the slice-to-volume correction [Andersson, 2017], which also corrects for inter- and intra-volume subject movement. The ability to correct for intra-volume subject movement was introduced in eddy v.5.0.11. This correction may be useful on our data, as the assumption that all participants remain still during the time is takes to acquire a DWI volume - in the case of SPINS, ~8.8 seconds – is likely not met. Intra-volume movement can be corrected as it does not cause signal loss: rather, signal merely needs to be relocated.
This comparison was run on SPINS participant SPN01_MRC_0042, a 48 year old male with SSD (age of onset = early 20s, number of hospitalizations = approximately 20). This participant was uniquely noted in Dashboard QC to have a pronounced zebra-pattern artifact, in which interleaved slices no longer stack up to a “valid” volume, likely the result of intra-volume movement.
Before comparing the two eddy methods, the data was preprocessed identically, using tools from FSL and MRtrix3, and in keeping with Nat’s recommendations (note, however, that no fieldmap/top-up/synthetic b0 correction was used):
#denoise
dwidenoise dwi.nii.gz dwi_denoise.nii.gz
#gibbs unringing
mrdegibbs dwi_denoise.nii.gz dwi_denoise_unring.nii.gz
#take first volume as B0
fslroi dwi_denoise_unring.nii.gz b0 0 1
#perform BET
bet b0 nodif_brain -m -R -f 0.3
Methods
Our regular eddy method does not include slice-to-volume correction. To preform eddy as usual, I ran the following command command in FSL/5.0.11:
eddy_cuda --imain=dwi_denoise.nii.gz \
--mask=nodif_brain_mask --acqp=acqparams.txt \
--index=index.txt --bvecs=bvec.txt --bvals=bval.txt \
--out=eddy_orig --repol --residuals --cnr_maps
The --repol parameter was introduced in v. 5.0.10. It detects slice dropout that results from motion coinciding in time with the diffusion encoding part of the sequence (visualized as a diagonal band artifact), and replaces them with Gaussian Process predictions. See the eddy UsersGuide for description of the other parameters.
The regular eddy correction looks better than the raw data, but the zebra artifact is still seen [again, slice 58, min: 0, max: 300]:
Slice-to-volume correction requires that eddy be provided with 5 new parameters. The first, --slspec, is a text file that specifies slice order. To obtain slice order, I ran the newest version of dcm2niix, which extracts slice time information from dicom headers and prints it to the ‘JSON sidecar’ file, with the ‘SliceTiming’ tag. The slice timing information for participant SPN01_MRC_0042 is as follows:
## [1] 4.4175 0.0000 4.5500 0.1350 4.6850 0.2675 4.8200 0.4025 4.9525 0.5350
## [11] 5.0875 0.6700 5.2200 0.8025 5.3550 0.9375 5.4875 1.0700 5.6225 1.2050
## [21] 5.7550 1.3375 5.8900 1.4725 6.0225 1.6075 6.1575 1.7400 6.2925 1.8750
## [31] 6.4250 2.0075 6.5600 2.1425 6.6925 2.2750 6.8275 2.4100 6.9600 2.5425
## [41] 7.0950 2.6775 7.2275 2.8100 7.3625 2.9450 7.4950 3.0800 7.6300 3.2125
## [51] 7.7625 3.3475 7.8975 3.4800 8.0325 3.6150 8.1650 3.7475 8.3000 3.8825
## [61] 8.4325 4.0150 8.5675 4.1500 8.7000 4.2825Then, I used MATLAB code provided by FSL to deduce slice acquisition order and create the --slspec text file. Though it appears slice timing varies minutely between participants, slice order appears constant, as follows:
## [1] 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45
## [24] 47 49 51 53 55 57 59 61 63 65 0 2 4 6 8 10 12 14 16 18 20 22 24
## [47] 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64Then, I built on the regular eddy command by additionally providing the slice order text file via --slspec, and used the defaults for 3 of the 4 other required arguments. Specifically, --s2v_niter is the number of iterations for estimating slice-to-volume movement [DEFAULT: 5]; --s2v_lambda is the strength of temporal regularisation of slice-to-volume movement [DEFAULT: 1], and --s2v_interp is the interpolation model for estimating slice-to-volume movement [DEFAULT: trilinear].
The --mporder parameter provides the temporal order of movement [DEFAULT: 0]. When performing slice-to-volume correction, however, this value must be set to an integer value greater than 0 but less than the number of excitations in a volume. The larger the value of --mporder, the more degrees of freedom for modelling movement. If --mporder is set to N-1, where N is the number of excitations in a volume, the location of each slice is individually estimated. However, the FSL group recommends against going that high, and in their tests, used values of N/4 – N/2.
An FSL mailing list post indicates that, for single-band acquisitions, the number of excitations (per volume) is the same as the number of slices. Our DWI data has 66 slices. I choose an --mporder value of 66/4 = 16.5 = ~16.
Thus, the command read:
eddy_cuda --imain=dwi_denoise.nii.gz \
--mask=nodif_brain_mask --acqp=acqparams.txt \
--index=index.txt --bvecs=bvec.txt --bvals=bval.txt \
--out=eddy_slspec16 --repol --residuals --cnr_maps \
--slspec=my_slspec.txt \
--mporder=16 --s2v_niter=5 --s2v_lambda=1 --s2v_interp=trilinear
The new eddy slice-to-volume correction looks better than before. The zebra artifact is largely corrected [again, slice 58, min: 0, max: 300]:
The difference is especially evident in a gif:
QC metrics
To quantitatively examine data quality, we ran dtifit on the outputs from the regular eddy method, as well as slice-to-motion corrected outputs, with the following command:
dtifit --data=eddy_slspec16.nii.gz \
--mask=nodif_brain_mask --bvecs=bvec.txt --bvals=bval.txt \
--out=dtifit_slspec16 --sse
We see that the tensor residuals (sse files) look better after slice-to-volume correction [x=64, y=62, z=32, min:0, max:1]:
Recent work in our DWI group has suggested that CSD residuals are also especially informing for some artifacts. More information about CSD residuals can be found on the MRtrix documentation. I ran the following code:
dwidenoise eddy_slspec16.nii.gz \
eddy_slspec16_CSD.nii.gz \
-noise eddy_slspec16_CSD_res.nii.gz
fslmaths eddy_slspec16_CSD_res.nii.gz \
-mas nodif_brain_mask.nii.gz \
eddy_slspec16_CSD_res_masked.nii.gz
However, we do not see much difference in CSD residuals after slice-to-volume correction [x=63, y=63, z=32, min:0, max:3]:
To quantitatively examine data quality, we ran eddy quad on the outputs from the regular eddy method, as well as slice-to-motion corrected outputs, with the following command:
eddy_quad eddy_slspec16 \
-idx index.txt -par acqparams.txt \
-m nodif_brain_mask.nii.gz -b bval.txt -s my_slspec.txt
Below is a summary of values of interest from the generated qc.pdf:
qc.pdf metric
|
regular eddy |
eddy with --slspec
|
|---|---|---|
| volume-to-volume motion | ||
| absolute motion (mm) | 2.09 | 1.88 |
| relative motion (mm) | 0.68 | 0.68 |
| x translation (mm) | -0.14 | -0.21 |
| y translation (mm) | -0.50 | -0.18 |
| z translation (mm) | -1.50 | -1.39 |
| x rotation (deg) | -0.92 | -0.67 |
| y rotation (deg) | 0.09 | 0.11 |
| z rotation (deg) | -0.08 | -0.02 |
| within volume motion | ||
| std. x translation (mm) | – | 0.18 |
| std. y translation (mm) | – | 0.41 |
| std. z translation (mm) | – | 0.27 |
| std. x rotation (deg) | – | 0.49 |
| std. y rotation (deg) | – | 0.14 |
| std. z rotation (deg) | – | 0.18 |
| outliers | ||
| total outliers (%) | 1.84 | 1.31 |
| outliers (b=1000 s/mm^2) | 1.84 | 1.31 |
| outliers (PE dir=[ 0. -1. 0.]) | 1.68 | 1.19 |
| SNR/CNR | ||
| average SNR (b=0 s/mm^2) | 12.29 | 16.34 |
| average CNR (b=1000 s/mm^2) | 1.35 | 1.50 |
We also see vastly different “mean slice difference” plots:
We also decided to compare the ratings of our original 5 raters using different methods, to see if their rating would change post --slspec correction. Results were as follows (note: all raters provide ratings apart from other QC, so ratings may not be reliable, and * indicates John’s algorithm results are likely not valid, as no MRC scans were included in his SPINS sample):
| method | rater | rating_eddy_old | rating_eddy_new |
|---|---|---|---|
| volume-by-volume | Hajer | INDETERMINATE/FAIL | PASS |
| tensor residuals (dtifit) | Grace | INDETERMINATE | INDETERMINATE |
| CSD residuals (MRTrix) | Nat | PASS | PASS |
| mean outliers (eddy) | Navona | FAIL | INDETERMINATE |
| algorithm | John (algorithm) | INDETERMINATE* | INDETERMINATE* |
Comparative analysis
Next, I ran ENIGMA TBSS code to register and skeletonize FA images, as described here, and extract FA values per JHU atlas ROI, described here, for both our original eddy and eddy with --slspec preprocessing.
We seen an increase in FA values (without exception) after the --slspec correction.
Below, the same data is sorted by the difference in FA value between correction with --slspec and our regular eddy correction.
Lastly, we arranged the JHU atlas regions from posterior to anterior, i.e., by the value of the Y co-ordinate of the centre of gravity (note, we only fewer regions here, as coordinates not available for all regions):
| ROI tag | ROI name |
|---|---|
| ACR-L | Anterior corona radiata left |
| ACR-R | Anterior corona radiata right |
| ALIC-L | Anterior limb of internal capsule left |
| ALIC-R | Anterior limb of internal capsule right |
| BCC | Body of corpus callosum |
| CGC-L | Cingulum (cingulate gyrus) left |
| CGC-R | Cingulum (cingulate gyrus) right |
| CGH-L | Cingulum (hippocampus) left |
| CGH-R | Cingulum (hippocampus) right |
| CP-L | Cerebral peduncle left |
| CP-R | Cerebral peduncle right |
| CST-L | Corticospinal tract left |
| CST-R | Corticospinal tract right |
| EC-L | External capsule left |
| EC-R | External capsule right |
| FX | Fornix (column and body of fornix) |
| FX/ST-L | Fornix (cres) / Stria terminalis (can not be resolved with current resolution) left |
| FX/ST-R | Fornix (cres) / Stria terminalis (can not be resolved with current resolution) right |
| GCC | Genu of corpus callosum |
| ICP-L | Inferior cerebellar peduncle left |
| ICP-R | Inferior cerebellar peduncle right |
| IFO-L | Inferior fronto-occipital fasciculus left |
| IFO-R | Inferior fronto-occipital fasciculus right |
| ML-L | Medial lemniscus left |
| ML-R | Medial lemniscus right |
| PCR-L | Posterior corona radiata left |
| PCR-R | Posterior corona radiata right |
| PLIC-L | Posterior limb of internal capsule left |
| PLIC-R | Posterior limb of internal capsule right |
| PTR-L | Posterior thalamic radiation (include optic radiation) left |
| PTR-R | Posterior thalamic radiation (include optic radiation) right |
| RLIC-L | Retrolenticular part of internal capsule left |
| RLIC-R | Retrolenticular part of internal capsule right |
| SCC | Splenium of corpus callosum |
| SCP-L | Superior cerebellar peduncle left |
| SCP-R | Superior cerebellar peduncle right |
| SCR-L | Superior corona radiata left |
| SCR-R | Superior corona radiata right |
| SFO-L | Superior fronto-occipital fasciculus (could be a part of anterior internal capsule) left |
| SFO-R | Superior fronto-occipital fasciculus (could be a part of anterior internal capsule) right |
| SLF-L | Superior longitudinal fasciculus left |
| SLF-R | Superior longitudinal fasciculus right |
| SS-L | Sagittal stratum (include inferior longitidinal fasciculus and inferior fronto-occipital fasciculus) left |
| SS-R | Sagittal stratum (include inferior longitidinal fasciculus and inferior fronto-occipital fasciculus) right |
| UNC-L | Uncinate fasciculus left |
| UNC-R | Uncinate fasciculus right |
- Slice-to-volume correction must be run on cuda (cf. openmp), as it is not (yet) implemented in the CPU version of eddy.
- All data can be found in
/projects/ncalarco/Kimel/projects/DWI_harmonization/eddyTest_sliceTiming/data/SPN01_MRC_0042_01 - eddy with
--slspectook approximately 1.5 hours for one participant on local system