function  [im_out]  =  CSR_Thresholding( im, opts, A, blk_arr, wei_arr, Tau1, Tau2, flag, PCA_D, PCA_idx, s_idx, seg)
[h w ch]   =   size(im);

b          =   opts.win;
b2         =   b*b*ch;
k          =   0;
s          =   opts.step;

N     =  h-b+1;
M     =  w-b+1;
L     =  N*M;
r     =  [1:s:N];
r     =  [r r(end)+1:N];
c     =  [1:s:M];
c     =  [c c(end)+1:M];
X     =  zeros(b*b,L,'single');
for i  = 1:b
    for j  = 1:b
        k    =  k+1;
        blk  =  im(i:h-b+i,j:w-b+j);
        blk  =  blk(:);
        X(k,:) =  blk';            
    end
end

X_m     =  zeros(length(r)*length(c),b*b,'single');
X0 = X';
for i = 1:opts.nblk
   v        =  wei_arr(:,i);
   X_m      =  X_m(:,:) + X0(blk_arr(:,i),:) .*v(:, ones(1,b2));
end
X_m=X_m';


ind         =   zeros(N,M);
ind(r,c)    =   1;
X           =   X(:, ind~=0);

N           =   length(r);
M           =   length(c);
L           =   N*M;
Y           =   zeros( b2, L );

idx          =   s_idx(seg(1)+1:seg(2));
tau1         =   opts.tau1;
tau2         =   opts.tau2;
if  flag==1
    tau1    =   Tau1(:, idx);
    tau2    =   Tau2(:, idx);
end

Y(:, idx)    =   A'*CSR_Shrinkage(A*X(:, idx), A*X_m(:, idx), tau1, tau2);

for   i  = 2:length(seg)-1   
    idx  =  s_idx(seg(i)+1:seg(i+1));    
    cls  =  PCA_idx(idx(1));
    P    =   reshape(PCA_D(:, cls), b2, b2);
    tau1         =   opts.tau1;
    tau2         =   opts.tau2;
    if  flag==1 
        tau1    =   Tau1(:, idx);
        tau2    =   Tau2(:, idx);        
    end
    
    Y(:, idx)  =  P'*CSR_Shrinkage(P*X(:, idx), P*X_m(:, idx), tau1, tau2);
end


% Output the processed image
im_out   =  zeros(h,w);
im_wei   =  zeros(h,w);
k        =  0;
for i  = 1:b
    for j  = 1:b
        k    =  k+1;
        im_out(r-1+i,c-1+j)  =  im_out(r-1+i,c-1+j) + reshape( Y(k,:)', [N M]);
        im_wei(r-1+i,c-1+j)  =  im_wei(r-1+i,c-1+j) + 1;       
    end
end

im_out  =  im_out./(im_wei+eps);
return;



function Y  = CSR_Shrinkage( x, mu, tau1, tau2 )

m1  =  (mu>=0).*mu;
m2  =  (mu<0).*mu;

y1   =  ( x < (-tau1-tau2) ).*( x+tau1+tau2 );
y3   =  ( x>(tau1-tau2) & x<(tau1-tau2+m1) ).*( x-tau1+tau2 );
y4   =  ( x>=(tau1-tau2+m1) & x<=(tau1+tau2+m1) ).*m1;
y5   =  ( x>(tau1+tau2+m1) ).*( x-tau1-tau2 );
Y1   =  (y1 + y3 + y4 + y5).*(mu>=0);

y1   =  ( x<(m2-tau1-tau2) ).*( x+tau1+tau2 );
y2   =  ( x>=(m2-tau1-tau2) & x<=(tau2+m2-tau1) ).*m2;
y3   =  ( x>(tau2+m2-tau1) & x<(tau2-tau1) ).*( x+tau1-tau2 );
y5   =  ( x>(tau1+tau2) ).*( x-tau1-tau2 );
Y2   =  (y1 + y2 + y3 + y5).*(mu<0);

Y    =  Y1 + Y2;
return;
