Convert TrackVis, MRtrix, DTI-TK, FSL data to DSI Studio format

Convert MRtrix files to DSI Studio format


Python code for converting mif to fib or fib to mif

     The Python code for converting between MRtrix and DSI Studio can be found here: https://github.com/mattcieslak/dmri_convert 
     

Matlab code for conversion


     MRtrix is a diffusion MRI analysis tool that uses spherical deconvolution to obtain the fiber orientation distribution. DSI Studio can also make use of the data output from MRtrix.

     Note that the format of the b-table file is different in MRtrix. In DSI Studio, the first column is b-value, whereas, in MRtrix, the last column is b-value.

Procedures:

1. Process the diffusion MRI data using MRtrix to get the CSD.mif file. An example is shown as follows:
    
mrconvert hardi160.nii hardi160.mif
dwi2tensor hardi160.mif -grad b_table.txt hardi160.dt.mif
average hardi160.mif -axis 3 - | threshold - - | median3D - - | median3D - mask.mif
tensor2FA hardi160.dt.mif - | mrmult - mask.mif hardi160.fa.mif
erode mask.mif -npass 1 - | mrmult hardi160.fa.mif - - | threshold - -abs 0.38 hardi160.sf.mif
estimate_response.exe hardi160.mif -grad b_table.txt -lmax 6 hardi160.sf.mif     hardi160.rs.txt
csdeconv hardi160.mif -grad b_table.txt hardi160.rs.txt -lmax 8 -mask mask.mif hardi160.CSD6.mif

    Note: you may use DSI Studio to convert src to a 4D nifti file. Load the src file, switch to the "Source images" tab, and click the "save 4d nifti button". You may need to supply the b_table.txt, which has a similar format as the one used in DSI Studio, except that the b_value is at the last column.


2. Use find_SH_peaks to get the fiber orientations

gendir 60 dir60.txt
find_SH_peaks hardi160.CSD6.mif dir60.txt peaks_CSD.mif

3. Use the following function to convert peaks_CSD.mif to peaks_CSD.mif.fib, which can be loaded into DSI Studio. This function makes use of the "read_mrtrix" function, which is provided in the MRtrix package.


function mif2fib(filename)
I = read_mrtrix(filename);
I.data(isnan(I.data)) = 0;
dimension = I.dim(1:3);
voxel_size = I.vox(1:3);
dir0 = I.data(:,:,:,1:3);
dir1 = I.data(:,:,:,4:6);
dir2 = I.data(:,:,:,7:9);
fa0 = sqrt(dir0(:,:,:,1).*dir0(:,:,:,1) + dir0(:,:,:,2).*dir0(:,:,:,2) + dir0(:,:,:,3).*dir0(:,:,:,3));
fa1 = sqrt(dir1(:,:,:,1).*dir1(:,:,:,1) + dir1(:,:,:,2).*dir1(:,:,:,2) + dir1(:,:,:,3).*dir1(:,:,:,3));
fa2 = sqrt(dir2(:,:,:,1).*dir2(:,:,:,1) + dir2(:,:,:,2).*dir2(:,:,:,2) + dir2(:,:,:,3).*dir2(:,:,:,3));

% to get rid of more noisy fibers, increase 0.01
mask0 = fa0 < 0.01*max(reshape(fa0,1,[])) | isnan(fa0); 
mask1 = fa1 < 0.01*max(reshape(fa0,1,[])) | isnan(fa1);  
mask2 = fa2 < 0.01*max(reshape(fa0,1,[])) | isnan(fa2);  
fa0(mask0) = 1.0;
fa1(mask0) = 1.0;
fa2(mask0) = 1.0;
dir0(:,:,:,1) = dir0(:,:,:,1)./fa0;
dir0(:,:,:,2) = dir0(:,:,:,2)./fa0;
dir0(:,:,:,3) = dir0(:,:,:,3)./fa0;
dir1(:,:,:,1) = dir1(:,:,:,1)./fa1;
dir1(:,:,:,2) = dir1(:,:,:,2)./fa1;
dir1(:,:,:,3) = dir1(:,:,:,3)./fa1;
dir2(:,:,:,1) = dir2(:,:,:,1)./fa2;
dir2(:,:,:,2) = dir2(:,:,:,2)./fa2;
dir2(:,:,:,3) = dir2(:,:,:,3)./fa2;
fa0(mask0) = 0.0;
fa1(mask0) = 0.0;
fa2(mask0) = 0.0;
fa0 = reshape(fa0,1,[]);
fa1 = reshape(fa1,1,[]);
fa2 = reshape(fa2,1,[]);
dir0 = permute(dir0,[4 1 2 3]);
dir1 = permute(dir1,[4 1 2 3]);
dir2 = permute(dir2,[4 1 2 3]);
dir0 = reshape(dir0,3,[]);
dir1 = reshape(dir1,3,[]);
dir2 = reshape(dir2,3,[]);
clear I;
clear mask0;
clear mask1;
clear mask2;
save(strcat(filename,'.fib'),'-v4');
gzip(strcat(filename,'.fib'));
delete(strcat(filename,'.fib'));
end

