tv_l2_decom.m.bak

纹理提取的matlab 源代码, 研究图像处理算法的可参考

% 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 

function [u,v,u_1] = TV_L2_decom(img_input_path)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT u0 = noisy image 
% OUTPUT: u(i,j) = denoised image 
% read the noisy data u0 in BMP format

u0=double(imread(img_input_path));
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.05;%028;

% space discretization 
h=1;

% number of iterations (depends on the image) 
IterMax=20;

% 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');

u_1 = u0-v;
shown(u+v); title('u+v');
shown(u0-(u+v)); title('f-u-v');
%figure,plot(Energy);title('Energy/Iterations');