Developed by
Fang-Cheng Yeh
Department of Biomedical Engineering, Carnegie Mellon University  

e-mail: frankyeh (at) cmu.edu

Documentation‎ > ‎

Reconstruction(DTI, QBI, DSI, GQI)

Introduction

The following manual introduces the reconstruction processes implemented on DSI Studio. The method to parse DICOM images parsing is detailed in parsing DICOM/Analyze/2dseq images, which generate ".src" files for further reconstruction. 

Open DSI Studio and press the "step 2:reconstruction" button in the main project window to select the .src file. DSI Studio will present a reconstruction window as shown in the figure to the right.
 

Setup a Mask


The purpose of the mask is to filter out the background region, increase the reconstruction efficacy, and facilitate further visualization. In the mask selection window, there is a track bar at the bottom for slice selection. The mask can be generated using several build-in functions provided by DSI Studio:
  
"Thresholding" generates an initial selection. 
"Smoothing" smooths contour. 
"Expansion" expands the current selection.
"Erosion" shrinks the contour.
"Defragment" filters out small fragments.
 
Recommendation steps include thresholding, smoothing, and defragment. 

You may open/save the mask to a txt file or nifti file.
 

Reconstruction Options


ODF Tessellation


This parameter determines the number of the sampling directions of the ODF. 8-fold results in 642 sampling directions, where 6-fold results 362. It is recommended to use the highest fold to ensure the best angular resolution.

In DSI Studio, the sampling directions are discretized into 162, 252, 362, or 642 directions, which are equally distributed on a unit sphere. The sampling directions of ODF are obtained by projecting a tessellated icosahedron to a sphere. 4-fold or 5-fold tessellated icosahedron is capable of generating 162 or 252 sampling directions. Higher fold of ODF is recommended for a more detailed presentation of fiber directions.


ODF Sharpening


This will conduct diffusion deconvolution to the ODFs generated from DSI, QBI, or GQI [10]. The parameter determines the error percentage defined in the paper [10]. Lower value results in more regularization and less sharpening effect but also more sensitive to noise.

If you have limited diffusion sampling e.g. less than 100, ODF sharpening may improve the angular resolution and assist fiber tracking.

Number of fibers resolved

This number determines the maximum number of fiber resolved in each voxel.

Multi-thread reconstruction

Enable this option to make use of multi-thread reconstruction, which can greatly reduced the computation time.

Record ODF in the output file

Enable this option to store all ODF in the .fib file and allow DSI Studio to visualize the ODFs for inspection. Note enabling this option will greatly increase the size of the .fib file.

Half-sphere scheme

Enable this option if you are using a half sphere DSI or HARDI scheme.



Reconstruction Methods


DTI
The DTI was proposed by Basser et.al. [1]. It is able to characterize the major diffusion direction of the fiber in human brains. The DTI implementation detail in DSI Studio is following the paper [2], the same as DTI Studio, which is a powerful tool developed by the Johns Hopkins MRI team. The reconstruction performs eigenanalysis on the calculated tensor, and the indices such as FA, MD, and three eigenvalues are also exported.
 

QBI (Funk-Radon)

. In the QBI reconstruction method proposed in the paper, two major parameters are defined before reconstruction. One is the width of the Guassian smoothing kernel, and another is the width of the interpolation kernel. 
 
Based on experiences, the recommended values for the parameters is 5 degree for interpolation kernel and 15 for the smoothing kernel. Moreover, it is recommended that the sampling directions of the ODF be less than the number of sampling points in ODFs; otherwise, the noise will be amplified during the convolusion/deconvolusion processes.
 

QBI (Spherical Harmonics)

DSI Studio is able to reconstruct QBI[3]. Three reconstruction approaches to reconstruct QBI using spherical harmonics were independently proposed by Anderson, Hess et. al. and Descoteaux et.al . DSI Studio adopted Descoteaux 's method [4].

A regularization parameter is required for the reconstruction. The recommended value is 0.006 [4].
 

DSI

