第五章 图像复原
一、图像退化、复原处理模型
在空间域,h(x,y)称为点扩散函数(PSF),对于任何种类的输入,让h(x,y)作用于光源的一个点来得到退化的特征。
二、噪声模型
模拟噪声的行为和影响的能力是图像复原的核心。
g = imnoise(f,type,parameters)
典型的函数语法形式:
g = imnoise(f,'gaussian',m,var) %将均值为m、方差为var的高斯噪声加到图像f上,默认值为均值是0,方差是0.01的噪声 g = imnoise(f,'salt & pepper',d) %用椒盐噪声污染图像f,其中d是噪声密度 g = imnoise(f,'poisson') %从数据中生成泊松噪声,而不是人工添加
2.2、使用指定的分布产生空间随机噪声
g = imnoise2(type,M,N,a,b)
与imnoise不同,下面的M函数产生一个大小为M*N的噪声数组R,它不以任何方式缩放;imnoise输出一个有噪声的图像,而imnoise2产生噪声模式本身。
注意:为了使这个数组有用,需要对它进行进一步的处理,若要用这个数组来污染一幅图像,则我们可使用函数find寻找R中所有值为0的坐标,并把图像中相应得坐标置为可能的最小灰度值(通常为0),同样我们可以使用函数find寻找R中所有值为1的坐标,并把图像中相应得坐标置为可能的最高灰度值(8bit通常为255)。
模拟产生各种噪声
clc clear r1=imnoise2('gaussian',100000,1,1);%imnoise2是产生噪声矩阵,这里是产生高斯噪声,矩阵大小为10000*1 bins=100; %均值为0,方差为1 figure,hist(r1,bins);%将r矩阵中的数值直方图表示出来,共100个bin title('guassian'); r2=imnoise2('uniform',1);%同上,此处产生的是0~1的均匀分布噪声,矩阵大小为100000*1 figure,hist(r2,bins); title('uniform'); r3=imnoise2('salt & pepper',1000,0.1,0.27);%注意这里的salt & pepper中间记得留空格。 figure,hist(r3,bins); title('salt & pepper'); r4=imnoise2('lognormal',1);%产生对数正态噪声,大小100000*1,偏移值默认为1,形状参数默认0.25 figure,hist(r4,bins); title('lognormal'); r5=imnoise2('rayleigh',1);%rayleigh噪声,该噪声好像是对整体矢量合成的一种分布 figure,hist(r5,bins); title('rayleigh'); r6=imnoise2('exponential',1);%指数分布的参数默认为1 figure,hist(r6,bins); title('exponential'); % figure,subplot(321),imshow(r1),title('高斯噪声'); % subplot(322),imshow(r2),title('均匀分布噪声'); % subplot(323),imshow(r3),title('椒盐噪声'); % subplot(324),imshow(r4),title('对数正态噪声'); % subplot(325),imshow(r5),title('瑞利噪声'); % subplot(326),imshow(r6),title('指数噪声');
模拟结果:
三、仅有噪声的复原:空间滤波
g(x,y)=f(x,y)+η(x,y)
3.1、空间噪声滤波器
%% 空间噪声滤波 clc clear f=imread('.\images\dipum_images_ch05\Fig0507(a)(checkeboard8_pixeldup_8).tif'); imshow(f),title('原始图像'); [M N]=size(f); R=imnoise2('salt & pepper',0);%仅仅产生0.1的椒噪声 c=find(R==0);%把R中产生的椒点找出来 gp=f; gp(c)=0;%在原始图像对应椒点的位置赋值为0 % figure,imshow(gp); % title('被概率为0.1的胡椒噪声污染的图像'); R=imnoise2('salt & pepper',0.1);%仅仅产生0.1的盐噪声 c=find(R==1);%把R中产生的盐点找出来 gs=f; gs(c)=255;%在原始图像对应盐点的位置赋值为1 % figure,imshow(gs); % title('被概率为0.1的盐噪声污染的图像'); fp=spfilt(gp,'chmean',3,1.5);%反调和滤波,尺寸大小为3*3,Q默认为1.5,正的Q过滤胡椒噪声 % figure,imshow(fp); % title('正Q反调和滤波器滤波的结果') fs=spfilt(gs,-1.5);%负的Q过滤盐粒噪声 % figure,imshow(fs); % title('负Q反调和滤波器滤波的结果') fpmax=spfilt(gp,'max',3); % figure,imshow(fpmax); % title('最大值滤波后的结果'); fsmin=spfilt(gs,'min',imshow(fsmin); % title('最小值滤波后的结果'); figure,imshow(gp);title('被概率为0.1的胡椒噪声污染的图像'); subplot(322),imshow(gs);title('被概率为0.1的盐噪声污染的图像'); subplot(323),imshow(fp);title('正Q反调和滤波器滤波的结果') subplot(324),imshow(fs);title('负Q反调和滤波器滤波的结果') subplot(325),imshow(fpmax);title('最大值滤波后的结果'); subplot(326),imshow(fsmin);title('最小值滤波后的结果');
实验结果:
3.2、自适应空间滤波器
f = adpmedian (g,Smax)
%% 自适应中值滤波器 clc clear f=imread('.\images\dipum_images_ch05\Fig0507(a)(checkeboard8_pixeldup_8).tif'); figure,subplot(221),imshow(f),title('原始图像'); g=imnoise(f,0.25);%噪声点有白有黑,因为这是用的imnoise,不是imnoise2和imnoise3 subplot(222),imshow(g); title('被概率为0.25的椒盐噪声污染过的图像'); f1=medfilt2(g,[7 7],'symmetric'); %边缘扩展模式为对称扩展 subplot(223),imshow(f1); title('用普通中值滤波器滤波后的图像'); f2=adpmedian(g,7); subplot(224),imshow(f2); title('用Smax=7的自适应中值滤波器滤波后的图像');
实验结果:
四、退化函数建模和维纳滤波
%% 模糊噪声图像建模 clc clear f=checkerboard(8);%直接生产黑白棋盘,每个小正方形一边的像素点为8个 figure,imshow(f); title('原始图像'); PSF=fspecial('motion',7,45);%创建一个motion滤波器,长度为7,仰角为45度。 gb=imfilter(f,PSF,'circular');%具体motion滤波器是什么还不是很清楚。 subplot(222),imshow(gb); title('使用motion滤波器模糊后'); noise=imnoise(zeros(size(f)),0.001);%产生高斯噪声,且后面要用到该噪声图像 subplot(223),imshow(noise,[]); % subplot(223),imshow(noise); %全是黑的 title('高斯纯噪声') g=gb+noise; subplot(224),imshow(g,[]);%模糊加噪声后的图像 title('模糊加噪声后的图像'); %========================================================================== %%使用deconvwnr函数复原模糊噪声图像 figure,[]);%模糊加噪声后的图像 title('模糊加噪声后的图像'); fr1=deconvwnr(g,PSF);%维纳滤波,去模糊化 subplot(222),imshow(fr1,[]); title('简单维纳滤波后的结果'); Sn=abs(fft2(noise)).^2;%噪声功率谱 nA=sum(Sn(:))/prod(size(noise));%平均噪声功率谱,prod是元素相乘的意思,这里就是noise的长和宽相乘 Sf=abs(fft2(f)).^2;%信号功率谱 fA=sum(Sf(:))/prod(size(f));%平均信号功率谱 R=nA/fA;%求得信噪比 fr2=deconvwnr(g,R);%信噪比为常数的参数维纳滤波器 subplot(223),imshow(fr2,[]); title('常数比率信噪比维纳滤波器'); NCORR=fftshift(real(ifft(Sn)));%噪声的自相关函数,why? ICORR=fftshift(real(ifft(Sf)));%信号的自相关函数,why? fr3=deconvwnr(g,NCORR,ICORR);%自相关后的维纳滤波 subplot(224),imshow(fr3,[]); title('使用自相关函数的维纳滤波后');
实验结果:
退化函数建模
维纳滤波
五、约束的最小二乘方(正则)滤波
fr = deconvreg(g,NOTSEPOWER,RANGE)
%(约束最小二乘方)正则滤波 clc clear f=checkerboard(8);%直接生产黑白棋盘,每个小正方形一边的像素点为8个 PSF=fspecial('motion','circular');%具体motion滤波器是什么还不是很清楚。 noise=imnoise(zeros(size(f)),0.001);%产生高斯噪声,且后面要用到该噪声图像 g=gb+noise; figure,subplot(131),[]);%模糊加噪声后的图像 title('模糊加噪声后的图像'); fr1=deconvreg(g,4); subplot(132),[]); title('噪声功率为4的正则滤波后结果'); fr2=deconvreg(g,0.4,[1e-7,1e7]); subplot(133),[]); title('带搜索范围的正则滤波后结果');
实验结果:
六、总结
本章描述了怎样用Matlab和IPT函数进行图像复原的方法,以及怎样将其做描述退化图像而生成模型的方法。函数imnoise2和imnoise3的开发对产生噪声的能力有了显著地提高,同样,函数spfilt可以实现的空间过滤。
附录:IPT复原代码