function [FiltImage,Slim,Res,Src]=filter_energy_dep(EvtTable,varargin)
% Filter a Poisson-noise image with energy dependent PSF and background.
% Package: ImUtil.X
% Description: Filter a Poisson-noise image with energy dependent PSF
%              background, and spectrum.
% Input  : - A clean event table - e.g., [X, Y, Energy] per photon.
%          * Arbitrary number of pairs of arguments: ...,keyword,value,...
%            where keyword are one of the followings:
%            'ColX' - Column index of X coordinate.
%            'ColY' - Column index of Y coordinate.
%            'ColE' - Column index of energy.
%            'Beta' - False alarm probability - the probability that the
%                     source above this threshold was generated by the
%                     noise. Default is 1e-4.
%            'EnergyEdges' - Vector of energy edges defines the bins in
%                     which to apply the filters. Each energy bin have its
%                     own filter.
%                     Default is [200, 500, 1000, 2000, 4000, 6000, 10000]
%            'PSF'  - PSF per energy bin. This is either a scalar (same PSF
%                     for all energy bins); Vector of sigmas of Gaussian
%                     PSF per energy bin; a cell array of PSFs; a ClassPSF
%                     object with number of elements equal to the number of
%                     energy bins.
%                     Default is [2    2    2     2     2     2].
%            'Spec' - Vector of spectrum (intensity) per energy bin. This
%                     will be used for optimal filtering.
%                     Default is [1    1    1     1     1     1].
%            'Back' - Vector of background levels (events per pixel) in
%                     each energy bin.
%                     Default is [0.01 0.01 0.01  0.01  0.01  0.01].
%            'CCDSEC'- CCDSEC region to filter [Xmin Xmax Ymin Ymax].
%                     Default is [0 1024 0 1024].
%            'PoissonFilter' - Use Poisson-noise filter.
%                     Default is true.
%            'KernelSize' - PSF size (if generated). Default is 15.
%            'PSFPar' - Cell array of additional arguments to pass to
%                     ClassPSF/getmpsf. Default is {}.
%            'PixSize' - Image pixel size. Default is 1.
%            'FilterFun' - Filter function: @filter2 | @filter2_fft |
%                     @filter2_fftc. Default is @filter2.
%            'ThreshPar' - Cell array of additional arguments to pass to
%                     ImUtil.Im.poiss_f_threshold. Default is {}.
%            'Conn' - imregionalmax connectivity parameter. Default is 4.
% Output : - Filtered image statistics (S).
%          - The threshold in S corresponding to the requested Beta.
%          - Structure of additional information.
%          - Structure of local maxima (sources) found in filtered image.
%            Coordinates are relative to the CCDSEC corner.
% License: GNU general public license version 3
%     By : Eran O. Ofek                    May 2017
%    URL : http://weizmann.ac.il/home/eofek/matlab/
% Example:
% [FiltImage,Slim,Res,Src]=ImUtil.X.filter_energy_dep([rand(1e4,2).*1024,rand(1e4,1).*1e4]);
% Reliable: 2
%--------------------------------------------------------------------------


error('problem with S_E summation!');


DefV.ColX                = 1;
DefV.ColY                = 2;
DefV.ColE                = 3;
DefV.Beta                = 1e-4;
DefV.EnergyEdges         = [200, 500, 1000, 2000, 4000, 6000, 10000];
DefV.PSF                 = [2    2    2     2     2     2];
DefV.Spec                = [1    1    1     1     1     1];
DefV.Back                = [0.01 0.01 0.01  0.01  0.01  0.01];
DefV.CCDSEC              = [0 1024 0 1024];
DefV.PoissonFilter       = true;
DefV.KernelSize          = 15;
DefV.PSFPar              = {}; % additional arguments to ClassPSF/getmpsf
DefV.PixSize             = 1;
DefV.FilterFun           = @filter2;  % @filter2 | @filter2_fft | @filter2_fftc
DefV.ThreshPar           = {};
DefV.Conn                = 4;

InPar = InArg.populate_keyval(DefV,varargin,mfilename);

% select events by position
FlagPos = EvtTable(:,InPar.ColX)>InPar.CCDSEC(1) & ...
          EvtTable(:,InPar.ColX)<InPar.CCDSEC(2) & ...
          EvtTable(:,InPar.ColY)>InPar.CCDSEC(3) & ...
          EvtTable(:,InPar.ColY)<InPar.CCDSEC(4);
      
EvtTable = EvtTable(FlagPos,:);

