%% NODDIx 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)
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');
A=mask.dat;
 %% load data
yd= nifti('data.nii');         
yd1= yd.dat;
%% Recording Results
f11 =zeros(size(A,1),size(A,2),size(A,3));
f21 = zeros(size(A,1),size(A,2),size(A,3));
f12 = zeros(size(A,1),size(A,2),size(A,3));
f22 = zeros(size(A,1),size(A,2),size(A,3));
f3 = zeros(size(A,1),size(A,2),size(A,3));
theta1 = zeros(size(A,1),size(A,2),size(A,3));
phi1 =zeros(size(A,1),size(A,2),size(A,3));
OD1 = zeros(size(A,1),size(A,2),size(A,3));
theta2 = zeros(size(A,1),size(A,2),size(A,3));
phi2 = zeros(size(A,1),size(A,2),size(A,3));
OD2 = zeros(size(A,1),size(A,2),size(A,3));
dir1 = zeros(size(A,1),size(A,2),size(A,3),3);
dir2 = zeros(size(A,1),size(A,2),size(A,3),3);
for p=  1:size(A,3)
    fitting_slice_z=p
    for n= 1:size(A,2)   
        parfor  m= 1:size(A,1)
             if A(m,n,p)>0
                  ydata1=squeeze(yd1(m,n,p,:));

                    [f11e, f21e, f12e, f22e, f3e, OD1e, theta1e, phi1e, OD2e, theta2e, phi2e] = nod_2O(xdata,Gamma,D_intra,D_iso,ydata1)

                 if (OD1e)>(OD2e)
                    f11(m,n,p)= f11e;
                    f21(m,n,p)= f21e;
                    f12(m,n,p)= f12e;
                    f22(m,n,p)= f22e;
                    f3(m,n,p)= f3e;
                    theta1(m,n,p)= theta1e;
                    phi1(m,n,p)= phi1e;
                    OD1(m,n,p)= OD1e;
                    theta2(m,n,p)= theta2e;
                    phi2(m,n,p)= phi2e;
                    OD2(m,n,p)= OD2e;  
                    n1e=[cos(phi1e)*sin(theta1e) sin(phi1e)*sin(theta1e) cos(theta1e)];
                    dir1(m,n,p,:)= n1e;
                    n2e=[cos(phi2e)*sin(theta2e) sin(phi2e)*sin(theta2e) cos(theta2e)];
                    dir2(m,n,p,:)= n2e;
                 else
                    f11(m,n,p)= f12e;
                    f21(m,n,p)= f22e;
                    f12(m,n,p)= f11e;
                    f22(m,n,p)= f21e;
                    f3(m,n,p)= f3e;
                    theta1(m,n,p)= theta2e;
                    phi1(m,n,p)= phi2e;
                    OD1(m,n,p)= OD2e;
                    theta2(m,n,p)= theta1e;
                    phi2(m,n,p)= phi1e;
                    OD2(m,n,p)= OD1e; 
                    n1e=[cos(phi1e)*sin(theta1e) sin(phi1e)*sin(theta1e) cos(theta1e)];
                    dir2(m,n,p,:)= n1e;
                    n2e=[cos(phi2e)*sin(theta2e) sin(phi2e)*sin(theta2e) cos(theta2e)];
                    dir1(m,n,p,:)= n2e;
                 end
             else
                 f11(m,n,p)= 0
                 f21(m,n,p)= 0;
                 f12(m,n,p)= 0;
                 f22(m,n,p)= 0;
                 f3(m,n,p)= 0;
                 theta1(m,n,p)= 0;
                 phi1(m,n,p)= 0;
                 OD1(m,n,p)= 0;
                 theta2(m,n,p)= 0;
                 phi2(m,n,p)= 0;
                 OD2(m,n,p)= 0; 
                 dir1(m,n,p,:)= zeros(3,1);
                 dir2(m,n,p,:)= zeros(3,1);
            end
        end
    end
end
 %% Save output nifti
save_nii( 'HCP_brain_mask', OD1, 'OD1_MIX.nii')
save_nii( 'HCP_brain_mask', OD2, 'OD2_MIX.nii')
save_nii( 'HCP_brain_mask', f11, 'v_ic1_MIX.nii')
save_nii( 'HCP_brain_mask', f12, 'v_ic2_MIX.nii')
save_nii( 'HCP_brain_mask', f3, 'v_iso_MIX.nii')
% © 2016 Regents of the University of Minnesota. All rights reserved.