%% ZCDx model parameter estimation with MIX using any 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 xdata/ scheme file 
mask = nifti('scheme.nii');    % has to be in STEJSKALTANNER format 
xdata=mask.dat;
D_intra = 1.7e3;              %(mircometer^2/sec for in vivo human)
D_iso= 3.0e3;
Gamma = 2.675987E8;
g=[xdata(:,1) xdata(:,2) xdata(:,3)];
G_mag= (xdata(:,4)./1000000); %(Tesla/micro meter)
Del=xdata(:,5);               %(in seconds)
del=xdata(:,6);               %(in seconds)
td=(Del-(del/3));             %(in seconds)
b1=(td).*((Gamma.*del.*G_mag).^2); %(sec/mircometer^2)
%% load mask file
mask = nifti('mask.nii');
M=mask.dat;
 %% load data
yd= nifti('data.nii');         
yd1= yd.dat;
%% Recording Results
 f11 =zeros(size(M,1),size(M,2),size(M,3));
 f21 = zeros(size(M,1),size(M,2),size(M,3));
 f12 = zeros(size(M,1),size(M,2),size(M,3));
 f22 = zeros(size(M,1),size(M,2),size(M,3));
 f13 = zeros(size(M,1),size(M,2),size(M,3));
 f23 = zeros(size(M,1),size(M,2),size(M,3));
 f3 = zeros(size(M,1),size(M,2),size(M,3));
 R1 = zeros(size(M,1),size(M,2),size(M,3));
 R2 = zeros(size(M,1),size(M,2),size(M,3));
 R3 = zeros(size(M,1),size(M,2),size(M,3));
 theta1 = zeros(size(M,1),size(M,2),size(M,3));
 phi1 = zeros(size(M,1),size(M,2),size(M,3));
 dir1=zeros(size(M,1),size(M,2),size(M,3),3);
 theta2 = zeros(size(M,1),size(M,2),size(M,3));
 phi2 = zeros(size(M,1),size(M,2),size(M,3));
 dir2=zeros(size(M,1),size(M,2),size(M,3),3);
 theta3 = zeros(size(M,1),size(M,2),size(M,3));
 phi3 = zeros(size(M,1),size(M,2),size(M,3));
 dir3=zeros(size(M,1),size(M,2),size(M,3),3);
 v1 = zeros(size(M,1),size(M,2),size(M,3));
 v2 = zeros(size(M,1),size(M,2),size(M,3));
 v3 = zeros(size(M,1),size(M,2),size(M,3));
%% Recording Results
for p=  1:size(M,3)
    fitting_slice_z=p
    for n= 1:size(M,2)   
        parfor  m= 1:size(M,1)
             if M(m,n,p)>0
                  ydata1=squeeze(yd1(m,n,p,:))
            
                 [f11e, f21e, f12e, f22e, f13e, f23e, f3e, theta1e, phi1e, R1e, theta2e, phi2e, R2e, theta3e, phi3e, R3e] = zcd_x_gen(xdata,Gamma,D_intra,ydata1);
                 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
    end
end
 %% Save output nifti      
 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')
% © 2016 Regents of the University of Minnesota. All rights reserved.