function [Flim,Slim,Status,Sind]=poiss_f_threshold(Back,PSF,Beta,varargin)
% Find the flux threshold for Poisson matched filter
% Package: ImUtil.X
% Description: Find the flux threshold for Poisson matched filter
% Input  : - Image background expectency.
%          - If scalar than this is a Gaussian PSF sigma.
%            Alternatively, this is a matrix or ClassPSF object containing
%            the PSF.
%          - False alarm probability - the probability that the source
%            above this threshold was generated by the noise.
%            Default is 0.01.
%          * Arbitrary number of pairs of arguments: ...,keyword,value,...
%            where keyword are one of the followings:
%            'NsimIm' - Number of simulated images. Default is 1.
%                       If beta is small, this may requires large number of
%                       simulated images.
%            'ConvAccuracy' - Convergence accuracy. Default is 0.03.
%            'ConvMaxIter' - Maximum number of iterations. Default is 10.
%            'ImageSize' - Simulated image size. Default is 2048.
%            'RadSigmaPSF' - Radius of PSF in units of Gaussian sigma
%                          generated when the PSF is a scalar.
%                          Default is 5.
%            'Radius'   - Radius in which 2nd moments are calculated.
%                         Default is 5.
%            'Sigma'    - Sigma of 2nd moment weight. Default is 1.5 pix.
%            'MomMaxIter' - Maximum number of 2nd moment iterations.
%                         Default is 5.
%            'Verbose'  - Default is false.
% Output : - Flux threshold. NaN if not converged.
%          - The limit in units of the matched filter statistics.
%          - true if converged, false if not.
%          - Vector of the simulated values of S from which the probability
%            distribution of S can be calculated.
% License: GNU general public license version 3
%     By : Eran O. Ofek                    Apr 2017
%    URL : http://weizmann.ac.il/home/eofek/matlab/
% Example: [Flim,Slim]=ImUtil.X.poiss_f_threshold(0.005,2,0.001);
% Reliable: 2
%--------------------------------------------------------------------------

Step       = 1;
MinNperbin = 100;

if (nargin<3)
    Beta = 0.01;
end

DefV.NsimIm              = 1;
DefV.ConvAccuracy        = 0.03;
DefV.ConvMaxIter         = 10;
DefV.ImageSize           = 2048;
DefV.RadSigmaPSF         = 5;
DefV.Radius              = 5;
DefV.Sigma               = 1.5;
DefV.MomMaxIter          = 5;
DefV.SubBack             = true;
DefV.Verbose             = false;
InPar = InArg.populate_keyval(DefV,varargin,mfilename);

% check if simulations are meaningful
Npix = InPar.ImageSize.^2.*InPar.NsimIm;
if (min(Beta)<MinNperbin./Npix)
    error('Not enough simulations in order to estimate the flux threshold for such a low beta');
end

Sigma = [];
if (numel(PSF)==1)
    % user provided PSF sigma
    Sigma = PSF;
    P = Kernel2.gauss(Sigma,Sigma,0,ceil(InPar.RadSigmaPSF.*Sigma),ceil(InPar.RadSigmaPSF.*Sigma));
else
    % user provided PSF stamp
    if (ClassPSF.isClassPSF(PSF))
        % PSF is ClassPSF object
        P = getpsf(PSF,[],{});
    else
        % PSF is a matrix
        P = PSF;
    end
end

% make sure the PSF is normalized to unity
P = P./sum(P(:));
        
if (isempty(Sigma))
    SizeP = size(P);
    [~,Mom2]=ImUtil.Im.im_moments(P,SizeP(2).*0.5,SizeP(1).*0.5,InPar.Radius,InPar.Sigma,InPar.MomMaxIter);
    Sigma = sqrt(Mom2.X2 + Mom2.Y2);
end


%FluxNormFun = @(Fth,PSF,Back) sum( PSF(:)./Back .* log(1+Fth.*PSF(:)./Back) );
FluxNormFun = @(Fth,PSF,Back) sum( PSF(:) .* log(1+Fth.*PSF(:)./Back) );
% FluxNormFun = @(Fth,PSF,Back) sum( PSF(:).^2);
  

