function [OutImageName,OutImageMatrix,TForm,RSamp,Report,CandList,BestInd]=imregs_1fits(RefIm,ListIm,OutImageName,RefCat,ListCat,varargin)
%--------------------------------------------------------------------------
% imregs_1fits function                                            ImBasic
% Description: Register (align) a single image relative to a reference
%              image using commom stars detected in the images.
%              See imregs_fits.m for the multiple input files version.
% Input  : - Reference image. This is either a string containing a
%            reference FITS image file name or a matrix containing the
%            image.
%          - An image to align in respect to the reference image. This is
%            either a string containing a FITS image file name or a
%            matrix containing the image.
%          - Output FITS/mat image name. Default is empty matrix.
%            If empty matrix then use 'OutPrefix', 'OutDir' and ListIm
%            to construct the output file name.
%          - Reference catalog. Default is empty matrix (i.e., []).
%            If empty than will generate a reference catalog using
%            SExtractor. Otherwise, a SExtractor reference catalog should
%            be provided either as a string containing a FITS binary table
%            name, or a matrix, or a structure in which one of the fields
%            is the catalog matrix.
%            In case this is a structure then the field name containing
%            the catalog name is taken from the optional paramter
%            'CatName', and a field name containing the column names is
%            taken from the optional parameter 'ColName'. If 'ColName'
%            field is not given than the default value of the Col
%            structure is taken from the 'Col' optional parameter.
%          - Image catalog (of the image to align). If empty matrix
%            then the catalog will be generated using SExtractor.
%            Otherwise, a SExtractor reference catalog should be provided.
%            See previous parameter for details.
%          * Arbitrary number of pairs of input parameters: keyword,value,...
%            Available keywords are:
%            '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., '').
%            'MaxN'    - The maximum number of stars to use in each list in
%                        the initial matching process. Default is 500.
%            'Frac'    - Minimum fraction of stars in the list (not reference list)
%                        requireed for a possible solution. Default is 0.2;
%            'ColSort' - The column in the list and the reference list by
%                        which to sort before selecting the first 'MaxN' stars.
%                        For example, if you match two lists of stars, it is
%                        recommended to sort by the magnitude.
%                        If empty matrix (e.g., []) then do not sort.
%                        Default is [].
%            'BinSize' - Search radius for the initial matching process
%                        in the same units in which the coordinates are
%                        given. Default is 1.
%            'CheckAll'- Check all possible shifts {'y'|'n'}. Default is
%                        'n'. 'n' will first calculate the histogram of DX
%                        calculate the local S/N per bin and only check
%                        shifts which have S/N>LimitSN.
%            'LimitSN' - Limit for S/N in which to choose DX values to check.
%                        Default is 5.
%            'HistSmooth'- Size of histogram median smoothing which is used
%                        to calculate the background in the S/N estimation.
%                        Default is 10.
%            'Refine'  - Refine initial solution {'y'|'n'} by trying to
%                        match the catalogs using each one of the possible
%                        solutions and calculate the matching statistics
%                        for each solution. Default is 'y'.
%            'Tol'     - Matching tolerance used after initial matching
%                        Default is 1 [in units of the X/Y dimensions].
%            'Fit'     - Fit method. Options are:
%                        'xy'     - shift in x and y coordinates (default).
%                        For more options see fit_lists2d.m
%            'RSamp'   - Transformation interpolation kernel sampling
%                        structure created by
%                        makeresampler.m. Default is empty matrix.
%                        If empty matrix then use
%                        makeresampler('cubic','bound').
%            'XData'    - Output image X-axis range in the refernce image
%                         reference frame (see imtransform.m for details).
%                         Default is [1 size(InImageMatrix,2)].
%            'YData'    - Output image Y-axis range in the refernce image
%                         reference frame (see imtransform.m for details).
%                         Default is [1 size(InImageMatrix,1)].
%            'MeanType'- Type of function {'Mean'|'Median'} by which
%                        to estimate the best shift. Default is 'mean'.
%            'Tran'    - Transformation type to apply to lists
%                        'shiftxy' - shift in X and Y (default).
%            'MaxIter' - Maximum number of sigma clipping iterations.
%                        Default is 1.
%                        If 0, then will use only initial shift estimate.
%            'ColX'    - Column index of the X coordinate in the list
%                        and reference list. Default is 1.
%            'ColY'    - Column index of the Y coordinate in the list
%                        and reference list. Default is 2.
%            'Method'  - The method of clipping. Available methods are
%                        {'StD','StdP','Constant','MinMax','Perc'},
%                        Default is 'StD'.
%                        See clip_resid.m for details.
%            'Mean'    - Method by which to calculate the "mean" of the
%                        sample in the sigma clipping process.
%                        Available methods are:
%                        'mean'   - mean of the sample.
%                        'median' - median of the sample (default).
%                        See clip_resid.m for details.
%            'Clip'    - Two elements vector containing the lower and upper
%                        values for the sigma clipping.
%                        This is [Lower, Upper] number of sigmas (positive)
%                        below/above the mean.
%                        See clip_resid.m for details.
%            'RemConf' - Remove potentialy confused sources {'y'|'n'},
%                        default is 'y'. If 'y' then program will use only
%                        entries which does not have any other entries
%                        within ConfRad distance from them.
%            'ConfRad' - Distance threshold for confused sources selection.
%                        Default is 20.
%            'InShift' - Optional prelimnary shift to apply to each image
%                        [X, Y]. If this parameter is specified than
%                        the program will skip the first step of finiding
%                        the shift between the images and go immidetly to
%                        refining the shift.
%                        Default is empty matrix. If empty than will
%                        attempt to measure the initial shift
%                        automatically.
%            'CopyHead' - Copy header from original image {'y' | 'n'}.
%                         Default is 'y'.
%            '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.
%            'DelDataSec'- Delete original DATASEC keyword from header
%                         {'y' | 'n'}, default is 'y'.
%            'FitsReadPar' - Additional parameters to send to the fitsread
%                         function. Default is {}.
%            'DataType' - Output data type (see fitswrite.m for options), 
%                         default is float32.
%            'Save'     - Save FITS image to disk {'y' | 'n'}.
%                         Default is 'y'.
%            'Col'      - Default columns in the catalog containing the
%                         XWIN_IMAGE, YWIN_IMAGE and MAG_AUTO columns.
%                         Default is [1 2 3].
%            'ColName'  - Column name in the catalog structure containing
%                         the column structure. Default is 'Col'.
%            'CatName'  - Catalog name in the catalog structure containing
%                         the catalog matrix. Default is 'Cat'.
% Output  : - Output image Name
%           - Output image in matrix form.
%           - Transformation structure.
%           - Interpolation kerenel structure.
%           - Structure containing the best fit matching and transformation
%             information.
%           - Array structure of all possible shift solutions.
%           - Index of best fit solution.
% Tested : Matlab 7.13
%     By : Eran O. Ofek                     April 2012
%    URL : http://wise-obs.tau.ac.il/~eran/matlab.html
% Example: 
%--------------------------------------------------------------------------

