数字图像处理6

数字图像处理第6章图像恢复

6.1退化模型及恢复技术基础

6.2空间域滤波恢复

6.3频率域滤波恢复6.4逆滤波

6.5最小均方误差滤波器-维纳滤波

前言

图像恢复和图像增强一样，都是为了改善图像视觉效果，以及便于后续处理。只是图像增强方法更偏向主观判断，而图像恢复则是根据图像畸变或退化原因，进行模型化处理。本章安排如下：首先介绍退化模型及恢复技术基础；随后分析空间域滤波恢复技术；最后介绍频率域图像恢复技术及应用。

6.1退化模型及恢复技术基础-

退化的原因

成像系统的象差、畸变、带宽有限等造成图像图像失真；由于成象器件拍摄姿态和扫描非线性引起的图像几何失真；运动模糊，成象传感器与被拍摄景物之间的相对运动，引起所成图像的运动模糊；灰度失真，光学系统或成象传感器本身特性不均匀，造成同样亮度景物成象灰度不同；辐射失真，由于场景能量传输通道中的介质特性如大气湍流效应、大气成分变化引起图像失真；图像在成象、数字化、采集和处理过程中引入的噪声等。6.1退化模型及恢复技术基础-退化模型

实际所得退化图像可简化为6.1退化模型及恢复技术基础-

恢复技术的概念及分类

定义：图像恢复是根据退化原因，建立相应的数学模型，从被污染或畸变的图像信号中提取所需要的信息，沿着使图像降质的逆过程恢复图像本来面貌。6.1退化模型及恢复技术基础-

恢复技术的概念及分类图像恢复技术的分类：(1)在给定退化模型条件下，分为无约束和有约束两大类;(2)根据是否需要外界干预，分为自动和交互两大类;(3)根据处理所在域，分为频域和空间域两大类。6.2噪声模型（NoiseModels）

数字图像的噪声的来源噪声的频率特性噪声的空间特性6.2.1噪声类型-高斯噪声定义：6.2.1噪声类型-均匀分布噪声

定义：6.2.1噪声类型-脉冲噪声（椒盐噪声）

定义例6.1A=imread('fig606a.jpg');%读取图像figure,imshow(A);%显示图像A=rgb2gray(A);figure,hist(double(A),10);%求出A的直方图并显示B=imnoise(A,'gaussian',0.05);%对A附加高斯噪声figure,imshow(B);%显示附加高斯噪声后的图像Bfigure,hist(double(B),10);%求出B的直方图并显示C=imnoise(A,'speckle',0.05);%对A附加均匀分布噪声figure,imshow(C);%显示附加均匀噪声后的图像Cfigure,hist(double(C),10);%求出C的直方图并显示D=imnoise(A,'salt&pepper',0.05);%对A附加椒盐噪声figure,imshow(D);%显示附加椒盐噪声后的图像Dfigure,hist(double(D),10);%求出D的直方图并显示6.2.1噪声类型-示例(a)原图(b)高斯噪声图6.2.1噪声类型-示例(c)均匀分布噪声(d)椒盐噪声6.3空间域滤波恢复定义：空间域滤波恢复即是在已知噪声模型的基础上，对噪声的空间域滤波6.3.1均值滤波器

采用均值滤波模板对图像噪声进行滤除6.2.2均值滤波-类型算术均值滤波器：几何均值滤波器

6.2.2均值滤波-类型谐波均值滤波器

逆谐波均值滤波器

6.2.2均值滤波-示例(a)输入图像；(b)高斯噪声污染图像；(c)用均值滤波结果6.2.2均值滤波-示例

(d)几何均值滤波（e）Q＝－1.5的逆谐波滤波(f)Q=1.5滤波的结果6.3.2顺序统计滤波

1.中值滤波

