function [flatfield, varargout] = BaSiC(IF,varargin)

% Estimation of flatfield for optical microscopy. Apply to a collection of 
% monochromatic images. Multi-channel images should be separated, and each 
% channel corrected separately.
% 
% Usage:  FLATFIELD = ShadingEst(D_original, ...)
%
% Input:  IF is an NROWSxNCOLZxNIMAG array containing 
%                the images to be corrected.
%
%                SLIC also supports the following optional arguments,
%                which override default parameters. 
%         
%          'darkfield'(default = 'false') whether you would like to estimate
%          darkfield, keep 'false' if the input images are brightfield
%          images or bright fluoresence images, set 'true' only if only if
%          the input images are fluorescence images are dark and have a strong
%          darkfield contribution.  
%
%         'lambda' (default value estimated from input images) high values (eg. 9.5) increase the
%                spatial regularization strength, yielding a more smooth
%                flatfield
%
%          'lambda_dark' (default value estimated from input images) high values (eg. 9.5) increase the
%                spatial regularization strength, yielding a more smooth
%                darkfield
%             
%          'optimization_tol' (default = 1e-5) specifies the tolerance of error in the optimization.
%
%          'max_iterations' (default = 500) specifies the maximum
%                number of iterations allowed in the optimization.
%
% Output: flatfield: estimated gain
%         darkfield: estimated offset
%        
% 

% 15.2.2015. Tingying Peng: equiry for usage addressed to
% tying84ster@gmail.com

addpath('dcttool');

% parse the input arguments, return a structure containing parameters
options = BaSiC_parseInputs(varargin);

% set default values for options that are not specified

if isempty(options.estimation_mode)
    options.estimation_mode = 'l0';
end
if isempty(options.max_iterations)
    options.max_iterations = 500;
end
if isempty(options.optimization_tol)
    options.optimization_tol = 1e-6;
end


if not(isempty(options.initial_flatfield)) % check size
   if size(options.initial_flatfield,1)==size(IF,1)&&size(options.initial_flatfield,2)==size(IF,2)
       flag_initialflatfield = true;
       initial_flatfield = imresize(options.initial_flatfield,[options.working_size options.working_size]);
   else
       Disp('Wrong size of input flatfield, skip the initial flatfield then...');      
       flag_initialflatfield = false;
   end
else
    flag_initialflatfield = false;
end
flag_segmentation = false;
if not(isempty(options.segmentation)) % check size
   if isequal(size(IF), size(options.segmentation))
       options.segmentation = imresize(options.segmentation,[options.working_size options.working_size]);
       flag_segmentation = true;
   else
       Disp('Wrong size of segmentation, skip segmentation then...');      
   end
end

if isempty(options.darkfield)
    if flag_initialflatfield
        options.darkfield = 'true';
    else
        options.darkfield = 'false';
    end
end

nrows = options.working_size; ncols = options.working_size; 
D = double(imresize(IF,[nrows ncols],'bilinear'));
%medianD = median(D,3);
meanD = mean(D,3);
meanD = meanD./mean(meanD(:)); % mean 1.0
%meanD = (meanD-1.0)./std(meanD(:))+1.0; % std 1.0
W_meanD = mirt_dct2(meanD);

if isempty(options.lambda)
    options.lambda = sum(abs(W_meanD(:)))./(400)*0.5;
end
if isempty(options.lambda_darkfield)
    if flag_initialflatfield
        options.lambda_darkfield = options.lambda*0.1;
    else
        options.lambda_darkfield = options.lambda*0.2;
    end
end
if ~flag_segmentation
%     D = sort(D,3);
%     options.segmentation = sort(options.segmentation,3);
    %clear DI;
    D = sort(D,3);
end

XAoffset = zeros(nrows,ncols);

weight = ones(size(D));
if flag_segmentation
    weight(options.segmentation) = options.eplson;
end
i = 0;
flag_reweighting = true(1,1);
flatfield_last = ones(nrows,ncols);
darkfield_last = randn(nrows,ncols);
 while (flag_reweighting)
       i = i+1;
       disp(['Reweighting Iteration' num2str(i)]); 
       if flag_initialflatfield
           [X_k_A,X_k_E,X_k_Aoffset] = inexact_alm_rspca_l1_intflat(D, initial_flatfield,options,'weight',weight);
       else
           [X_k_A,X_k_E,X_k_Aoffset] = inexact_alm_rspca_l1(D, options,'weight',weight);
       end
       
       XA = reshape(X_k_A,nrows,ncols,[]);
       XE = reshape(X_k_E,nrows,ncols,[]);
       XAoffset = reshape(X_k_Aoffset,nrows,ncols);
       XE_norm = XE./(repmat(mean(mean(XA)),nrows,ncols)+1e-6);
       weight = 1./(abs(XE_norm)+options.eplson);
       if flag_segmentation
         weight(options.segmentation)= weight(options.segmentation).*options.eplson;
       end
       weight = weight.*numel(weight)./sum(weight(:));
       temp = mean(XA,3)-XAoffset;
       flatfield_current = temp./mean(temp(:));
       darkfield_current = XAoffset;
       mad_flatfield = sum(abs(flatfield_current(:)-flatfield_last(:)))./sum(abs(flatfield_last(:)));
       temp_diff = sum(abs(darkfield_current(:)-darkfield_last(:)));
       if (temp_diff<1e-7)
          mad_darkfield = 0;
       else
          mad_darkfield =temp_diff./(max(sum(abs(darkfield_last(:))),1e-6));
       end
       flatfield_last = flatfield_current;
       darkfield_last = darkfield_current;
       if or(max(mad_flatfield,mad_darkfield)<=options.reweight_tol,i>options.max_reweightiterations)
          flag_reweighting = 0;
       end
end
  

shading =  mean(XA,3)-XAoffset;
%XAoffset = XAoffset+B_offset.*shading;
flatfield = imresize(shading,[size(IF,1) size(IF,2)]);
flatfield = flatfield./mean(flatfield(:));
if strcmpi(options.darkfield,'true')
    darkfield = imresize(reshape(XAoffset,nrows,ncols,[]),[size(IF,1) size(IF,2)]);
else
    darkfield = [];
end
varargout{1} = darkfield;


