Documentation‎ > ‎

Introduction to reconstruction methods

DSI Studio supports several reconstruction methods, but what is their difference? To better understand the pros and cons of each reconstruction method, we can categorize them into either model-based methods or model-free methods. This categorization is similar to the classification of the parametric and non-parametric methods in statistics. The parametric or model-based approaches assume a known distribution/model (e.g. Gaussian) to obtain inference, whereas non-parametric or model-free approaches assume no underlying distribution/model and obtain inference using empirical distribution. Both model-based and model-free methods have their strength and weakness:

Model-based methods

Model-based methods include DTI, ball-and-sticks model, NODDI as well as more complicated model like CHARM and AxCaliber. Model based methods assume a particular diffusion distribution pattern/function, and the parameters are calculated by fitting diffusion signals with the model. For example, DTI assumes that velocity of water diffusion follows a 3D Gaussian distribution, and the tensor calculated is exactly the covariance matrix of the Gaussian. Ball-and-sticks models is a kind of multiple tensor model, whereas the ball is the isotropic Gaussian, and the sticks is a purely anisotropic Gaussian. The strength of model-based methods, similar to the parametric methods in statistics, is that they only requires few samples to get the whole distribution. However, the results of model-based methods are limited by the model, and it is common that the diffusion pattern does not follows the assumption. A complicated model may also have overfitting problem. For example, past studies have used bi-exponential model to fit intra-cellular and extra-cellular diffusion, but it turned out that they failed to reveal such biophysics.

There are still other methods which are both model-based and model-free. The spherical deconvolution methods (including CSD and the derived approaches) can be viewed as model-based because they assume a underlying distribution for a fiber population (i.e. response function) and used it for deconvolution. It can be also viewed as model-free because they also get an ODF for the fiber distribution (also termed fiber orientation distribution, FOD). The advantage of these methods is that they have both the benefit of model-based and model-free method (i.e. less sampling to get a sharp ODF). The disadvantage is that they also have flaws from both (model violated, sensitive to outliers, prone to create false fibers). CSD can give rise to false crossing geometry if there is a model mismatch (Parker, G. D., Neuroimage 65: 433-448., 2013, Schilling, K., Neuroimage 129: 185-197. 2016) and the tracking results can be noisy.

Model-free methods

Model-free methods estimate the empirical distribution of the water diffusion, and there is no assumption on the distribution. The methods include diffusion spectrum imaging (DSI), q-ball imaging (QBI), and generalized q-sampling imaging (GQI). DSI use Fourier transformation and numerical integration to calculate the orientation distribution function (ODF, which is the empirical distribution of water diffusion at different orientations) of water diffusion. The Fourier transform requires a specific grid  diffusion sampling scheme (multiple b-values multiple directions). QBI uses either Funk-Randon transform or spherical harmonics to calculate the ODF. The calculation requires a shell-like diffusion sampling scheme (i.e. HARDI acquisition, single b-value, multiple directions). GQI provides an analytical relation between diffusion signals and the spin distribution function (SDF) and can be applied to any diffusion sampling scheme. SDF is the density of diffusing water at different direction and is a kind of diffusion ODF. This relation allows GQI to directly compute SDF and thus free from the error in the numerical estimationThe advantage of model-free methods is that they are not limited by a model. They do not assume a particular diffusion structure, and there is no risk for violation of the model. It does not have the overfitting problem in the model-based method. The calculation does not require complicated optimization or fitting and thus is less affected by outliers in comparison with model-based method. The down side of model-free methods is that they often need more diffusion samplings, at least 60 to get a more robust estimation (In comparison,DTI only needs 6 sampling in addition to b0).

If you are using DSI or QBI reconstruction, consider placing them with GQI for the following reasons: DSI and QBI offers only a numerical estimation of the diffusion ODF while GQI offers a direct analytical relation for the diffusion ODF (termed SDF here). In DSI, the numerical estimation includes a Fourier transform, followed by a filter to remove noise and a radial integration with numerical interpolation. In QBI, the numerical estimation includes re-sample diffusion at the radial directions before Funk-Randon transform or use spherical decomposition to estimate spherical harmonics for calculating the transform before converting it back to the ODF. In GQI, since an analytical relation is provided, the calculation is not a numerical estimation but a direct calculation using a simple matrix multiplication to convert diffusion signals to SDF. There is no error due to numerical estimation, and GQI is more accurate in terms of its mathematical form. GQI's analytical relation also bring additional benefit: DSI only works on grid data, and QBI only works on shell data, but GQI can be applied to grid, shell, multi-shell, and non-grid-non-shell data (e.g. body center cubic or face center cubic). GQI can also be coupled with linear and nonlinear transformation to reconstruction data in the T1W space or MNI space. This derived approach is called q-space diffeomorphic reconstruction (QSDR), a method that reconstructs GQI diffusion pattern directly in the MNI space. This makes group comparison and regression studies much easier. QSDR enables template construction, connectome fingerprinting, and connectometry analysis. The SDF calculated from GQI also provides density-based measurements such as QA, ISO, and RDI. It is different from an ODF, which is a probability density function on a unit sphere.

