MrDiffusion

From VISTA LAB WIKI

Jump to: navigation, search

mrDiffusion is a diffusion imaging analysis and visualization package. It is available as part of Stanford VISTA Lab's open-source and freely distributed mrVista package. To run mrDiffusion you should download the entire mrVista package from our software downloads page.

In addition to mrDiffusion, we frequently use the methods in MRtrix and LiFE to analyze tractography data. These are described in their own sections of the wiki. After publication, we will begin distributing an additional method, Ensemble Tractography developed by H. Takemura. Finally, we note that Pestilli has made enormous progress in speeding up the LiFE calculation and shrinking its

Contents

[edit] Introduction

mrDiffusion is useful for processing diffusion-weighted data and measuring and visualizing fractional anisotropy, mean diffusivity, axial and radial diffusivity, and RGB fiber direction maps. It is also useful for performing deterministic tractography. It will soon support probabilistic tractography through the ConTrack tools that Anthony Sherbony developed and it is integrated with CINCH that David Akers developed.

The mrVista code, mrDiffusion included, is under very active development, so you should try the latest version before reporting a bug. As things stabilize, we will begin to put together stable releases from time to time.

[edit] Software Set-Up

  • Before you run mrDiffusion you will need to download the VISTASOFT package - details can be found on the Software page. You should also follow the steps outlined in our Matlab page to ensure that your Matlab environment is properly set-up and that you can run our software.

Some parts of mrDiffusion depend on routines in SPM (SPM2 or SPM8), so you should download that as well - instructions for how that is done can be found on the SPM page. After you've added all the files from mrDiffusion, mrVista and spm8 to your Matlab path (eg., with "addpath(genpath('/path/to/packages'))") you should be able to run dtiFiberUI, our main DTI analysis front-end.

Sample Data:

  • If you would like to see how mrDiffusion works, get the sample dataset and load it using 'File -> Open tensor (dt6) file...'.

[edit] Starting mrDiffusion

dtiFiberUI is our main DTI analysis front-end. This program provides the interface through which DTI data is visualized and analyzed. Among the many things you can do with dtiFiberUI are:

  • Create ROIs
  • Import ROIs from mrVista
  • Track fibers passing through, or ending in, a given ROI
  • Visualize 3D fiber groups (with mrMesh)
  • Compute ROI statistics
  • Compute fiber group density

To start the mrDiffusion package you type at the Matlab prompt

>> dtiFiberUI

Here are some screenshots to give you a sense of the mrDiffusion visualization features:

[edit] A Quick Tour

This section is a quick tour to allow new users to become oriented to the software and be exposed to some of the most common features and uses. The best way to "learn" the software is to explore it for yourself. You should review this section and use what you learn to work with the sample dataset.

[edit] Opening a Dataset (dt6)

There are a couple of ways to open data in dtiFiberUI. Because you will probably be loading and saving Fibers and ROIs in your subject's folder it is best to change directory in MATLAB to your subject's directory and run the dtiFiberUI command from the directory that contains the dt6.mat file you will be using.

  • Open the dtiFiberUI: Assuming that dtiFiberUI is in your path...
In MATLAB type >> dtiFiberUI
  • Once dtiFiberUI is open you can point it to your dt6.mat file:
File > Open Tensor (dt6) file...
  • Alternatively, once you are in the directory with your dt6.mat file you can type:
>> dtiFiberUI YourDt6FileName.mat

[edit] Backgrounds

The following table provides some information about the different backgrounds you may come across in your DTI analysis.

