www.pudn.com > texture_extract.rar > TV_L2_decombak, change:2006-02-14,size:6192b


% Matlab code to denoise an image by the Rudin-Osher-Fatemi model 
% Min TV(u) + lambda*(L^2[(f-u)])^2 
% lambda = coefficient of the fidelity term  
% lambda needs to be tunned to obtain a sharp denoised image 
% h = space discretization step  
% The stationary problem (no time-regularization) is used here 
% A fixed point iteration is used  
% The numerical scheme is implicit in the max norm and unconditionally stable 
% A stopping criteria can be included, or the "IterMax" must be given.  
 
% LAST MODIFIED: 2005-04-19  
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% INPUT u0 = noisy image  
% OUTPUT: u(i,j) = denoised image  
% read the noisy data u0 in BMP format 
 
u0=double(imread('barbara512.bmp')); 
u0=double(u0);%(1:200,313:512)); 
%  
% std_n=20; % Gaussian noise standard deviation 
% In = randn(size(u0))*std_n; % White Gaussian noise 
% u0 = u0 + In;  % noisy input image 
 
[M N]=size(u0); 
 
% visualize the image u0 in Matlab (rescaled)  
%imagesc(u0); axis image; axis off; colormap(gray); 
%shown(u0); 
 
%--------------------------------------------------------------------------- 
% initialize u by u0 (not necessarily) or by a constant, or by a random image 
%--------------------------------------------------------------------------- 
 
u=u0; 
[M,N]=size(u); 
 
%-----------------------------------------------  
%           PARAMETERS   
%----------------------------------------------- 
 
% coefficient of the TV norm (needs to be adapted for each image)  
lambda=0.1;%028; 
 
% space discretization  
h=1; 
 
% number of iterations (depends on the image)  
IterMax=10; 
 
% needed to regularize TV at the origin  
eps=1;%0.000001; 
 
% norm 
p=1; 
mu=0.01; 
 
[fx,fy] = gradient(u0); 
g1 = -1/(2*lambda) * fx./sqrt(eps*eps+fx.*fx+fy.*fy); 
g2 = -1/(2*lambda) * fy./sqrt(eps*eps+fx.*fx+fy.*fy); 
 
