%% ActiveAx ex-vivo 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 protocol file
load('ActiveAx_xdata.txt'); % has to be in STEJSKALTANNER format 
xdata=ActiveAx_xdata;
D_intra= 0.6e3;             %(mircometer^2/sec fixed for ex vivo)
D_iso= 2.0e3;               %(mircometer^2/sec fixed for ex vivo)
Gamma = 2.675987E8;
%% load mask file
mask = nifti('synth_mask.nii');
A = mask.dat; 
%% Load synthetic data 
ydata  = nifti('Aax_synth_data.nii'); 
Y = ydata.dat; 
%% Recording Results
DI = zeros(size(A,1),size(A,2),size(A,3));
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 = 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(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_synth_Aax.nii')
 save_nii( 'synth_mask.nii', DI, 'DI_MIX_synth_Aax.nii')
% save_nii( 'synth_mask.nii', f1, 'f1_MIX_synth_Aax.nii')
% save_nii( 'synth_mask.nii', f2, 'f2_MIX_synth_Aax.nii')
% save_nii( 'synth_mask.nii', f3, 'f3_MIX_synth_Aax.nii')
% save_nii( 'synth_mask.nii', f4, 'f4_MIX_synth_Aax.nii')
% save_nii( 'synth_mask.nii', dir, 'dir_MIX_synth_Aax.nii')
%% load true data 
nii = nifti('R_true.nii');R_true=nii.dat; R_true=double(squeeze(R_true)); 
nii = nifti('DI_true.nii');DI_true=nii.dat; DI_true=double(squeeze(DI_true)); 
DI_true = reshape(DI_true,[1 100]); DI_MIX= reshape(DI,[1 100]);
R_true = reshape(R_true,[1 100]); R_MIX= reshape(R,[1 100]);
%% Plot comparison
figure(1)
    h1=plot(R_true,R_true,'b');
    hold on
    h2=scatter(R_true,R_MIX,'r*');
    title('ActiveAx Radius Index Estimates Comparison at SNR = 1000')
    legend('True Values','MIX')
    xlabel('R'); ylabel('R estimated')
    hold off
    
figure(2)
    h1=plot(DI_true,DI_true,'b');
    hold on
    h2=scatter(DI_true,DI_MIX,'r*');
    title('ActiveAx Density Index Estimates Comparison at SNR = 1000')
    legend('True Values','MIX')
    xlabel('DI'); ylabel('DI estimated')
    hold off 
% © 2016 Regents of the University of Minnesota. All rights reserved.