Which method to use?

The choice of the method really depends on the applications as well as the SNR and the sampling scheme of the image data. If the data quality is excellent (good SNR), complicated model-based methods such as ball-and-stick model, NODDI, or CSD can work very well (CSD only works on shell data). If the data are noisy, a better choice would be DTI or GQI. A recent comparison study have shown that the methods that offered the highest validation connections are GQI and DTI tractography (https://www.nature.com/articles/s41467-017-01285-x). This is because the competition data have poor SNR.

My personal method of choice is GQI, QSDR, and DTI. The reason follows my personal preference in statistics: If there is an unknown distribution, the first thing I would do is fitting it with a Gaussian and see its mean and standard deviation. This is the same as applying DTI to diffusion MRI data (What DTI does is fitting data with a 3D Guassian model). A simple Gaussian fitting will be robust against noise and will allow us to see the overall distribution pattern even if there are limited data samples. If I have enough data samples, then I will compute the empirical distribution to observe the overall pattern of the unknown distribution. This is the same as applying GQI to diffusion MRI to see the empirical distribution of the diffusion. The empirical distribution will faithfully reveal the diffusion pattern and not be biased toward the model assumption or give false estimation just due to model mismatch (e.g. fit a Poisson or Beta with a Gaussian). The recent advance in multi-band acquisition allows us to acquire numerous diffusion samplings within a short scanning time. Model-free methods were once criticized by its long acquisition but now the new development in fast imaging will favor non parametric methods and model-free method in the near future.

Data reconstruction

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.

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 built-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.

STEP2: Select a Reconstruction Method

Diffusion Tensor Imaging (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 (in 10-3 mm2/s), and three eigenvalues are also exported.

Ball-and-Sticks Model

The ball-and-sticks model is provided by FSL Diffusion Toolkit (bedpostx) to estimate fiber orientations. At the bottom of this webpage are the Matlab codes (please download bs_recon.m and read_src.m) to reconstruct the DSI Studio src file using this model. This code uses trust-region-reflective to optimize the ball-and-sticks model.

The "fa0" and "fa1" variables from bs_recon.m are volume fractions for the first and second fiber population (note that the ball-stick model gives volume fraction, not "FA"). The volume fraction for the isotropic ball component is thus 1-fa0-fa1.

You may also use FSL's bedpostx to process the diffusion data. The resulting files can be converted to DSI Studio readable format, as detailed in here.

The QBI reconstruction method proposed in the Tuch's MRM paper has two major parameters. One is the width of the Gaussian smoothing kernel, and another is the width of the interpolation kernel. 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 convolution/deconvolution processes.

In the newer version of DSI Studio, Tuch's version of QBI reconstruction is no longer supported since the spherical harmonics (below) can provide much better reconstruction results. Users are recommended to use SH-based QBI or GQI instead.

Q-ball imaging (QBI, Spherical Harmonics)

QBI can also be reconstructed by using spherical harmonics based transformation [3]. There are three different versions proposed by Anderson, Hess et. al. and Descoteaux et.al, respectively. DSI Studio adopted Descoteaux 's method [4], since it is much easier to implement. A regularization parameter is required for the reconstruction. The recommended value is 0.006 [4].

One should note that QBI can be applied to multiple b-value datasets, and to make use of multiple b-value, users are recommended to use GQI.

Diffusion Spectrum Imaging (DSI)

DSI reconstruction method [5] is based on the Fourier transform relation between q-space signals and the diffusion PDF (or the "averaged propagator"). The detailed processing flow is listed 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 discrete 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 q table is n*3.
% Note that your q table loaded here is the same as the mgh_dsi_q.txt.

% 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 with DSI reconstruction in DSI Studio is an index called quantitative index (QA), which is defined in GQI paper [9].

DSI reconstruction does not consider the "image orientation" information in the DICOM header, and the tilted slice may result in incorrect reconstruction. Here is the solution to correct slice tilting for DSI.

1. After loading DICOM images, replace the b-table by the table you used to acquire the DICOM.
2. perform DSI reconstruction to get the fib file.
3. uncompress the .fib.gz file to .fib file. rename it to .mat file and load it in Matlab.
4. use dicominfo command to get slice orientation matrix from one of your DICOM files. The group element number is (0x0020,0x0037)
5. The 6 values in (0x0020,0x0037) are the first and second rows of the orientation matrix, and the third row can be calculated by the cross product of the first and second rows. If the sum of the third row is negative, you may need to multiply -1 to it.
6. Multiple the slice orientation matrix to the odf_vertices matrix (perform reorientation here)
7. save the fib file use save xxx.fib -v4

Limitation and drawbacks
In addition to the re-orientation problem, DSI also suffers from the interpolation and discrete estimation error in radial projection. This error can be better handled by GQI, which provides an analytical relation that converts DWI data to diffusion ODF directly. Furthermore, the effect of the Hanning filter is not fully investigated, whereas, in GQI, it is replaced by a well-defined smoothing parameter called the diffusion sampling distance.

DSI reconstruction is limited to DWI acquired by grid sampling scheme. It cannot be applied to HARDI or multi-shell sampling. This limitation is overcome by GQI.

Implementation difference from TrackVis
DSI Studio uses the minimum non-zero b-value as the reference of the grid distance, and the grid coordinate of other DWI are then calculated by floor(sqrt(b_value/min_b_value)*gradient_direction+0.5).The final grid size is always 16x16x16. TracVis (or more precisely the diffusion Toolkit), however, uses its built-in q vector table and does not use the b-table in the DICOM header.

Generalized Q-sampling Imaging (GQI)

Generalized q-sampling imaging (GQI) is a model-free reconstruction method that quantifies the density of diffusing water at different orientations. This measurement, termed spin distribution function (SDF), is an orientation distribution function of diffusing spins. Studies have shown its greater sensitivity and specificity to white matter characteristics and pathology.

GQI can calculate SDF from a variety of diffusion data sets, including DSI dataset, HARDI, multiple-shell, combined DTI dataset or even body center cubic (BCC) dataset [6][7][9]. GQI provides an analytical relation to compute SDF, and the reconstruction requires only a simple matrix multiplication:

The value 0.018 is the 6D, where D is the self-diffusion coefficient of water at the body temperature (3.00*10^-3 mm^2/s), [Nature Reviews Neuroscience 4, 469-480, 2003 ]. In DSI Studio, the GQI implementation uses 0.015 for 6D, and thus the diffusion length ratio appears longer. (e.g., a ratio 1.25 is in fact 1.145).

The optimal value for "diffusion sampling length ratio" depends on the signal to noise ratio and the b-table, and it is highly recommended to optimize the ratio using a preliminary scan following these steps:

a) Reconstruct FIB files using ratios of 0.3, 0.4, 0.5, ...2.0 etc.
b) Open each file and locate non-crossing regions (e.g. mid corpus callosum) and crossing regions (e.g. lateral corpus callosum).
c) The optimal value should be the largest ratio that resolves crossing patterns in crossing regions and still give single fiber direction at mid corpus callosum.

