function xs  = pwls_bh2_os_pcm_monochromatic(x, Gb, yi,bi, ri, R, niter, mask)  
 
%   Penalized weighted least squares reconstruction
%	using separable paraboloidal surrogates algorithm with ordered subsets
%	
%	
%       poly-energetic surrogate
%       curvature is second derivative of hi, evaluated at 
%				log (bin/(yi - bin bar{yi}))
%       uses pre-computed denominator (effective energy).
%       new model allows for mixed pixels.
%
%  Input
%		x	        	old estimate
%		Gb	   		    system matrix' 
%		yi		        sinogram (tran. data)
%       bi              Flood image
%		ri              Scatter
%     	R               Regularization variable
%       niter           Number of iteration
%       mask            Field of view mask
%       Ths             Thesholds for segmentation functions
%     
%	Output
%		xs			new estimate
%
% Copyright November 2001,  Idris Elbakri, University of Michigan 
% Updated March 2002.

%[lamLMAtab,xaxisMA, yaxisMA]=deal(lamLdata{:});


[nb, na] = size(yi);


   gi = yi.*(reshape(sum(Gb'), size(yi)));
   denoms  = Gb'*reshape(gi,nb*na,1);


% ------------ penalty term

C1 = R.C1; % instantiate object
Rdercurv = R.dercurv;

nblock = Gb.nblock;  %1;%
[starts, nsubset] = subset_start(nblock);

x=x(mask);
xs = zeros(length(x(:)), niter+1);
xs(:,1) = x(:);

for it=2:niter+1
	disp(sprintf('PWL-OS iteration %d', it-1))

	%
	%	loop over subsets
	%
  for iset=1:nblock
     iblock = starts(iset);
     ia      = iblock:nblock:na;
     hx      = x;


   te = (Gb{iblock}*hx); 

     
        F=reshape(te,size(yi(:,ia)));   
        ei = exp(-F);                                         
        mi = bi .* ei + ri;                                    

        Ni = (yi(:,ia)./mi - 1).*bi .* ei;  

        Njs = reshape(Ni,nb*length(ia),1);

        Nj = Gb{iblock}'*(Njs);

%--------------------     
      %penalty part   
       C1x = C1 * x;
       [pderiv pcurv] = feval(Rdercurv, R, C1x);
       pgrad = C1' * pderiv; 
       pcurv= (C1.^2)'*pcurv;  

%------------------------
     num = -nsubset*Nj - pgrad;   
     den = denoms+pcurv; 

     update=num./den;
     x = x + update ;		



     x = max(x,0);			% enforce nonnegativity
      
    end  %subset 

   xs(:,it) = x;
end
xs = embed(xs,mask);