Def.OutImageName   = [];
Def.RefCat  = [];
Def.ListCat = [];
if (nargin==2),
    OutImageName     = Def.OutIm;
    RefCat    = Def.RefCat;
    ListCat   = Def.ListCat;
elseif (nargin==3),
    RefCat    = Def.RefCat;
    ListCat   = Def.ListCat;
elseif (nargin==4),
    ListCat   = Def.ListCat;
else
    % do nothing
end


DefV.XData       = [];
DefV.YData       = [];
DefV.RSamp       = [];
DefV.OutPrefix   = '';
DefV.OutDir      = '';
DefV.CopyHead    = 'y';
DefV.AddHead     = [];
DefV.DataType    = 'float32';
DefV.Save        = 'y';
DefV.ColName     = 'Col';
DefV.CatName     = 'Cat';
DefV.Col         = [1 2 3];
DefV.Fit        = 'xy';
DefV.HistSmooth = 10;
DefV.LimitSN    = 5;
DefV.CheckAll   = 'n';
DefV.MaxN      = 500;
DefV.Frac      = 0.2;
DefV.ColSort   = [];
DefV.BinSize   = 1;
DefV.Refine    = 'y';
DefV.Tol       = 1;
DefV.MaxIter   = 1;
DefV.MeanType  = 'mean';
DefV.Method    = 'StD';
DefV.Mean      = 'median';
DefV.Clip      = [];
DefV.ColX      = 1;
DefV.ColY      = 2;
DefV.RemConf   = 'y';
DefV.ConfRad   = 20;
DefV.InShift   = [];
DefV.DelDataSec = 'y';
DefV.FitsReadPar = {};

InPar = set_varargin_keyval(DefV,'y','use',varargin{:});