The following is the Matlab code for GQI reconstruction. The reconstruction codes require users to input the preferred "diffusion sampling length ratio" (e.g. 1.2). This length ratio defines the radius of the diffusion spins included in the ODF estimation. For example, one may need to calculate the number of spins within a displacement radius of 40 microns, but this radius should be considered with respect to the diffusion time (short diffusion time => smaller radius, longer diffusion =>larger radius). A way to consider diffusion time is dividing this radius by the mean diffusion distance of free water diffusion. The resulting is this "diffusion sampling length ratio". Thus, a ratio of 1.2 means that GQI detects diffusion spins within a displacement radius of 1.2 mean diffusion distance of free water.

`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.018);`
`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 available in the attachment file provided at the bottom of this web page.

GQI reconstruction offers quantitative anisotropy (QA) to replace FA. The definition for QA is documented in GQI paper [9].

QA at u = Z(ODF(u)-iso(ODF))

where u is the fiber orientation. iso(ODF) calculates the isotropic component of the ODF. In the implementation, the minimum value of an ODF function is used as the isotropic component. Z is a scaling constant that makes the maximum of all iso(ODF) in the space equal to one. More discussions about this ODF scaling can be found in [9] and [11]. In QSDR reconstruction, Z is estimated using two ventricle voxels that contain free water diffusion.

Normalized QA is calculated by normalizing the maximum QA value to one.

CSF calibration uses spatial normalization to identify the location of CSF and use it as the free water diffusion to unify the amount of diffusion with respect to it.

Restricted Diffusion Imaging (RDI)

RDI[13] is a model-free method to quantify the density of restricted diffusion given a diffusion displacement range (e.g. 10 microns). Its calculation uses different "diffusion sensitization strengths" (or b-values) to separate non-restricted diffusion from restricted diffusion according to their sensitivity to diffusion gradient. This feature allows for selective quantification of restricted diffusion while ignoring non-restricted diffusion (or vice, versa for nRDI), leading to higher sensitivity and specificity to the structural change due to pathological condition.

RDI can be used to detect increased cellularity caused by inflammation or a tumor mass [13]. It can be used to study any pathological condition that leads to a change in restricted diffusion.

RDI also provides its counterpart: non-restricted diffusion imaging (nRDI), which allows for quantifying diffusion with displacement greater than a specified length (e.g. edema, free water). nRDI can be used to visualize edematous tissues.

*To calculate RDI, first select GQI as the reconstruction method and in the advanced options, click on a check box labeled with "restricted diffusion imaging". This will calculate RDI indices along with GQI reconstruction. The indices include rdi02L, rdi04L, ...as well as nrdi04L, nrdi 04L. "rdi02L" quantified the restricted diffusion within "0.2 L", where L is the diffusion distance. "nrdi02L" quantifies non-restricted diffusion with displacement greater than "0.2 L".

The calculation is a simple linear combination of the DWI signals acquired by long diffusion time:

rho(L) is the density of diffusing spins restricted with a displacement distance of L. In DSI Studio, the L value is a ratio with respect to diffusion distance (squared root of 6Dt). For example, rdi02L quantifies the diffusion restricted within 0.2*diffusion distance.

To acquire RDI, the DWI should be ideally acquired using multiple b-values. The optimal setting is using a grid sampling scheme (e.g., DSI scheme) with a long diffusion time. Stimulated echo is the best choice to achieve a long diffusion time without increase the echo time. Nonetheless, this does not mean that RDI cannot be applied to DTI or HARDI, which can be viewed as a partial sampling in the q-space. The only concern is that RDI calculated from DTI, HARDI, or multishell scheme may be biased toward fast or slow diffusion due to parital sampling in the q-space.

Q-Space Diffeomorphic Reconstruction (QSDR)

Q-Space diffeomorphic reconstruction (QSDR) [11] is the generalization of GQI that allows users to construct spin distribution functions (SDFs, a kind of ODF) in any given template space (e.g. MNI space). QSDR can be applied to DTI data, multi-shell data, DSI data, none-shell-none-grid data, or a combination of the above-mentioned data sets.

The reconstruction equation for QSDR is the following:

By reconstructing SDFs in the template space, QSDR provides a direct way to analyze the group difference (see Group connectometry analysis). QSDR can potentially work with any spatial registration method such as SPM normalization, fnirt, ..., etc. However, most normalization methods concerns only point-to-point matching and allow for large deformation that brings extreme values to the Jacobian matrix. Consequently, the fiber orientations can be heavily distorted and not suitable for fiber tracking. An eligible method should enforce smoothness by using a well-behaved basis function or penalize on the extreme values in the Jacobian matrix.

There are two major normalization approaches implemented in QSDR to obtain the diffeomorphic mapping. One set of it is an enhanced version of the SPM-like normalization, including SPM 7-9-7, SPM 14-18-14, and SPM 21-27-21. SPM normalization uses Fourier basis, which has built-in smoothness. The parameter 7-9-7 means that there are 7, 9, and 7 transformation parameters (Fourier basis) used for deforming x, y, and z axis. Higher numbers give better accuracy, but it comes with a dramatically increased demand on computation time and memory. The original SPM routine cannot use a larger number of parameters, but the enhanced version used in DSI Studio breaks this limit to 21-27-21. Another set of normalization, including CDM (constrained diffeomorphic mapping) and T1W-CDM uses a diffeomorphic field with parameter space equal to the image volume. This approach maximizes the intensity correlation and penalizes the second derivative of the diffeomorphic field so that the Jacobian matrix is well-posed for QSDR.

T1W-CDM uses the T1-weighted images to guide the normalization so that QSDR can be applied to a diffusion set with atypical anisotropy map obtained from low b-value or high b-value acquisitions. Please note that the BET-processed T1W cannot be used here because DSI Studio expects original T1W with skull. The DWI should be corrected for gradient distortion using FSL TOPUP so that there will be no mismatch between subject's T1W and DWI. T1W-CDM is recommended if (1) the original T1W is available and (2) DWI is corrected for gradient distortion using FSL TOPUP. If any of the requirement is not meet, CDM should be used instead.

The bounding box information for QSDR output:

``x = -78  to 78``
``y = -112 to 76``
``z = -50  to 85``

To convert the reconstruction coordinate (2mm output) to MNI coordinate:

mni_x =  78.0-2.0*x
mni_y =  76.0-2.0*y
mni_z = -50.0+2.0*z

The 2mm output image has the following sform information:

`srow_x = -2  0  0  76`
`srow_y =  0 -2  0  78`
`srow_z =  0  0  2 -50`

Users can assign arbitrary output resolution (integer) and perform reconstruction without resorting to external normalization tools. DSI Studio first calculates the quantitative anisotropy (QA) mapping in the native space and then normalizes it to the MNI QA map. QSDR also records the R-squared value between the subject QA and MNI QA map in the filename (e.g. .R72.fib.gz means a R-squared value of 0.72). A value greater than 0.6 suggests good registration results, whereas a low value may indicate possible error in the registration. The most common cause for low R2 value is a flipping of the slice order at the Z direction, which causes the brain volume to be placed up side down. This can be corrected using the [Flip z] function provided in the reconstruction window.

You may warp subject T1W/T2W together with QSDR. To add these images, click on the "Add T1WT2W" button at the right side of the output resolution box. If you have T1W-based ROI to be transformed along with the T1W, you may add the T1W first and then the ROI. DSI Studio will prompt you to apply for the previous transform. Once the QSDR reconstruction is done, the added image volume will be stored in the FIB file. You can open the FIB file in DSI Studio and switch the imaging modality from "qa" to whatever images you have added in QSDR reconstruction.

QSDR can work with any template (mice, rodents, macaques...etc.). To add your template, copy the template nifti file to the "template" folder in the DSI Studio package (For Mac version, right click on the app and "show package content" to find the template folder). The nifti file should have a correct pixdim and srow values. If QSDR does not work on your template, please contact Frank.

Diffusion Difference Imaging (DDI)

Diffusion Difference Imaging (DDI) combines quantitative MRI and GQI to calculate difference in diffusion ODF between two dMRI scans in the native space. The reconstruction can look at either increased or decreased diffusion, and the anisotropy value calculated is the anisotropy of the difference in diffusion ODF.

Reconstruction Using MATLAB

Occasionally we may encounter huge diffusion datasets, e.g. a high spatial resolution DSI dataset, and it is impossible for DSI Studio to process the images in the memory since the image volume may have a size of several gigabytes. To reconstruct such images, here I offer the following MATLAB codes that conduct direct GQI reconstruction from a huge image file. This skips the src file creation and directly obtains fib file that can be used to conduct fiber tracking.

The following codes assume that the image pixels are arranged in an order of x - y - z - d, where d represents the serial number of the diffusion weighted image. You may also need to download find_peak.m and odf8.mat (both provided at the bottom of this page) to run this code.

`function gqi_reco(filename,file_type,pixel_size,b_table,mean_diffusion_distance_ratio)`
`% Direct GQI reconstruction from huge image data`
`% You may need to include find_peak.m to run these codes.`
`% parameters:`
`% ` `filename: the filename of the image volume`
`% ` `file_type: the pixel format of the image, can be 'int8', 'int16', 'int32', 'single', or 'double'`
`% ` `pixel_size: the size of the pixel in bytes.  `
`% ` `b_table: the b-table matrix with size of 4-by-d, where d is the number of the diffusion weighted images`
`%          b(1,:) stores the b-value, whereas b(2:4,:) stores the grandient vector`
`% ` `mean_diffusion_distance_ratio: check out GQI reconstruction for detail. Recommended value=1.2  `
`% `
`% example:`
`% ` `gqi_reco('2dseq','int32',4,b_table,1.0);`

`load odf8.mat;`

`% you may need to change the dimension, number of diffusion images, and voxel size`
`dim = [128 128 128];`
`dif = 515;`
`voxel_size = [6/128 6/128 6/128];`

`fa0 = zeros(dim);`
`fa1 = zeros(dim);`
`fa2 = zeros(dim);`
`index0 = zeros(dim);`
`index1 = zeros(dim);`
`index2 = zeros(dim);`

`reco_temp = zeros(dim(1),dim(2),dif);`
`plane_size = dim(1)*dim(2);`

`% GQI reconstruciton matrix A`
`l_values = sqrt(b_table(1,:)*0.01506);`
`b_vector = b_table(2:4,:).*repmat(l_values,3,1);`
`A = sinc(odf_vertices'*b_vector*mean_diffusion_distance_ratio/pi);`

`f =fopen(filename);`
`max_dif = 0;`
`for z = 1:dim(3)`
`    for d = 1:dif`
`        fseek(f,((z-1)*plane_size+(d-1)*dim(1)*dim(2)*dim(3))*pixel_size,'bof');`
`        reco_temp(:,:,d) = reshape(fread(f,plane_size,file_type),dim(1),dim(2));`
`    end`
`    for x = 1:dim(1)`
`        for y = 1:dim(2)`
`            ODF=A*reshape(reco_temp(x,y,:),[],1);`
`            p = find_peak(ODF,odf_faces);`
`            max_dif = max(max_dif,mean(ODF));`
`            min_odf = min(ODF); `
`            fa0(x,y,z) = ODF(p(1))-min_odf;`
`            index0(x,y,z) = p(1)-1;`
`            if length(p) > 1`
`                fa1(x,y,z) = ODF(p(2))-min_odf;`
`                index1(x,y,z) = p(2)-1;`
`            end`
`            if length(p) > 2`
`                fa2(x,y,z) = ODF(p(3))-min_odf;`
`                index2(x,y,z) = p(3)-1;`
`            end`
`        end`
`    end`
`end`
`fa0 = fa0/max_dif;`
`fa1 = fa1/max_dif;`
`fa2 = fa2/max_dif;`
`fa0 = reshape(fa0,1,[]);`
`fa1 = reshape(fa1,1,[]);`
`fa2 = reshape(fa2,1,[]);`
`index0 = reshape(index0,1,[]);`
`index1 = reshape(index1,1,[]);`
`index2 = reshape(index2,1,[]);`
`dimension = dim;`
`save('result.fib','dimension','voxel_size','fa0','fa1','fa2','index0','index1','index2','odf_vertices','odf_faces','-v4');`
`fclose(f);`
`end`

DSI Studio detects fiber directions by searching for local maximum on an ODF. This is implemented in DSI Studio.

The following is the Matlab codes for getting the 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);

