function [OutImageCell,OutMatrixCell,FlagedPixCell]=imcr_fits(InImages,OutImages,varargin);
%--------------------------------------------------------------------------
% imcr_fits function                                               ImBasic
% Description: Search for radiation hits ("cosmic rays") in a set of
%              images using different techniques. The cosmic ray
%              positions are rturned as a list of positions and a
%              FITS mask image.
% Input  : - List of input FITS images (see create_list.m for details).
%            Alternatively, this can be a cell array in which each cell
%            contains the image in a matrix form, or a matrix.
%          - Lits of output mask FITS images (see create_list.m for details).
%            If empty matrix will use the input images, along with the
%            'OutPrefix' and 'OutDir' parameters to determine
%            the output image names.
%            If empty matrix than the input images must be in FITS format.
%          * Arbitrary number of argument pairs: ...,keyword,value,...
%            Where each keyword is followd by a value.
%            The following keywords are available:
%            'Method'    - Method by which to search for cosmic rays:    
%                          Default is 'Del2Ps'.
%                          'Del2C' - Thresholding the discrete Laplacian of
%                                    the image.
%                                    Extrem value are stated in terms of
%                                    a consatnt value above below
%                                    to mark a pixel as cosmic ray
%                                    (see 'Del2Const').
%                          'Del2P' - Thresholding the discrete Laplacian of
%                                    the image.
%                                    Use percentile above/below to mark the
%                                    pixels as cosmic rays ('Del2Perc').
%                          'Del2Ps'- Like 'Del2Ps' but for the discrete Laplacian
%                                    of the image divided by the square root of
%                                    the image. Note that in this case the
%                                    image must be positive (i.e., not
%                                    background subtracted).
%            'Del2Const' - Constant values [below above] to mark as cosmic rays.
%                          Default is [-100 100].
%            'Del2Perc'  - Percentil below/above to mark as cosmic rays.
%                          Default is 0.999.
%            'BitPos'    - The bit position in the mask which to change.
%                          Default is 1.
%            'BitVal'    - Bit value (0|1), default is 1.
%            'Mask'      - Optional list of input masks which to edit in
%                          order to mark the cosmic cays
%                          This can be list of files (see create_list.m for details).
%                          Alternatively, this can be a cell array in which each cell
%                          contains the mask image in a matrix form, or a matrix.
%                          Default is []. of [], then will generate a new
%                          mask file which initial value is zero.
%                          Mask file type is determinde by the FileType
%                          parameter.
%            'OutPrefix' - Add prefix before output image names,
%                          default is empty string (i.e., '').
%            'OutDir'    - Directory in which to write the output images,
%                          default is empty string (i.e., '').
%            'AddHead'   - Cell array with 3 columns containing additional
%                          keywords to be add to the header.
%                          See cell_fitshead_addkey.m for header structure
%                          information. Default is empty matrix.
%            'CopyHead'  - Copy header from input image {'y'|'n'}, default is 'y'.
%            'DataType'  - Output data type (see fitswrite.m for options), 
%                          default is int16.
%            'DelDataSec'- Delete DATASEC keyword from output image header
%                          {'y'|'n'}, default is 'n'.
%            'Save'      - Save output images {'y'|'n'|'fits','mat'},
%                          default is 'y', which saves a fits file.
%                          If 'n' then will return only a matrix images.
%            
%            'FileType'  - Input file type {'FITS'|'mat'}, default is 'FITS'.
%                          This keyword is overrid if input is a matrix.
%            'FitsReadPar'-A cell array of additional parameters to pass to
%                          fitsread.m (e.g., {2}; which means read the image
%                          from the second extension.
%                          Default is {}.
% Output : -
% Tested : Matlab 7.11
%     By : Eran O. Ofek                   June 2011
%    URL : http://wise-obs.tau.ac.il/~eran/matlab.html
%--------------------------------------------------------------------------


% All im*_fits functions
DefV.OutPrefix   = '';
DefV.OutDir      = '';
DefV.DelDataSec  = 'n';
DefV.CopyHead    = 'y';
DefV.AddHead     = [];
DefV.DataType    = 'int16';
DefV.Save        = 'y';      % {'n' | 'y' | 'FITS' | 'mat'} 
DefV.FileType    = 'FITS';   % {'FITS' | 'mat'}
DefV.FitsReadPar = {};
% specific parameters for imtrans_fits.m
DefV.Method      = 'Del2Ps';
DefV.Del2Const   = [-100 100];
DefV.Del2Perc    = 0.999;
DefV.BitPos      = 1;
DefV.BitVal      = 1;
DefV.Mask        = [];
InPar = set_varargin_keyval(DefV,'y','use',varargin{:});


   
[InImagesCell,InMatCell,IsNumeric,Nim] = imget_fits(InImages,...
                                                    'FileType',InPar.FileType,...
                                                    'ReadAll','n');

%--- Deal with specific parameters ---                                         
[MaskInImagesCell,MaskInMatCell,MaskIsNumeric,MaskNim] = imget_fits(InPar.Mask,...
                                                            'FileType',InPar.FileType,...
                                                            'ReadAll','n');
if ~(MaskNim==1 || MaskNim==Nim),
    error('List of mask images have illegal size');
end


if (IsNumeric==1 && isempty(OutImages) && ~strcmpi(InPar.Save,'n'))
    error('Save=y option while input is a numeric matrix is possible only when OutImages are specified');
end
switch lower(InPar.Save)
    case {'y','fits','mat'}
        [OutImageCell] = prep_output_images_list(InImages,OutImages,varargin{:});
    otherwise
        OutImageCell = [];
end
        
HeaderInfo = [];
if (nargout>1),
   OutMatrixCell = cell(Nim,1);
end

FlagedPixCell = cell(Nim,1);
OutMatrixCell = cell(Nim,1);
for Iim=1:1:Nim,
   HeaderInfo = [];

   %--- read image ---
   if (IsNumeric==1),
      InImage = InMatCell{Iim};
   else
      [~,InImage] = imget_fits(InImagesCell{Iim},...
                               'FileType',InPar.FileType,...
                               'FitsReadPar',InPar.FitsReadPar,...
                               'ReadAll','y',...
                               'Cell','n');
      switch lower(InPar.FileType)
       case 'fits'
          ImHead = fitsinfo(InImagesCell{Iim});
          HeaderInfo = ImHead.PrimaryData.Keywords;
       otherwise
          HeaderInfo = [];
      end                   
   end
   
   %---------------------
   %--- work on image ---
   %---------------------
   switch lower(InPar.Method)
       case {'del2p','del2c','del2ps'}
           Del2 = del2(InImage);
           switch lower(InPar.Method)
               case 'del2ps'
                  Del2 = Del2./sqrt(InImage);
                  CL = err_cl(Del2(:),InPar.Del2Perc);
                  Ibad = find(Del2>CL(2) | Del2<CL(1));
               case 'del2p'
                  CL = err_cl(Del2(:),InPar.Del2Perc);
                  Ibad = find(Del2>CL(2) | Del2<CL(1));
               case 'del2c'
                  Ibad = find(Del2>InPar.Del2Const(2) | Del2<InPar.Del2Const(1));
               otherwise
                  error('Unknown Method option');
           end
   
       otherwise
           error('Unknown Method option');
   end
   
   [Ybad,Xbad] = ind2sub(size(InImage),Ibad);

   FlagedPixCell{Iim} = [Xbad, Ybad];
   
   Imaski = min(MaskNim,Iim);
   if (MaskIsNumeric),
       if (isempty(MaskInMatCell{Imaski})),
           % use zero intial mask
           MaskImage = zeros(size(InImage));
       else
          MaskImage = MaskInMatCell{Imaski};
       end
   else
       [~,MaskImage] = imget_fits(MaskInImagesCell{Iim},...
                               'FileType',InPar.FileType,...
                               'FitsReadPar',InPar.FitsReadPar,...          
                               'ReadAll','y',...
                               'Cell','n');
   end
   OutImage = MaskImage;
  
   OutImage(Ibad) = bitset(MaskImage(Ibad),InPar.BitPos,InPar.BitVal);
   
 
   %--- write image ---
   % assume output matrix is in: OutImage
   % and outoput image name is in: OutImageCell{Iim}
   
   if (nargout>1),
       OutMatrixCell{Iim} = OutImage;
   end
   
   switch lower(InPar.Save)
       case {'y','fits'}
           %------------------------
           %--- write FITS image ---
           %------------------------

           %--- Add to header comments regarding file creation ---
           [HeaderInfo] = cell_fitshead_addkey(HeaderInfo,...
                                       Inf,'COMMENT','','Created by imcr_fits.m written by Eran Ofek');


           %--- Modify WCS keywords if needed ---
           % TBD!!!!!                                   
                                   
           if (~isempty(InPar.AddHead)),
              %--- Add additional header keywords ---
              HeaderInfo = [HeaderInfo; InPar.AddHead];
           end

           switch lower(InPar.DelDataSec)
              case 'n'
                  % do nothing
              case 'y'
                  % delete DATASEC keyword from header
                  [HeaderInfo] = cell_fitshead_delkey(HeaderInfo,'DATASEC');
              otherwise
                   error('Unknown DelDataSec option');
           end

           %--- Write fits file ---
           fitswrite(OutImage,OutImageCell{Iim},HeaderInfo,InPar.DataType);
        
       case 'mat'
           %----------------------
           %--- wtite mat file ---
           %----------------------
           save(OutImageCell{Iim},'OutImage');
   
      
       case 'n'
           % do not write output
   
       otherwise
           error('Unknown Save option');
     
   end
end
   
   
