xlili.m

一个简单的图象去噪的算法

clear
%
  n=64
u0 = zeros(n,n);
    p1=round(44*n/64); p2=round(54*n/64); p3=max(round(15*n/64),1);
    p4=round(49*n/64);
    p0=mean([p1 p2]);
    q0=mean([p3 p4]);
    [x,y]=meshgrid(p1:p2,p3:p4);
  u0(p1:p2,p3:p4) = 99*ones(size(u0(p1:p2,p3:p4)));
 for i=12:round(n/2+1)
     for j=2:n
      if j>=i & j<=(n-i+1)
        u0(i,j) = 90;
      end
     end
 end
  randn('seed',0); h=1/(n-1); sigma=0.2*max(max(u0));
  noise = randn(size(u0));
  noise = noise/(norm(noise,'fro')*h);
  z = u0 + sigma*noise;  

  alpha=40;
  beta=10^(-4);
  nit=100;
  tol=10^(-4);
    
  tic;
   [u2,iter2,res2] =cgm_2d2(z, alpha,beta, nit, tol);
  t1=toc;
  fprintf('t= %.8e\n', t1);
  fprintf(' iter2= %3d\n',iter2);
  
  z=reshape(z,n,n);
   u2=reshape(u2,n,n);
  

 j0 = j2(z,z,alpha);
  jj2 = j2(u2,z,alpha); 
  
  figure (1)
  imagesc(z)
  ylabel(['J = ' num2str(j0) ]), title('Original data u=z'), a=axis;
  figure (3)
  imagesc(u2)
  ylabel(['J = ' num2str(jj2) ]), a=axis; 
  a=title(['cgm\_2d Step ' num2str(iter2) ', \alpha =' num2str(alpha)]);
  s=get(a,'Fontsize'); set(a,'Fontsize',s); set(a,'Color','r')
  
%   b1=sum(sum(z))/(n^2);b2=sum(sum(z-u0))/(n^2);
%   b3=sum(sum((z-b1).^2));b4=sum(sum((z-u0-b2).^2));
%   b5=b3/b4;snr=10*log(10*b5)
%   a1=sum(sum(u2))/(n^2);a2=sum(sum(u2-u0))/(n^2);
%   a3=sum(sum((u2-a1).^2));a4=sum(sum((u2-u0-a2).^2));
%   a5=a3/a4;snr1=10*log(10*a5)
%   rmse=sqrt(sum(sum((z-u0).^2)))/(n^2)
%   rmse1=sqrt(sum(sum((u2-u0).^2)))/(n^2)

  s=sum(sum(u0));  %suo you yuan su xiang jia
   s=s/(n*n);
   for i=1:n
   for j=1:n
   A(i,j)=u0(i,j)-s;
   end
   end
   A=A.^2;
   s1=sum(sum(A));
 B=sigma*noise;
   B=B.^2;
   s2=sum(sum(B));

   snr=s1/s2  
  aa=sum(sum(u2));
aa=aa/(n^2);
for i=1:n
   for j=1:n
       AA(i,j)=u2(i,j)-aa;
   end
end
AA=AA.^2;
ss1=sum(sum(AA));
BB=u2-u0;
s2=sum(sum(BB.^2));
snr2=ss1/s2
fprintf('noise= %.8e\n',norm(BB,'fro'))