%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [FEFA, FA, FE, SE, eigval,Cl,Cp,Cs,MD] = View3DDTMaps(DT, mask)

h = waitbar(0,'Memory allocation...');

FE = repmat(double(0), [size(DT,1) size(DT,2) size(DT,3) 3]);
SE = repmat(double(0), [size(DT,1) size(DT,2) size(DT,3) 3]);
eigval = repmat(single(0), [size(DT,1) size(DT,2) size(DT,3) 3]);
Cl= repmat(single(0), [size(DT,1) size(DT,2) size(DT,3) ]);
Cp= repmat(single(0), [size(DT,1) size(DT,2) size(DT,3) ]);
Cs= repmat(single(0), [size(DT,1) size(DT,2) size(DT,3) ]);
MD= repmat(single(0), [size(DT,1) size(DT,2) size(DT,3) ]);
FEFA = repmat(double(0), [size(DT,1) size(DT,2) size(DT,3) 3]);
FA = repmat(single(0), [size(DT,1) size(DT,2) size(DT,3)]);
close(h);pause(0.05)


h = waitbar(0,'Calculating FEFA, FA, MD, and FE maps...');
for i=1:size(DT,1)
    for j=1:size(DT,2)
        for k=1:size(DT,3)
            
            if mask(i,j,k)
                [V, E] = AbsoluteEigenvalueDT(double(DT(i,j,k,:)));
                
                FE(i,j,k,:) = double(V(:,1));
                SE(i,j,k,:) = double(V(:,2));
                lan1=single(E(1,1));
                lan2=single(E(2,2));
                lan3=single(E(3,3));
                eigval(i,j,k,1) = lan1;
                eigval(i,j,k,2) = lan2;
                eigval(i,j,k,3) = lan3;
                MD(i,j,k)=1/3*(lan1+lan2+lan3);
                Cl(i,j,k) = (lan1 - lan2 )/( lan1+lan2+lan3 );
                Cp(i,j,k) = 2*(lan2 - lan3 )/( lan1+lan2+lan3 );
                Cs(i,j,k) = 3*lan3/( lan1+lan2+lan3 );
                
            end
            
        end
    end
    waitbar(i/size(DT,1));
end
close(h);pause(0.05)


h = waitbar(0,'Further processing of FEFA, FA, and FE maps...');
% uint(0)
warning off MATLAB:divideByZero

for k=1:size(DT,3)
    
    M = mask(:,:,k);
    fm = find(M);
    
    L1 = double(eigval(:,:,k,1));
    L2 = double(eigval(:,:,k,2));
    L3 = double(eigval(:,:,k,3));
    FE1 = abs(double(FE(:,:,k,1)));
    FE2 = abs(double(FE(:,:,k,2)));
    FE3 = abs(double(FE(:,:,k,3)));
    
    
    fa = repmat(single(0), [size(DT,1) size(DT,2)]);
    fa(fm) = sqrt(3*((L1(fm)-L2(fm)).^2 + (L2(fm) - L3(fm)).^2 +...
        (L3(fm)-L1(fm)).^2))./sqrt(2*(L1(fm).^2 + L2(fm).^2 + L3(fm).^2));
    fa(fm(find(isnan(fa(fm))))) = 0;
    FA(:,:,k) = single(fa);
    
    FE1(fm) = double((fa(fm).*FE1(fm))./sqrt(3));
    FE2(fm) = double((fa(fm).*FE2(fm))./sqrt(3));
    FE3(fm) = double((fa(fm).*FE3(fm))./sqrt(3));
    
    FEFA(:,:,k,1) = FE1;
    FEFA(:,:,k,2) = FE2;
    FEFA(:,:,k,3) = FE3;
    
    waitbar(k/size(DT,3));
    
end
close(h);pause(0.05)

warning on MATLAB:divideByZero
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [VV, EE] = AbsoluteEigenvalueDT(DT)

try
    [V, E] = eig(reshape(DT,3,3));
