图像复原MATLAB实现

前言：本篇博客先介绍滤波器滤除噪声，再介绍滤波器复原，侧重于程序的实现。

一：三种常见的噪声

二：空间域滤波

空间域滤波复原是在已知噪声模型的基础上，对噪声的空间域进行滤波。
空间域滤波复原方法主要包括：

　　均值滤波器

　　　　算术均值滤波器

　　　　几何均值滤波器
　　　　谐波均值滤波器
　　　　逆谐波均值滤波器
　　顺序统计滤波器
　　　　中值滤波器
　　　　最大值/最小值滤波器

2.1算数均值滤波器

1 img=imread('D:/picture/ZiXia.jpg');
2 img=rgb2gray(img);
3 figure,imshow(img);//原图
4 img_noise=double(imnoise(img,'gaussian',0.06));
5 figure,imshow(img_noise,[]);//含有高斯噪声的图
6 img_mean=imfilter(img_noise,fspecial('average',3));//滤波后的图
7 figure;imshow(img_mean,[]);

2.2几何均值滤波器

1 img=imread('cameraman.tif');
2 img=rgb2gray(img);
3 figure,imshow(img);
4 img_noise=double(imnoise(img,'gaussian',0.06));
5 figure,imshow(img_noise,[]);
6 img_mean=exp(imfilter(log(img_noise+1),fspecial('average',3)));
7 figure;imshow(img_mean,[]);

2.3谐波均值滤波器

2.4逆谐波均值滤波器

采用逆谐波均值滤波器对附加胡椒噪声图像进行滤波的matlab程序如下：

1 img=imread('cameraman.tif');  figure,imshow(img);
2 [M,N]=size(img);R=imnoise2('salt & pepper',M,N,0.1,0);
3 img_noise=img;img_noise(R==0)=0;
4 img_noise=double(img_noise); figure,imshow(img_noise,[]);
5 Q=1.5;
6 img_mean=imfilter(img_noise.^(Q+1),fspecial('average',3))./imfilter(img_noise.^Q,fspecial('average',3));
7 figure;imshow(img_mean,[]);

采用逆谐波均值滤波器对附加盐噪声图像进行滤波的matlab程序如下：

1 img=imread('csboard.tif');figure,imshow(img);
2 [M,N]=size(img);R=imnoise2('salt & pepper',M,N,0,0.1);
3 img_noise=img;img_noise(R==1)=255;
4 img_noise=double(img_noise); figure,imshow(img_noise,[]);
5 Q=-1.5;
6 img_mean=imfilter(img_noise.^(Q+1),fspecial('average',3))./imfilter(img_noise.^Q,fspecial('average',3));
7 figure;imshow(img_mean,[]);

2.5中值滤波器

1 img=imread('cameraman.tif');
2 img_noise=double(imnoise(img,'salt & pepper',0.06));
3 img_mean=imfilter(img_noise,fspecial('average',5));
4 img_median=medfilt2(img_noise);%一次中值滤波
5 img_median2=medfilt2(img_median);%二次中值滤波

2.6最大值，最小值滤波器

利用最大值滤波器消除胡椒噪声污染图像的matlab程序如下。

1 img=imread('csboard.tif');
2 [M,N]=size(img);
3 R=imnoise2('salt & pepper',M,N,0.1,0);
4 img_noise=img;
5 img_noise(R==0)=0;
6 img_noise=double(img_noise);
7 imwrite(uint8(img_noise),'csbord_pepper.jpg');
8 img_max=imdilate(img_noise,ones(3,3));
9 imwrite(uint8(img_max),'cameraman_saltpepper_max.jpg');

利用最小值滤波器消除盐噪声污染图像的matlab程序如下。

1 img=imread('csboard.tif');
2 [M,N]=size(img);
3 R=imnoise2('salt & pepper',M,N,0,0.1);
4 img_noise=img;
5 img_noise(R==1)=255;
6 img_noise=double(img_noise);
7 imwrite(uint8(img_noise),'csbord_salt.jpg');
8 img_min=imerode(img_noise,ones(3,3));
9 imwrite(uint8(img_min),'cameraman_saltpepper_min.jpg');

2.7带阻滤波器

 1 I=imread('pout.tif');
 2 [m,n]=size(I);
 3 J=I;
 4 for i=1:m
 5     for j=1:n
 6         J(i,j)=I(i,j)+20*sin(20*i)+20*sin(20*j);%增加周期性噪声
 7     end
 8 end
 9 IF=fftshift(fft2(I));
