function so = pix2pix_red(si)
%SAXS/pix2pix_red Reduces SAXS data on a pixel-to-pixel basis.
%   pix2pix_red(S) returns a SAXS image object S with image
%   pixels reduced according to the instructions resident in the SAXS
%   object. 
%     
%   In other words, pix2pix_red returns the reduced image.
%

if nargin ~= 1
    error('The routine pix2pix_red must have a SAXS object as argument')
end

stemp=si;
%if beamstop mask is not desired...define it as 1 everywhere
if str2num(si.reduction_state(12))==0 % if the mask is not to be used
    stemp.bsmask_pixels=1;
end
%if mask is not desired...define it as 1 everywhere
if str2num(si.reduction_state(8))==0 % if the mask is not to be used
    stemp.mask_pixels=1;
end

%if flatfield  is not desired...define it as 1 everywhere
if str2num(si.reduction_state(6))==0 % if the flatfield is not to be used
    stemp.flatfield_pixels=1;
    flatfield=stemp.flatfield_pixels;
else
    maxflatfield=ffvalue(si,3); % maxflatfield is actually not the maximum but rather an "average" defined in ffvalue.m
    %3 means to calculate it as average of log values
    flatfield=stemp.flatfield_pixels/maxflatfield; 
end
stemp.pixels=stemp.pixels./flatfield.*stemp.mask_pixels.*stemp.bsmask_pixels;

warning off MATLAB:divideByZero
%if darkcurrent  is not desired...define it as 0 everywhere
if str2num(si.reduction_state(3))==0 % if the darkcurrent is not to be used
    stemp.darkcurrent_pixels=0;
    stemp.darkcurrent_livetime=stemp.livetime;
else
    stemp.darkcurrent_pixels=(si.darkcurrent_pixels.*stemp.mask_pixels)./flatfield;
    stemp.darkcurrent_pixels=prepare_pixels(stemp.darkcurrent_pixels,stemp); %maybe do readout_subtraction and zinger_removal
end


%if empty is not desired...define it as 0 everywhere
if str2num(si.reduction_state(4))==0 % if the empty is not to be used
    empty_corr=0;
    stemp.empty_pixels=0;
    stemp.empty_livetime=stemp.livetime;
    stemp.empty_transfact=1;
else
    empty_corr=1;
    stemp.empty_pixels=(si.empty_pixels.*stemp.mask_pixels)./flatfield;
    stemp.empty_pixels=prepare_pixels(stemp.empty_pixels,stemp); %maybe do readout_subtraction and zinger_removal
end

%if sample trans is not desired...define it as 1
if str2num(si.reduction_state(5))==0 % if the sample trans is not to be used
    stemp.transfact=1;
end

%if sample thickness is not desired...define it as 1 everywhere

if str2num(si.reduction_state(1))==0 % if sample thickness is not to be used
    stemp.sample_thickness=1;
end
        
if str2num(si.reduction_state(2))==0 % if absolute intensity scaling factor is not to be used
    stemp.abs_int_fact=1;
end

warning on MATLAB:divideByZero

% now we are ready look at some data.
pixels_red=(stemp.pixels-(stemp.livetime/stemp.darkcurrent_livetime)*stemp.darkcurrent_pixels)/stemp.transfact...
    -empty_corr*(stemp.empty_pixels-(stemp.empty_livetime/stemp.darkcurrent_livetime)*stemp.darkcurrent_pixels)/stemp.empty_transfact*...
    (stemp.livetime/stemp.empty_livetime);
pixels_red=pixels_red*stemp.abs_int_fact/stemp.sample_thickness./stemp.mask_pixels;

so=si;
so.pixels=pixels_red;

function pixelsout=prepare_pixels(pixelsin,si)
% this routine makes sure that the data to be operated on is corrected for
% zingers
% prior to later application
if str2num(si.reduction_state(9))==1    %this means that readoutnoise should be subtracted
    pixelsin=pixelsin-si.readoutnoise;
end
if str2num(si.reduction_state(10))==1    %this means that data should be "zinger removed"
    pixelsin=zinger_removal(pixelsin,si.error_norm_fact);
end
pixelsout=pixelsin;