Other Functions

Image Rotation

DSI Studio can apply rotation to the source DWI volume. The functions are under the  menu (see figure to the right). The images can be flipped in x, y, or z-direction, and the b-table will be rotated accordingly.

The DWI can be rotated to the T1W or T2W space using [To T1W/T2W space...]. DSI Studio will ask for the T1W or T2W image in the nifti format for the registration. B-table will also be rotated based on the registration matrix.

You can also manually rotate the image volume using [manual rotation]. DSI Studio will provide an interface for manual registration.

[Trim image] removes the background based on the mask you assigned.

Export images and b-table

The images volume in the SRC file can be exported to a 4D nifti file. The function is under the [File] Menu.

ODF Sharpening

Turning on "ODF sharpening" to increase the angular resolution of the resolved fiber. "Deconvolution" conducts diffusion deconvolution to the diffusion ODFs generated from DSI, QBI, or GQI [10]. The regularization parameter controls the smoothness of the generated fiber ODF, and the optimal value may vary from case to case. It may require several trial-and-error to figure out the best setting. "Decomposition"[12] offers a sparse solution of the fODF. The parameter, m, assigns the maximum fiber population for each ODF, and delta assigns the decomposition volume.

Output in the fib file

This option stores extra information to the .fib file. The complete ODF include all ODF information and allows DSI Studio to visualize the ODFs for inspection. Note that enabling this option will greatly increase the size of the .fib file.