Convert a DTI-TK 4D Tensor Nifti file to DSI Studio format




% You may need to download a nifti file loader to run this function.
% http://www.mathworks.com/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze-image

I = load_nii('tensor.nii.gz');
image_size = size(I.img);
fib.dimension = image_size(1:3);
fib.voxel_size = [2.5 2.5 2.5];
fib.dir0 = zeros([3 fib.dimension]);
fib.fa0 = zeros(fib.dimension);
for z = 1:image_size(3)
z
for y = 1:image_size(2)
for x = 1:image_size(1)
      tensor = zeros(3,3);
      tensor(1,1) = I.img(x,y,z,1,1);
      tensor(1,2) = I.img(x,y,z,1,2);
      tensor(2,1) = I.img(x,y,z,1,2);
      tensor(1,3) = I.img(x,y,z,1,4);
      tensor(3,1) = I.img(x,y,z,1,4);
      tensor(2,2) = I.img(x,y,z,1,3);
      tensor(2,3) = I.img(x,y,z,1,5);
      tensor(3,2) = I.img(x,y,z,1,5);
      tensor(3,3) = I.img(x,y,z,1,6);
      [V D] = eig(tensor);
      if D(3,3) == 0
          continue;
      end
      l1 = D(3,3);
      if(l1 < 0)
          continue;
      end
      l2 = D(2,2);
      l3 = D(1,1);
      if(l2 < 0)
          l2 = 0;
      end
      if(l3 < 0)
          l3 = 0;
      end
      ll = (l1+l2+l3)/3;
      ll1 = l1-ll;
      ll2 = l2-ll;
      ll3 = l3-ll;
      fib.fa0(x,y,z) = sqrt(1.5*(ll1*ll1+ll2*ll2+ll3*ll3)/(l1*l1+l2*l2+l3*l3));
      V(1,3) = -V(1,3); 
      fib.dir0(:,x,y,z) = V(:,3);
end
end
end

% enable the following codes if the image need to flip x and y
% fib.fa0 = fib.fa0(image_size(1):-1:1,image_size(2):-1:1,:);
% fib.dir0 = fib.dir0(:,image_size(1):-1:1,image_size(2):-1:1,:);
% fib.dir0(3,:,:,:) = -fib.dir0(3,:,:,:);

fib.fa0 = reshape(fib.fa0,1,[]);
fib.dir0 = reshape(fib.dir0,3,[]);
save('tensor.fib','-struct','fib','-v4');
gzip('tensor.fib');



Convert FSL bedpostx files to DSI Studio format


The bedpostx routine in FSL uses ball-and-sticks model to fit diffusion data. The output files include mean_th1samples.nii.gz, mean_ph1samples.nii.gz, .. etc. 

Here is the Matlab code for n=2. It converts the FSL  files to a fib file, which can be opened in DSI Studio. 

% You may need to download a nifti file loader to run this function.
% http://www.mathworks.com/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze-image
% This example shows converting a maximum of 2 fiber populations. 
% For populations more than 2, you may need to modify the codes.