catch
    disp('nan or inf encountered...')
    V = [1 0 0; 0 1 0; 0 0 1];
    E = [eps 0 0; 0 eps 0; 0 0 eps];
end
E = abs(E);

eigvals = [E(1,1) E(2,2) E(3,3)];
[seigvals,I] = sort(eigvals);

VV = zeros(3,3);
EE = zeros(3,3);

EE(1,1) = seigvals(3);
EE(2,2) = seigvals(2);
EE(3,3) = seigvals(1);

VV(:,1) = V(:,I(3));
VV(:,2) = V(:,I(2));
VV(:,3) = V(:,I(1));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [X, chi_sq, flag, O_M] = RESTORE_func(dwi,b)

% flag = 1 --> unphysical DWIs (b0<DWIs)
% flag = 2 --> negative DWIs
% flag = 3 --> negative DWIs that where also flagged 1
% flag = 4 --> outliers detected
% flag = 5 --> outliers detected that where also flagged 1
% flag = 6 --> outliers detected that where also flagged 2
% flag = 7 --> outliers detected that where also flagged 3


dwi = double(dwi);

MaxNoIter = 10;
RelPrec = 10^-2;

O_M = (2:length(dwi))';

flag = 0;

if any(dwi(2:end)>dwi(1))
    flag = 1;
    if mean(double(any(dwi(2:end)>dwi(1))))>0.1
        dwi(1) = max(dwi(:))+eps;
    else
        dwi(dwi>dwi(1)) = dwi(1)-eps;
    end
end

if any(dwi<0)
    if flag==0
        flag = 2;
    else
        flag = 3;
    end
end

dwi(dwi<=0) = 1;

A = log(dwi(1)./dwi(2:end));
% w = (dwi(2:end).^2)./(stdev_DWI(2:end).^2);
w = ones(size(dwi(2:end)));
X = lscov(b,A,w);

fit = dwi(1)*exp(-b*X);
resid = fit - dwi(2:end);

err_weighted = (resid.^2).*w;
chi_sq = sum(err_weighted(:));

mean_err = mean(resid);
std_err = std(resid);

% figure;
% plot(dwi(2:end),'.b');
% hold on;
% plot(fit,'.r');
% hold off;

logi = or(resid<(mean_err-3*std_err),resid>(mean_err+3*std_err));

if any(logi)
    if flag==0
        flag=4;
    elseif flag==1
        flag=5;
    elseif flag==2
        flag=6;
    elseif flag==3
        flag=7;
    end
end

weight = w;

while any(logi)
    
    dwi([logical(0);logi])=[];
    b(logi,:)=[];
    w(logi)=[];
    A(logi)=[];
    weight(logi)=[];
    O_M(logi)=[];
    
    
    X = lscov(b,A,weight);
    
    fit = dwi(1)*exp(-b*X);
    resid = fit - dwi(2:end);
    
    err_weighted = (resid.^2).*weight;
    chi_sq = sum(err_weighted(:));
    
    cn = 0;
    prec = 1;
    
    while cn<MaxNoIter && prec > RelPrec
        
        cn = cn+1;
        chi_sq_T = chi_sq;
        X = lscov(b,A,weight);
        
        fit = dwi(1)*exp(-b*X);
        resid = fit - dwi(2:end);
        m_err = median(resid);
        C = 1.4826*median(abs(resid-m_err));
        if C<eps
            chi_sq=single(nan);
            return;
        else
            weight = 1./(C^2+resid.^2);
            weight = weight./mean(weight(:));
            err_weighted = (resid.^2).*weight;
            chi_sq = sum(err_weighted(:));
        end
        
        prec = abs(chi_sq_T-chi_sq)/chi_sq_T;
        
    end
    
    mean_err = mean(resid);
    std_err = std(resid);
    
    logi = or(resid<(mean_err-3*std_err),resid>(mean_err+3*std_err));
    
end

X = lscov(b,A,w);

fit = dwi(1)*exp(-b*X);
resid = fit - dwi(2:end);

err_weighted = (resid.^2).*w;
chi_sq = sum(err_weighted(:));

end