%------------------------------------------------------------------------
%--- If no catalog is provided then generate catalog using SExtractor ---
%------------------------------------------------------------------------
ParSex = {'XWIN_IMAGE','YWIN_IMAGE','MAG_AUTO'};
if (isempty(RefCat)),
    [TmpRefCat,~,RefCol] = run_sextractor(RefIm,[],ParSex);
    RefTable    = TmpRefCat{1};
else
    if (ischar(RefCat)),
        % read RefCat from FITS binary table
        [~,~,~,RefCol,RefTable]=get_fitstable_col(RefCat);
    else
        if (isstruct(RefCat))
            FN = fieldnames(RefCat);
            if (isempty(find(strcmp(InPar.ColName,FN))));
                % no column structure - use default.
                RefCol.XWIN_IMAGE = InPar.Col(1);
                RefCol.YWIN_IMAGE = InPar.Col(2);
                RefCol.MAG_AUTO   = InPar.Col(3);
            else
                RefCol = RefCat.(InPar.ColName);
            end
            if (isempty(find(strcmp(InPar.CatName,FN))));
                error('RefCat structure must contain a Cat field');
            else
                RefTable = RefCat.(InPar.CatName);
            end
        else
            RefTable = RefCat;
            % no column structure - use default.
            RefCol.XWIN_IMAGE = InPar.Col(1);
            RefCol.YWIN_IMAGE = InPar.Col(2);
            RefCol.MAG_AUTO   = InPar.Col(3);
        end
    end
end



%------------------------------------------------------------------------
%--- If no catalog is provided then generate catalog using SExtractor ---
%------------------------------------------------------------------------
if (isempty(ListCat)),
    [TmpCat,~,Col] = run_sextractor(ListIm,[],ParSex);
    Table          = TmpCat{1};
else
    if (ischar(ListCat)),
        % read RefCat from FITS binary table
        [~,~,~,Col,Table]=get_fitstable_col(ListCat);
    else
        if (isstruct(ListCat))
            FN = fieldnames(ListCat);
            if (isempty(find(strcmp(InPar.ColName,FN))));
                % no column structure - use default.
                Col.XWIN_IMAGE = InPar.Col(1);
                Col.YWIN_IMAGE = InPar.Col(2);
                Col.MAG_AUTO   = InPar.Col(3);
            else
                Col = ListCat.(InPar.ColName);
            end
            if (isempty(find(strcmp(InPar.CatName,FN))));
                error('RefCat structure must contain a Cat field');
            else
                Table = ListCat.(InPar.CatName);
            end
        else
            Table = ListCat;
            % no column structure - use default.
            Col.XWIN_IMAGE = InPar.Col(1);
            Col.YWIN_IMAGE = InPar.Col(2);
            Col.MAG_AUTO   = InPar.Col(3);
        end
    end
end


% match two tables using match2s_lists.m  
[CandList,BestInd]=match2s_lists(RefTable,Table,'ColX',Col.XWIN_IMAGE',...
                                    'ColY',Col.YWIN_IMAGE,...
                                    'ColSort',Col.MAG_AUTO,...
                                    varargin{:});
                                
                                %    'MaxN',100,'Frac',0.1,'BinSize',3,'Fit','xy');
                                           
                                
   save CandList.mat CandList

% prepare transformation:
if (isempty(CandList)),
    TForm          = [];
    RSamp          = [];
    Report         = [];
    OutImageMatrix = [];
else
   Report = CandList(BestInd);
   TForm = maketform('custom',2,2,...
                     CandList(BestInd).Forward_tf,...
                     CandList(BestInd).Inverse_tf,...
                     CandList(BestInd).Par);
   %TForm.Matching = CandList(BestInd);
    
    
   %switch lower(InPar.Save)
   %    case 'y'
   %        [OutImageCell] = prep_output_images_list(ListIm,OutIm,varargin{:});
   %        OutImageName   = OutImageCell{1};
   %    otherwise
   %        OutImageName   = [];
   %end

   [OutImageName,OutImageMatrix,RSamp] = imtrans_1fits(ListIm,OutImageName,TForm,...
                                              'RSamp',InPar.RSamp,...
                                              'OutPrefix',InPar.OutPrefix,...
                                              'OutDir',InPar.OutDir,...
                                              'CopyHead',InPar.CopyHead,...
                                              'AddHead',InPar.AddHead,...
                                              'DelDataSec',InPar.DelDataSec,...
                                              'FitsReadPar',InPar.FitsReadPar,...
                                              'DataType',InPar.DataType,...
                                              'Save',InPar.Save);
                                          
end