-
Notifications
You must be signed in to change notification settings - Fork 141
DWI Initialization
This page describes the steps involved in preprocessing DTI data. It includes sections for dealing with GE data from Stanford as well as general information for those who have data from Siemens and Phillips scanners. It also describes the preprocessing pipeline.
The VISTASOFT preprocessing pipeline for DWI data is described in detail in the sections below. Breifly, the steps are:
- Process a high-resolution T1-weighted anatomical image for DWI alignment: This step is optional. If you have raw scanner dicoms that you would like to process please visit the Anatomical-Processing page. If you would like to align your DWI data to the MNI-EPI template you can do so by passing in the string 'MNI' instead of a filename when initializing your DWI data - more on that later.
- Convert raw DWI dicom data to 4-D NIFTI: This is the first step for most users.
- Bvals and Bvecs: Many methods exist for doing this, you must assure that whatever method you choose has provided you with two text files (bvals and bvecs) which have information regarding the bvalue and the gradient vectors. Some users will have these files from whomever set up their DWI scan-sequence.
- Initialize the DWI data: This is done with a single call to the function dtiInit. Many parameters can be set (with a call to dtiInitParams) for particular needs and to call different algorithms for tensor fitting.
- Wait: Depending on your parameters the pipeline can take between ~30 min and ~3 hrs. Relax and be patient....
At the end of this section you will have a dt6.mat file that you can load into mrDiffusion. Note: This pipeline was written for use with Matlab 2009 and later.
The first step in the processing of diffusion data is to go from the raw dicom files that the scanner outputs to a 4-D NIFTI file that can be fed into our pipeline. Depending on your scanner vendor, and your location, this method of going from raw dicoms to NIFTI will differ. Below we have details that will help you get from dicom to NIFTI.
Before you begin: You will want an anatomical image to align to: Skip to the section below if you need help with processing your anatomical data.
The directions in this section are specific to the VISTA LAB computers and data collected at CNI. Please see the sections below for more general directions.
-
1. Download the data (.tgz file) from NIMS.
2. Go to the folder in with the NIMS session.
$ cd /biac2/wandell2/data/diffusion/<subject></subject></3.>
$ tar -xvzf YYYYMMDD_XXXX.tgz
3. Set the folder containing the raw Diffusion files.
$ mv -v YYYYMMDD_XXXX raw
4. Perform field map correction. (in Matlab, about 2 hours)
>> cd /biac2/wandell2/data/diffusion/<subject>/YYYMMDD_XXXX/ >> preprocessing_script_YYYYMMDD_XXXX.m
</ol></subject>
The directions in this section are specific to the VISTA LAB computers. Please see the sections below for more general directions.
If you have collected your diffusion data from the CNI, then you already have the required files available for download from the CNI database. However, if you have downloaded just the dicoms, and you prefer to start from there then you can convert those dicoms to a 4-D NIFTI file in the following way.
-
1. Zip the DTI directory which contains the RAW data from the scanner.
- This directory is usually something like subjectCode/raw/dti_data.
- We like to keep the raw dicom files in a zip archive to save space and keep the large numbers of individual files in one container. If your DTI dicom files are not already zipped up, you can zip the directory with the following command in the terminal (replacing the folder names with yours):
- The niftify shell-script is installed on all the vision group linux machines (in /usr/local/) and works with a .zip file that contains the images from the scanner. 'niftify' will unzip the raw data to a tmp directory (no need to unzip them yourself) and build the nifti file using dinifti.
- This is done in a terminal with the command:
- Then gzip the nifti file with the command:
- This will result in a file named dti_data.nii.gz
The directions in this section are meant for users outside the VISTA Lab.
- Use dinifti to create the nifti file from the raw data.
- This is done in a terminal with the command:
$ dinifti -g -s <# of slices> <input></input> <out>
- Note: the -g creates a .nii.gz file and -s indicates that the next file is the number of slices.
- This process will return a file with the .nii.gz extension.
- Now change the file name to reflect the b-value and dwepi (grad code) in the file name. This will allow dtiRawPreProcess to choose the correct files for processing your data.
Example: If dwepi = 87 and b-value = 900 your file name should be: '''dti_g87_b900.nii.gz'''
2. Zip the DTI directory which contains the RAW data from the scanner.
- This is an optional step. We like to keep the raw dicom files in a zip archive to save space and keep the large numbers of individual files in one container. If your DTI dicom files are not already zipped up, you can zip the directory with the following command in the terminal (replacing the folder names with yours):
$ zip -r dti_data.zip dti_data
</ol></out>
This section is for those who have collected diffusion MR data on a Siemens or Phillips scanner, or on a GE scanner using the GE-supplied diffusion sequence (rather than Roland Bammer's DTI sequence at the Stanford Lucas center).
- Note: For Siemens data, we need to use newer version of micron (newer than 11/11)
-
1. Create a NIFTI from your raw DTI data:
dcm2nii to create the nifti file from the raw data.
- This is done in a terminal with the command: $ dcm2nii -g -d fullPath/*.dcm
- Note: the -g creates a .nii.gz file and -d places the date in the name of the nifti file.
- This process will return three files: (1) A file with the .nii.gz extension (2) a file with a .bval extension and (3) a file with a .bvec extension.
- NOTE: Don't change the file names - this will cause problems.
The .bval and .bvec files should be created from the raw dicom headers by whichever method you used to create the raw dwi NIFTI file in the previous section. Some methods (e.g., dinifti, dicomToNifti.m) do not create these files for you. If you don't have those files (and you can't create them yourself (no easy task)) you should talk to the person who programmed your DWI sequence and get the files that are specific to your use case.
- For simplicity place your bvals and bvecs files in your raw data directory where your raw DWI NIFTI file is (hopefully in a folder called 'raw') with the same prefix as your raw DWI NIFTI.
[Step] - dtiInit
dtiInit is our main workhorse initialization code. dtiInit acts as a sophisticated wrapper for the dtiRaw preprocessing pipeline (see the dti Raw Preprocessing Pipeline section below).
- Make sure you have your directory structure set up such that you have:
<subdir></subdir> _____|_____ | | t1.nii.gz raw | dwi.nii.gz dwi.bval dwi.bvec
Specific parameters can be set using the function dtiInitParams. A call to this function will return the default parameters (as shown below). The user can change these params using 'varargin' or by simply assigning them to a structure and redefining the fields in the structure. For information about what each of these options (and the up to the second defaults) see the help for the function by clicking the function name --> dtiInitParams.m, or by typing help dtiInitParams in Matlab.
>> dwParams = dtiInitParams dwParams = bvalue: [] gradDirsCode: [] clobber: 0 dt6BaseName: '' flipLrApFlag: 0 numBootStrapSamples: 500 fitMethod: 'ls' nStep: 50 eddyCorrect: 1 excludeVols: [] bsplineInterpFlag: 0 phaseEncodeDir: [] dwOutMm: [2] rotateBvecsWithRx: 0 rotateBvecsWithCanXform: 0
- Note: For Siemens data, we need to change the parameter as follows:
dwParams.rotateBvecsWithCanXform = 1
Running dtiInit is very simple and it can be called from Matlab without any input arguments - provided you wish to use the default parameters during initialization (see the previous heading for details on parameters). If you don't provide inputs dialog boxes will pop and ask you to select the raw DWI NIFTI data as well as the T1 NIFTI. If there is more than one bvecs and/or bvals file in your raw directory another dialog box will pop and ask you to select those you wish to use.
Example Usage:
Using default params:
[dt6FileName,] = dtiInit <or></or> [dt6FileName,] = dtiInit('rawDti.nii.gz','t1.nii.gz')
Using previously set parameters:
[dt6FileName,] = dtiInit([],[],dwParams);
Using varargin to set specific params (arbitrary example using clobber and phaseEncodeDir):
[dt6FileName,] = dwParams = dtiInitParams('clobber',1,'phaseEncodeDir',2); [dt6FileName,] = dtiInit('rawDti.nii.gz','t1.nii.gz', dwParams)
... and you're done! You should have a directory structure that looks something like the following.
______________ <subdir></subdir> ____________ _____|_____ | | | | t1.nii.gz raw dti40trilin | ____________|_____________ dwi.nii.gz | | | 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
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 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. In matlab, run:
>> 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 MRIConvert from the University of Oregon, dinifti from the NYU CBI, and 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.
(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 NAMIC wiki page on DTI.
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 FSL pipeline 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 FSL pipeline 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.
<caption>Motion/Eddy Correction</caption>Original | Corrected |
---|---|
uncorrected | eddy/motion corrected |
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, 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 the DICOM-to-NIFTI section above). Most likely, though, you'll need to edit dtiRawPreprocess to suit your needs.
Since you'll likely have to edit 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.
====dtiRawRohdeEstimateEddyMotion==== 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 (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.
====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 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 Rohde et. al. method).
====dtiRawReorientBvecs==== 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 this report from Mori's lab.
To find the voxels in your data that have negative eigenvalues:
Something that might improve tensor fits is a 'robust' tensor fitting algorithm that ignores outliers. We have an experimental implementation of the Chang, Jones and Pierpaoli 2005 RESTORE algorithm in dtiRawFitTensor- just pass in 'rt' for the fitMethod parameter. Unlike 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.
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 Schwartzman et. al., MRM 2005 for details). The default bootstrap algorithm now uses a residual bootstrap method (see Chung et al., NeuroImage 2006 and DK Jones, IEEE Trans. med. Imaging 2008). But, you can also do a traditional repetition-based bootstrap if you have enough repeated measurements.
dtiRawFitTensorMex is a matlab wrapper around the compiled mex function dtiFitTensor. The core of this matlab wrapper does something like this:
bn = '/path/to/raw/rawDti' dwRaw = readFileNifti([bn]); bvecs = dlmread([bn]); bvals = dlmread([bn]); 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)]; 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/pi*180); 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.
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.
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.
This page describes an alternative method for computing tensors from NIFTI data using the FSL tool-kit.
NOTE: see the dtiRaw code 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 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) 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