## Introduction to diffusion MRI 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 a 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 the 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 require 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 follow the assumption. A complicated model may also have an overfitting problem. For example, past studies have used a 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 an 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 the 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 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 uses 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 in a 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 the model-based methods. The downside 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 samplings in addition to b0).

If you are using DSI or QBI reconstruction, consider placing them with GQI for the following reasons: DSI and QBI offer 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 transforms 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 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 the 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 a 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 the 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 has 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 Gaussian 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 nonparametric methods and model-free method in the near future.

In terms of performance, CSD was shown to outperformance model-free methods and most of the model-based methods in simulation studies and capillary tube studies. In histology studies, a recent one using histology validation has raised questions to the accuracy of CSD (Schilling, Kurt G., Neuroimage 165 (2018): 200-221). This could be due to its assumption that all axonal shared the same diffusion pattern. GQI was shown to be consistent with histology (Gangolli, Neuroimage 153, 152-167 and Modo, Hum Brain Mapp 37, 780-795.).

## Step T2: 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 T2 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. ### Step T2a: 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 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.

### Step T2a: additional image rotation for small animal scans (optional)

Most animal scans are likely to be acquired at an orientation different from the animal template in DSI Studio. It is recommended that users rotate the image volume to match the template so that the atlas function can be used in the following steps.

To facilitate this additional procedure, DSI Studio provides rotation functions under the top  menu. You can use [Image Flip] or [Image Swap] to make sure your image volume matches that of the template.

The following is the template orientation for mouse/rat scans. Moving the sliding bar to the "right" should go to the "top" of the brain. ### Step T2b(1): Select a Reconstruction Method

#### Diffusion Tensor Imaging (DTI)

The DTI was proposed by Basser et.al. . 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 , 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.

#### 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 . 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 in [Step T3 Fiber Tracking]
c) Switch to the coronal view [View Coronal View] and adjust [Tracking Threshold] option at the right upper corner under [Step T3c: Options][Tracking Parameters] so that the value covers only the white matter: d) Locate non-crossing regions (e.g. mid corpus callosum marked by light red circle) and crossing regions (e.g. lateral corpus callosum marked by light blue). The optimal value should resolves crossing patterns in crossing regions (light blue) and still have clean fiber direction at mid corpus callosum (light red). GQI reconstruction offers quantitative anisotropy (QA) to replace FA. The definition for QA is documented in GQI paper .

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  and . 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 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 . 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)  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). By reconstructing SDFs in the template space, QSDR provides a direct way to analyze the group difference (see Group connectometry analysis). 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. In QSDR, 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 a 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 upside down. This can be corrected using the [Flip z] function provided in the reconstruction window.

The QSDR reconstruction requires the assignment of a template (e.g. human, monkey, rat, mouse). Please make sure that your data match the template. If QSDR does not work on the template, please contact Frank.

TIPS:

1. For animal studies, it is very common that the image orientation is not the same as the template. You will need to check the template orientation (use O4: View images button to open template file under the template folder (For Mac users, the template folder is inside the app package. Right click on dsi_Studio.app to show package content). To change the image orientation, in the reconstruction windows, use [Swap xxx] or [Flip xxx] to make sure that your image orientation match that of the template.
2. For animal studies, you can try [remove background] button in Step T2a and [Trim Images] to reduce the influence of signals outside the brain tissue.
3. You may warp subject other images modalities together with QSDR. To add these images, click on the [Attach Images...] 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.

[ODFs] will export all ODF information and allows DSI Studio to run connectometry analysis or visualize the ODFs for inspection. Note that enabling this option will greatly increase the size of the .fib file.

[Restricted diffusion imaging] will export RDI measures (see RDI above for details)

[Spatial mapping] will export voxel to voxel mapping used in the QSDR reconstruction. It output the mapping function (_x, _y, _z) in the fib file.

[Jacobian determinant] will export Jacobian determinant used in the QSDR reconstruction.

[Tensor matrix] will export the entire tensor matrix.

[DTI measures] will output FA and diffusivities.

[Helix angle] will output helix angle.

#### 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 and correcting whether the b-table is flipped and/or swapped in x, y or z direction. The function is realized by using a fiber coherence index to check which settings a better result .
The checkbox for "check b-table" is located on the bottom of the reconstruction dialog after selecting the reconstruction method (e.g. DTI, GQI, ...etc.). If the b-table is flipped, the generated fib file will have ".fx", ".fy", or ".fz" added to the file extension. If the b-table is rotated, then the dimension order will be added to the file extension (e.g. 012 is the original x-y-z order, and 021 is x-z-y).

*This function only works if the imaging voxel has isotropic resolution. The new version of DSI Studio (after 10/24/2019) will check whether the slice thickness to determine whether this function will be utilized.
*If the diffusion signals are corrupted, this b-table checking may result in a random configuration that can be easily identified from the file name of the FIB file.
*The FSL bval and bvec often results in an "fy" extension.

### 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 FIB 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 swap, 210 x-z swap, and 021 y-z- swap.

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

## Optional: Reconstruction Using MATLAB

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.

The following is MATLAB code that conducts direct GQI reconstruction. This skips the src file creation and directly obtains fib file that can be used to conduct fiber trackingThe code assumes 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`

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

DSI reconstruction method  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
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 .

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.

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);

## Optional: Reconstruct Diffusion Data in T1W or T2W space 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 [Step T2 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.

DSI Studio can also transform DWI data to the T1W or T2W space (also correct the b-table accordingly). To use this function, open the SRC file using [Step T2 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. 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.

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

References

 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.
 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.
 Tuch D.S. Q-ball imaging. Magn. Reson. Med. 2004;52:1358–1372.
 Descoteaux, M., Angelino, E., Fitzgibbons, S., Deriche, R., 2007. Regularized, fast, and robust analytical Q-ball imaging. Magn Reson Med 58, 497-510.
 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
 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.
 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)
 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)
 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]
 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)
 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]
 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).
 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]
 Schilling, Kurt G., et al. "A fiber coherence index for quality control of B-table orientation in diffusion MRI scans." Magnetic resonance imaging (2019).

ċ
Fang-Cheng Yeh,
Feb 16, 2015, 9:21 PM
ċ
data.mat
(2445k)
Fang-Cheng Yeh,
Feb 3, 2010, 9:23 PM
ċ
Fang-Cheng Yeh,
Jul 3, 2012, 12:42 PM
ċ
Fang-Cheng Yeh,
Jul 3, 2012, 12:42 PM
ċ
Fang-Cheng Yeh,
Jun 4, 2012, 10:33 PM
ċ
Fang-Cheng Yeh,
Jun 4, 2012, 10:33 PM
ċ
odf8.mat
(10k)
Fang-Cheng Yeh,
Jun 4, 2012, 10:33 PM
ċ
Fang-Cheng Yeh,
Feb 16, 2015, 9:21 PM