%% ActiveAx ex-vivo model parameter estimation with MIX, using exvivo data, used in the original ActveAx study
% © 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 CC
% download from http://camino.cs.ucl.ac.uk/index.php?n=Tutorials.ActiveAx
mask = nifti('MidSagCC.hdr');
A=mask.dat;      % mask file has to be of the same size as data file             
%% load data 
% read instructions at http://camino.cs.ucl.ac.uk/index.php?n=Tutorials.ActiveAx 
% or at https://github.com/daducci/AMICO/tree/master/matlab/doc/demos/ActiveAx
fid= fopen('All.Bfloat','r','b');               
ydata_camino=fread(fid,'float');
fclose(fid);
Y = reshape(ydata_camino,[372 128 256 3]);
%% 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=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_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.