%sinogramFilename='barbie-20121210-eyes-3.txt';

% 360 deg rotation but parallel sinogram data is from 0-180 deg!

%close all

sinogramFilename='/Users/mamo/Documents/KTH/CTrecon/todaystest.txt';
%sinogramFilename='/Users/mamo/Documents/KTH/CTrecon/sethdany.txt';
beamProfileFilename='/Users/mamo/Documents/KTH/CTrecon/beamProfile-Cu-48kV-78uA-20121213.txt';

% iradon reconstruction parameters
filter='Hann';
interpolation='linear';
filterCutoff=1.0;

% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% !!! DO NOT MODIFY LINES BELOW !!!
% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

% load sinogram and beam profile
fanbeamdata=load(sinogramFilename);
beamProfile=load(beamProfileFilename);

%first two numbers are the valid and invalid events in the profile
validEventsInBeamProfile=beamProfile(1);
invalidEventsInBeamProfile=beamProfile(2);
beamProfile=beamProfile(3:end)';

% first row of fanbeam data contains the camera parameters, the other rows
% are the sinogram

parameters=fanbeamdata(1,:);
fanbeamdata=fanbeamdata(2:end,:);

cameraParameters.detectorWidth=parameters(1);
cameraParameters.detectorHeight=parameters(2);
cameraParameters.detectorPosition=[parameters(3) parameters(4) 1];
cameraParameters.XRayTubePosition=[parameters(5) parameters(6) 1];
cameraParameters.centerOfRotation=[parameters(7) parameters(8) 1];
cameraParameters.numberOfSamples=parameters(9);
cameraParameters.acquisitionTime=parameters(10);
cameraParameters.verticalSteps=parameters(11);
cameraParameters.verticalRange=parameters(12);
cameraParameters.rotationSteps=parameters(13);
cameraParameters.extraRotationSteps=parameters(14);
cameraParameters.numberOfDetectorPixels=128;

cameraParameters.centerOfRotation=[0 0 1];

% Convert data to integral of attenuation coefficients using beam profile
% and measured profiles

correctedSinogram=zeros(size(fanbeamdata));
% Ratio between the samples in the beam pofile and measured sinogram
sampleNumberRatio=validEventsInBeamProfile/cameraParameters.numberOfSamples;
 for i=1:size(correctedSinogram,1)
     correctedSinogram(i,:)=-log((squeeze(fanbeamdata(i,:))./beamProfile)*sampleNumberRatio);
 end;
% noise rejection
correctedSinogram(correctedSinogram<0)=0;
correctedSinogram(correctedSinogram==Inf)=0;
% Swap the X-Y axis of the sinogram 90 deg because iradon uses this
% orientation
correctedSinogram=correctedSinogram';

% Reconstruct data slice by slice and put into a volume
% number of projections per slice
numberOfFanbeamProjections=ceil((cameraParameters.extraRotationSteps*2+180.0/cameraParameters.rotationSteps));
%number of slices stored in the file
nSlices=floor(size(correctedSinogram, 2)/numberOfFanbeamProjections);

numberOfParallelbeamProjections=ceil(180.0/cameraParameters.rotationSteps);

% reconstructed 3D volume contains axial slices
reconVolume=zeros(nSlices, cameraParameters.numberOfDetectorPixels, cameraParameters.numberOfDetectorPixels);

theta=[0:cameraParameters.rotationSteps:(numberOfParallelbeamProjections-1)*cameraParameters.rotationSteps];
originalRotationCenter1=cameraParameters.centerOfRotation(1);
for sliceNo=1:nSlices
    sinogram=correctedSinogram(:,(sliceNo-1)*numberOfFanbeamProjections+1:sliceNo*numberOfFanbeamProjections);
    
    % figure;imagesc( sinogram )
    
    parallelSino=genParalellSinogram(sinogram, cameraParameters);
      figure;imagesc( parallelSino )
    % reconstruct the parallel sinogram
    img=iradon(parallelSino, theta, interpolation, filter, filterCutoff, cameraParameters.numberOfDetectorPixels);
    
    % remove statistical error of the corrections in air
    img(img<0)=0;
    
    figure;imagesc( img );
    
    % put the reconstructed slice into a volume
    reconVolume(sliceNo, :,:)=img;
end