Image Description
B0
B0
The "b0" map is the non-diffusion-weighted image. It is called "b0" because it the image acquired with a b-value (diffusion weighting strength) of 0.
MD
Mean Diffusivity
Mean diffusivity (MD) is the mean of the eigenvalues of the diffusion tensor or, equivalently, the trace of the diffusion tensor divided by 3. This parameter quantifies the size of the diffusion ellipsoid.
The MD measure is invariant with respect to the orientation of the diffusion tensor and tends to vary little across brain tissue. (Ref.)
FA
Fractional Anisotropy
Fractional anisotropy (FA) is a measure of the degree of anisotropy in a tensor model of diffusion at a given voxel. FA values are bounded between zero (a perfect shape) and one (an infinitely long cigar shape).
In the FA map FA values closer to one appear brighter. (Ref.)
vRGB
Vector RGB
The Vector RGB color map is created from FA values and the three vector elements. The three colors in this map correspond to the direction of the largest eigenvalue of the diffusion tensor (principle diffusion direction) with proportional intensity according to the FA.
RED = left/right, BLUE = superior/inferior, and GREEN = anterior/posterior. (Ref.)
T1
T1
The T1 weighted image provides structural contrast for anatomical images of the brain. The relative signal intensity of voxels within the image depends upon the T1 value of the tissue in that voxel.
Tissues that have a smaller T1 value appear brighter than those with larger values:
(White matter T1 value = ~600ms. Gray matter T1 value = ~900ms. Cerebrospinal fluid T1 value = ~4000ms) - all values for Field strength of 1.5T (Ref.)

[edit] Creating Regions Of Interest (ROIs)

ROIs can be created many different ways. In the ROIs drop-down menu you will find all the options for creating and working with ROIs.

