function [f,D]   =  BMLSVDTV( g, Ima_ori, Par)


f           = g;                                                        % Estimated Image



[Height Width]  = size(f);   
TotalPatNum     = (Height-Par.patsize+1)*(Width-Par.patsize+1);                 %Total Patch Number in the image
Dim             = Par.patsize*Par.patsize;  


[Neighbor_arr Num_arr Self_arr Keyr Keyc] =	NeighborIndex(g, Par);                  % PreCompute the all the patch index in the searching window 
            NL_mat              =   zeros(Par.patnum,length(Num_arr));          % NL Patch index matrix
            Pat_g_hat              =	zeros( Dim, TotalPatNum );
            Sigma_arr           =   zeros( 1, TotalPatNum);            
            TTalpha                =   zeros( size(Pat_g_hat) );     
            PTP                   =   zeros( size(Pat_g_hat) ); 
           Lambda          =   zeros(Dim, Par.patnum,TotalPatNum);
           
  f_old=f;         
            
for k = 1 : Par.Iter        
    
    g_hat=g*1./(1+Par.eta)+f*(Par.eta/(1+Par.eta));  % to produce U^j and V^j.
 
    [Pat_g_hat,Pat_f, Sigma_arr]	=	Im2Patch( g_hat, g, f, Par);                      % image to patch and estimate local noise variance            
    
    if (mod(k-1,Par.Innerloop)==0)
      %  Par.patnum = Par.patnum-2;                                             % Lower Noise level, less NL patches
        Lambda          =   Lambda(:,1:Par.patnum,:);
        fprintf('Block matching...  ');
      tic;
        NL_mat=  Block_matching(Pat_g_hat, Par, Neighbor_arr, Num_arr, Self_arr);% Caculate Non-local similar patches for each 
       t0=toc;
       fprintf('Computation time: %f \n',t0);
        if(k==1)
            Sigma_arr = Par.sigma * ones(size(Sigma_arr));                       % First Iteration use the input noise parameter
        end
    end    
    fprintf('Sparsity and basis updating...  ');
     tic;
     [TTalpha, PTP, Lambda,D]  =  Sparsityreg( NL_mat, Self_arr, Sigma_arr, Pat_g_hat, Par ,Pat_f,Lambda);   % Sparsity regularization
     t1=toc;
    fprintf('Computation time: %f \n',t1);
     
    fprintf('TV...  ');
     tic;
     [R, Rv] = TVpreprocess( TTalpha, PTP,Par, Height, Width );  
     f=TVreg(Rv,R,Par.mu*R,0.001);  %TV regularization   %100 0.5; 50 0.2;
     t2=toc;
    fprintf('Computation time: %f \n',t2);

     error=sum((f(:)-f_old(:)).^2)/sum(f_old(:).^2);
     f_old=f;       
     PSNR(k)  = computePSNR(f(1:end-Par.patsize+1,1:end-Par.patsize+1),Ima_ori); 
    fprintf( 'k = %d, error=%2.7f , PSNR = %2.2f \n', k, error, PSNR(k));
    
    figure(1),
    subplot(1,2,2),imshow(f(1:end-Par.patsize+1,1:end-Par.patsize+1),[]);
    title(num2str(k));drawnow;  
     if error<1e-6 
           break;
     end   
    if k>2 & PSNR(k)-PSNR(k-1)<0.001
         break;
     end
end
return;