The DSI reconstruction method [5] is based on the Fourier transform relation between q-space signals and PDF of the diffusion pattern. The detailed processing flow of our standard DSI reconstruction is as follows:
  1. The dimension of the q-space is set to 16x16x16. (3-dimensional double floating-point array)
  2. To reduce the truncation artifact, the q-space data are applied by a Hanning filter, 0.5*(1+cos(2*pi*q/16)), where q is the q-space distance to the origin in the unit of the grid point, and 16 is the filter bandwidth.
  3. Perform 3D Fourier transform on the q-space data to obtain PDFs.
  4. The ODFs are estimated by integrating several discrete points along each sampling direction in the PDFs (from 2.1 to 6.0, step = 0.2). Tri-Iinear interpolation is used to obtain the estimated values from the discrete data of the PDFs. Before integration, the value of each sampling point is weighted by the square of the distance to the center and integrated to obtain the distribution value.
  5. The ODF obtained can be used to determine fiber directions by searching the local maximums or applying the ODF decomposition method 
The matlab code for DSI reconstruction is the follows:

% Thanks Jiaying Zhang for sharing this DSI code
% the dimension of data used here is n*1, and n refers to the number of 
% images, including the image of b0 .
% the dimension of odf_vertices used here is 642*3.
% the dimension of q table is n*3.
% Note that your q table loaded here is the same as the mgh_dsi_q.txt.

q=load('mgh_dsi_q.txt');
% The download link for this file in at the bottom of this page

% Hanning filter
hanning_filter=zeros(size(b_table,1),1);% b_table is 203*4,515*4 etc
for i=1:n
    hanning_filter(i,1)=0.5*(1.0+cos(2.0*pi*sqrt(sum(q(i,:).*q(i,:)))/16));
end
% the signal after hanning filter
value=s.*hanning_filter;   %s, n*1, the modulus of the signal

% to get the q space index
q=q+9;

% To value sq
sq=zeros(16,16,16);
for i=1:n
    sq(q(i,1),q(i,2),q(i,3))=value(i);
end

%To calculate PDF of each voxel
Pr=fftshift(abs(real(fftn(fftshift(sq),[16,16,16]))));

% To calculate the ODF of each voxel
odf=zeros(length(odf_vertices)/2,1);
for m=1:length(odf_vertices)/2
   origin=9;
   r=2.1:0.2:6.0;
   xi=origin+r*odf_vertices(1,m);
   yi=origin+r*odf_vertices(2,m);
   zi=origin+r*odf_vertices(3,m);
   PrI=interp3(Pr,xi,yi,zi,'linear');   
   for i=1:20
       odf(m)=odf(m)+PrI(i)*r(i)^2;
   end
end

A parameter for Hanning filter is required. The recommended value is 16 for the 16x16x16 window. The obtained ODF can be visualized, as shown in here.

One should note that the "fa" generated is in fact a new index called quantitative index (QA), which is defined in GQI paper [9].

GQI


Generalized q-sampling imaging (GQI) is a q-space reconstruction method that can reconstruct ODFs from a veriety of diffusion datasets, including DSI dataset, HARDI, multiple-shell, combined DTI dataset or even body center cubic (BCC) dataset [6][7][9]. It is shown that QBI is a special case of GQI, and DSI reconstuction is equivalent to GQI [9]. The reconstruction requires only a matrix multiplication, and the implementation is simple.

The reconstruction requires users to input the preferred diffusion sampling length ratio (e.g. 1.2) The following is the matlab code for GQI reconstruction.

function ODF=gqi_l(S,odf_vertices,b_table,mean_diffusion_distance_ratio);
% S is an Nx1 vector, storing the diffusion MR signals of a voxel. N is the number of the diffusion images acquired.
% odf_vertices is a 3x362 matrix storing the ODF directions
% b_table is a 4-by-N matrix, where N is the number of diffusion images. For each diffusion image, the b-value, bx, by, and bz are the four elements for the matrix. Note that [bx by bz] should be a unit vecter.
% mean_diffusion_distance_ratio a constant parameter, typically around 1.0 ~ 1,3


