close all

%xRaydataFilename='barbiex-xray-37kV113uA.txt';
xRaydataFilename='sethdany.txt';
beamProfileFilename='beamProfile-Cu-48kV-78uA-20121213.txt';


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

% load sinogram and beam profile
fanbeamdata=load(xRaydataFilename);
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);


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

correctedImg=zeros(size(fanbeamdata));
% Ratio between the samples in the beam pofile and measured sinogram
sampleNumberRatio=validEventsInBeamProfile/cameraParameters.numberOfSamples;
for i=1:size(correctedImg,1)
    correctedImg(i,:)=-log((squeeze(fanbeamdata(i,:))./beamProfile)*sampleNumberRatio);
end;
% noise rejection
correctedImg(correctedImg<0)=0;
correctedImg(correctedImg==Inf)=0;

figure;imagesc( correctedImg );
% pixle size of the detector in mm
detectorPixelSize=0.4;
daspect([cameraParameters.verticalSteps/detectorPixelSize 1 1]);
