function DifT = dti_recon(DTIdata, b0_thr)
% This function will perform DTI calculations on a certain
% DTI dataset.
%
% [DifT]=DTI(DTIdata,parameters)
%
% DTIdata:  A struct containing DTIdata(i).VoxelData, DTIdata(i).Gradient
%           and DTIdata(i).Bvalue of all the DTI datasets
%
% DifT: A 3D matrix with all  Diffusion tensors [Dxx,Dxy,Dxz,Dyy,Dyz,Dzz]
%
% Example,
%  see DTI_example.m
%
% This function is written by D.Kroon University of Twente (August 2008)
    
% Make a 4D matrix to store the different gradient voxel volumes
% (The constant 6 is just a minimum number of volumedatasets,
% this matrix will automaticaly expand if more volumedatasets are used.)
S=zeros([size(DTIdata(1).VoxelData) 6],'single');
% Make a 3D matrix to store the zero gradient voxel volume(s)
S0=zeros(size(DTIdata(1).VoxelData),'single');
% Make a matrix to store the gradients
H=zeros(6,3);
% Make a vector to store the different B-values (timing)
Bvalue=zeros(6,1);

% Read the input data (DTIdata), and seperate the zero gradient
% and other gradients into different matrices.
voxel0=0; voxelg=0;
for i=1:length(DTIdata)
    if(nnz(DTIdata(i).Gradient(:)==[0;0;0])==3)
        voxel0=voxel0+1;
        S0=S0+single(DTIdata(i).VoxelData);
    else
        voxelg=voxelg+1;
        S(:,:,:,voxelg)=single(DTIdata(i).VoxelData);
        H(voxelg,:)=single(DTIdata(i).Gradient);
        Bvalue(voxelg) = single(DTIdata(i).Bvalue);
    end
end
% The zero gradient matrix is the mean of all zero gradient datasets
S0=S0/voxel0;

% Create the b matrices
% (http://www.meteoreservice.com/PDFs/Mattiello97.pdf)
b=zeros([3 3 size(H,1)]);
for i=1:size(H,1),
    b(:,:,i)=Bvalue(i)*H(i,:)'*H(i,:);
end
% Convert measurement intentensity into absorption (log)
Slog=zeros(size(S),'single');
for i=1:size(H,1),
    Slog(:,:,:,i)=log((S(:,:,:,i)./S0)+eps);
end

% Sort all b matrices in to a vector Bv=[Bxx,2*Bxy,2*Bxz,Byy,2*Byz,Bzz];
Bv=squeeze([b(1,1,:),2*b(1,2,:),2*b(1,3,:),b(2,2,:),2*b(2,3,:),b(3,3,:)])';

% Create a matrix to store the Diffusion tensor of every voxel
% [Dxx,Dxy,Dxz,Dyy,Dyz,Dzz]
DifT=zeros([size(S0) 6],'single');

for x=1:size(S0,1),
    for y=1:size(S0,2),
        for z=1:size(S0,3),
            
            % Only process a pixel if it isn't background
            if(S0(x,y,z)>b0_thr)
                
                % Calculate the Diffusion tensor [Dxx,Dxy,Dxz,Dyy,Dyz,Dzz],
                % with a simple matrix inverse.
                Z=-squeeze(Slog(x,y,z,:));
                M=Bv\Z;
               
                DifT(x,y,z,:)=M(:);

            end
        end
    end
end

end