function [fi_base,flat_amp] = BasefluorEst(IF,flatfield,darkfield,varargin) %,segmentation)

nrows = 128; ncols = 128;%nframe = size(IF,3);
D = double(imresize(IF,[nrows ncols],'bilinear'));
D = reshape(D,nrows*ncols,[]);
flatfield= double(imresize(flatfield,[nrows ncols],'bilinear'));
if isempty(darkfield)
    darkfield = zeros(size(flatfield));
else
    darkfield = double(imresize(darkfield,[nrows ncols],'bilinear'));
end
if ~isempty(varargin)
    D_seg = varargin{1};
    D_seg = imresize(D_seg,[nrows,ncols],'nearest');
    D_seg = reshape(D_seg,nrows*ncols,[]);
else
    D_seg = [];
end

eplson = 0.1;
tol = 1e-6;
maxIter = 500;
weight = ones(size(D));
if ~isempty(D_seg)
   weight(D_seg)=eplson;
end
for reweighting_iter = 1:2
    W_idct_hat = flatfield(:);
    A_offset = darkfield(:);
    A1_coeff = mean(D);
    %C1_coeff = ones(1,length(A1_coeff));
    % main iteration loop starts
    temp = svd(D, 'econ');
    norm_two = temp(1);
    mu = 12.5/norm_two;% this one can be tuned
    mu_bar = mu * 1e7;
    rho = 1.5 ;         % this one can be tuned
    d_norm = norm(D, 'fro');
    ent1 = 1;
    iter = 0;
    total_svd = 0;
    converged = false;
    A1_hat = zeros(size(D));
    E1_hat = zeros(size(D));
    Y1 = 0;
    while ~converged       
        iter = iter + 1;
        A1_hat = bsxfun(@plus, W_idct_hat*A1_coeff,A_offset);
        %A1_hat = bsxfun(@plus, bsxfun(@times, (1+(W_idct_hat-1.0)*C1_coeff),A1_coeff),A_offset);
        % update E1 using l0 norm
        E1_hat = E1_hat + (D - A1_hat - E1_hat + (1/mu)*Y1)./ent1;
        E1_hat = max(E1_hat - weight/(ent1*mu), 0)+min(E1_hat + weight/(ent1*mu), 0);
        % update A1_coeff, A2_coeff and A_offset
        %if coeff_flag

            R1 = bsxfun(@minus,D-E1_hat,A_offset); 
            %if isempty(D_seg)
                A1_coeff = mean(R1);
            %else
               % A1_coeff = maskmean(R1,~D_seg,1);
            %end
               
%             R1 = D-E1_hat;
%             A1_coeff = mean(R1)-mean(A_offset);
            %[A1_mincoeff,A1_mincoefftile]= min(A1_coeff);
            %A1_back(A1_coeff<1) = 1;
            A1_coeff(A1_coeff<0) = 0;
            %C1_coeff = mean(R1./(W_idct_hat*A1_coeff));
            %C1_coeff = C1_coeff./mean(C1_coeff(:));

            %A_offset= max(min(A_offset,A_uplimit),0);
         Z1 = D - A1_hat - E1_hat;
        %Z2 = D - A2_hat - E2_hat;

         Y1 = Y1 + mu*Z1;
        %Y2 = Y2 + mu*Z2; 

        mu = min(mu*rho, mu_bar);

        %% stop Criterion  
        stopCriterion = norm(Z1, 'fro')/ d_norm;
        %stopCriterion = max(norm(Z1, 'fro')/ d_norm, norm(Z2,'fro')/d_norm);
        if stopCriterion < tol
            converged = true;
        end

        if mod( total_svd, 10) == 0
            disp(['Iteration ' num2str(iter) ...
                ' |E1|_0 ' num2str(sum((abs(E1_hat(:))>0)))...
                ' stopCriterion ' num2str(stopCriterion)]);
        end  
    end
    % %% fine tuning flat_amp
    A2_coeff = zeros(1,length(A1_coeff));
    %C1_coeff = ones(1,length(A1_coeff));
    % main iteration loop starts
    temp = svd(E1_hat, 'econ');
    norm_two = temp(1);
    mu = 12.5/norm_two;% this one can be tuned
    mu_bar = mu * 1e7;
    rho = 1.5 ;         % this one can be tune
    d_norm = norm(E1_hat, 'fro');
    ent1 = 1;
    iter = 0;
    total_svd = 0;
    converged = false;
    %A2_hat = zeros(size(E1_hat));
    E2_hat = zeros(size(E1_hat));
    Y1 = 0;
    while ~converged 
        iter = iter + 1;
        A2_hat = (W_idct_hat-1.0)*A2_coeff;
        % update E1 using l0 norm
        E2_hat = E2_hat + (E1_hat - A2_hat - E2_hat + (1/mu)*Y1)./ent1;
        E2_hat = max(E2_hat - weight/(ent1*mu), 0)+min(E2_hat + weight/(ent1*mu), 0);
        %E2_hat(E2_hat<0) = 0;
        % update A1_coeff, A2_coeff and A_offset
        %if coeff_flag

            R1 = bsxfun(@rdivide,(E1_hat-E2_hat),(W_idct_hat-1.0));
            A2_coeff = median(R1);
            %C1_coeff = C1_coeff./mean(C1_coeff(:));

            %A_offset= max(min(A_offset,A_uplimit),0);
         Z1 = E1_hat - A2_hat - E2_hat;
        %Z2 = D - A2_hat - E2_hat;

         Y1 = Y1 + mu*Z1;
        %Y2 = Y2 + mu*Z2; 

        mu = min(mu*rho, mu_bar);
        
        %% stop Criterion 
        stopCriterion = norm(Z1, 'fro')/ d_norm;

        if stopCriterion < tol
            converged = true;
        end

        if mod( total_svd, 10) == 0
            disp(['Iteration ' num2str(iter) ...
                ' |E2|_0 ' num2str(sum((abs(E2_hat(:))>0)))...
                ' stopCriterion ' num2str(stopCriterion)]);
        end    


        if ~converged && iter >= maxIter
            disp('Maximum iterations reached') ;
            converged = 1 ;       
        end
    end
    % update weight
    XE_norm = bsxfun(@rdivide, E2_hat, A1_coeff);
    %XE_norm = E2_hat;
    weight = 1./(abs(XE_norm)+eplson);
    if ~isempty(D_seg)
        weight(D_seg)=weight(D_seg).*eplson;
    end
     weight = weight.*numel(weight)./sum(weight(:)); 
end
fi_base = A1_coeff;flat_amp = 1.0+A2_coeff./A1_coeff;
