紫馨86
2009-04-16, 15:19
clear all;
close all;
I=imread('lena.bmp');
X=rgb2gray(I);
[height width]=size(X);
X=im2double(X);
A=imnoise(X,'gaussian',0,0.01);
wname='sym4';
n=3;
[c,s]=wavedec2(A,n,wname);
for i=1:3
step(i)=s((i+1),1)*s((i+1),2);
end
num(1,1)=s(1,1)*s(1,2)+1; %获取各层各高频分量在c向量中的坐标
num(1,2)=num(1,1)+s(2,1)*s(2,2);
num(1,3)=num(1,2)+s(2,1)*s(2,2);
num(2,1)=num(1,3)+s(2,1)*s(2,2);
num(2,2)=num(2,1)+s(3,1)*s(3,2);
num(2,3)=num(2,2)+s(3,1)*s(3,2);
num(3,1)=num(2,3)+s(3,1)*s(3,2);
num(3,2)=num(3,1)+s(4,1)*s(4,2);
num(3,3)=num(3,2)+s(4,1)*s(4,2);
m=0.02;
C=c;
for i=1:3
[H,V,D]=detcoef2('a',c,s,i);%提取第i层各高频系数
B=[H V D];
[L,T]=size(B);
for k=1:L
for w=1:T
sigma=median(abs(B(k,w)))/0.6745;
end
end
th=sigma*sqrt(2*log10(L*T));%计算阈值
ch=c(1,num(4-i,1):num(4-i,3)+step(4-i)-1);%对各高频系数进行阈值处理
chl=length(ch);
for j=1:chl
if abs(ch(j))>=th
ch(j)=sign(ch(j))*(abs(ch(j))-2*th/(1+exp(m*(ch(j)^2))));%阈值处理函数
else
ch(j)=0;
end
end
C(1,num(4-i,1):num(4-i,3)+step(4-i)-1)=ch(1,1:chl);
end
X0=waverec2(C,s,wname);
disp('原图像信噪比为：'),sny=10*log10(sum(sum(double(X).^2))/sum(sum((double(X)-double(A)).^2)))
disp('去噪后图像信噪比为：'),sndz=10*log10(sum(sum(double(X).^2))/sum(sum((double(X)-double(X0)).^2)))
figure(1)
imshow(X)
title('原图像')
figure(2)
imshow(A)
title('添加噪声后的图像')
figure(3)
imshow(X0,[])
title('去噪后图像')
close all;
I=imread('lena.bmp');
X=rgb2gray(I);
[height width]=size(X);
X=im2double(X);
A=imnoise(X,'gaussian',0,0.01);
wname='sym4';
n=3;
[c,s]=wavedec2(A,n,wname);
for i=1:3
step(i)=s((i+1),1)*s((i+1),2);
end
num(1,1)=s(1,1)*s(1,2)+1; %获取各层各高频分量在c向量中的坐标
num(1,2)=num(1,1)+s(2,1)*s(2,2);
num(1,3)=num(1,2)+s(2,1)*s(2,2);
num(2,1)=num(1,3)+s(2,1)*s(2,2);
num(2,2)=num(2,1)+s(3,1)*s(3,2);
num(2,3)=num(2,2)+s(3,1)*s(3,2);
num(3,1)=num(2,3)+s(3,1)*s(3,2);
num(3,2)=num(3,1)+s(4,1)*s(4,2);
num(3,3)=num(3,2)+s(4,1)*s(4,2);
m=0.02;
C=c;
for i=1:3
[H,V,D]=detcoef2('a',c,s,i);%提取第i层各高频系数
B=[H V D];
[L,T]=size(B);
for k=1:L
for w=1:T
sigma=median(abs(B(k,w)))/0.6745;
end
end
th=sigma*sqrt(2*log10(L*T));%计算阈值
ch=c(1,num(4-i,1):num(4-i,3)+step(4-i)-1);%对各高频系数进行阈值处理
chl=length(ch);
for j=1:chl
if abs(ch(j))>=th
ch(j)=sign(ch(j))*(abs(ch(j))-2*th/(1+exp(m*(ch(j)^2))));%阈值处理函数
else
ch(j)=0;
end
end
C(1,num(4-i,1):num(4-i,3)+step(4-i)-1)=ch(1,1:chl);
end
X0=waverec2(C,s,wname);
disp('原图像信噪比为：'),sny=10*log10(sum(sum(double(X).^2))/sum(sum((double(X)-double(A)).^2)))
disp('去噪后图像信噪比为：'),sndz=10*log10(sum(sum(double(X).^2))/sum(sum((double(X)-double(X0)).^2)))
figure(1)
imshow(X)
title('原图像')
figure(2)
imshow(A)
title('添加噪声后的图像')
figure(3)
imshow(X0,[])
title('去噪后图像')