function [Norm,StD,Nnn,Flat,Mat3D]=optimal_flatnorm_fits(Input,varargin)
%-----------------------------------------------------------------------------
% optimal_flatnorm_fits function                                      ImBasic
% Description:
% Input  : - Set of input images for which to find the optimal global
%            normalization.
%            This can be either:
%            A string containing FITS images name (including wild cards -
%            see create_list.m for options).
%            A cell array in which each cell contains a FITS image name.
%            A cell array in which each cell contains a 2D matrix.
%            A 3D matrix, in which the third dimension is the image index.
%            All the images should have the same size
%            (see create_mat3d_images.m).
%          * Arbitrary number of pairs of arguments: ...,keyword,value,...
%            Optional keywords are:
%            'MaxIter'    - Maximum number of iterations, default is 10.
%            'RatioThresh'- Convergance criterion for ratio between
%                           sucessive normalizations, default is 0.001.
% Output : - Vector of normalization - one per image.
%          - StD image of all normalized flats.
%          - Matrix in which each pixel contains the number of images
%            in which this pixel is not NaN.
%          - Flat field image. Note that this image should not be used
%            since it is better to use a median combined image.
%          - 3D matrix, of normalized images, in which the third dimension
%            is the image index.
% Tested : Matlab 7.11
%     By : Eran O. Ofek                     March 2011
%    URL : http://wise-obs.tau.ac.il/~eran/matlab.html
%-----------------------------------------------------------------------------

DefV.MaxIter     = 10;
DefV.RatioThresh = 0.0001;

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

Mat3D=create_mat3d_images(Input);




%--- debuging ---
%Mat3D = randn(2000,1000,10)+100;
%Mat3D(:,:,1) = Mat3D(:,:,1).*1.1;
%Mat3D(:,:,2) = Mat3D(:,:,2).*1.2;
%Mat3D(:,:,3) = Mat3D(:,:,3).*0.6;
%Mat3D(1,1,1:5)=NaN;


Nim   = size(Mat3D,3);
Ny    = size(Mat3D,1);
Nx    = size(Mat3D,2);
Npix  = Nx.*Ny;
Mat2D = reshape(Mat3D,[Npix Nim]);   % 2D representation: f_{ij}

% make sure that Mat2D (i.e., f_{ij}) doesnt contain zeros
Mat2D(find(Mat2D==0))=NaN;
 
% guess Alpha (nomalization), based on median value
Alpha = nanmedian(Mat2D);

OldAlpha = ones(1,Nim).*Inf;
MaxRatio = max(abs(Alpha./OldAlpha-1));
Iter     = 0;
while (Iter<InPar.MaxIter & MaxRatio>InPar.RatioThresh)
   Iter = Iter + 1;

NormMat2D = bsxfun(@times,Mat2D,Alpha);
NormMat2D = bsxfun(@times,NormMat2D,1./nanmedian(NormMat2D,1));
StD  = nanstd(NormMat2D,[],2);
nanmedian(StD)
MaxRatio

   OldAlpha = Alpha;
   % MatAlpha is a matrix of Alpha for each pixel,
   % however it is set to zeros in cases the pixel is NaN
   MatAlpha     = bsxfun(@times,~isnan(Mat2D),Alpha);
   %F_j   = sum(bsxfun(@times,Mat2D,Alpha),2)./sum(Alpha);

   % mean-based calculation for F_j
   %   F_j   = nansum(bsxfun(@times,Mat2D,Alpha),2)./sum(MatAlpha,2);
   % median-based calculation for F_j
   F_j   = nanmedian(bsxfun(@times,Mat2D,Alpha),2);

   F_j   = F_j./nanmedian(F_j(:));
   %nanmean(F_j(:))
   % mean-based calculation for Alpha
   %Alpha = nansum(bsxfun(@rdivide,F_j,Mat2D),1);
   % median-based calculation for Alpha
   Alpha = nanmedian(bsxfun(@rdivide,F_j,Mat2D),1).*Nim;

% deal with outliers
%   NormMat2D    = bsxfun(@times,Mat2D,Alpha);
%   MedNormMat2D = bsxfun(@times,NormMat2D,1./nanmedian(Mat2D,1));
%   StD2D        = repmat(nanstd(MedNormMat2D,[],2),1,Nim);
%   Nsig = 3;
%   Ioutl = find(abs(MedNormMat2D-1)>StD2D.*Nsig);
%   Mat2D(Ioutl) = NaN;

   MaxRatio = max(abs(squeeze(Alpha)./squeeze(OldAlpha)-1));
end

Norm = 1./(squeeze(Alpha)./Npix);
%Norm = 1./(squeeze(Alpha));

%[nanmedian(Mat2D,1).', nanmean(Mat2D,1).', Norm.', nanmedian(Mat2D,1).'./Norm.', nanmean(Mat2D,1).'./Norm.', nanmedian(bsxfun(@times,Mat2D,1./Norm),1)']

% StD image
NormMat2D = bsxfun(@times,Mat2D,1./Norm);
StD  = nanstd(NormMat2D,[],2);
StD = reshape(StD,[Ny Nx]);

if (nargout>2),
   Nnn = sum(~isnan(NormMat2D),2);
   Nnn = reshape(Nnn,[Ny Nx]);
end

if (nargout>3),
   Flat = reshape(F_j,[Ny Nx]);
end

if (nargout>4),
   % normalize images
   Mat3D = reshape(NormMat2D,[Ny Nx Nim]);
end
