function Rec = Split_Bregman_demosaicking_solver(Y_matrix, Rec, alpha, beta, mu, Mask, tol, iter)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  This function solves the following sub-problem:
%  argmin 1/2*\|Phi*X -Y\|^2_2 + alpha*\|DCX\|_1 + beta*sum(\|F_iX - Li\|^2_F)
%  Author: Kan Chang, Guangxi University
%  Date: 1/14/2015
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initialization
[Height, Width, ~] = size(Rec);
[D, Dt, DtD] = DiffOper_d(Height, Width);
n = Height*Width;
   
Y = double(Y_matrix(:));
X = Y;

d = zeros(6*n,1);   % 2n for R/G/B
b = zeros(6*n,1);    

% bregman iteration
for i = 1:iter
    X_p = X;
    
    %% Step 1: solving X subproblem
    H_L = @(x) Hfun_L(x, DtD, mu, beta, Height, Width, Mask);
    H_R = Hfun_R(Y, mu, beta, b, d, Dt, Rec(:));
    [X,~] = cgs(H_L, H_R, 1e-5, 100);
    
    %% Step 2: Solving d subproblem
    R_G = X(1:n) - X(n+1:2*n);
    G_B = X(n+1:2*n) - X(2*n+1:3*n);
    B_R = X(2*n+1:3*n) -  X(1:n);
    DCX = [D*R_G
           D*G_B
           D*B_R];

    d = shrinkage(DCX + b, alpha/mu);
    
    %% Step 3: update b
    b = b + DCX -d;
    
    %% Early termination
    if norm(X(:)-X_p(:))/norm(X_p(:)) < tol
%            fprintf(1,'Loop num of Bregman is %d\n', i);
        break; 
    end   
end

Rec(:,:,1) = reshape(X(1:n), Height, Width);
Rec(:,:,2) = reshape(X(n+1:2*n), Height, Width);
Rec(:,:,3) = reshape(X(2*n+1:3*n), Height, Width);

end


%% Aux function
function y = Hfun_L(x, DtD, mu, beta, Height, Width, Mask)
    
    Mask_vec = Mask(:);
    x_sampled = x.*Mask_vec;
    
    x_matrix = reshape(x, Height, Width, 3);
    R_G = x_matrix(:,:,1) - x_matrix(:,:,2);  % R-G
    G_B = x_matrix(:,:,2) - x_matrix(:,:,3);  % G-B   
    B_R = x_matrix(:,:,3) - x_matrix(:,:,1);  % B-R
    
    R_G = R_G(:);
    G_B = G_B(:);
    B_R = B_R(:);
    DtD_RG =  DtD*R_G;
    DtD_GB =  DtD*G_B;
    DtD_BR =  DtD*B_R;
    CtDtDC_x = [ DtD_RG          - DtD_BR
                -DtD_RG + DtD_GB
                         -DtD_GB + DtD_BR];
    
    y = x_sampled + mu*CtDtDC_x + beta*x;        
end 

function Out = Hfun_R(Y, mu, beta, b, d, Dt, Li)
    n = size(d,1);
    temp_db = d-b;
    
    Dtdb_1 = Dt*temp_db(         1:  n/3);
    Dtdb_2 = Dt*temp_db(  n/3  + 1:2*n/3);
    Dtdb_3 = Dt*temp_db(2*n/3  + 1:  n);
    CtDtdb = [ Dtdb_1          - Dtdb_3
              -Dtdb_1 + Dtdb_2
                      - Dtdb_2 + Dtdb_3];
    
   Out = Y + mu*CtDtdb + beta*Li;
end

function y = shrinkage(x,b)
    y = max(abs(x)-b,0).*sign(x);
end