function [im_out PSNR SSIM ]   =  CSR_Denoising( par )
time0         =   clock;
nim           =   par.nim;
[h  w ch]     =   size(nim);

par.step      =   1;
par.h         =   h;
par.w         =   w;
par.cls_num   =   64;

dim         =   uint8(zeros(h, w, ch));
ori_im      =   zeros(h,w);

% RGB->YUV
%
if  ch == 3
    n_im           =   rgb2ycbcr( uint8(nim) );
    dim(:,:,2)     =   n_im(:,:,2);
    dim(:,:,3)     =   n_im(:,:,3);
    n_im           =   double( n_im(:,:,1));
    
    if isfield(par, 'I')
        ori_im         =   rgb2ycbcr( uint8(par.I) );
        ori_im         =   double( ori_im(:,:,1));
    end
else
    n_im           =   nim;
    
    if isfield(par, 'I')
        ori_im             =   par.I;
    end
end
disp(sprintf('The PSNR of the blurred image = %f \n', csnr(n_im(1:h,1:w), ori_im, 0, 0) ));

[d_im]   =  Denoising(n_im, par, ori_im);

if isfield(par,'I')
   [h w ch]  =  size(par.I);
   PSNR      =  csnr( d_im(1:h,1:w), ori_im, 0, 0 );
   SSIM      =  cal_ssim( d_im(1:h,1:w), ori_im, 0, 0 );
end

if ch==3
    dim(:,:,1)   =  uint8(d_im);
    im_out       =  double(ycbcr2rgb( dim ));
else
    im_out  =  d_im;
end    
disp(sprintf('Total elapsed time = %f min\n', (etime(clock,time0)/60) ));
return;


function  [d_im ]    =   Denoising(n_im, par, ori_im)

[h1 w1]    =   size(ori_im);
b2         =   par.win*par.win;

par.tau1    =   0.1;
par.tau2    =   0.2;
d_im        =   n_im;
lamada      =   0.02;
v           =   par.nSig;
cnt         =   1;

for k    =  1:3
    
    if k>2
        par.step =  2;
    end

    [PCA_D, cls_idx, s_idx, seg, D0]     =   Clustering_PCA_New( d_im, par, 70 );
    D0       =   dctmtx(b2);
    [blk_arr wei_arr]    =   find_blks_fast( d_im, par);
        
    for i  =  1 : 3
        
        d_im       =   d_im + lamada*(n_im - d_im);
        
        dif=d_im-n_im;
        vd=v^2-(mean(mean(dif.^2)));        
        
        if (i ==1 && k==1)
            par.nSig=sqrt(abs(vd));            
        else
            par.nSig=sqrt(abs(vd))*par.lamada;
        end
    
        [Tau1, Tau2]    =   Cal_Sparsity_Parameters( d_im, cls_idx, PCA_D, par, s_idx, seg, D0, blk_arr, wei_arr );        
        
        d_im            =   CSR_Thresholding(d_im, par, D0, blk_arr, wei_arr, Tau1, Tau2, 1, PCA_D, cls_idx, s_idx, seg);
        
        PSNR            =   csnr( d_im(1:h1,1:w1), ori_im, 0, 0 );
        fprintf( 'Preprocessing, Iter %d : PSNR = %f,   nsig = %3.2f\n', cnt, PSNR, par.nSig );
        cnt   =  cnt + 1;
        imwrite(d_im./255, 'Results\tmp.tif');                
    end    
    
end


%--------------------------------------------------------------------------
% Compute the local sparsity parameters
%
%--------------------------------------------------------------------------
function  [Tau1, Tau2]   =  Cal_Sparsity_Parameters( im, PCA_idx, PCA_D, par, s_idx, seg, A, blk_arr, wei_arr )

b        =  par.win;
s        =  par.step;
b2       =  b*b;
[h  w]   =  size(im);

N       =  h-b+1;
M       =  w-b+1;
r       =  [1:s:N];
r       =  [r r(end)+1:N];
c       =  [1:s:M];
c       =  [c c(end)+1:M];
X0      =  zeros(b*b, N*M);
X_m     =  zeros(b*b,length(r)*length(c),'single');

N       =  length(r);
M       =  length(c);
L       =  N*M;

% For the Y component
k    =  0;
for i  = 1:b
    for j  = 1:b
        k        =  k+1;        
        blk      =  im(i:end-b+i,j:end-b+j);
        X0(k,:)  =  blk(:)';
    end
end

idx      =   s_idx(seg(1)+1:seg(2));
set      =   1:size(X_m,2);
set(idx) =   [];

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

C0          =   zeros(b2, L, 'single' );
Cu0         =   zeros(b2, L, 'single' );
   
idx         =   s_idx(seg(1)+1:seg(2));
    
for  k  =  1 : length(idx)
    i           =   idx(k);
    coe         =   A*X0(:, blk_arr(i, :));    
    C0(:,i)     =   mean(coe.^2, 2);
    
    coe        =   A*(X0(:, blk_arr(i, :)) - repmat(X_m(:, i), 1, par.nblk) );
    Cu0(:,i)   =   mean(coe.^2, 2);        
end

set        =   1:L;
set(idx)   =   [];
L          =   length(set);

for  k  =  1 : L
    i         =   set(k);    
    cls   =   PCA_idx(i);
    P     =   reshape(PCA_D(:, cls), b2, b2);
    
    coe       =   P*X0(:, blk_arr(i,:));
    C0(:,i)   =   mean(coe.^2, 2);
    
    coe       =   P*(X0(:, blk_arr(i,:)) - repmat(X_m(:, i), 1, par.nblk) );
    Cu0(:,i)  =   mean(coe.^2, 2);    
end

C0       =   max(0, C0-par.nSig^2);
Cu0      =   max(0, Cu0-par.nSig^2);
Tau1     =   (par.c1*sqrt(2)*par.nSig^2)./(sqrt(C0) + eps);
Tau2     =   (par.c2*sqrt(2)*par.nSig^2)./(sqrt(Cu0) + eps);
return;