Output mapping only applies only to QSDR reconstruction. It output the mapping function (_x, _y, _z) in the fib file.

Output Jacobian determinant applies only to QSDR reconstruction.

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 folds of ODF is recommended for a more detailed presentation of fiber directions.

Number of fibers resolved

This is the maximum number of peaks resolved on an ODF.

Half-sphere scheme

Enable this option if you are using a half sphere DSI scheme. In full sphere DSI scheme, each DWI acquired with (bx, by, bz) has its corresponding DWI acquired by (-bx, -by, -bz), whereas the b0 image does not have a corresponding image. In the half-sphere scheme, we acquire only positive bz (or bx, by, depending on the setting), and to ensure consistent ODF value, the b0 signal should be divided by two because b0 image does not have a corresponding image.

This option will affect QA estimation.

Check b-table

DSI Studio allows for checking whether the b-table is flipped in x, y or z direction. The function is realized by using DSI Studio uses a fiber coherence index (how well the fiber directions connect to each other) to check whether flip x, flip y, and flip z gives a better result.

The checkbox in the reconstruction window enables this checking problem. If the b-table is flipped, the generated fib file will have ".fx", ".fy", or ".fz" added to the file extension. The default setting is having this function enabled. Users can disable this b-table checking by unchecking this checkbox.