Note that mrDiffusion ROIs are implemented as sets of points in AC-PC space, not the voxel space of any specific image. These points do get mapped to a specific voxel grid when you do something with them (like compute the ROI stats). They are mapped by taking all voxels of the underlying image that contain at least one point. Thus, for the same ROI, you might get a different volume estimates when you map the ROI to images with different voxel sizes (e.g., by calling the ROI stats command from the "Analyze' menu on a high-res anatomy image and on a lower res b0 image). These differences are due to the different quantization error for the different voxel sizes. In other words, mrDiffusion ROIs have no intrinsic volume (points are infinitely small). To compute an ROI volume, mrDiffusion maps the points to the voxels of the underlying image. In general, the higher the resolution of the image, the smaller the volume will be. Imagine the extreme with just one big voxel. Any number of points at any location will produce the same- very large- volume estimate.

Placing an ROI Seed

  • You can click the + or - buttons above a given image to move through slices in that image. Note how the coordinates change under Current Position (you can also type in the coordinates for the seed of your ROI here). Using the drop-down list you can choose image coordinates (default) or talairach (if loaded).
  • ROIs can be created against the background of your choice (see the drop-down list under the Z-image). Most of our ROIs are seeded using the Fractional Anisotropy (FA) background, where white matter fascicles appear brighter. You can also choose to draw your ROI using the directional information given in the Vector RGB map (red = left/right, blue = superior/inferior, and green = anterior/posterior), or using the anatomical data in the T1 image.
  • Left Click to set the seed for your ROI.

Create a Spherical ROI (3D)

  1. Click in any image to identify where you would like the center of the ROI to be.
  2. ROIs > New (sphere)
  3. Enter the radius in mm - click OK.

Create a Rectangular ROI (3D)

  1. Click in any image to identify where you would like the center of the ROI to be.
  2. ROIs > New (rectangular)
  3. In the Build Rectangle dialog box enter the coordinates which you would like the rectangle to span in each dimension - click OK.

Draw a Polygon ROI (2D)

  1. Navigate to the slice in which you would like to draw your ROI. Note the image you want to use to draw the ROI - you have to choose it in the next step.
  2. ROIs > New (polygon) > X, Y, or Z Image
  3. Left click to set points. Right click to close the polygon.

Other Types of ROIs

  • From the ROIs > New Other menu you can create an ROI from the image mask, create an automatic ROI that captures the corpus callosum, create an ROI using the endpoints of the current fiber group and create an ROI from the selected fiber group.

Saving ROIs

You can choose to save one ROI or a group of ROIs
  1. File > ROIs > Save (current/ALL) ROI/s
  2. Choose those ROIs you wish to save and click OK
  • You will be prompted for a directory and a name for the ROI.

Deleting ROIs

You can delete a single ROI, or a group of ROIs
  1. ROIs > Delete some ROIs
  2. Highlight those ROIs you wish to delete - Click OK
  • Note that this will not permanently remove the ROI, it will just remove it from your dtiFiberUI session.

Loading ROIs

You can load one, or multiple ROIs that you previously saved
  1. File > ROIs > Load (many) ROI(s)
  2. Highlight those ROIs you wish to load- Click OK

[edit] Transforming ROIs from mrVista

Often, you will define an ROI in the course of fMRI analyses and want to use the ROI in your DTI analyses. In order to 'see' such functionally-defined ROIs in dtiFiberUI, the ROIs need to be transformed.

  1. Load the ROI of interest in mrVista (Inplane-view)
  2. Open the Volume 3-view window (Inplane 1: Window -> Open volume 3-view window).
  3. In the new volume window, transform the ROI (Volume 1: Xform -> Inplane to volume -> ROIs). It will now appear in the ROI menu of the volume window.
  4. Save the ROI (Volume 1: File -> ROI -> Save ROI). It will be stored in subject/Volume/ROIs
  5. Start dtiFiberUI, and compute a mrVista transform (dtiFiberUI: Xform -> Compute mrVista xform). This process requests a vAnatomy file and requires the spm2 toolbox. It will take a few minutes and spit numbers to the command window like this:
    Xform.png
  6. Import the ROI from mrVista (dtiFiberUI: Xform -> import mrVista volume ROI)
  7. Save the ROI in desired DTI/ROIs directory
  • Note: try not to close any of the mrVista, volume, or dtiFiberUI windows until you are done transforming between them.

This process can be scripted with the function dtiXformMrVistaVolROIs.m (located in VISTASOFT/mrDiffusion/analysisScripts).

[edit] Editing ROIs

Below are some tips on editing the shape, size, name and color of ROIs. Tip: To find your ROI go to ROIs > Find Current ROI. This will bring the ROI to the center of each image. You can also Click the Zoom button and click the area you wish enlarge. This will allow for a better view of your ROI for editing (turn off zoom before manual editing).

Editing the Shape/Size of an ROI

Manual Editing: You can add or remove individual voxels from your ROI.
  1. Choose the ROI you want to edit in the pull-down list in the Current ROI box
  2. Check the box "editable"
  3. To get a better view of your ROI use the Zoom On/Off button. You can click (or drag to create a selection box) in any image view to enlarge the area. Be sure that Zoom is off to edit.
  4. Choose an Edit Shape (default = cube 1)
  5. Click the voxel you wish to add/remove
  • Right Click - Adds to the ROI in the shape and size chosen in the Edit Shape pull-down list
  • Wheel Click - Removes voxels of the size and shape chosen in the Edit Shape pull-down list
  • Left Click - Moves the position of the cross hairs (this will change the view in the other images)
Adding a Predefined Shape: You can also edit an ROI in a predefined way: Add/Remove Sphere, Add/Remove Polygon etc.
  1. ROIs > Edit Shape > Choose action/shape
  2. Enter value for Size (in mm)
  3. Click OK
2D or 3D Floodfilling: You can automatically add to an ROI by using the Grow Similar option. For each option you will see a dialog box that will prompt you to enter the tolerance (smaller values = more strict), the amount of smoothing, and the max size of the ROI (in voxels).
  1. ROIs > Edit Shape > Grow Similar(2D or 3D)
  2. Enter values for tolerance, smoothing, and max radius
  3. Click OK

Editing the Name and Color of an ROI

  1. Choose the ROI you wish to edit in the drop-down list in the Current ROI box at the bottom left of the dtiFiberUI.
  2. Click the edit button to the left of the ROI drop-down list.
  3. A dialog box will appear in which you can enter the name of the ROI and choose the color.
  4. Click OK


[edit] Tracking Fibers

There are many options for tracking fibers. Fibers can be tracked from the Current ROI, with the current ROI, not with the current ROI, and with the current ROI, etc. You can choose which option to use based on experience and your research question. We will provide a couple of examples, but first a word on tracking parameters...

Fiber Tracking Parameters
The Fiber Tracking Parameters dialog window. Values shown are default values.

Once you click on a tracking option in the Fibers menu (e.g., Track from current ROI...) the dialog box shown at right will appear. Below are brief descriptions of how each parameter effects fiber tracking:

  • Step Size - Sets the size (in mm) of each tracking step when generating fibers (default 1mm).
  • Angle Threshold - Sets the threshold for "turns" (in degrees) for a fiber. Fibers making a "turn" above the angle threshold are not generated (default 30 degrees).
  • FA Threshold - This value determines the FA value below which fibers will not be generated. Fibers will end when the FA value is less than what is entered here (default 0.15).
  • Length Threshold - Fibers shorter than this value (given in mm) are discarded. You can also enter two numbers (e.g., 20 50): The first number will be the lower bound (fibers shorter than this value are discarded); the second number will be the upper bound (fibers longer than this value are discarded) (default is 20mm).
  • Puncture Coefficient - See tensor lines paper [Citation Here] (default 0.2).
  • Algorithm Type - Select algorithm used to track fibers (default 1=STT RK4 algorithm).
  • Interpolation Type - Enter 0 for Nearest Neighbor, 1 for Linear.
  • Seed Voxel Offset - This parameter sets the number of seed points per voxel and thus determines how many fibers will be tracked in each voxel. For example, a Seed Voxel Offset of 0.5 will give one seed point per voxel and will track only one fiber per voxel. The default values of 0.25 and 0.75 set 8 seed points per voxel and thus track 8 fibers per voxel.
  • Name - Enter the name for the fiber group (this can also be edited later).


Tracking Fibers From the Current ROI
To get started we'll track fibers from one ROI:

  • Assuming you have already defined your ROI...
  1. In the Current ROI box select the ROI from which you would like to track fibers
  2. Fibers > Track from current ROI...
  3. Click OK - Your fiber group (with the name you entered) will appear in the Current Fiber Group pull-down menu in the dtiFiberUI. The name, visible thickness of the fibers, visibility (0 off, 1 on), and color of the fiber group can be edited by clicking the edit button to the left of the fiber group name.


Tracking Fibers AND/NOT/SPLIT with/by Current ROI
You can use these option to track fibers that are related to your current ROI in some way, either passing through, not passing through, ending in (endpoints) etc. The steps below can be used generally for these options.

  • Assuming you have already tracked fibers from one ROI...
  1. Select a fiber group in Current Fiber Group
  2. Select the ROI of interest in Current ROI
  3. Select the tracking option Fibers > AND/NOT/SPLIT etc...
  4. Enter tracking parameters, or keep defaults, and click OK


====Visualizing Fibers: mrMesh====
mrMesh: Partial Corpus Callosum ROI. Fibers tracked with default parameters. T1 Background.

Once you have tracked fibers you can visualize them in 3D using mrMesh, our 3D DTI visualization tool. Note: Fibers can also be visualized and further manipulated using CINCH.

  1. Select the fiber group(s) you wish to visualize in the Current Fiber Group box.
  2. Check the box labled mrMesh
  3. Select the views you wish to see (radial buttons: Axial, Cor, Sag)
  4. Select the background you want from the pull-down menu
  5. Click Update


Note for 64-bit users: mrMesh uses mrMeshSrv. Users on 64-bit platforms may experience problems with dtiFiberUI starting mrMeshSrv, expressed by a Error connecting to server message in MatLab. This problem is addressed below in the dtiFiberUI FAQ section. I will say here that this problem is most often solved by running mrMeshSrv.exe or mrMeshSrv.glx manually before you try to visualize fibers. The mrMeshSrv executable for Windows and Linux can be found in \VISTASOFT\Anatomy\mrMesh. You can start the server via the command line in Matlab on Linux machines with mrmStart.
If you're having further problems with mrMeshSrv crashing you proabably need to compile it to run on your system. For help with this please see the Mesh Compile page.

====Visualizing fiber endpoints on a 3D mesh====
3D mesh: fiber endpoints density map projected to 3d mesh surface. Colored regions indicate that nearby white matter voxels contain fiber endpoints.

Once you obtained the tractography data, you could visualize the fascicle endpoints in cortical mesh representation.

1. Save tractography data (either whole-brain connectome or tract of interest) into .pdb or .mat (mrDiffusion) format

2. Making the nifti file representing fascicle endpoints (for example, using dtiFiberendpointNifti.m).

3. Loading the nifti on cortical mesh.

We have two mesh tools for visualizing fascicle endpoints. One tool is AFQ Cortical Meshes, which is developed by Jason Yeatman. See AFQ page for instruction on AFQ tool. Hereafter, we will describe the methods for using mrVista tool.

If you have fMRI data from identical subjects, you could visualize the fascicle endpoint in mrVista mesh tools (see Figure). The strong advantage of this methods is that it is easy to compare with Visual Field Maps.

For visualizing the fascicle endpoints saved as nifti in mrVista window, you could use

Xform --> Read Nifti

in Gray Window in mrVista.

This could load the fascicle endpoint on gray matter representation. The requirement is that nifti file of fascicle endpoint should be saved in same resolution as the T1 anatomy file used in mrVista.

[edit] dtiFiberUI FAQ

I was able to download your DTI software and load the sample data set that you have to download, but I need help to do fiber tracking. Do you have a manual or more instructions for that?

  1. Start Matlab
  2. Type dtiFiberUI (assuming dtiFiberUI is in your path)
  3. In the dtiFiberUI GUI: click File->Open Tensor (dt6) file ->select the sample data
  4. Define region of interest from which to trace:
    1. click on Set for either the Z,Y,or X images and click where you want to set your seed (ROI).
    2. ROI->New(sphere) select size (play with new polygon later). Sphere should appear on image in GUI.
  5. Fibers->Track from current ROI (use defaults for now).
  6. click on mrMesh at bottom of GUI.
  7. click on Update to view and rotate fibers.
  8. Note that the brain image displayed will depend on which image is selected under Background on the main GUI and whether you have axial,Cor,and/or Sag selected (next to update).

Tracking works but I couldn't present fibers with mrMesh. I get multiple messages: 'error connecting to server'. How can I fix this?

If it's a 64-bit machine (windows or linux), you may need to manually launch the mrMesh server by finding the executable (in VISTASOFT/Anatomy/mrMesh/mrMeshSrv.glx for linux or VISTASOFT/Anatomy/mrMesh/mrMeshSrv.exe for windows) and running it (double-click it in a file browser or use the command line). Once the server launches, try visualizing something.

  • In Windows, you may need a couple of extra .dlls, msvcr70.dll and MSVCP70.dll, to run mrMeshSrv.exe. If you get this error message when you double click mrMeshSrv.exe, go ahead and download those .dlls, unzip and then save them in C:\Windows\system32\ (if it asks you if you want to overwrite one or both click yes). Those .dll files should also be automatically copied in C:\Windows\wow64. This may fix your problem and mrMesh will then open fine through dtiFiberUI without independently running it. But...
  • If mrMesh does not run in dtiFiberUI without running mrMeshSrv.exe first you may want to create a shortcut to mrMeshSrv.exe on your desktop (right-click mrMeshSrv.exe --> Send to --> desktop (create shortcut)). You will need to double-click that shortcut to run mrMeshSrv.exe before you try to use mrMesh in dtiFiberUI. * Hint: While using dtiFiberUI, don't close the mrMesh Server window. If you do you will have to run mrMeshSvr.exe to connect to the server each time you want to use mrMesh.


The manual ROI drawing tools in dtiFiberUI are somewhat primitive and not ideal for precise drawing of ROIs. Just out of curiosity how do you guys place ROIs? Another software package?

We usually directly import ROIs from our mrVista fMRI analysis package, which is linked to dtiFiberUI using a variety of commands in the Xform menu. We also spend a lot of time defining ROIs using fiber tracking (eg, see this poster). And we sometimes define ROIs by 3d flood-fill ('ROIs->Edit shape->Grow similar (3D)') of hot-spots in functional activation maps that we align and import into dtiFiberUI with 'File->add NIFTI/Analyze image'.

We also sometimes import MNI-space atlases, such as the Amunts cytoarchitechtonic maps and the AAL maps, and use these to create ROIs (by 3d flood-fill and/or 'ROIs->Edit Shape->Restrict to image values'). You can import an MNI map with the 'File->Add normalized map' menu item. Of course, you first need to have computed the MNI spatial normalization parameters. My dtiMakeDt6 does this automatically, but if you used dtiMakeDt6FromFsl, you won't have it. To compute it, be sure that you have spm2 or spm5 in your path and that you are viewing the t1 image. Then select 'Xform->Compute spatial normalization' and select the SPM MNI T1 template and be sure to name the normalization 'MNI'. You can also try normalizing your b=0 image to the MNI EPI template, but I've found this to often produce a pretty bad match. Be sure to check the normlaization by 'File->Add normalized map' and selecting the template that you just normalized to. It should get warped by the inverse deformation to match your native subject space.

When I compute FA for a fiber group using the GUI it is very different than if I convert that fiber tract to an ROI and compute the FA. Can you advise which method is the preferable way to calculate FA?

The two results represent two fundamentally different approaches to summarizing the properties of a fiber group.

  • Weighted-mean approach: When you use the fiber properties GUI, you are calling dtiFiberSummary.m, which converts the fiber group into a set of points, interpolates the tensor value at each point, computes diffusion properties from these many interpolated tensors, and summarizes them with a simple mean. (See also dtiFiberProperties.m.) This approach heavily over-samples voxels in the core of the fiber tract. In effect, this gives you a weighted average of the voxels that the fibers pass through, with more weight going to voxels at the core (with many fibers) and less weight o voxels at the fringe 9with fewer fibers). It's the approach that we used in the 2007 PNAS paper, and we describe and justify it a bit in there.
  • Unweighted-mean approach: When you first convet to an ROI and take the stats, you are using dtiGetRoiStats.m, which takes the mean of the unique voxels in the ROI. When you conver a fiber group to an ROI, any voxel that has at least one fiber passing though it is included in the ROI. So, this approach is a non-weighted average of the voxels that the fibers pass through. It also does the interpolation differently, but this is a negligable effect compater to the weighted vs. non-weighted average.

