-
Notifications
You must be signed in to change notification settings - Fork 141
DWI Initialization
This page describes diffusion weighted imaging (DWI) preprocessing that is carried out by the Matlab function dtiInit.
The dtiInit begins with the diffusion data in a nifti file, a T1-anatomical in a nifti file, and the bvecs and bvals files. These are typically provided in the NIMS or Flywheel data system. We expect the T1-anatomical is assumed to be AC-PC aligned. For more information about how to do that see [ACPC alignment](ACPC alignment).
A wide range of parameters are set in dtiInit. The collection of parameters, their description, and their defaults are set in '''dtiInitParams'''.
Depending on how you set these parameters and the size of your data, the dtiInit can take anywhere from 30 minutes to 3 hours. Eddy current correction is a particularly time-consuming step.
dtiInit produces a series of files that are used in other software: e.g. mrtrix, AFQ, LiFE, and mrDiffusion. A critical file is the dt6.mat file. We describe the full set below.
Using NIMS, you can download an entire Session, which will appear as a zipped tar file. When you unzip the data the session unpacks into a directory tree. One of the folders contains the DWI data (nii.gz format) and the relevant bvec and bval file.
Using Flywheel, the individual files can be downloaded.
We are using the Flywheel Gears (based on docker) to upload dicom files from other vendors and convert them to NIfTI files. The core technology for the conversion is dcm2niix.
Please arrange your directory tree like this
<subdir>
_____|___________
| |
t1.nii.gz <ROIS>
dwi.nii.gz
dwi.bval
dwi.bvec
The ROIS directory is optional. Some people make a sub-directory () and put the dwi files in there.
Specific parameters can be set using the function dtiInitParams. A call to this function will return the default parameters (as shown below). For information about the meaning of these parameters and how to set them, go to the DWI Init Parameters page.
dwParams = dtiInitParams
dwParams =
bvalue: []
gradDirsCode: []
clobber: 0
dt6BaseName: ''
flipLrApFlag: 0
numBootStrapSamples: 500
fitMethod: 'ls'
nStep: 50
eddyCorrect: 1
excludeVols: []
bsplineInterpFlag: 0
phaseEncodeDir: []
dwOutMm: [2 2 2]
rotateBvecsWithRx: 0
rotateBvecsWithCanXform: 0
bvecsFile: ''
bvalsFile: ''
noiseCalcMethod: 'b0'
outDir: ''
To run dtiInit without any parameters, accepting the defaults, or you can simply specific the raw diffusion data and the T1 anatomical
[dt6FileName, outBaseDir] = dtiInit('rawDWI.nii.gz','t1.nii.gz')
This is equivalent to
dwParams = dtiInitParams;
[dt6FileName, outBaseDir] = dtiInit('raw/dwi.nii.gz','t1.nii.gz', dwParams)
You can set parameters using param/value pairs with dtiInitParams, as in this example
dwParams = dtiInitParams('clobber',1,'phaseEncodeDir',2,'eddyCorrect',-1);
[dt6FileName, outBaseDir] = dtiInit('rawDWI.nii.gz','t1.nii.gz', dwParams)
Or just set dwParams by hand. E.g.,
dwParams = dtiInitParams;
dwParams.clobber = 1;
dwParams.phaseEncodeDir = 2;
dwParams.eddyCorrect = -1;
You should have a directory structure that looks something like the following.
<subdir>
_____________________|_____________________________
| | |
t1.nii.gz raw <Directory: dtiNumDirsInterpMethod>
dti_aligned.nii.gz | ____________|_____________
dti_aligned.bvecs dwi.nii.gz | | |
dti_aligned.bvals dwi.bval bin dt6.mat t1pdd.png
dwi.bvec |
tensors.nii.gz
b0.nii.gz
brainMask.nii.gz
vectorRGB.nii.gz
wmMask.nii.gz
wmProb.nii.gz
faStd.nii.gz
mdStd.nii.gz
pddDispersion.nii.gz
This DWI Files page describes the files themselves.
LMP to update with dcm2niix comments. Gears.
If you cannot use 'niftify', and neither of the sections above worked for you, here are many other ways to convert DICOM files to NIFTI. For DTI data, it's especially important that your DICOM to NIFTI converter set the NIFTI qform correctly using the volume orientation information in the DICOM header, since this rigid-body transform is used to reorient your gradient directions from the scanner coordinate space to image space. Currently, we have a robust but slow Matlab script ('''niftiFromDicom''') that works reliably for our GE data. But for Siemens or Philips data, Chris Rorden's [http://www.sph.sc.edu/comd/rorden/mricron/dcm2nii.html dcm2nii] works much better and generates proper bvecs/bvals files for you.)
niftiFromDicom(dicomDir, [outDir=dicomDir/..], [studyId=''], [sortByFilenameFlag=false], [makeAxialFlag=true])
If you acquired your dti in separate series, put both directories inside another directory and give that directory as input to niftiFromDicom. It will recursively search the directory tree and try to use all the files that it finds, correctly ordering them based on dicom header info (filenames are ignored).
Alternatively, you can generate 2 niftis and then merge them.
ni1 = readFileNifti('first.nii.gz');
ni2 = readFileNifti('second.nii.gz');
ni1.data = cat(4, ni1.data, ni2.data);
ni1.fname = 'merged.nii.gz';
writeFileNifti(ni1);
Alternatives to our code include [http://lcni.uoregon.edu/~jolinda/MRIConvert MRIConvert] from the University of Oregon, [http://www.cbi.nyu.edu/public/software/dinifti/ dinifti] from the NYU CBI, and [http://www.sph.sc.edu/comd/rorden/mricron/dcm2nii.html dcm2nii] from Chris Rorden at U. South Carolina. All of these have a simple command-line interface that is easy to script and they seem to handle our dti data correctly for the most part (but in some cases come out flipped or with the wrong xForm, so we use niftiFromDicom as our fallback). dinifti doesn't correctly guess the number of slices for DTI data, so you specify that manually. E.g., for a 50-slice dataset:
$ dinifti -s 50 -g dicomFileDir/ dti_raw
Unlike MRIConvert, dinifti won't build your bvals and reoriented bvecs file for you. This is not a problem with data from Stanford, which will be handled correctly by dtiRawPreprocess.m. The rest of this section concerns data acquired elsewhere.
If you have Siemens data, then dcm2nii should be able to generate the bvecs and bvals for you. Or, you can build these using our dtiRawBuildBvecs.m function (in the 'preprocess' directory of the mrDiffusion package). To use this function, you need to know your b-value. If you don't know your b-value, it should be available in the dicom header. For our GE data using the Hedehus/Bammer DTI sequence, it's in the (0019, 10B0) DICOM slot. You will also need to know the gradient directions. dtiReorientBvecs.m assumes that these are stored in the Hedehus/Bammer dwepi.NN.grads format, where NN is an integer code specifying the particular gradient scheme file. The integer code for the dwepi grads file is also in the DICOM header, in (0019, 10B2).
To get the bval, grads or Phase Encode Direction file code, run this in Matlab:
'''info = dicominfo('/path/to/your/dicomFile');''' bval = info.Private_0019_10b0; gradsFileCode = info.Private_0019_10b2; phaseEncode = info.InPlanePhaseEncodingDirection; %'COL' = A/P, 'ROW' = L/R te = info.EchoTime; cv24 = info.Private_0019_10df; % should correspond to the UserCV called "duration of diffusion gradient" in units of microseconds
The gradient dir files that we commonly use are copied on our system under /usr/local/dti/diffusion_grads. If your file isn't there, you'll have to get it from the scanner (ask Roland if you don't know where these are stored). An example from our old 6 direction workhorse 'dwepi.13.grads':
0.0 0.0 0.0 1.0 1.0 0.0 -1.0 -1.0 0.0 0.0 1.0 1.0 0.0 -1.0 -1.0 1.0 0.0 1.0 -1.0 0.0 -1.0 -1.0 1.0 0.0 1.0 -1.0 0.0 0.0 -1.0 1.0 0.0 1.0 -1.0 1.0 0.0 -1.0 -1.0 0.0 1.0
After you run dinifiti and dtiRawBuildBvecs, you should have three files that can be used as inputs to the pipeline for doing motion and eddy-current correction and then estimating tensors (below). These three files are: one big NIFTI file with all the raw diffusion data and two little text files- bvals and bvecs. These two little files specify the diffusion weighting params for each volume in the NIFTI file and are in a simple text format. Note that the gradient directions in the bvecs file might need to be adjusted for the orientation of the imaged volume with respect to the scanner coordinate frame. This is because some DTI sequences specify the directions in the physical scanner reference frame. Thus, if the scan prescription is oblique, the image reference frame is rotated with respect to the scanner reference frame, and that rotation must be applied to the bvecs so that they are in the image reference frame. If this hasn't already been done for you by your DICOM-to-NIFTI converter, you can specify the appropriate rotation in the 'xform' param to dtiRawBuildBvecs. For our sequence at Stanford, the gradient directions are specified in the image reference frame (sometimes called 'logical' frame), so we don't need to rotate our bvecs.
====Tips for Manual Data Inspection====
(a) In Matlab, you can use the "DICOMViewer" command to inspect the original DICOM files. (b) You can use the "fslview" command in FSL to inspect the NIFTI output.
For more info about getting bvecs/bvals from various vendor-specific diffusion sequences, see the [http://wiki.na-mic.org/Wiki/index.php/NAMIC_Wiki:DTI:DICOM_for_DWI_and_DTI NAMIC wiki page on DTI].
= [http://vistalab.stanford.edu/trac/vistasoft/browser/trunk/mrDiffusion/preprocess '''dtiRaw'''] Preprocessing Pipeline =
The tools in mrDiffusion/preprocessing form a new (as of May 2007) DTI processing pipeline to replace our old pipeline. All these function names begin with 'dtiRaw' to help you find them with Matlab's tab-completion. The old processing pipeline was based on tensorcalc- a GE-specific, closed-source eddy-current correction and tensor-fitting program from Maj Hedehus and Roland Bammer. The #Alternative_Method_for_Computing_Tensors_from_NIFTI_Data_Using_FSL was an alternative pipeline that we recommended for external users who could not use tensorcalc. However, we are now switching over to the new open-source pipeline described here. We recommend that off-site users of mrDiffusion switch over as well. That said, the #Alternative_Method_for_Computing_Tensors_from_NIFTI_Data_Using_FSL is fine, so if you have other reasons to use the FSL tools try that method. Note, however, that we may eventually drop support for making a dt6 session file from the FSL output.
{|align="right" class="gallery" width="20%" |+'''Motion/Eddy Correction''' ! style="background:#efefef;" | Original ! style="background:#efefef;" | Corrected |- |align="center"|image:rawDti_orig.gif || image:rawDti_warp.gif |- |align="left" colspan="2"|A single slice from 10 repeats of a 12 direction scan (+1 non-DWI = 13 X 10 = 130 image volumes). On the left is the uncorrected data. The data on the right are corrected using the mrDiffusion implementation of the Rohde et. al. algorithm. (Gif movies made using dtiRawMovie.m.) |}
If the gods currently favor you, you can just type
'''dtiRawPreprocess.m dtiRawPreprocess]'''
point it to your raw dti data (4-d NIFTI format) and an ac-pc aligned t1-weighted NIFTI from the same subject (see Anatomical_Methods), and a few hours later you will have a dt6 file to load into dtiFiberUI. However, this currently only works for our GE data from the Lucas center. If you have data from elsewhere, just make sure that you have your bavls and bvecs files already built, or create a Lucas-center-style dwepi.XX.grads file to pass to dtiRawPreprocess when it asks (see DTI#Converting_DICOM_data_to_NIFTI_format). Most likely, though, you'll need to edit [http://vistalab.stanford.edu/trac/vistasoft/browser/trunk/mrDiffusion/preprocess/dtiRawPreprocess.m dtiRawPreprocess] to suit your needs.
Since you'll likely have to edit [http://vistalab.stanford.edu/trac/vistasoft/browser/trunk/mrDiffusion/preprocess/dtiRawPreprocess.m dtiRawPreprocess], we thought a little tour of the major processing steps would be useful. Note that the only ''really'' slow stage is the eddy-current correction estimation, which will take several hours. The next slowest step is tensor fitting. We use boot-strap variance estimates for our probabilistic tractography (ConTrack), so the default tensor fitting will do 500 tensor fits per voxel. This can take an hour or so for a typical dataset. If you don't the bootstrap variance estimates, you can set the boostrap samples to 0 and your tensor fitting stage will then only take a few minutes. All remaining steps of the preprocessing take at most a few minutes.
dtiRawComputeMeanB0 The eddy-current/motion correction and t1-alignment require a good reference image from the dMRI sequence. This function simply finds all your non-DW images (using the supplied bvals file), aligns them to the first one to correct for motion, and averages them all. The result is saved as a NIFTI file.
** RFD says we should change to "eddy" instead of this and he will help us. That is the form inserted in dipy. **
The first major DTI processing step is eddy-current correction. We have implemented the algorithm described in Rohde et. al. (2004) Comprehensive approach for correction of motion and distortion in diffusion-weighted MRI. MRM, 51(1):103-114 ([http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?db=pubmed&cmd=Retrieve&list_uids=14705050 pubmed]). This algorithm combines a rigid-body 3d motion correction (6 parameters) with a constrained non-linear warping (8 parameters) based on a model of the expected eddy-current distortions. Our implementation builds upon tools from SPM5 to estimate the Rohde model parameters that maximize the normalized mutual information between each diffusion-weighted image (DWI) and the mean of the non-diffusion weighted images (computed in the previous step). We've built the core computational routines as compiled mex functions (C-code) for speed. Even so, estimating the 14 Rohde parameters on a typical dataset of 100 dMRI volumes will take 3-4 hours on a 2GHz Opteron. Note that the image data are not resampled at this stage- the transform parameters are computed and saved.
====[http://vistalab.stanford.edu/trac/vistasoft/browser/trunk/mrDiffusion/preprocess/dtiRawAlignToT1.m dtiRawAlignToT1]==== We like our dMRI data to be oriented in a standard way (AC-PC aligned) and in-register with other scans (like fMRI) that we may have done on this same subject, so we align the dMRI data (using the mean non-DW image) to the t1-weighted high-res volume that is already AC-PC aligned. If you don't have an ac-pc aligned t1, we do have experimental support to allow alignment to a standard template (the MNI EPI image). Simply pass in 'MNI' for the t1 filename.
dtiRawResample We've computed two transforms on the dMRI data- an eddy-current/motion correction and a rigid-body alignment to the anatomical reference image (usally the AC-PC aligned t1-weighted image from the same subject). At this stage, we apply these transforms to the raw dMRI data. The default parameters will resample to 2x2x2 mm voxels with a 7th order b-spline interpolation using the very nice interpolation tools from SPM (spm_bsplins). You can also tweak a flag to get tri-linear interpolation if you prefer that. However, if you use tri-linear, you will likely get banding patterns in your variance estimates of the tensor fits (see [http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?db=pubmed&cmd=Retrieve&list_uids=15955477 this other Rohde et. al. paper]). We've found that by using an interpolation method with a large support (like the 7th order b-spline), these variance errors are below our visual detection threshold, so we don't worry about correcting for them. But, if you choose to use tri-linear interpolation, you will need to deal with this problem (perhaps by implementing the [http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?db=pubmed&cmd=Retrieve&list_uids=15955477 Rohde et. al. method]).
dtiRawReorientBvecs.m Since the dMRI data are rotated somewhat from their original orientation by the t1 alignment and motion correction stages, the diffusion-weighting directions need to be adjusted appropriately before computing the tensors. This function combines the rotation matrix from the t1-alignment with the rotation matrix from the rigid-body component of the Rohde transform and applies it to the bvecs. The result is a new set of bvecs that are correctly oriented with respect to the resampled dMRI data computed in the previous step.
dtiRawFitTensorMex This step fits the tensor model to the resampled dMRI data using the bvals and the reoriented bvecs. We currently use a simple least-squares fit. This does not guarantee a positive-definite 3x3 matrix, so you can end up with some odd things like negative diffusivites and FA values > 1, both caused by negative eigenvalues. Most of the voxels with negative eigenvalues are in fact noise or artifacts and are outside the white matter- usually at the brain edge or in the large blood vessels. But we've found that a precious few (maybe 10 voxels or so for a typical dataset) are valid white matter voxels with data that might be worth looking at. In our experience, these savable voxels are always found in white matter regions of very high anisotropy, usually in the corpus callosum. It is clear that the 2nd and 3rd eigenvalues should be small for such voxels, and sometimes the 3rd eigenvalue dips below zero due to noise in the data. So, you might want to try to fix the tensors in such voxels so that the eigenvalues are positive. The easiest fix is to clip the eigenvalues to zero (try ''dtiFixTensor''). But, if you are a purist, you might fix them by replacing them with a tensor that is the average of their neighbors. Note that the fiber direction (PDD) is unaffected by the negative eigenvalues, and is just as accurate and reliable as any other highly anisotropic voxel. So, fiber tracking is unaffected by this issue. A simulation of the effects of clipping eigenvalues (among other hacks) is described in [http://sirl.stanford.edu/~bob/pdf/DTI/Methods/Landman_Mori_Prince_tensor_fit_ISMRM_2007.pdf this report] from Mori's lab.
To find the voxels in your data that have negative eigenvalues:
dt = dtiLoadDt6('path/to/dt6.mat');
[vec,val] = dtiEig(dt.dt6);
badData = any(val<0, 4);
% To exclude non-white matter voxels, continue with this code
wmProb=dtiFindWhiteMatter(dt.dt6,dt.b0,dt.xformToAcpc);
badData(wmProb<0.8)=0;
numberVoxelsWNegEig=sum(badData(:));
showMontage(badData);
Something that might improve tensor fits is a 'robust' tensor fitting algorithm that ignores outliers. We have an experimental implementation of the [http://www.ncbi.nlm.nih.gov/pubmed/15844157 Chang, Jones and Pierpaoli 2005 RESTORE] algorithm in [http://vistalab.stanford.edu/trac/vistasoft/browser/trunk/mrDiffusion/preprocess/dtiRawFitTensor.m dtiRawFitTensor]- just pass in 'rt' for the fitMethod parameter. Unlike [http://vistalab.stanford.edu/trac/vistasoft/browser/trunk/mrDiffusion/preprocess/dtiRawFitTensorMex.m dtiRawFitTensorMex], this function is not mexified and thus slow, especially when doing the RESTORE fit method. So, it might take several hours to do a single brain.
[http://vistalab.stanford.edu/trac/vistasoft/browser/trunk/mrDiffusion/preprocess/dtiRawFitTensorMex.m dtiRawFitTensorMex] also includes a bootstrap tensor variance estimate option. Currently, this option generates FA and Mean Diffusivity standard deviation maps as well as a principal diffusion direction dispersion map (where dispersion is the direction variance parameter in the Watson distribution). These dispersion maps are in radians, with 0.945 (i.e. cos(1/3) or 54 degrees) representing maximal dispersion on the sphere (see [http://www.ncbi.nlm.nih.gov/pubmed/15906307 Schwartzman et. al., MRM 2005] for details). The default bootstrap algorithm now uses a residual bootstrap method (see [http://www.ncbi.nlm.nih.gov/pubmed/16938472 Chung et al., NeuroImage 2006] and [http://www.ncbi.nlm.nih.gov/pubmed/18779066 DK Jones, IEEE Trans. med. Imaging 2008]). But, you can also do a traditional repetition-based bootstrap if you have enough repeated measurements.
[http://vistalab.stanford.edu/trac/vistasoft/browser/trunk/mrDiffusion/preprocess/dtiRawFitTensorMex.m dtiRawFitTensorMex] is a matlab wrapper around the compiled mex function [http://vistalab.stanford.edu/trac/vistasoft/browser/trunk/mrDiffusion/src/dtiFitTensor.c dtiFitTensor]. The core of this matlab wrapper does something like this: bn = '/path/to/raw/rawDti' dwRaw = readFileNifti([bn '.nii.gz']); bvecs = dlmread([bn '.bvecs']); bvals = dlmread([bn '.bvals']); d = double(dwRaw.data); % a simple brain mask based on thresholding the mean b0 mnB0 = mean(d(:,:,:,bvals==0),4); mask = uint8(dtiCleanImageMask(mrAnatHistogramClip(mnB0,0.4,0.99)>0.25)); % Set up the X matrix. The least-squares fit of the diffusion tensor is computed % by multiplying the log of the diffusion data by inv(X). q = [bvecs.sqrt(repmat(bvals,3,1))]'; X = [ones(size(q,1),1) -q(:,1).^2 -q(:,2).^2 -q(:,3).^2 -2q(:,1).q(:,2) -2q(:,1).q(:,3) -2q(:,2).q(:,3)]; permutations = dtiBootGetPermutations(length(bvals),300,0); [dt,pdd,mdStd,faStd,pddDisp] = dtiFitTensor(d,X,0,permutations,mask); makeMontage3(abs(pdd)); % look at the bootstrap variance maps [fa,md] = dtiComputeFA(dt(:,:,:,2:7)); showMontage(pddDisp/pi180); inds = find(mask); figure;scatter(fa(inds),pddDisp(inds)/pi*180,1); xlabel('FA'); ylabel('PDD dispersion (deg)');
- You will have a set of newly created files in the same directory that your ''dwRawFileName'' was located. As of 09/12/2008, these include a set of raw files (e.g., dti_g865_b900.nii.gz, dti_g865_b900.bvals, dti_g865_b900.bvecs) and aligned files (e.g., all of the above with "_aligned" added to the file names). You will also have an xform computed by the eddy current motion estimation function (e.g., dti_g865_b900_ecXform.mat).
- You will also have a newly created dt6 directory (usually at the top-level of your subject's directory) with a dt6.mat file inside and a bin directory. Inside the bin directory are several niftis (b0, tensors, brainMask), including the results of the bootstrap analyses.
= DTI Preprocessing User Manual = Clicking on the above link takes you to a currently updated account for user best practices, including descriptions of checks to be conducting in/around the preprocessing stage, scripting tips, and more.
= Reprocessing and Debugging Tips = You might find that you want to re-preprocess a dataset with different parameters. If so, it can be worthwhile to understand the ''clobber'' and ''numBootStrapSamples'' input arguments, as these will influence greatly the running time for the function. *''clobber'': if clobber is true, you will recompute and write-over all the files created in the "raw directory" (aka, the same directory that your ''dwRawFileName'' was located, including the "_aligned" files). If all you want to do is recompute the tensor fits (perhaps chaning the number of bootstrap samples), it's okay to set clobber to false and skip recomputing/writing-over these files. *''numBootStrapSamples'': If what you want is a new dt6.mat to do deterministic tracking in dtiFiberUI as soon as possible, you can skip computing the bootstraps, which is very time-consuming. To do this, set this argument to 0.
= Preprocessing DTI Data Using FSL =
Preprocessing DTI Data Using FSL describes an alternative method for computing tensors from NIFTI data using the [http://www.fmrib.ox.ac.uk/fsl/ FSL] tool-kit.
''NOTE:'' see the DTI_Preprocessing#dtiRaw_Preprocessing_Pipeline as an alternative to the FSL preprocessing.
Prior to running the preprocessing pipleline on your diffusion data you want to process a high-resolution t1 anatomical image that you can align the diffusion data to - thus, we begin here...
- In order to align the dti data to the anatomical data you need to provide a t1-weighted image. If you have not previously done this this should be your first step.
- If you need help with this see the [http://white.stanford.edu/newlm/index.php/Anatomical-Processing Anatomical-Processing page].
- We usually place a soft-link to the t1.nii.gz at the top-level of the subject's directory like this: $ ln -s path/to/t1.nii.gz subCode/t1.nii.gz
- '''Note:''' If you have already processed your T1 files through another pipeline you could have difficulties aligning the dti data to your T1. This is due to the header information being different from what is expected by our software. This can be fixed by running (in matlab) '''[http://white.stanford.edu/newlm/index.php/Anatomical-Processing#3._Resample_and_Align_to_AcPc mrAnatAverageAcpcNifti]'''.
- Vistasoft
- Getting Started
- mrVista Overview
- Anatomy
- Functional MRI
- mrVista
- Retinotopy tutorial
- Population RF methods also prf Model, prf_tutorial, prf tutorial
- Diffusion weighted MRI
- Visualization
- Tractography
- Tutorials
- Software overview