theta1 = load_untouch_nii('mean_th1samples.nii.gz');
theta2 = load_untouch_nii('mean_th2samples.nii.gz');
phi1 = load_untouch_nii('mean_ph1samples.nii.gz');
phi2 = load_untouch_nii('mean_ph2samples.nii.gz');
f1 = load_untouch_nii('mean_f1samples.nii.gz');
f2 = load_untouch_nii('mean_f2samples.nii.gz');
fib.dimension = size(f1.img);
fib.voxel_size = f1.hdr.dime.pixdim(2:4);
fib.fa0 = double(reshape(f1.img,1,[]));
fib.fa1 = double(reshape(f2.img,1,[]));
fib.dir0(1,:,:,:) = sin(theta1.img).*cos(phi1.img);
fib.dir0(2,:,:,:) = sin(theta1.img).*sin(phi1.img);
fib.dir0(3,:,:,:) = cos(theta1.img);
fib.dir1(1,:,:,:) = sin(theta2.img).*cos(phi2.img);
fib.dir1(2,:,:,:) = sin(theta2.img).*sin(phi2.img);
fib.dir1(3,:,:,:) = cos(theta2.img);
fib.dir0 = double(reshape(fib.dir0,3,[]));
fib.dir1 = double(reshape(fib.dir1,3,[]));

% flip xy: you may need to make sure that this orientation is correct

fib.fa0 = reshape(fib.fa0,fib.dimension);
fib.fa0 = fib.fa0(fib.dimension(1):-1:1,fib.dimension(2):-1:1,:);
fib.fa0 = reshape(fib.fa0,1,[]);

fib.fa1 = reshape(fib.fa1,fib.dimension);
fib.fa1 = fib.fa1(fib.dimension(1):-1:1,fib.dimension(2):-1:1,:);
fib.fa1 = reshape(fib.fa1,1,[]);


fib.dir0 = reshape(fib.dir0,[3 fib.dimension]);
fib.dir0 = fib.dir0(:,fib.dimension(1):-1:1,fib.dimension(2):-1:1,:);
fib.dir0(3,:,:,:) = -fib.dir0(3,:,:,:);
fib.dir0 = reshape(fib.dir0,3,[]);

fib.dir1 = reshape(fib.dir1,[3 fib.dimension]);
fib.dir1 = fib.dir1(:,fib.dimension(1):-1:1,fib.dimension(2):-1:1,:);
fib.dir1(3,:,:,:) = -fib.dir1(3,:,:,:);
fib.dir1 = reshape(fib.dir1,3,[]);

save('out.fib', '-struct','fib','-v4');

Here is the code for n=3

% You may need to download a nifti file loader to run this function.
% http://www.mathworks.com/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze-image
% This example shows converting a maximum of 2 fiber populations. 
% For populations more than 2, you may need to modify the codes.theta1 = load_untouch_nii('mean_th1samples.nii.gz');
theta1 = load_untouch_nii('mean_th1samples.nii.gz');
theta2 = load_untouch_nii('mean_th2samples.nii.gz');
theta3 = load_untouch_nii('mean_th3samples.nii.gz');
phi1 = load_untouch_nii('mean_ph1samples.nii.gz');
phi2 = load_untouch_nii('mean_ph2samples.nii.gz');
phi3 = load_untouch_nii('mean_ph3samples.nii.gz');
f1 = load_untouch_nii('mean_f1samples.nii.gz');
f2 = load_untouch_nii('mean_f2samples.nii.gz');
f3 = load_untouch_nii('mean_f3samples.nii.gz');
fib.dimension = size(f1.img);
fib.voxel_size = f1.hdr.dime.pixdim(2:4);
fib.fa0 = double(reshape(f1.img,1,[]));
fib.fa1 = double(reshape(f2.img,1,[]));
fib.fa2 = double(reshape(f3.img,1,[]));
fib.dir0(1,:,:,:) = sin(theta1.img).*cos(phi1.img);
fib.dir0(2,:,:,:) = sin(theta1.img).*sin(phi1.img);
fib.dir0(3,:,:,:) = cos(theta1.img);
fib.dir1(1,:,:,:) = sin(theta2.img).*cos(phi2.img);
fib.dir1(2,:,:,:) = sin(theta2.img).*sin(phi2.img);
fib.dir1(3,:,:,:) = cos(theta2.img);
fib.dir2(1,:,:,:) = sin(theta3.img).*cos(phi3.img);
fib.dir2(2,:,:,:) = sin(theta3.img).*sin(phi3.img);
fib.dir2(3,:,:,:) = cos(theta3.img);
fib.dir0 = double(reshape(fib.dir0,3,[]));
fib.dir1 = double(reshape(fib.dir1,3,[]));
fib.dir2 = double(reshape(fib.dir2,3,[]));