NenegyBin = numel(InPar.EnergyEdges) - 1;

% prepare PSFs
if (isnumeric(InPar.PSF) && numel(InPar.PSF)==NenegyBin)
    % InPar.PSF is a vector of Gaussian sigmas per energy bin
    CellPSF = cell(NenegyBin,1);
    for Ieb=1:1:NenegyBin
        CellPSF{Ieb} = Kernel2.gauss(InPar.PSF(Ieb),InPar.PSF(Ieb),0,InPar.KernelSize,InPar.KernelSize);
    end
    
elseif (isnumeric(InPar.PSF) && numel(InPar.PSF)==1)
    % InPar.PSF is a scalar of Gaussian sigmas for all energy bins
    CellPSF = cell(NenegyBin,1);
    CellPSF = Kernel2.gauss(InPar.PSF(1),InPar.PSF(1),0,InPar.KernelSize,InPar.KernelSize);
    for Ieb=2:1:NenegyBin
        CellPSF{Ieb} = CellPSF{1};
    end
    
elseif (ClassPSF.isClassPSF(InPar.PSF))
    % InPar.PSF is a ClassPSF array
    
elseif (iscell(InPar.PSF))
    % InPar.PSF is a cell array of PSFs
    CellPSF   = InPar.PSF;
    [CellPSF] = getmpsf(InPar.PSF,InPar.PSFPar{:});
    
else
    error('Unknown InPar.PSF option');
end

% Prepare spectrum (weights)
if (numel(InPar.Spec)==NenegyBin)
    SpecW = InPar.Spec./sum(InPar.Spec);
else
    error('Illegal format for InPar.Spec'); 
end

% prepare pixel edges
XE = (InPar.CCDSEC(1):InPar.PixSize:InPar.CCDSEC(2));
YE = (InPar.CCDSEC(3):InPar.PixSize:InPar.CCDSEC(4));


% Filter for each energy bin
for Ieb=1:1:NenegyBin
    % Energy edges for current bin
    E1 = InPar.EnergyEdges(Ieb);
    E2 = InPar.EnergyEdges(Ieb+1);
 
    % select events by energy
    FlagE = EvtTable(:,InPar.ColE)>=E1 & EvtTable(:,InPar.ColE)<E2;
    ET    = EvtTable(FlagE,:);
    
    % create 2-D image by 2D histogram
    Image = histcounts2(ET(:,InPar.ColX),ET(:,InPar.ColY),XE,YE);
    
    if (Ieb==1)
        ImageSize = size(Image);
        FiltCube = zeros(NenegyBin,ImageSize(1),ImageSize(2));
    end
    
    % Detection threshold for specific energy bin
    [Res.Flim(Ieb),Res.Slim(Ieb),Res.Status(Ieb),Sval(Ieb).S]=ImUtil.Im.poiss_f_threshold(InPar.Back(Ieb),CellPSF{Ieb},InPar.Beta,InPar.ThreshPar{:});
    % Apply weight also for S simulations
    Sval(Ieb).WS = Sval(Ieb).S;    %.*InPar.Spec(Ieb)./InPar.Back(Ieb);
    
    if (InPar.PoissonFilter)
        % calculate the PMF for the PSF
        % Ofek & Zackay
       
        CellPSF{Ieb} = log(1 + CellPSF{Ieb}.*Res.Flim(Ieb)./InPar.Back(Ieb));
    end
    
    % Filter the NH image with the appropriate filter
    switch func2str(InPar.FilterFun)
        case 'filter2'
             FiltCube(Ieb,:,:) = InPar.FilterFun(CellPSF{Ieb},Image,'same');
        otherwise
             FiltCube(Ieb,:,:) = InPar.FilterFun(Image,CellPSF{Ieb});
    end
    
    
    
end

% Calculate the joint threshold in S
Slim = quantile([Sval.WS],1-InPar.Beta);
    
% Join the filtered images
Weight    = InPar.Spec./InPar.Back;
FiltImage = bsxfun(@times,FiltCube,Weight(:));

% coadd images
FiltImage = squeeze(sum(FiltImage,1));

if (nargout>3)
    % Threshold the image and look for local maxima
    % set value below threshold to zero
    ThImage = FiltImage;
    ThImage(FiltImage<Slim) = 0;
    % look for local maxim
    BW = imregionalmax(ThImage,InPar.Conn);
    % Find sources
    Isrc  = find(BW>0);
    [Src.Y,Src.X] = ind2sub(size(BW), Isrc);
end

    
    