%% ActiveAx ex-vivo model parameters estimation with MIX
% © 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 protocol file
mask = nifti('scheme.nii');    % has to be in STEJSKALTANNER format 
xdata=mask.dat;
D_intra= 0.6e3;             %(mircometer^2/sec fixed for ex vivo, for invivo use 1.7e3)
D_iso= 2.0e3;               %(mircometer^2/sec fixed for ex vivo, for invivo use 3e3)
Gamma = 2.675987E8;
%% load mask file
mask = nifti('mask.nii');
A = mask.dat; 
%% Load synthetic data 
ydata  = nifti('data.nii'); 
Y = ydata.dat; 
%% Recording Results
R = zeros(size(A,1),size(A,2),size(A,3));
f1 = zeros(size(A,1),size(A,2),size(A,3));
f2 = zeros(size(A,1),size(A,2),size(A,3));
f3 = zeros(size(A,1),size(A,2),size(A,3));
f4 = zeros(size(A,1),size(A,2),size(A,3));
dir = zeros(size(A,1),size(A,2),size(A,3),3);
tic
for p = 1:size(A,3)
     for n = 1:size(A,2)
         parfor m = 1:size(A,1)
               if A(m,n,p)>0
                  Y2=double(squeeze(Y(m,n,p,:)));
                  [f1e, f2e, f3e, f4e, thetae, phie, Re] = Aax_exvivo_gen(xdata,Gamma,D_intra,D_iso,Y2);
                  R(m,n,p)= Re;
                  f1(m,n,p)= f1e; 
                  f2(m,n,p)= f2e;
                  f3(m,n,p)= f3e; 
                  f4(m,n,p)= f4e;
                  f12e=f1e+f2e;  ve=f1e/f12e; DIe=ve/(pi*(Re.^2))
                  DI(m,n,p)= DIe;
                  ne=[cos(phie)*sin(thetae) sin(phie)*sin(thetae) cos(thetae)];
                  dir(m,n,p,:)= ne;
               else
                  R(m,n,p)= 0;
                  f1(m,n,p)= 0; 
                  f2(m,n,p)= 0;
                  f3(m,n,p)= 0; 
                  f4(m,n,p)= 0;
                  DI(m,n,p)= 0;
                  dir(m,n,p,:)= zeros(1,3);
               end
        end
    end
end
toc
%% Saving results
% DI(DI>0.05)=0.05;                  
 save_nii( 'synth_mask.nii', R, 'R_MIX_Aax.nii')
 save_nii( 'synth_mask.nii', DI, 'DI_MIX_Aax.nii')
 save_nii( 'synth_mask.nii', f1, 'f1_MIX_Aax.nii')
 save_nii( 'synth_mask.nii', f2, 'f2_MIX_Aax.nii')
 save_nii( 'synth_mask.nii', f3, 'f3_MIX_Aax.nii')
 save_nii( 'synth_mask.nii', f4, 'f4_MIX_Aax.nii')
 save_nii( 'synth_mask.nii', dir, 'dir_MIX_Aax.nii')
% © 2016 Regents of the University of Minnesota. All rights reserved.