% flip xy: you may need to make sure that this orientation is correct

fib.fa0 = reshape(fib.fa0,fib.dimension);
fib.fa0 = fib.fa0(fib.dimension(1):-1:1,fib.dimension(2):-1:1,:);
fib.fa0 = reshape(fib.fa0,1,[]);

fib.fa1 = reshape(fib.fa1,fib.dimension);
fib.fa1 = fib.fa1(fib.dimension(1):-1:1,fib.dimension(2):-1:1,:);
fib.fa1 = reshape(fib.fa1,1,[]);

fib.fa2 = reshape(fib.fa2,fib.dimension);
fib.fa2 = fib.fa2(fib.dimension(1):-1:1,fib.dimension(2):-1:1,:);
fib.fa2 = reshape(fib.fa2,1,[]);


fib.dir0 = reshape(fib.dir0,[3 fib.dimension]);
fib.dir0 = fib.dir0(:,fib.dimension(1):-1:1,fib.dimension(2):-1:1,:);
fib.dir0(3,:,:,:) = -fib.dir0(3,:,:,:);
fib.dir0 = reshape(fib.dir0,3,[]);

fib.dir1 = reshape(fib.dir1,[3 fib.dimension]);
fib.dir1 = fib.dir1(:,fib.dimension(1):-1:1,fib.dimension(2):-1:1,:);
fib.dir1(3,:,:,:) = -fib.dir1(3,:,:,:);
fib.dir1 = reshape(fib.dir1,3,[]);

fib.dir2 = reshape(fib.dir2,[3 fib.dimension]);
fib.dir2 = fib.dir2(:,fib.dimension(1):-1:1,fib.dimension(2):-1:1,:);
fib.dir2(3,:,:,:) = -fib.dir2(3,:,:,:);
fib.dir2 = reshape(fib.dir2,3,[]);

save('out.fib', '-struct','fib','-v4');


Convert FSL qboot files to DSI Studio format


The qboot routine in FSL outputs files such as mean_f1samples.nii.gz, mean_f2samples.nii.gz, .. etc. 

Here is the Matlab code for converting these files to the fib file, which can be opened in DSI Studio. 


% You may need to download a nifti file loader to run this function.
% http://www.mathworks.com/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze-image
% This example shows converting a maximum of 2 fiber populations. 
% For populations more than 2, you may need to modify the codes.
theta1 = load_untouch_nii('merged_th1samples.nii.gz');
theta2 = load_untouch_nii('merged_th2samples.nii.gz');
theta3 = load_untouch_nii('merged_th3samples.nii.gz');
theta1.img = mean(theta1.img,4);
theta2.img = mean(theta2.img,4);
theta3.img = mean(theta3.img,4);
phi1 = load_untouch_nii('merged_ph1samples.nii.gz');
phi2 = load_untouch_nii('merged_ph2samples.nii.gz');
phi3 = load_untouch_nii('merged_ph3samples.nii.gz');
phi1.img = mean(phi1.img,4);
phi2.img = mean(phi2.img,4);
phi3.img = mean(phi3.img,4);
f1 = load_untouch_nii('mean_f1samples.nii.gz');
f2 = load_untouch_nii('mean_f2samples.nii.gz');
f3 = load_untouch_nii('mean_f3samples.nii.gz');
fib.dimension = size(f1.img);
fib.voxel_size = f1.hdr.dime.pixdim(2:4);
fib.fa0 = double(reshape(f1.img,1,[]));
fib.fa1 = double(reshape(f2.img,1,[]));
fib.fa2 = double(reshape(f3.img,1,[]));
fib.dir0(1,:,:,:) = sin(theta1.img).*cos(phi1.img);
fib.dir0(2,:,:,:) = sin(theta1.img).*sin(phi1.img);
fib.dir0(3,:,:,:) = cos(theta1.img);
fib.dir1(1,:,:,:) = sin(theta2.img).*cos(phi2.img);
fib.dir1(2,:,:,:) = sin(theta2.img).*sin(phi2.img);
fib.dir1(3,:,:,:) = cos(theta2.img);
fib.dir2(1,:,:,:) = sin(theta3.img).*cos(phi3.img);
fib.dir2(2,:,:,:) = sin(theta3.img).*sin(phi3.img);
fib.dir2(3,:,:,:) = cos(theta3.img);
fib.dir0 = double(reshape(fib.dir0,3,[]));
fib.dir1 = double(reshape(fib.dir1,3,[]));
fib.dir2 = double(reshape(fib.dir2,3,[]));

