function [] = gated_fdk( volsize, voxsize, f_projdata, f_matrices, f_phases, f_output, gating )
%GATED_FDK CAVAREV 1.0 (CArdiac VAsculature Reconstruction EValuation)
% example reconstruction in Matlab.
%
%   This is an example Matlab function which is based on the data of the CAVAREV platform.
%	It implements an ECG-gated FDK reconstruction as described on the website.
%	It requires you to work with the post-processed projection data, the projection matrices
%	and the heart phase data.
%	
%	The related data can be downloaded on our website: www.cavarev.com!
%
%	Contact:
%	Christopher Rohkohl
%	Chair of Pattern Recognition, Department of Computer Science
%	Friedrich-Alexander University Erlangen-Nuremberg
%	Martensstr. 3, 91058 Erlangen, Germany.
%
%	Email: contact@cavarev.com
%
% USAGE:
% GATED_FDK(VOLSIZE, VOXSIZE, F_PROJDATA, F_MATRICES, F_PHASES, F_OUTPUT, GATING)
%   performs an ECG-gated FDK reconstruction. The parameters are as
%   follows:
%       VOLSIZE         3-D vector containing the requested reconstruction size
%       VOXSIZE         Voxel size in [mm]
%       F_PROJDATA      Path to the post-processed projection data file
%       F_MATRICES      Path to the projetion matrices file
%       F_PHASES        Path to the heart phase data file
%       F_OUTPUT        Path to where the result is stored
%       GATING          2-D vector containing the gating configuration.
%                       First element is the target reconstruction heart phase.
%                       Second element is the gating function window width
%
%   The output of the reconstruction F_OUTPUT can be read in via the
%   function CAVAREV_READVOLUME.


% fixed data properties for CAVAREV
N = 133;    % number of projection images
Sx = 960;   % projection image width
Sy = 960;   % projection image height

% read in the projection matrices
A = readProjectionMatrices(f_matrices, N);

% read in the heart phases
H = readHeartPhases(f_phases, N);

% open handle for projection images
fid = fopen(f_projdata, 'rb');

% allocate memory for the reconstruction
VOL = zeros(volsize);

% Determine the world coordinates for the first volume voxel.
% It is set such that the iso-center (0,0,0) in world coordinates
% is located in the center of the volume
ORIG = -voxsize.*0.5.*(volsize - 1);

% Compute the world coordinates for each voxel of the grid
[WX WY WZ] = meshgrid( ORIG(1)+(0:(volsize(1)-1)).*voxsize, ORIG(2)+(0:(volsize(2)-1)).*voxsize, ORIG(3)+(0:(volsize(3)-1)).*voxsize  );

% iterate over the projection images
for i=1:N
    
    disp(sprintf('Processing image %d of %d', i, N));
    
   % read the post-processed projection image
   IPP = fread(fid, [Sx Sy], 'single')';
   
   % get projection matrix
   A_i = A(:,:,i);
   
   % get gating weighting value
   lambda = getGatingWeight(H(i), gating);

   % skip the projection if it has a zero-weight
   if (lambda == 0)
      continue; 
   end
   
   % weight the projection image by the gating function value
   IPP = IPP .* lambda;
   
   % project each of the voxels onto the projection image
   UVW = A_i * [WX(:) WY(:) WZ(:) ones(numel(WX),1)]';

   % interpolate linearly in the projection image and apply distance weighting
   IMVALS = interp2(IPP, UVW(1,:)./UVW(3,:)+1, UVW(2,:)./UVW(3,:)+1, 'linear', 0) ./ (UVW(3,:).^2);
   
   % reshape the values to the volume format
   IMVALS = reshape(IMVALS, size(WX));
   
   % update the volume
   VOL = VOL + IMVALS;
end

% close projection file handle
fclose(fid);

% scale to 8-bit using min/max scaling
minval = min(VOL(:));
maxval = max(VOL(:));
VOL = round(255 .* (VOL-minval) ./ (maxval-minval));

% save the result
fid = fopen(f_output, 'wb');

%%%%%%%%% save the volume in the CAVAREV 1.0 format. All data is expected in little-endian format.
% 3*float						3*4 bytes				first voxel x=(x0,x1,x2) origin in world coordinates
% 3*unsigned int				3*4 bytes				reconstruction size (Sx0, Sx1, Sx2) in voxels
% 1*float						1*4 bytes				voxel size in mm
% Sx0*Sx1*Sx2*unsigned char     Sx0*Sx1*Sx2*1 bytes		reconstructed volume in row-major format
fwrite(fid, ORIG, 'single');
fwrite(fid, volsize, 'uint32');
fwrite(fid, voxsize, 'single');
for i=1:size(VOL,3),
   % the slice needs to be transposed for row-major format
   S = VOL(:,:,i)';
   fwrite(fid, S, 'uint8');
end
fclose(fid);

end


function [A] = readProjectionMatrices(f_matrices, N)
% read in the projection matrices
    fid = fopen(f_matrices, 'rb');
    
    A = zeros(3,4, N);
    for i=1:N
       A(:,:,i) = fread(fid, [4 3], 'single')'; 
    end
    
    fclose(fid);
end

function [H] = readHeartPhases(f_phases, N)
% read in the heart phase file
% if the file is invalid, all heart phases are defaulted to 0
    fid = fopen(f_phases, 'rt');
    
    if (fid <= 0)
        H = zeros(N, 1);
        return;
    end
    
    H = fscanf(fid, '%f', N);
    
    if (numel(H)~= N)
        H = zeros(N, 1);
    end
    
    fclose(fid);
end

function [lambda] = getGatingWeight(phase, gprops)
% determine the gating function value for a specific heart phase
    d = min(abs((phase-gprops(1)) + [0 -1 +1]));
    
    if (d >= 0.5*gprops(2))
       lambda = 0;
       return;
    end
    
    lambda = cos(d / gprops(2) * pi);
    lambda = lambda * lambda;
end