l_values = sqrt(b_table(1,:)*0.01506);
b_vector = b_table(2:4,:).*repmat(l_values,3,1);
ODF = sinc(odf_vertices'*b_vector*mean_diffusion_distance_ratio/pi)*S;


Use Case Example:

ODF12 = gqi_l(S,odf_vertices,b_table,1.25);
 
To show the reconstructed odf, use the example code in ODF visualization page:
     
The test example is also avaliable in the attachment file provided at the bottom of this web page.

One should note that the "fa" generates from GQI is in fact an index called quantitative index (QA), which is defined in GQI paper [9].


Recommended scanning parameters for GQI

Description

High angular resolution diffusion MRI sequence for GQI analysis

Acquisition time

14min 21s

Type

2D

Slices

55

Slice gap

none

Slice thickness

2.5 mm

Slice order

Interleaved

FOV

240mm x 240mm

Matrix

96 x 96

Resolution

2.5 mm x 2.5 mm

TR

8200 ms

TE

69 ms

Diff directions

101, b-table, b_table for SIEMENS

Diff weightings

multiple

max b-value

4000 s/mm2

flip angle

90 deg

Bandwidth

2004 Hz/Px

Parallel imaging

GRAPPA 
Acceleration factor PE: 3 
Reference lines PE: 40


Q-Space Diffeomorphic Reconstruction


Q-Space diffeomorphic reconstruction (QSDR) [11] is a further generalization of GQI that allows users to construct ODF in any template space (e.g. MNI space). QSDR provides a direct approach to analyze the group difference and also facilitates the comparison using fiber tracking. QSDR can be used with any spatial registration method (e.g. SPM normalization, LDDMM, ...etc), and it allows users to combine other registration method with QSDR to conduct ODF reconstruction in a template space. 

In current version, DSI Studio users the transformation and normalization matrix provided by the SPM5 to do QSDR. The steps are the following:

1. Perform linear mapping (SPM coregister method) between subject's b0 image and structural image (T1/T2 images)

    1.1 [Coregister]->[Coreg:Estimate]
    1.2 Assign b0 image in the source image
    1.3 Assign structural image in the reference image
    1.4 Export the transformation matrix
         In the spm_coreg.m, find out the following codes which shows the coregisteration result.

text(0,0.5, sprintf('X1 = %0.3f*X %+0.3f*Y %+0.3f*Z %+0.3f',Q(1,:)),'Parent',ax);
text(0,0.3, sprintf('Y1 = %0.3f*X %+0.3f*Y %+0.3f*Z %+0.3f',Q(2,:)),'Parent',ax);
text(0,0.1, sprintf('Z1 = %0.3f*X %+0.3f*Y %+0.3f*Z %+0.3f',Q(3,:)),'Parent',ax);

         Before or after this code segment, insert the command save('Q.mat','Q') to export the "Q" matrix. You may need to run the program again to get the matrix.

2. Perform spatial normalization between subject's structural image and MNI template

    2.1 [Normalize]->[Normalise:Estimate]
    2.2 Assign the subject's structural image to the [Data]->[Subject]->[Source Image]
    2.3 Assign the corresponding MNI template to the [Template Image]
    2.4 Perform normalization and the obtain the "sn" matrix (e.g. s07080_-0008-00001-000001-00_sn.mat)

3. Combine the transformation and normalization result using the following matlab code

    You need to assign the path for the xxxxxx_sn.mat, Q.mat, and the output tr.mat files. The b0_matrix_size is the image dimension of your b0 images. For example, a b0 image of 128x128 has b0_matrix_size = 128.

function spm_convert_matrix(sn,Q,output,b0_matrix_size)
% spm_convert_matrix('sn.mat','Q.mat','tr.mat');
sn = load(sn);
Q = load(Q);
TrX = reshape(sn.Tr(:,:,:,1),7,63);
TrY = reshape(sn.Tr(:,:,:,2),7,63);
TrZ = reshape(sn.Tr(:,:,:,3),7,63);
Affine = [1,0,0,0;
          0,-1,0,b0_matrix_size+1;
          0,0,1,0;
          0,0,0,1]*inv(Q.Q)*sn.Affine;
save(output,'TrX','TrY','TrZ','Affine','-v4');

return

    This gives the tr.mat which can be used by DSI Studio to do the QSDR.

4. Perform QSDR

    Load .src file in the DSI Studio and select q-space diffeomorphic reconstruction. Assign the tr.mat matrix and finish the reconstruction.

Detect Fiber Orientations

The popular way to determine the fiber directions from the ODFs is to search for local maximums.

The following is the example code for detecting local maximum:

function p = find_peak(odf,odf_faces)
is_peak = odf;
odf_faces = odf_faces + 1;
odf_faces = odf_faces - (odf_faces > length(odf))*length(odf);
is_peak(odf_faces(1,odf(odf_faces(2,:)) >= odf(odf_faces(1,:)) | ...
    odf(odf_faces(3,:)) >= odf(odf_faces(1,:)))) = 0;
is_peak(odf_faces(2,odf(odf_faces(1,:)) >= odf(odf_faces(2,:)) | ...
    odf(odf_faces(3,:)) >= odf(odf_faces(2,:)))) = 0;
is_peak(odf_faces(3,odf(odf_faces(2,:)) >= odf(odf_faces(3,:)) | ...
    odf(odf_faces(1,:)) >= odf(odf_faces(3,:)))) = 0;
[values,ordering] = sort(-is_peak);
p = ordering(values < 0);
end

The return vector p stores the "index" of the peak orientations. To access the orientation, use odf_vertices(:,p);
 
  
References
 
[1] P.J. Basser, J. Mattiello and D. LeBihan, Estimation of the effective self-diffusion tensor from the NMR spin-echo, J Magn Reson 1994;B 103(3):247-254.
[2] Jiang H, van Zijl PC, Kim J, Pearlson GD, Mori S. DtiStudio: resource program for diffusion tensor computation and fiber bundle tracking. Comput Methods Programs Biomed. 2006;81:106–116.
[3] Tuch D.S. Q-ball imaging. Magn. Reson. Med. 2004;52:1358–1372.
[4] Descoteaux, M., Angelino, E., Fitzgibbons, S., Deriche, R., 2007. Regularized, fast, and robust analytical Q-ball imaging. Magn Reson Med 58, 497-510.
[5] Van J. Wedeen, Patric Hagmann, Wen-Yih Isaac Tseng, Timothy G. Reese, Robert M. Weisskoff, Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging, MRM 2005 54:1377–1386
[6] F-C. Yeh, V. J. Wedeen, and W-Y. I. Tseng "Dataset-independent reconstruction of high angular resolution diffusion sampling schemes by generalized q-space imaging" Electronic poster, proc: 17th Scientific Meeting & Exhibition, Honolulu, Hawai'i, USA 18-24 April 2009.
[7] F-C. Yeh, V. J. Wedeen, and W-Y. I. Tseng "Practical crossing fiber imaging with combined DTI datasets and generalized reconstruction algorithm" Oral presentation, proc: 17th Scientific Meeting & Exhibition, Honolulu, Hawai'i, USA, 18-24 April 2009. (link)
[8] Fang-Cheng Yeh, Van J. Wedeen, Wen-Yih Isaac Tseng. “A recursive algorithm to decompose orientation distribution function and resolve intra-voxel fiber directions." Proc 16th Scientific Meeting & Exhibition, Toronto, Ontario, Canada, May 3-9, 2008. (link)
[9] Fang-Cheng Yeh, Van J. Wedeen, Wen-Yih Isaac Tseng. “Generalized Q-Sampling Imaging." IEEE Trans Med Imaging 27, 1415-1424, 2010.
[10] Yeh, F.C., Wedeen, V.J., Tseng, W.Y., 2011. Estimation of fiber orientation and spin density distribution by diffusion deconvolution. Neuroimage 55, 1054-1062.
[11] Yeh, F.C., Tseng, W.Y., 2011. NTU-90: A high angular resolution brain atlas constructed by q-space diffeomorphic reconstruction. Neuroimage 58, 91-99.
 
 
Č
ċ
ď
DiffusionVectors_DSI_HALF_101.txt
(6k)
Fang-Cheng Yeh,
Feb 3, 2012 9:51 AM
ċ
ď
b_table101.txt
(6k)
Fang-Cheng Yeh,
Feb 3, 2012 9:50 AM
ċ
ď
data.mat
(2445k)
Fang-Cheng Yeh,
Feb 3, 2010 9:23 PM
ċ
ď
mgh_dsi_q.txt
(4k)
Fang-Cheng Yeh,
Jun 22, 2011 11:58 AM