function [ fbpC   ] = f_correctionJS( A,B,yi,ftab1,G,f)
%	
%	Joseph and Spital method
%	
%   Input
%       A,B: Parameters for beam-hardening correction to approximate the 2D
%       beam-hardening function by using the 1D function correspoding to water
%		ftab	    Beam-hardening function of the water
%		yi	        Polychromatic sinogram
%       bi          Flood image
%       G           System Matrix
%		ri          Scatter
%       f           Struct for the filtered backprojection
%
%   Output
%       fbpC        Corrected filtered backprojection
%
% Jeff Fessler, University of Michigan
% Cristobal Martinez and Monica Abella, Universidad Carlos III de Madrid
% Updated April 2019.

% FBPU
fbpu = fbp2(yi, f.fbp, 'window', 'hanning,0.8');


% Thresholds for the segmentation functions
Th0=0.14;
ThST=0.2;
Th=0.3;
ThBone=0.5;

xx=[Th0-0.05 Th0 ThST ThST+0.05];
yy=[0 0 1 1];
[xData, yData] = prepareCurveData(xx,yy);

% Set up fittype and options.
ft = fittype( 'poly3' );
opts = fitoptions( 'Method', 'LinearLeastSquares' );
opts.Robust = 'LAR';

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
coeff=coeffvalues(fitresult);
pST=[coeff(1) coeff(2) coeff(3) coeff(4)];
pSTdot=[3*coeff(1) 2*coeff(2) coeff(3)];
pSTddot=[6*coeff(1) 2*coeff(2)];

xx=[Th-0.05 Th ThBone ThBone+0.05];
yy=[1 1 0 0];
[xData, yData] = prepareCurveData( xx, yy );

% Set up fittype and options.
ft = fittype( 'poly3' );
opts = fitoptions( 'Method', 'LinearLeastSquares' );
opts.Robust = 'LAR';

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
coeff=coeffvalues(fitresult);
p=[coeff(1) coeff(2) coeff(3) coeff(4)];
pdot=[3*coeff(1) 2*coeff(2) coeff(3)];
pddot=[6*coeff(1) 2*coeff(2)];



% Bone segmentation
fwj= polyval_STinterv(p,pST,fbpu,Th0,ThST,Th,ThBone);
fwj2=polyval_Boneinterv(p,fbpu,Th, ThBone);
fwj2= min(1-fwj2,1);
hxt = fbpu.*(fwj2);
hxt(hxt>0)=1;

% Projection of the bone segmentation
ti = (G*hxt);        

% Correction factor
correctedProj=(ftab1.inv1.fun(yi))- A*ti+B*ti.^2;

% Filtered Backprojection
fbpC = fbp2(correctedProj, f.fbp, 'window', 'hanning,0.8');



end

