function [OutImageCell,OutMatrixCell]=pixinterp_fits(InImages,OutImages,PixList,varargin);
%--------------------------------------------------------------------------
% pixinterp_fits function                                            ImBasic
% Description: Interpolate over a set of pixels in a set of images.
% 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 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.
%          - List of pixels over which to interpolate or a mask image.
%            This can be a single list to apply to all input images
%            or a cell array of lists in which
%            each list is applied to its corresponding input image.
%            Single or  multiple lists can be either strings
%            or cell arrays (see create_list.m for details).
%            or matrix or cell array of matrices.
%            A "list" maybe either a coordinates list [X, Y],
%            a matrix, mat file containing matrix, or a FITS file.
%            If the list has two columns (i.e., [X, Y]) then 
%            will interpolate over the pixels identified by X,Y.
%            If the list is some kind of image or matrix
%            (having the same size as the input image) then
%            it will be treated as a mask image.
%            In this case the program
%            will interpolate over pixels which vale range is
%            specified by the 'Range' parameter keyword.
%            In addition pixels which are set to NaN will be interpolated.
%            The input lits can be empty - in this case the action
%            is according to the 'NaN' keyword value.
%          * Arbitrary number of argument pairs: ...,keyword,value,...
%            Where each keyword is followd by a value.
%            The following keywords are available:
%            'Range'     - If the list of pixels to interpolate is given
%                          in a mask image format then the program will
%                          search the mask for pixels that their value
%                          is within Range(1) and Range(2),
%                          and interpolate over this pixels.
%                          In addition, pixels that their value in the mask image
%                          is NaN will be interpolated.
%                          If empty matrix, then will search only for NaN values.
%                          Default is [].
%            'NaN'       - Specify if to interpolate also over pixels which their
%                          value is NaN in the input image {'y'|'n'}, default is 'n'.
%            'Xin'       - Vector of X coordinates corresponding to the pixels
%                          of the input image. Default is [1:1:size(InImage,2)].
%            'Yin'       - Vector of Y coordinates corresponding to the pixels
%                          of the input image. Default is [1:1:size(InImage,1)].
%            'InterpMethod'-Interpolation method. See interp2.m for options.
%                          Default is 'cubic'.
%            '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 float32.
%            '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 : - A cell array containing the names of the output images.
%          - A cell array containing the output images in a mtrix format.
% 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    = 'float32';
DefV.Save        = 'y';      % {'n' | 'y' | 'FITS' | 'mat'} 
DefV.FileType    = 'FITS';   % {'FITS' | 'mat'}
DefV.FitsReadPar = {};
% specific parameters for imtrans_fits.m
DefV.Range       = [];
DefV.NaN         = 'n';
DefV.Xin         = [];
DefV.Yin         = [];
DefV.InterpMethod= 'cubic';
InPar = set_varargin_keyval(DefV,'y','use',varargin{:});

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

%--- Deal with specific parameters ---                                         
%  PixList
[PixImagesCell,PixMatCell,PixIsNumeric,PixNim] = imget_fits(PixList,...
  					                    'FileType',InPar.FileType,...
						            'ReadAll','n');


if (PixNim==1 || PixNim==Nim),
   % ok
else
   error('PixList has illegal size');
end

if (PixNim==1),
   % read PixList
   if (~PixIsNumeric),
      [~,PixImage] = imget_fits(PixImagesCell{1},...
                            'FileType',InPar.FileType,...
                            'FitsReadPar',InPar.FitsReadPar,...
                            'ReadAll','y',...
                            'Cell','n');
      PixIsNumeric = 1;
      PixMatCell{1} = PixImage;
      clear PixImage;
   end
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

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

   %--------------------
   %--- Read PixList ---
   %--------------------   
   if (IsNumeric==1),
      PixImage = PixMatCell{min(PixNim,Iim)};
   else
      [~,PixImage] = imget_fits(PixImagesCell{Iim},...
                               'FileType',InPar.FileType,...
                               'FitsReadPar',InPar.FitsReadPar,...
                               'ReadAll','y',...
                               'Cell','n');
   end

   %------------------------------------
   %--- select pixels to interpolate ---
   %------------------------------------

   switch lower(InPar.NaN)
    case 'y'
       [Yi, Xi] = find(isnan(InImage));
    case 'n'
       Yi = [];
       Xi = [];
    otherwise
       error('Unknown NaN option');
   end

   if (size(PixImage,2)==2),
      % assume PixList is a list of [X,Y] coordinates
      Xi = [Xi; PixImage(:,1)];
      Yi = [Yi; PixImage(:,2)];
   else
      if (~isempty(InPar.Range)),
         [Yii,Xii] = find((PixImage>InPar.Range(1) & PixImage<InPar.Range(2)) | isnan(PixImage));
         Xi = [Xi; Xii];
         Yi = [Yi; Yii];
      end
   end

   
   %---------------------
   %--- work on image ---
   %---------------------
   % interpolate
   if (isempty(InPar.Xin)),
      InPar.Xin = [1:1:size(InImage,2)]';
   end
   if (isempty(InPar.Yin)),
      InPar.Yin = [1:1:size(InImage,1)]';
   end
   OutImage = InImage;
    OutImage(Yi,Xi) = interp2(InPar.Xin,InPar.Yin,InImage,Xi,Yi,InPar.InterpMethod);
not good because interp2 uses also the value pre-interpolation....


   %--- 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 pixinterp_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
   
   