其中，其中，g为输入图像，sxy为滤波窗口。修正后的阿尔法均值滤波器img=imread('cameraman.tif');imshow(img);%显示图像img_noise=double(imnoise(img,'gaussian',0.06));%对图像附加高斯噪声figure,imshow(img_noise,[]);%显示加噪图像img_mean=imfilter(img_noise,fspecial('average',3));%对附加有高斯噪声的图像实行算术均值滤波figure;imshow(img_mean,[]);%显示算术均值滤波后的图像img_mean=exp(imfilter(log(img_noise),fspecial('average',3)));%对附加有高斯噪声的图像实行几何均值滤波figure;imshow(img_mean,[]);%显示几何均值滤波后的图像Q=-1.5;%对高斯噪声图像实行Q取负数的逆谐波滤波img_mean=imfilter(img_noise.^(Q+1),fspecial('average',3))./imfilter(img_noise.^Q,fspecial('average',3));figure;imshow(img_mean,[]);%显示逆谐波滤波后的图像Q=1.5;%对高斯噪声图像实行Q取正数的逆谐波滤波img_mean=imfilter(img_noise.^(Q+1),fspecial('average',3))./imfilter(img_noise.^Q,fspecial('average',3));figure;imshow(img_mean,[]);%显示逆谐波滤波后的图像例6.3img=rgb2gray(imread('football.jpg'));%读取图像并转换成灰度图像figure;imshow(img);%显示图像img_noise=double(imnoise(img,'salt&pepper',0.06));%加椒盐噪声figure,imshow(img_noise,[]);%显示附加椒盐噪声后的图像img_mean=imfilter(img_noise,fspecial('average',5));%对附加有椒盐噪声的图像实行算术均值滤波figure;imshow(img_mean,[]);%显示算术均值滤波后的图像img_median=medfilt2(img_noise);%对附加有椒盐噪声的图像实行中值滤波figure;imshow(img_median,[]);%显示中值滤波后的图像img_median2=medfilt2(img_median);%对中值滤波处理后的图像再次实行中值滤波figure;imshow(img_median2,[]);%显示再次中值滤波后的图像中值滤波示例(a)椒盐噪声污染的图像(b)均值滤波结果；中值滤波示例(续）

(c)中值滤波结果(d)对c图再次中值滤波最大/最小滤波2.最大/最小滤波1）最大值滤波器为：

2）最小值滤波器为：最大/最小滤波示例(a)噪声图像(b)最大滤波结果(c)最小滤波结果6.4频率域滤波恢复

原理：时域卷积相当于频域乘积。因此可以在频率域中直接设计滤波器，对信号进行恢复处理。分类：常用的图像恢复方法有带阻滤波器、带通滤波器、陷波滤波器等6.4.1带阻滤波器

1.理想带阻滤波器2.巴特沃思带阻滤波器6.4.1带阻滤波器高斯带阻滤波器(a)理想带阻滤波器;(b)巴特沃思带阻滤波;(c)高斯带阻滤波器6.4.1带阻滤波器示例(a)(b)(c)(d)(a)被正弦噪声污染的图像；(b)图(a)的频谱；(c)巴特沃思带阻滤波器；(d)滤波效果图

I=imread('woman1.bmp');%读取图像I=rgb2gray(I);%转换成灰度图像[M,N]=size(I);%得到图像的高度和宽度P=I;fori=1:Mforj=1:NP(i,j)=P(i,j)+20*sin(20*i)+20*sin(20*j);%添加周期噪声endendfigure,imshow(I);%显示原图像figure,imshow(P);%显示加噪图像IF=fftshift(fft2(I));%对原图像作傅里叶变换，并将原点移至中心IFV=log(1+abs(IF));%原图像的频谱PF=fftshift(fft2(P));%对加噪图像作傅里叶变换，并将原点移至中心PFV=log(1+abs(PF));%加噪图像的频谱figure,imshow(IFV,[]);%显示原图像的频谱figure,imshow(PFV,[]);%显示加噪图像的频谱freq=50;%设置带阻滤波器中心频率width=5;%设置带阻滤波器频带宽度ff=ones(M,N);fori=1:Mforj=1:Nff(i,j)=1-exp(-0.5*((((i-M/2)^2+(j-N/2)^2)-freq^2)/(sqrt(i.^2+j.^2)*width))^2);%高斯带阻滤波器endendfigure,imshow(ff,[]);%显示高斯带阻滤波器out=PF.*ff;%矩阵点乘实现频域滤波out=ifftshift(out);%原点移回左上角out=ifft2(out);%傅里叶逆变换out=abs(out);%取绝对值out=out/max(out(:));%归一化figure,imshow(out,[]);%显示滤波结果6.4.2带通滤波器

带通滤波器执行与带阻滤波器相反的操作可用全通滤波器减去带阻滤波器来实现带通滤波器6.4.3陷波滤波器

陷波滤波器被用于阻止(或通过)事先定义的中心频率领域内的频率由于傅立叶变换时对称的,因此陷波滤波器必须以关于原点对称的形式出现。6.3.3陷波滤波器(a)(b)(c)(a)理想陷波滤波器;(b)巴特沃思陷波滤波器;(c)高斯陷波滤波器6.5估计退化函数

退化模型:(1)如果:则系统H是一个线性系统.则系统H称为位置不变系统(或空间不变系统).(2)如果退化模型为线性和位置不变的,其可表示为:即:6.5估计退化函数许多退化类型可以近似表示为线性的位置不变过程.非线性的与位置有关的技术难以求解.由于退化模型为卷积的结果,且图像复原需要滤波器,应此术语”图像去卷积”常用于表示线性图像复原,而用于复原处理的滤波器称为”去卷积滤波器”.6.5估计退化函数退化函数通常未知,因此在复原之前需要估计退化函数.估计退化函数的方法:(1)观察法(2)实验法(3)数学建模法6.5估计退化函数(1)观察法收集图像自身的信息来估计退化函数.例如:对于模糊图像,选择一小部分图像,强信号区,减少噪声影响.并构建一个不退化的图像6.5估计退化函数(2)试验估计法使用与获取退化图像的设备相似的装置,得到准确的退化估计.小亮点成像系统H由于冲激的傅立叶变换为常数A,可得:实验估计模型如下:6.5估计退化函数冲激特性的退化估计一个亮脉冲图像化的(退化的)冲激6.5估计退化函数(3)模型估计法建立退化模型,模型要把引起退化的环境因素考虑在内.例如退化模型就是基于大气湍流的物理特性而提出来的,其中k为常数,与湍流特性相关.6.5估计退化函数大气湍流模型的解释可忽略的湍流剧烈湍流,k=0.0025中等湍流,k=0.001轻微湍流,k=0.00025另外也可以从基本原理开始推导出退化模型.如匀速直线运动造成的模糊就可以运用数学推导出其退化函数.2.运动模糊模型

当成像传感器与被摄景物之间足够快的相对运动时，所摄取的图像就会出现“运动模糊”，运动模糊是场景能量在传感器拍摄瞬间内在像平面上的非正常积累。例6.7I=imread('i_camera.bmp');%读取图像I=rgb2gray(I);%转换为灰度图像figure,imshow(I);%显示图像LEN=25;%设置线性运动位移THETA=11;%设置旋转角度PSF=fspecial('motion',LEN,THETA);%图像线性运动Blurred=imfilter(I,PSF,'circular','conv');%图像被线性运动模糊figure,imshow(Blurred);%显示运动模糊后的图像6.6逆滤波

逆滤波原理：图像退化模型：

傅立叶变换逆滤波恢复逆滤波原理分析：

1.实验证明，当退化图像的噪声较小，即轻度降质时，采用逆滤波恢复的方法可以获得较好的结果。通常，在离频率平面原点较远的地方数值较小或为零，因此图像恢复在原点周围的有限区域内进行，即将退化图像的傅立叶谱限制在没出现零点而且数值又不是太小的有限范围内。

2.当噪声作用范围很大时，逆滤波不能从噪声中提取图像。逆滤波示例(a)原图(b)退化图像(c)逆滤波结果6.7最小均方误差滤波器-维纳滤波

逆滤波恢复方法对噪声极为敏感，要求信噪比较高，通常不满足该条件。为了解决高噪声情况下的图像恢复问题，可采用最小均方滤波器来解决，其中，用得最多的是维纳滤波器6.7最小均方误差滤波器-维纳滤波目标函数:采用拉格朗日乘数法，在有噪声条件下，从退化图像g(x,y)复原出f(x,y)的估计值，该估计值符合一定的准则。维纳滤波器结果上式称为维纳滤波，括号中的项组成的滤波器通常称为最小均方误差滤波器,或最小二乘方误差滤波器。

6.7最小均方误差滤波器-维纳滤波示例(a)运动模糊退化图像(b)7次循环(c)15次循环例6.8I=imread(i_camera.bmp');%读取图像figure(1);imshow(I,[]);%显示图像PSF=fspecial('motion',25,11);%运动模糊函数，运动位移是25像素，角度是11Blurred=imfilter(I,PSF,'conv','circular');%对图像运动模糊处理Noise=0.1*randn(size(I));%正态分布的随机噪声BlurredNoisy=imadd(Blurred,im2uint8(Noise));%对退化后的图像附加噪声figure(2);imshow(BlurredNoisy,[]);%显示运动模糊且加噪声后图像WI1=deconvwnr(BlurredNoisy,PSF);%不带参数的维纳滤波（逆滤波）复原figure(3);imshow(WI1,[]);%显示逆滤波复原结果NSR=sum(Noise(:).^2)/sum(im2double(I(:)).^2);%计算噪信比WI2=deconvwnr(BlurredNoisy,PSF,NSR);%带噪信比参数的维纳滤波复原figure(4);imshow(WI2,[]);%显示带噪信比参数维纳滤波复原结果NP=abs(fftn(Noise)).^2;NCORR=real(ifftn(NP));%计算噪声的自相关函数IP=abs(fftn(im2double(I))).^2;ICORR=real(ifftn(IP));%计算信号的自相关函数WI3=deconvwnr(BlurredNoisy,PSF,NCORR,ICORR);%带自相关函数的维纳滤波复原figure(5);imshow(WI3,[]);6.8几何失真校正在图像的获取或显示过程中往往会产生几何失真，例如成像系统有一定的几何非线性。这主要是由于视像管摄像机及阴极射线管显示器的扫描偏转系统有一定的非线性，因此会造成如图6.20所示的枕形失真或桶形失真。

6.8.1空间变换（SpatialTransformation）

假设一幅图像为，经过几何失真变成了，这里的表示失真图像的坐标，它已不是原图像的坐标了。上述变化可表示为：这里，和是空间变换，产生了几何失真图像。

1．已知和条件下的几何校正（1）直接法。先由

推出，然后依次计算每个像素的校正坐标值，

（2）间接法。设恢复的图像像素在基准坐标系统为等距网格的交叉点，从网格交叉点的坐标(x,y)出发算出在已知畸变图像上的坐标(u,v)，即：2．和未知条件下的几何校正假定基准图像像素的空间坐标(x,y)和被校正图像对应像素的空间坐标(u,v)之间的关系用二元多项式来表示。式中，N为多项式的次数，aij和bij为各项待定系数。对于线性失真：对于一般的（非线性）二次失真：利用“连接点”建立失真图像与校正图像之间其他像素空间位置的对应关系，而这些“连接点”在输入（失真）图像和输出（校正）图像中的位置是精确已知的。

例6.9f=imread('fig620.jpg');%读取图像figure(1),imshow(f);%显示原图像k=0.7;%变换拉伸系数theta=pi/6;%变换旋转角度T=[k*cos(theta)k*sin(theta)0-k*sin(theta)k*cos(theta)0001];%变换矩阵tform1=maketform('affine',T);%几何变换结构g1=imtransform(f,tform1,'nearest');%以最近邻插值进行仿射变换figure(2),imshow(g1);%显示变换结果g2=imtransform(f,tform1,'bilinear');%以双线性插值进行仿射变换figure(3),imshow(g2);%显示变换结果例6.10f=imrea