%----------------------------------------------------------- 
%     BEGIN ITERATIONS IN ITER  
%----------------------------------------------------------- 
for Iter=1:IterMax,      
    Iter     
    for i=2:M-1, 
        for j=2:N-1,             
            %-----------computation of coefficients co1,co2,co3,co4---------             
            ux=(u(i+1,j)-u(i,j))/h; 
            uy=(u(i,j+1)-u(i,j-1))/(2*h); 
            Gradu=sqrt(eps*eps+ux*ux+uy*uy); 
            co1=1./Gradu;             
             
            ux=(u(i,j)-u(i-1,j))/h; 
            uy=(u(i-1,j+1)-u(i-1,j-1))/(2*h); 
            Gradu=sqrt(eps*eps+ux*ux+uy*uy); 
            co2=1./Gradu; 
             
            ux=(u(i+1,j)-u(i-1,j))/(2*h); 
            uy=(u(i,j+1)-u(i,j))/h; 
            Gradu=sqrt(eps*eps+ux*ux+uy*uy); 
            co3=1./Gradu; 
             
            ux=(u(i+1,j-1)-u(i-1,j-1))/(2*h); 
            uy=(u(i,j)-u(i,j-1))/h; 
            Gradu=sqrt(eps*eps+ux*ux+uy*uy); 
            co4=1./Gradu;             
             
            co=1.+(1/(2*lambda*h*h))*(co1+co2+co3+co4); 
            div = co1*u(i+1,j)+co2*u(i-1,j)+co3*u(i,j+1)+co4*u(i,j-1);                            
             
            g1x = (g1(i+1,j) - g1(i-1,j))/(2*h); 
            g2y = (g2(i,j+1) - g2(i,j-1))/(2*h);                                                                        
            u(i,j)=(1./co)*(u0(i,j)-g1x-g2y+(1/(2*lambda*h*h))*div);             
             
            Hg = norm(sqrt(g1(i,j)*g1(i,j) + g2(i,j)*g2(i,j)),p).^(1-p) .* (sqrt(g1(i,j)*g1(i,j) + g2(i,j)*g2(i,j))).^(p-2);                         
            gc = 2*lambda/(mu*Hg+4*lambda/(h*h));             
            g1xy = 1/(2*h*h)*(2*g1(i,j)+g1(i-1,j-1)+g1(i+1,j+1)-g1(i,j-1)-g1(i-1,j)-g1(i+1,j)-g1(i,j+1));            
            g2xy = 1/(2*h*h)*(2*g2(i,j)+g2(i-1,j-1)+g2(i+1,j+1)-g2(i,j-1)-g2(i-1,j)-g2(i+1,j)-g2(i,j+1));                                                                        
            g1(i,j) = gc*( (u(i+1,j)-u(i-1,j))/(2*h) - (u0(i+1,j)-u0(i-1,j))/(2*h) + (g1(i+1,j)+g1(i-1,j))/(h*h) + g2xy ); 
            g2(i,j) = gc*( (u(i,j+1)-u(i,j-1))/(2*h) - (u0(i,j+1)-u0(i,j-1))/(2*h) + (g2(i,j+1)+g2(i,j-1))/(h*h) + g1xy );                                                             
        end 
    end     
     
    %------------ FREE BOUNDARY CONDITIONS IN u ------------------- 
    for i=2:M-1, 
        u(i,1)=u(i,2); 
        u(i,N)=u(i,N-1); 
         
        g1(i,1)=g1(i,2); 
        g1(i,N)=g1(i,N-1);         
 
        g2(i,1)=g2(i,2); 
        g2(i,N)=g2(i,N-1);                 
    end 
     
    for j=2:N-1, 
        u(1,j)=u(2,j); 
        u(M,j)=u(M-1,j); 
 
        g1(1,j)=g1(2,j); 
        g1(M,j)=g1(M-1,j); 
         
        g2(1,j)=g2(2,j); 
        g2(M,j)=g2(M-1,j);         
    end 
     
    u(1,1)=u(2,2); 
    u(1,N)=u(2,N-1);  
    u(M,1)=u(M-1,2); 
    u(M,N)=u(M-1,N-1); 
     
    g1(1,1)=g1(2,2); 
    g1(1,N)=g1(2,N-1);  
    g1(M,1)=g1(M-1,2); 
    g1(M,N)=g1(M-1,N-1); 
 
    g2(1,1)=g2(2,2); 
    g2(1,N)=g2(2,N-1);  
    g2(M,1)=g2(M-1,2); 
    g2(M,N)=g2(M-1,N-1);     
    %---------------------------------------------------------------------- 
     
    %---------------- END ITERATIONS IN Iter ------------------------------ 
     
    %%% Compute the discrete energy at each iteration 
%     en=0.0;   
%     for i=2:M-1, 
%         for j=2:N-1, 
%             ux=(u(i+1,j)-u(i,j))/h; 
%             uy=(u(i,j+1)-u(i,j))/h; 
%             fid=(u0(i,j)-u(i,j))*(u0(i,j)-u(i,j)); 
%             en=en+sqrt(eps*eps+ux*ux+uy*uy)+lambda*fid; 
%         end 
%     end 
    %%% END computation of energy  
     
%     Energy(Iter)=en;    
     
    if mod(Iter,10)==0 
        name=['Iter_L2_',num2str(Iter),'.mat']; 
        save(name); 
    end 
%     if Iter>=2 
%         if Energy(Iter)>Energy(Iter-1) 
%             break; 
%         end 
%     end 
end  
 
 
% visualize the image in Matlab (re-scaled) 
%figure  
%imagesc(u); axis image; axis off; colormap(gray); 
shown(u0); title('u0'); 
shown(u); title('u'); 
 
% use export to save the image in ps or eps format  
% Plot the Energy versus iterations  
[g1x,g1y] = gradient(g1); 
[g2x,g2y] = gradient(g2); 
v = g1y+g2x; 
shown(v); title('v'); 
shown(u+v); title('u+v'); 
shown(u0-(u+v)); title('f-u-v'); 
%figure,plot(Energy);title('Energy/Iterations');