Flip or rotate DWI volume before reconstruction

DSI Studio provides functions to flip or rotate DWI data and also corrects the b-table accordingly. To use this function, open the SRC file using "STEP2 Reconstruction". In the main menu on the top, click on  and select the function. DSI Studio will then transform the DWI accordingly and also correct the b-table. You can then choose a reconstruction method and proceed with the processing.

Reconstruct Diffusion Data in T1W or T2W space

DSI Studio provides a function to linearly transform DWI data to the T1W or T2W space (also correct the b-table accordingly). To use this function, open the SRC file using "STEP2 Reconstruction". In the main menu on the top, click on [To T1W/T2W space] and select the T1W or T2W nifti file. DSI Studio will then transform the DWI to the T1W or T2W space and also correct the b-table. You can then choose a reconstruction method and proceed with the processing.

In the registration dialog, users are allowed to use either rigid body transform or affine transform. Rigid body transform assumes that the shape of the targets are identical, and thus there is only translocation and rotation. Affine allows for deformation such as scaling and shearing, which allows for different shapes of targets to be transformed to match each other.

If the T1W you assigned is from the same subject, use rigid body. If the T1W is from a template, use affine.

Data requirement for reconstruction

DTI: For DTI, at least 6 diffusion images with different sampling directions and a b0 image are needed. Images with multiple b-value are also supported.

