%% ZCDx model parameter estimation using synthetic data
% © Regents of the University of Minnesota.
% This software is licensed under GPL v3.0.  
% For an alternative license please contact the Office for Technology Commercialization at the University of Minnesota. UMN Reference #20170140.
% How to cite: If you use this software for your work, please use the following citation: 
% Farooq, H. et al. Microstructure Imaging of Crossing (MIX) White Matter Fibers from diffusion MRI. Sci. Rep. 6, 38927; doi: 10.1038/srep38927 (2016).

close all
clear all
clc
%% load scheme file
% Scheme file can be downloaded from http://cmic.cs.ucl.ac.uk/wmmchallenge/
% but here we provide the text file with the toolbox
load('scheme_wmm.txt');  %has to be in STEJSKALTANNER format (G(T/m) Del(sec) del(sec)) 
xdata=scheme_wmm;
D_intra= 1.7e3;          %(mircometer^2/sec for in vivo)
Gamma = 2.675987E8;
%% Load Data 
yd=nifti('ZCDx_synth_data.nii'); 
yd1=yd.dat;
%% Recording Results
 f11 = zeros(10,10,1);
 f21 = zeros(10,10,1);
 f12 = zeros(10,10,1);
 f22 = zeros(10,10,1);
 f13 = zeros(10,10,1);
 f23 = zeros(10,10,1);
 f3 = zeros(10,10,1);
 R1 = zeros(10,10,1);
 R2 = zeros(10,10,1);
 R3 = zeros(10,10,1);
 theta1 = zeros(10,10,1);
 phi1 = zeros(10,10,1);
 dir1=zeros(10,10,1,3);
 theta2 = zeros(10,10,1);
 phi2 = zeros(10,10,1);
 dir2=zeros(10,10,1,3);
 theta3 = zeros(10,10,1);
 phi3 = zeros(10,10,1);
 dir3=zeros(10,10,1,3);
 v1 = zeros(10,10,1);
 v2 = zeros(10,10,1);
 v3 = zeros(10,10,1);
 p= 1;
 tic
for  n= 1:10
    parfor  m= 1:10
                Y2=squeeze(yd1(m,n,p,:));
            
                 [f11e, f21e, f12e, f22e, f13e, f23e, f3e, theta1e, phi1e, R1e, theta2e, phi2e, R2e, theta3e, phi3e, R3e] = zcd_x(xdata,Gamma,D_intra,Y2);
                 As=[R1e theta1e phi1e f11e f21e; R2e theta2e phi2e f12e f22e; R3e theta3e phi3e f13e f23e]';
                 [Y,I]=sort(As(1,:));
                    A=As(:,I)
                 
                    f11(m,n,p)= A(4,3);
                    f21(m,n,p)= A(5,3);
                    f12(m,n,p)= A(4,2);
                    f22(m,n,p)= A(5,2);
                    f13(m,n,p)= A(4,1);
                    f23(m,n,p)= A(5,1);
                    f3(m,n,p)= f3e;
                    theta1(m,n,p)= A(2,3);
                    phi1(m,n,p)= A(3,3);
                    R1(m,n,p)=A(1,3);
                    n1e=[cos(A(3,3))*sin(A(2,3)) sin(A(3,3))*sin(A(2,3)) cos(A(2,3))];
                    dir1(m,n,p,:)= n1e;
                    v1e=(A(4,3)+A(5,3));
                    v1(m,n,p,:)= v1e;
                    theta2(m,n,p)= A(2,2);
                    phi2(m,n,p)= A(3,2);
                    R2(m,n,p)= A(1,2);
                    n2e=[cos(A(3,2))*sin(A(2,2)) sin(A(3,2))*sin(A(2,2)) cos(A(2,2))];
                    dir2(m,n,p,:)= n2e;
                    v2e=(A(4,2)+A(5,2));
                    v2(m,n,p,:)= v2e;
                    theta3(m,n,p)= A(2,1);
                    phi3(m,n,p)= A(3,1);
                    R3(m,n,p)= A(1,1); 
                    n3e=[cos(A(3,1))*sin(A(2,1)) sin(A(3,1))*sin(A(2,1)) cos(A(2,1))];
                    dir3(m,n,p,:)= n3e;
                    v3e=(A(4,1)+A(5,1));
                    v3(m,n,p,:)= v3e;
            
    end
end
toc
%% Saving results               
 save_nii( 'synth_mask.nii', R1, 'R1_ZCDx.nii')
 save_nii( 'synth_mask.nii', R2, 'R2_ZCDx.nii')
 save_nii( 'synth_mask.nii', R3, 'R3_ZCDx.nii')
%% Plot estimated and true radii in all the three directions 
nii = nifti('R1_true.nii');R1_true=nii.dat; R1_true=double(squeeze(R1_true)); 
R1_true = reshape(R1_true,[1 100]); R1= reshape(R1,[1 100]);
figure(1)
    h1=plot(R1_true,R1_true,'b');
    hold on
    h2=scatter(R1_true,R1,'r*');
    title('ZCDx Radius Index Estimates Comparison at SNR = 1000')
    legend('True Values','MIX')
    xlabel('R1'); ylabel('R1 estimated')
    hold off
% Radius along second orientation
nii = nifti('R2_true.nii');R2_true=nii.dat; R2_true=double(squeeze(R2_true)); 
R2_true = reshape(R2_true,[1 100]); R2= reshape(R2,[1 100]);
figure(2)
    h1=plot(R2_true,R2_true,'b');
    hold on
    h2=scatter(R2_true,R2,'r*');
    title('ZCDx Radius Index Estimates Comparison at SNR = 1000')
    legend('True Values','MIX')
    xlabel('R2'); ylabel('R2 estimated')
    hold off
% Radius along third orientation    
nii = nifti('R3_true.nii');R3_true=nii.dat; R3_true=double(squeeze(R3_true)); 
R3_true = reshape(R3_true,[1 100]); R3= reshape(R3,[1 100]);
figure(3)
    h1=plot(R3_true,R3_true,'b');
    hold on
    h2=scatter(R3_true,R3,'r*');
    title('ZCDx Radius Index Estimates Comparison at SNR = 1000')
    legend('True Values','MIX')
    xlabel('R3'); ylabel('R3 estimated')
    hold off
% © 2016 Regents of the University of Minnesota. All rights reserved.