10 JF=fftshift(fft2(J));
11 IF=log(1+abs(IF));
12 JF=log(1+abs(JF));
13 
14 subplot(221)%显示频谱
15 imshow(IF,[])
16 subplot(222)
17 imshow(JF,[])
18 %高斯带阻滤波器构造
19 fbrf=ones(m,n);
20 for i=1:m
21     for j=1:n
22         fbrf(i,j)=1-exp(-0.5*(((i-m/2)^2+(j-n/2)^2-55^2)/(sqrt(i.^2+j.^2)*5))^2);%20为带阻中心，5为带宽
23     end
24 end
25  H=fbrf;
26 %频率域滤波
27 f=fftshift(fft2(J));
28 out=f.*H;%频率域滤波
29 out=ifft2(ifftshift(out));
30 out=abs(out);
31 out=out/max(out(:));%归一化【0,1】
32 subplot(223)
33 imshow(out,[]);
34 subplot(224)
35 imshow(J,[])

2.8带通滤波器

2.9陷波滤波器

 1 clear;
 2 close all;
 3 src = im2double(imread('D:/picture/ZiXia.jpg'));
 4 src = rgb2gray(src);
 5 subplot(221)
 6 imshow(src);
 7 title('原始图像');
 8 [w h] = size(src);
 9 srcf = fft2(src);
10 srcf = fftshift(srcf);
11 subplot(222)
12 imshow(srcf);
13 %  低通滤波
14 % flt = zeros(size(src));
15 % rx1 = w/2;
16 % ry1 = h/2;
17 % r = min(w,h)/3;
18 % for i = 1:w
19 %     for j = 1:h
20 %         if(rx1-i)^2 +(ry1 - j)^2 <= r*r
21 %             flt(i,j) = 1;
22 %         end
23 %     end
24 % end
25 % 陷波滤波
26 flt = ones(size(src));
27 r = min(w,h)/12;
28 rx1 = r
29 ry1 =h/2
30 for i = 1:w
31     for j = 1:h
32         if(rx1-i)^2 +(ry1 - j)^2 <= r*r
33             flt(i,j) = 0;
34         end
35         if(w-rx1-i)^2 +(h-ry1 - j)^2 <= r*r
36             flt(i,j) = 0;
37         end
38     end
39 end
40 subplot(223)
41 imshow(flt);
42 title('滤波器图像');
43 dfimg = srcf.*flt;
44 dfimg = ifftshift(dfimg);
45 dimg = ifft2(dfimg,'symmetric');
46 subplot(224)
47 imshow(dimg):title('滤波后');

2.9逆滤波

 1 image_o=imread('D:/picture/lenagray.jpg');
 2 subplot(1,3,1);
 3 imshow(image_o);
 4 title('原图像');
 5 %频率域退化图像，退化函数H(u,v)=exp(-0.0025*( (u-M/2).^2+(v-N/2).^2).^(5/6) )
 6 %傅里叶变换
 7 f=im2double(image_o);
 8 F=fft2(f);
 9 F=fftshift(F);
10 %执行退化
11 [M,N]=size(F);
12 [u,v]=meshgrid(1:M,1:N);%生成二维坐标系
13 H=exp(-0.0025* ( (u-M/2).^2+(v-N/2).^2).^(5/6) );
14 F=F.*H;
15 %傅里叶反变换
16 X=ifftshift(F);
17 x=ifft2(X);
18 x=uint8(abs(x)*256);
19 subplot(1,3,2);
20 imshow(x);
21 %
22 title('退化图像');
23 
24 image_d=x;
25 %直接逆滤波图像复原
26 
27 ff=im2double(image_d);%将图像灰度值归一化到0-1之间
28 
29 % 傅里叶变换
30 f_Id=fft2(ff);
31 f_Id=fftshift(f_Id);
32 fH_Id=f_Id;
33 [M,N]=size(fH_Id);
34 % 逆滤波
35 threshold=78;
36 if threshold>M/2
37         %全滤波
38         fH_Id=fH_Id./(H+eps);
39 else
40         %对一定半径范围内进行滤波
41         for i=1:M
42             for j=1:N
43                 if sqrt((i-M/2).^2+(j-N/2).^2)<threshold
44                     fH_Id(i,j)=fH_Id(i,j)./(H(i,j)+eps);
45                 end
46             end
47         end
48 end
49 
50 % 执行傅立叶逆变换
51 fH_Id1=ifftshift(fH_Id);
52 f_new=ifft2(fH_Id1);
53 f_new=uint8(abs(f_new)*255);
54 subplot(1,3,3);
55 imshow(f_new);
56 title('滤波半径=78的逆滤波复原图像');

2.10维纳滤波

直接截图了，没敲