Which approach you should use depends on what you want to do with the results. For simply summarizing the diffusion properties of a fiber group, I prefer the weighted average approach. Note, however, that variance estimates from this approach are not valid- they are underestimated due to the oversampling. However, if you simply want to collapse a fiber group in each subject to a simple summary and then compute variance across subjects (e.g., what we did in the PNAS paper), then this it fine.

For more on this topic, see the section on tensor summary measures.

[edit] mrDiffusion command-line tips

dtiFiberUI 'dt6file.mat' : to load a particular dt6file when opening dtiFiberUI

[edit] Loading stuff

% load the tensors into dt6 format
[dt6, xformToAcpc, mmPerVoxel] = dtiLoadTensorsFromNifti(niFile);
% load an ROI
roi = dtiReadRoi(roiName);
% load some fibers
fg = dtiReadFibers(fiberName);

[edit] Creating and Saving

% create a new ROI struct
sphere = dtiNewRoi(sphereName, 'r', coords);
% save the ROI
dtiWriteRoi(roiVar, desiredPath&Name);
% save the fiber group
dtiWriteFiberGroup(fiberVar, desiredPath&Name);
% build a spherical ROI with a center at the coordinates given, and of specified radius:
sphere.coords = dtiBuildSphereCoords(centerCoords, radius);
% track fibers from a particular ROI:
fg=dtiFiberTrack(many parameters, including opts, see parameters below);

