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');
end Convert a DTI-TK 4D Tensor Nifti file to DSI Studio format
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 formatThe 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
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 formatThe 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,[]));
save('out.fib', '-struct','fib','-v4'); Convert DSI Studio src files to FSL 4D niftiThe 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 directionThe 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) 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 |