QBI: Images should be acquired using the same b-value. b0 images will be discarded.

DSI: DSI reconstruction can be applied to "grid" scheme only. The exemplary b-tables for DSI can be found here.

GQI: GQI can be applied to any diffusion sampling schemes, including single-shell (e.g. HARDI), multiple-shell, or grid scheme (DSI).

Decoding the file extension

The FIB file generated during the reconstruction will include several extension. Here is a list of the explanation

odf8: An 8-fold ODF tessellation was used

f5: For each voxel, a maximum of 5 fiber directions were resolved

rec: ODF information was output in the FIB file

csfc: quantitative diffusion MRI was conducted using CSF location as the free water calibration

hs: A half sphere scheme was used to acquired DSI

reg0: the spatial normalization was conducted using (reg0: 7-9-7 reg1: 14-18-14 reg2: 21-27-21)

bal: The diffusion scheme was resampled to ensure balance in the 3D space

fx, fy, fz: The b-table was automatically flipped by DSI Studio in x-, y-, or z- direction. 012 means the order of the x-y-z coordinates is the same, whereas 102 indicates x-y flip, 210 x-z flip, and 021 y-z- flip.

rdi: The restricted diffusioin imaging metrics were calculated

de: deconvolution was used to sharpen the ODF

dec: decomposition was used to sharpen the ODF

gqi: The images were reconstructed using generalized q-sampling imaging

qsdr: The images were reconstructed using q-space diffeomorphic reconstruction

R72: The goodness-of-fit between the subject's data and the template has a R-squared value of 0.72

Diffusion Indices

A common question about diffusion indices is that how do they differ in term of their biophysical meaning? Here is a brief explanation. First of all, all tensor derived measurements, including FA, AD, RD, are based on "diffusivity", which by definition measures how fast water diffuses. By contrast, measurements derived from GQI or QSDR, such as QA, iso, RDI, and local connectome fingerprint are based on "density", which by definition measure how much water diffuses in a particular direction. It is like two different approaches to measure the "traffic". One way is to quantify traffic by measuring how fast vehicles travels on the highway (diffusivity) or count how many vehicles are traveling (density). These measurement may be related, but entirely different.

In application, diffusivity measurements are more sensitive to pathological conditions, whereas density measurements are more sensitive to individual/physiological difference (see Yeh et al. PLoS Comput Biol 12(11): e1005203, where local connectome fingerprint is a type of density measurement). To better understand this difference, we can compare axons to water pipes. If the pipes are in good condition, they will have same water transfusion rate (diffusivity), even though the amount of water (density) being transfused can vary a lot. This indicates that diffusivity is good for detecting whether the structural is still intact, whereas the density measurement is good for quantifying the "connectivity" because it quantifies the total quantity of the diffusing water. This is shown in Yeh et al. PLoS Comput Biol 12(11): e1005203, where the individuality can be only revealed by density-based measurement, not diffusivity measurement.

The following is a brief explanation for the diffusion indices provided by DSI Studio:

FA: fractional anisotropy AD: axial diffusivity RD: radial diffusivity (the average of two radial eigenvalues) MD: mean diffusivity (the average of all eigen values). The diffusivity (either RD, MD, or AD) calculated in DSI Studio has a unit of 10-3 mm2/s