[edit] Modify/analyze commands

% modify an existing fiber group relative to an ROI, creating a new fiber group:
fgAndRoi=dtiIntersectFibersWithRoi(many parameters, including and/or/endpts);
% find the center of an ROI:
centerCoords = round(mean(roi.coords,1)*10)/10
% prints summary of fiber group to the screen, including (1) mean length, (2) number of fibers, (3) mean FA 
% (higher = more spherical?), (4) mean diffusivity (lower = less diffusion scatter?)
txt = dtiFiberSummary(handles); disp(txt);
% script prints the fiber density volume to the screen -- units are number of voxels (?) 
% that are included in the fiber group (?)
dti_OTS_fiberDensityVolume

[edit] Useful code snippets

[edit] Save a mrDiffusion ROI as a binary NIFTI image
% Save a mrDiffusion ROI as a binary NIFTI image
dtiExportRoiToNifti('path/to/roi.mat', 'path/to/bin/b0.nii.gz', 'path/to/new/roi.nii.gz');
[edit] Save a fractional anisotropy (fa) and mean diffusivity (md) map
% Save a fractional anisotropy (fa) and mean diffusivity (md) map
% Note that tensors.nii.gz is in the dtiNN/bin folder, where NN is the # of directions. 
[dt6, xformToAcpc, mmPerVox] = dtiLoadTensorsFromNifti('path/to/tensors.nii.gz');
[fa,md] = dtiComputeFA(dt6);
% Save fa and md in NIFTI format:
dtiWriteNiftiWrapper(fa, xformToAcpc, 'fa.nii.gz', 1, 'FA');
dtiWriteNiftiWrapper(md, xformToAcpc, 'md.nii.gz', 1, 'MD');
% Or, to save them in Analyze format (useful for spm2), use this:
[f,p] = uigetfile({'*.nii.gz';'*.*'},'Select a NIFTI tensor file...');
[dt6, xformToAcpc, mmPerVox] = dtiLoadTensorsFromNifti(fullfile(p,f)); 
b0Nifti = readFileNifti(fullfile(p,'b0.nii.gz'));
[fa,md] = dtiComputeFA(dt6);
origin = mrAnatXformCoords(inv(xformToAcpc),[0 0 0]);
analyzeWrite(fa, 'NN_fa', mmPerVox, 'FA', origin);
analyzeWrite(md, 'NN_md', mmPerVox, 'MD', origin);
analyzeWrite(b0Nifti.data, 'NN_b0', mmPerVox, 'b0', origin);
[edit] Save FA and MD maps directly from the raw DW data
% Save FA and MD maps directly from the raw DW data
bn = '/path/to/raw/rawDti'
dwRaw = readFileNifti([bn '.nii.gz']);
bvecs = dlmread([bn '.bvecs']);
bvals = dlmread([bn '.bvals']);
sz = size(dwRaw.data);
d = double(dwRaw.data);
xform = dwRaw.qto_xyz;
% Make a simple brain mask
b0 = mean(d(:,:,:,bvals==0),4);
% Tidy up the data a bit- replace any dw value > the mean b0 with the mean b0.
% Such data must be artifacts and fixing them reduces the # of non P-D tensors.
for(ii=find(bvals>0)) tmp=d(:,:,:,ii); bv=tmp>b0; tmp(bv)=b0(bv); d(:,:,:,ii)=tmp; end 
[b0,clipVal] = mrAnatHistogramClip(b0, 0.4, 0.99, false);
% Mask out the bottom 20% or b0 intensities
maskThresh = 0.2*diff(clipVal)+clipVal(1);
mask = uint8(dtiCleanImageMask(b0>maskThresh));
q = [bvecs.*sqrt(repmat(bvals,3,1))]';
X = [ones(size(q,1),1) -q(:,1).^2 -q(:,2).^2 -q(:,3).^2 -2*q(:,1).*q(:,2) -2*q(:,1).*q(:,3) -2*q(:,2).*q(:,3)];
[dt6,pdd] = dtiFitTensor(d,X,[],[],mask);
% log(b0) is stored in dt6(:,:,:,1)
dt6 = dt6(:,:,:,2:7);
makeMontage3(abs(pdd));
[eigVec,eigVal] = dtiEig(dt6);
v1 = squeeze(eigVec(:,:,:,[1 2 3],1));
[fa,md,rd] = dtiComputeFA(eigVal);
Personal tools