if (numel(P)<InPar.ImageSize.^2)
    FilterMethod = 'filter2';
else
    FilterMethod = 'fft';
end

% first guess for Flim
Flim = poissinv(1-Beta,4.*pi.*Sigma.*Sigma.*Back);
% make sure Flim is not 0.
Flim = max(Flim,1);

Converge = false;
Iter = 0;
while ~Converge
    Iter = Iter + 1;
    
    Filter = log(1 + Flim.*P./Back);
    
    % simulate multiple images
    for Iim=1:1:InPar.NsimIm
        MB = poissrnd(Back.*ones(InPar.ImageSize,InPar.ImageSize));
        if (InPar.SubBack)
            MB = MB - Back;
        end
        
        switch lower(FilterMethod)
            case 'fft'
                S = ifft2(fft2(MB).*conj(fft2(Filter)));
            case 'filter2'
                S = filter2(Filter,MB,'same');
            otherwise
                error('Unknown FilterMethod option');
        end
        Tmp = S(1:Step:end,1:Step:end);
        Im(Iim).S = Tmp(:)';
    end
    Sind = [Im.S];

    Slim = quantile(Sind(:),1-Beta);
    % Need to be careful here
    % Background was subtracted from image so this will effect
    % Flim
    if (InPar.SubBack)
        Flim1 = Slim./FluxNormFun(Flim,P,Back) + Back;
    else
        Flim1 = Slim./FluxNormFun(Flim,P,Back);
    end
    
    if (Slim==min(Sind(:)))
        % Slim is found in the first "bin"
        % this means that the estimate for Slim is unreiable
        % set to 0
        Slim = 0;
        Flim1 = 0;
    end
    
%     if (Slim<=0)
%         % In this case the 1-Beta is consistent with zero flux
%         % so for this falde alarm probabilty no consistent solution can be
%         % found
%         Slim = 0;
%         Flim1 = 0;
%     else
%         Flim1 = Slim./FluxNormFun(Flim,P,Back);
%     end

    Converge = abs(Flim1 - Flim)<InPar.ConvAccuracy || Iter>InPar.ConvMaxIter;
    Flim  = Flim1;
    
    if (InPar.Verbose)
        fprintf('Iteration %d: Flim=%f\n',Iter,Flim);
    end
    
end

if (Iter>InPar.ConvMaxIter)
    % not converged
    Status = false;
else
    Status = true;
end

if (Flim==0)
    % Flim and Slim are unrelaible
    % set to NaN
    Flim = NaN;
    Slim = NaN;
end

% plots for the paper
% save -v7.3 Sind.mat Sind
% 
% Vec = (0:0.01:20)';
% N=histc(Sind(:),Vec);
% semilogy(Vec,N./sum(N),'k-','LineWidth',1)
% axis([0 15 1e-9 1e-1])
% H=xlabel('S');
% H.FontSize=18;
% H.Interpreter='latex';
% H=ylabel('$p_{S}$ (per 0.01 bin)');
% H.FontSize=18;
% H.Interpreter='latex';
% print p_S_simulation_sigma2_beta0.001_B0.005.eps -depsc2


% semilogy(Vec,1-cumsum(N)./sum(N),'k-','LineWidth',2)
% LogCN = log10(1-cumsum(N)./sum(N));
% PolyPar = polyfit(Vec(LogCN>-Inf & Vec>3), LogCN(LogCN>-Inf & Vec>3),1)
% Val     = polyval(PolyPar,Vec);
% hold on;
% semilogy(Vec(Vec>3),10.^Val(Vec>3),'r-','Color',[0.7 0.7 0.7])
% H=xlabel('S');
% H.FontSize=18;
% H.Interpreter='latex';
% %H=ylabel('$1-\int_{S}^{\infty}{p_{S}dS}$ (per 0.01 bin)');
% H=ylabel('False Alarm Prob. (per 0.01 bin)');
% H.FontSize=18;    
% H.Interpreter='latex';
% legend('Simulations','power-law fit')
% axis([0 17 1e-9 1e0])
% print Cummul_p_S_simulation_sigma2_beta0.001_B0.005.eps -depsc2