% flip xy: you may need to make sure that this orientation is correct

fib.fa0 = reshape(fib.fa0,fib.dimension);
fib.fa0 = fib.fa0(fib.dimension(1):-1:1,fib.dimension(2):-1:1,:);
fib.fa0 = reshape(fib.fa0,1,[]);

fib.fa1 = reshape(fib.fa1,fib.dimension);
fib.fa1 = fib.fa1(fib.dimension(1):-1:1,fib.dimension(2):-1:1,:);
fib.fa1 = reshape(fib.fa1,1,[]);

fib.fa2 = reshape(fib.fa2,fib.dimension);
fib.fa2 = fib.fa2(fib.dimension(1):-1:1,fib.dimension(2):-1:1,:);
fib.fa2 = reshape(fib.fa2,1,[]);


fib.dir0 = reshape(fib.dir0,[3 fib.dimension]);
fib.dir0 = fib.dir0(:,fib.dimension(1):-1:1,fib.dimension(2):-1:1,:);
fib.dir0(3,:,:,:) = -fib.dir0(3,:,:,:);
fib.dir0 = reshape(fib.dir0,3,[]);

fib.dir1 = reshape(fib.dir1,[3 fib.dimension]);
fib.dir1 = fib.dir1(:,fib.dimension(1):-1:1,fib.dimension(2):-1:1,:);
fib.dir1(3,:,:,:) = -fib.dir1(3,:,:,:);
fib.dir1 = reshape(fib.dir1,3,[]);

fib.dir2 = reshape(fib.dir2,[3 fib.dimension]);
fib.dir2 = fib.dir2(:,fib.dimension(1):-1:1,fib.dimension(2):-1:1,:);
fib.dir2(3,:,:,:) = -fib.dir2(3,:,:,:);
fib.dir2 = reshape(fib.dir2,3,[]);

save('out.fib', '-struct','fib','-v4');

Convert DSI Studio src files to FSL 4D nifti


The bedpostx routine in FSL takes 4d nifti files as the input data, and DSI Studio can convert its src files to this format. This is done by first opening the src file using the "STEP2 : Reconstruction" button. In the reconstruction window, switch to the first tab to review the diffusion MRI data. Click on a button labeled "save 4d nifti" to export diffusion data to 4d nifti format.

FSL uses a different format of b-table. You may get this table from the src file by loading it in Matlab (check out here to know how to open src files in Matlab). You may copy the b values and save them to a file named "bvals", b vectors to a file named "bvecs"

To process the data in bedpostx, you need to create a folder and put in 4D nifti files, bvals, bvecs, and a brain mask. The brain mask can be created from the b0 image (use MRIcron to save the first DWI of the 4D nifti).


Flip images in X and Y direction


The data from different tool may have a different coordinate system, and they may require correction to present correct visualization. You may use the following steps to flip the image volume in a FIB file.

(1) gunzip the fib file and rename it as a *.mat file
(2) load the fib file
(3) reshape all fa0,fa1, and dir0,dir1,.. matrices using the following codes. You may apply reshape to all fa and dir matrix

fa0 = reshape(fa0,dimension);
fa1 = reshape(fa1,dimension);
...
dir0 = reshape(dir0,[3 dimension]);
dir1 = reshape(dir1,[3 dimension]);
...

(4) flip all fa0,fa1, and dir0, dir1,...etc matrices. The following codes flip matrix fa0 and dir0. You may need to add codes to flip other matrices.

fib.fa0 = fib.fa0(image_size(1):-1:1,image_size(2):-1:1,:);
fib.dir0 = fib.dir0(:,image_size(1):-1:1,image_size(2):-1:1,:);
fib.dir0(3,:,:,:) = -fib.dir0(3,:,:,:);
% flip x y in dir0 is equal to flip z because (-x,-y,z) is the same as (x,y,-z)

(5) reshape fa and dir matrices back

fa0 = reshape(fa0,1,[]);
fa1 = reshape(fa1,1,[]);
...
dir0 = reshape(dir0,3,[]);
dir1 = reshape(dir1,3,[]);
...

(5) Save the fib file using  

save you_file.fib -v4

Comments