function Rec = Low_Rank_denosing_RGB(Init, Patch_size, patch_pos_arr, Par, h, w)

% patchset = patch_size * number of patch 
patchset  =  Im2Patch_RGB(Init, Patch_size);
        
b            =   sqrt(size(patchset,1));
N            =   h-b+1;
M            =   w-b+1;
r            =   1:N;
c            =   1:M; 

X            =   patchset;                              % 3 dimensions
Ys_R         =   zeros( size(X,1), size(X,2));   
Ys_G = Ys_R; 
Ys_B = Ys_R;
W_R          =   zeros( size(X,1), size(X,2));
W_G = W_R;
W_B = W_R;
L            =   size(patch_pos_arr,2);

% loop over every group
for  i  =  1 : L   
    B         =   [ X(:, patch_pos_arr(:, i),1)
                    X(:, patch_pos_arr(:, i),2)
                    X(:, patch_pos_arr(:, i),3)];       % the whole group
    
    [tmp_y]   =   W_SVT( double(B), Par.lambda );

    tmp_w = ones( size(tmp_y) );
    
    tmp_y_R = tmp_y(1:b^2,:);
    tmp_y_G = tmp_y(b^2 + 1:2*b^2,:);
    tmp_y_B = tmp_y(2*b^2+1:3*b^2,:);
    
    tmp_w_R = tmp_w(1:b^2,:);
    tmp_w_G = tmp_w(b^2 + 1:2*b^2,:);
    tmp_w_B = tmp_w(2*b^2+1:3*b^2,:);
    
    Ys_R(:, patch_pos_arr(:,i))   =   Ys_R(:, patch_pos_arr(:,i)) + tmp_y_R;
    Ys_G(:, patch_pos_arr(:,i))   =   Ys_G(:, patch_pos_arr(:,i)) + tmp_y_G;
    Ys_B(:, patch_pos_arr(:,i))   =   Ys_B(:, patch_pos_arr(:,i)) + tmp_y_B;
    
    W_R(:, patch_pos_arr(:,i))    =   W_R(:, patch_pos_arr(:,i)) + tmp_w_R;
    W_G(:, patch_pos_arr(:,i))    =   W_G(:, patch_pos_arr(:,i)) + tmp_w_G;
    W_B(:, patch_pos_arr(:,i))    =   W_B(:, patch_pos_arr(:,i)) + tmp_w_B;
end

dim     =  zeros(h,w,3);
wei   =  zeros(h,w,3);
Rec = zeros(h,w,3);

% weighting the output result
k       =   0;
for i  = 1:b
    for j  = 1:b
        k    =  k+1;
        dim(r-1+i,c-1+j,1)  =  dim(r-1+i,c-1+j,1) + reshape( Ys_R(k,:)', [N M]);
        dim(r-1+i,c-1+j,2)  =  dim(r-1+i,c-1+j,2) + reshape( Ys_G(k,:)', [N M]);
        dim(r-1+i,c-1+j,3)  =  dim(r-1+i,c-1+j,3) + reshape( Ys_B(k,:)', [N M]);
        wei(r-1+i,c-1+j,1)  =  wei(r-1+i,c-1+j,1) + reshape( W_R(k,:)', [N M]);
        wei(r-1+i,c-1+j,2)  =  wei(r-1+i,c-1+j,2) + reshape( W_G(k,:)', [N M]);
        wei(r-1+i,c-1+j,3)  =  wei(r-1+i,c-1+j,3) + reshape( W_B(k,:)', [N M]);
    end
end
ind = find(wei==0);          % pixels which have not been proceed
n_ind = find(wei~=0);
Rec(n_ind) =  (dim(n_ind))./(wei(n_ind));
Rec(ind) = Init(ind);
end


function  X   =   W_SVT( Y, lambda )

[U0,Sigma0,V0]    =   svd(full(Y),'econ');
Sigma0            =   diag(Sigma0);
   
% logdet approximation
S                 =   max( Sigma0.^2/size(Y, 2), 0 );   
thr               =   lambda./ ( sqrt(S) + eps );
S                 =   shrinkage(Sigma0, thr);

r                 =   sum( S>0 );
U                 =   U0(:,1:r);                      
V                 =   V0(:,1:r);
X                 =   U*diag(S(1:r))*V';

end

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