GFA: generalized fractional anisotropy. GFA is calculated from an ODF function. The definition is documented in the q-ball imaging paper [3]. GFA has a high correlation with FA (Fritzsche, K. H., et al. (2010). Neuroimage 51(1): 242-251.), and it also suffers from partial volume effect. The value deceases in fiber crossing or voxels with CSF partial volume. GFA values may not be comparable across different reconstruction methods since the sharpness of ODF may differ due to the reconstruction method used. It may not be comparable if different b-value and diffusion sampling scheme are used. Studies using GFA should be conducted in a careful control manner to make sure that the results are not due to other confounding factors.

QA: quantitative anisotropy (QA). QA is calculated from the peak orientations on a spin distribution function (SDF). Each peak orientation defines a QA value. The definition for QA is defined in the generalized q-sampling imaging paper [9]. One should note that QA is defined for each fiber orientation, whereas FA and GFA are defined for each voxel. This forms a big difference in fiber tracking, since QA can be used to filter out false fibers in crossing fiber scenario. In DSI Studio, qa0 represents the QA value for the most prominent fiber orientations, and qa1 the second, ...etc.

Another difference between QA and FA, GFA is that QA scales with spin density. These features make it less susceptible to partial volume effect (see a comparison in Yeh F-C et al. PLoS ONE 8(11): e80713.2013). However, since QA scales with spin density, it is affected by T2-shine through, receiver gain, and B1 inhomogeneity. The QA value may not be comparable across subject if different TE and diffusion sampling scheme are used.

QA scales with diffusion signals and thus can be affected by the gain of the receiver coil. To ensure a better consistency/reproducibility, DSI Studio uses a voxel with free water diffusion to calibrate QA. The "CSF calibration" option use spatial normalization to search for voxels in the 3rd ventricle and use it to calibrate QA and achieve better reproducibility.

NQA is the normalized QA. Normalizing QA is another solution to ensure QA consistency across subjects. DSI Studio has adopts strategies to calibrate QA (see the above description under the QA subtitle), but sometimes the calibration may not be accurate enough, causing a large variability in QA. A solution to this problem is to scales the maximum QA value of a subject to 1 so that QA may be more comparable across the subject. This assumes that all subjects share the same compactness of the white matter bundle.

ISO: the isotropic value of the ODF. The definition is documented in [9]. ISO is the minimum distribution value of an ODF, and thus it represent background isotropic diffusion.

RDI: an index quantifying the density of restricted diffusion given a displacement distance (L). nRDI quantifies non-restricted diffusion. See the "restricted diffusion imaging" reconstruction above for detail.

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] Yeh, Fang-Cheng, Van Jay Wedeen, and Wen-Yih Isaac Tseng, "Generalized-sampling imaging."Medical Imaging, IEEE Transactions on 29.9 (2010): 1626-1635. [pdf]
[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. (pdf)
[11] Yeh, Fang-Cheng, and Wen-Yih Isaac Tseng, "NTU-90: a high angular resolution brain atlas constructed by q-space diffeomorphic reconstruction." Neuroimage 58.1 (2011): 91-99. [pdf]
[12] Yeh, F.C., Tseng, W.Y.Sparse Solution of Fiber Orientation Distribution Function by Diffusion Decomposition”, PLoS One. 2013 Oct 11;8(10):e75747. doi: 10.1371/journal.pone.0075747. (link).
[13] Yeh, Fang-Cheng, Li Liu, T. Kevin Hitchens, and Yijen L. Wu, "Mapping Immune Cell Infiltration Using Restricted Diffusion MRI", Magn Reson Med. accepted, (2016) [pdf]

ċ
bs_recon.m
(3k)
Fang-Cheng Yeh,
Feb 16, 2015, 9:21 PM
ċ
data.mat
(2445k)
Fang-Cheng Yeh,
Feb 3, 2010, 9:23 PM
ċ
dsi_q_vector_203.txt
(2k)
Fang-Cheng Yeh,
Jul 3, 2012, 12:42 PM
ċ
dsi_q_vector_515.txt
(4k)
Fang-Cheng Yeh,
Jul 3, 2012, 12:42 PM
ċ
find_peak.m
(1k)
Fang-Cheng Yeh,
Jun 4, 2012, 10:33 PM
ċ
gqi_reco.m
(3k)
Fang-Cheng Yeh,
Jun 4, 2012, 10:33 PM
ċ
odf8.mat
(10k)
Fang-Cheng Yeh,
Jun 4, 2012, 10:33 PM
ċ