第4章-图像复原

1、数字数字图像处理图像处理第四章第四章 图像图像复原复原4.1 图像退化 / 复原过程的模型4.2 噪声模型4.3 只存在噪声的空间滤波复原4.4 频域滤波削减周期噪声4.5 线性、位置不变的退化模型4.6 离散图像退化的数学模型4.7 估计退化函数4.8 逆滤波复原法4.9 最小均方误差复原法 -Wiener滤波复原4.10 几何均值滤波复原4.11 几何失真校正产生原因产生原因几何几何畸变畸变图像质量不佳原因：图像质量不佳原因：光学系统光学系统的像差的像差传感器非传感器非线性畸变线性畸变光学系统光学系统中的衍射中的衍射摄影胶片摄影胶片的非线性的非线性大气流的大气流的扰动效应扰动效应图像运动造

2、图像运动造成的模糊成的模糊图像复原图像复原: Image Restoration 也称图像恢复，是指去除或减轻在图像获取过程中发生的图像质量的下降。它是图像处理的主要内容之一。 图像复原与图像增强图像复原与图像增强相同之处： 改进给定图像的视觉质量。 不同之处： 图像增强是借助人的视觉系统特性，以取得较好的视觉结果（主观过程，不考虑退化原因） 图像复原是利用先验知识、根据相应的退化模型来重建或恢复被退化的图像（客观过程，考虑退化原因） 实际的复原过程是设计一个滤波器，使其能从降质图像 g(x, y)中计算得到真实图像的估值 使其根据预先规定的误差准则，最大程度地接近真实图像 f(x, y) 。

3、，yxf, 图像复原是根据退化原因，建立相应的数学模型，从被污染或畸变的图像信号中提取所需要的信息，沿着使图像降质的逆过程恢复图象本来面貌。 广义上讲，图像复原是一个求逆问题，逆问题经常存在非唯一解，有时甚至无解。为了得到逆问题的有用解，需要有先验知识以及对解的附加约束条件。技术：无约束和有约束 策略：自动和交互处理所在域：频域和空域 从广义的角度上来看： 几何失真退化 、校正（恢复 ） 投影（退化 ）重建（恢复 ）图像复原方法分类图像复原方法分类退化函数H4.1 图像退化图像退化 / 复原过程的模型复原过程的模型复原滤波f (x,y)g(x,y)n(x,y)噪声退化复原输入图像f(x,y)被

4、线性操作 h(x,y) (退化函数的空间描述) 所模糊，并叠加上噪声 n(x,y)，构成了退化后的图像 g(x,y)。),(yxf退化后的图像与复原滤波器卷积得到复原的近似图像),(yxfh 和 n 的信息知道得越多，复原图像 就越接近原始图像 f (x,y)!),(yxf降质通道 f (x,y)(原始图像）（实际图像）图像复原),(yxf(估计图像）技术关键： （1）降质通道建模 (退化函数的选择)（2）复原算法图像复原的任务描述： g(x,y) g(x,y)h 和 n 的信息知道得越多，复原图像 就越接近原始图像 f (x,y)!),(yxf图像复原中的“病态性”：(1) 进行逆运算时，复

5、原近似图像问题不一定有解 - 奇异性问题。(2) 逆问题可能存在多个解。假设 n(x, y) = 0，则 g(x, y)= H f(x, y) 退化模型：g(x, y) = H f(x, y) + n(x, y)，退化系统H是综合所有退化因素的函数。 ),( ),( ),(),( 22112211yxfHkyxfHkyxfkyxfkH(1) 线性 退化系统 H 的性质：其中 f1(x, y)和 f2(x, y) 是任意两幅输入图像，k1 和k2 是比例常数。若则系统 H 是一个线性系统。),( ),( 1111yxfHkyxfkH),(),( byaxgbyaxfH(3) 一致性 (均匀性)

6、( f2(x, y) = 0) (4) 位置 (空间) 不变性 任何与常数相乘的输入的响应，等于该输入响应乘以相同的常数。对任意的 f (x,y)，a 和 b，如果 则一个具有输入输出关系 g(x, y) = H f(x, y) 的系统称为位置不变系统或空间不变系统。),( ),( ),(),( 2121yxfHyxfHyxfyxfH(2) 相加性 (k1 = k2 = 1)若系统 H 是一个线性系统，则两个输入之和的响应等于两个响应之和。图像中任意一点的响应只取决于该点的输入值，与该点位置无关。 不考虑加性噪声：g(x, y)= f(x , y)* h(x, y) 使用线性、位置不变系统的原

7、因： 尽管实际非线性和位移可变的情况能更加准确而普遍地反映图像复原问题的本质，但求解困难。 当退化不太严重时，一般有较好的复原结果。 很多退化都可以用线性位置不变模型来近似，可以借助数学工具求解图像复原问题。 线性、位置不变系统 H 的退化模型： 考虑加性噪声：g(x, y)= f(x, y)* h(x, y) + n(x, y) 空域卷积等价于频域内乘积：G(u, v) = F(u, v) H(u, v) + N(u, v) （其中 f 和 F, g 和 G, H 和 h 分别为傅里叶对）4.2 噪声模型噪声模型 数字图像的噪声主要来源于图像的获取(数字化过程) 和传输过程。例如：CCD摄像

8、机，光照程度和传感器温度是生成图像中产生大量噪声的主要因素。在传输过程中主要是由于所用传输信道被干扰而受到噪声污染。 噪声(Noise) 妨碍人的视觉器官或者系统传感器对所接收图像源信息进行理解或分析的各种因素； 它是最常见的图像退化因素之一；图像中不希望有的部分； 噪声是不可预测的随机信号，只能用概率统计的方法去认识； 噪声影响图像处理的输入、采集、处理和输出等全过程。4.2.1 噪声的特征噪声的特征噪声的空间特性非常复杂，我们只讨论加性噪声独立于空间坐标的情况，且它与图像本身无关联，即噪声分量与像素之间无相关。噪声的频率特性是指噪声在傅里叶域的频率内容。 当噪声的傅里叶谱是常量时，噪声通常

9、称为白噪声。噪声可用 n (x, y) 来表示，它是随机性的，只能用统计量（均值、方差等量）来描述。加性噪声： g (x, y) = f (x, y) + n (x, y) 乘性噪声： g (x, y) = f (x, y) 1+ n (x, y) = f (x, y) + f (x, y) n (x, y)形成波形是噪声和信号的叠加，特点是噪声与信号无关。形成波形是噪声和信号的叠加，特点是噪声受信号的调制。4.2.2 一些重要噪声的概率密度函数一些重要噪声的概率密度函数 (PDF)(1) 高斯噪声（正态噪声）：高斯噪声（正态噪声）：z26070.21)(zp222)(e21)(zzp噪声灰度

10、随机变量用概率密度函数（噪声灰度随机变量用概率密度函数（PDF）来表示。）来表示。高斯噪声的概率密度函数其值有70%落在范围(-),(+)之内，且有95%落在范围落在(-2),(+2)内。 z 表示灰度值， 表示 z 的平均值或期望值， 表示标准差， 2为方差。数学上易于处理，实践中经常使用。z0a2bab26070.)(zp(2) 瑞利噪声：瑞利噪声：azazazbzpbaz0e2)(2)()(概率密度的均值：概率密度的方差：4/ba442/ )( b距原点的位移和其密度图形的基本形状向右变形。瑞利分布密度对于近似偏移的直方图十分适用。Reyleighza(3) 伽马伽马(爱尔兰爱尔兰)噪声

11、：噪声：000e1)(1zzbzazpazbb)!(概率密度的均值： = b / a概率密度的方差：22ab/其中a 0, b为正整数ab/ )(1K)(zpz)()!()(1111bbebbaKGamma严格地说，只有当分母为伽马函数时才是正确的。当分母如此表达式所示时，该密度称为爱尔兰密度。噪声举例噪声举例 (理想情况理想情况)原始图像直方图高斯噪声瑞利噪声伽马噪声附加样本噪声图像及其直方图噪声图像的直方图和它们的概率密度函数曲线对应相似。(4) 指数分布噪声：指数分布噪声：000e)(zzazpaz概率密度的均值： = 1 / a其中a 0概率密度的方差：221 a/它是当b=1时的伽马

12、(爱尔兰)概率密度分布的特殊情况。000e1)(1zzbzazpazbb)!(a)(zpzExponential(5) 均匀分布噪声：均匀分布噪声：otherwisebzaifabzp01)(概率密度的均值： = (a+b) / 2概率密度的方差：ab1)(zpzUniform1222/)(abab(6) 脉冲噪声脉冲噪声(椒盐噪声椒盐噪声)：othersbzforpazforpzpba0)(ap)(zpzImpulseSalt & Pepperabbp如果 b a , 则灰度值 b 在图像中将显示为一个亮点，a 则显示为一个暗点。若 Pa 或 Pb 为零，则称为单极脉冲。若 Pa 或

13、 Pb 均不为零，且它们近似相等时，脉冲噪声值类似于随机分布在图像上的胡椒 (黑点) 和盐 (白点) 颗粒，所以双极脉冲噪声也称为椒盐噪声。附加指数、均匀和椒盐噪声后的图像及其直方图指数噪声均匀噪声椒盐噪声4.2.3 周期噪声周期噪声周期噪声是在图像获取时从电气、电力干扰中产生的。周期噪声的模型是二维正弦波：周期噪声通过频域滤波可以显著地减少。/ )(/ )(sin),(NByvMBxAyxryx0022 对应的MN DFT为：),(),(),(/00N2j002j002jvvuuevvuueAvuRyxBvMB 这是一对分别位于 ( u + u0, v + v0 ) 和 (u - u0, v

14、 - v0 ) 的复共轭冲击。4.2.4 噪声参数的估计噪声参数的估计周期噪声趋向于产生频率尖峰，其参数可以通过检测图像的傅里叶谱来进行估计。周期噪声污染带阻滤波器?噪声PDF参数一般可以从传感器的技术说明中得知，但对于特殊的成像装置，常常有必要估计这些参数。当仅有通过传感器产生的图像可以利用时，常常可以从合理的恒定恒定灰度值灰度值的一小部分估计PDF参数 (小条带，见下图)。?噪声种类高斯噪声瑞利噪声均匀噪声直方图的形状可以指出最接近的PDF匹配。确定PDF对应的噪声种类后，就可以计算灰度值的均值和方差。Sziiizpz)(Sziiizpz)(22）（对所取的小条带S（子图像）方差：均值：(

15、其中zi 值是S中像素的灰度， p(zi ) 是相应的归一化直方图值 )均值和方差求得后，就可以得到PDF中的参数 a 和 b。思考题：怎样获得脉冲噪声的PDF参数？4.2.5 噪声添加噪声添加g = imnoise( f , type, parameters)使用方法：其中 f 是输入图像，g 代表输出图像，type：噪声类型，parameters：噪声参数。例如：g = imnoise( f , gaussian, 0, 0.02);g = imnoise( f , salt & pepper, 0.02);%高斯噪声的均值为0，方差为0.02 在Matlab中，使用imnoise

16、语句来污染一副图像。% 椒盐噪声的噪声密度为 2%g = imnoise( f, speckle, 0.04); %均匀分布随机噪声(乘性噪声) Adds multiplicative noise to the image I, using the equation J = I+n*I, where n is uniformly distributed random noise with mean 0 and variance v. The default for v is 0.04.% MatLab code for 添加噪声: Noise.mI=imread(lena.bmp);figur

17、e(1),imshow(I);title(Original picture); %J=imnoise(I,type,parameters); 添加噪声% type: 噪声种类, parameters: 与噪声种类相关的参数Gauss1=imnoise(I,gaussian,0,0.01); % 高斯噪声Gauss2=imnoise(I,gaussian,0,0.10); Impulse1=imnoise(I, salt & pepper,0.01); % 椒盐噪声Impulse2=imnoise(I,salt & pepper,0.10); Speckle1=imnoise(I

18、, speckle, 0.01); % 均匀分布随机噪声Speckle2=imnoise(I, speckle, 0.10);figure(2); imshow(Gauss1);title(Gaussian,Sigma=0.01);figure(3); imshow(Gauss2);title(Gaussian,Sigma=0.10);figure(4),imshow(Impulse1);title(Salt & Pepper, Sigma=0.01);figure(5),imshow(Impulse2);title(Salt & Pepper, Sigma=0.01);fig

19、ure(6),imshow(Speckle1);title(Speckle, Sigma=0.01);figure(7),imshow(Speckle2);title(Speckle, Sigma=0.01);高斯噪声均值为0，方差为0.01和0.10椒盐噪声密度为 1%和10%Speckle噪声均值为0，方差为0.01和0.10例题：假设 w 是在区间 (0,1) 内的随机变量，试用其来产生具有下列分布的随机数 z ：azazzpbaz0e1)(2)(解：为了得到 z , 需解方程wbaz2e1)()ln(wbaz1这样就可以产生参数是 a 和 b 的随机分布变量 z 。在MatLab中可用

20、简单语句来得到一个MN的随机数组R： R = a + sqrt ( b*log (1-rand(M, N);注意：在Matlab中，rand产生的是0-1之间的均匀分布随机数 randn产生均值为 0 方差为 1 的正态分布随机数 rand(n)或randn(n)生成nn的随机数矩阵 rand(m,n)或randn(m,n)生成mn的随机数矩阵 可用Matlab中的Random函数生成均匀随机数 - 并产生噪声。4.3 只存在噪声的空间滤波复原只存在噪声的空间滤波复原g(x, y)= f (x, y)* h(x, y) + n(x, y) g(x, y) = f (x, y) + n(x, y

21、) G(u, v) = F(u, v) H(u, v) + N(u, v) G(u, v) = F(u, v) + N(u, v) 当一副图像中唯一存在的退化是噪声时，噪声项是未知的，一般无法从g(x, y) 或者G(u, v) 中直接减去它们。当仅有加性噪声存在时，可以选择空间滤波方法去掉噪声。这与图像增强中的滤波方法一样，图像的增强和复原几乎是不可区分的。自适应局部噪声消除滤波器、自适应中值滤波器图像复原的空间滤波器 均值滤波器算术均值滤波器、几何均值滤波器、谐波均值滤波器、逆谐波均值滤波器 顺序统计滤波器中值滤波器、最大值滤波器、最小值滤波器、中点滤波器、修正后的Alpha均值滤波器 自

22、适应滤波器 算术均值滤波器算术均值滤波器4.3.1 均值滤波器均值滤波器xyStstsgmnyxf),(),(1),(Sxy表示中心在(x, y)、尺寸为 mn 的矩形子图像窗口的坐标组。效果：平滑了一幅图像的局部变化，使图像模糊，同时减少了噪声。白条宽7像素，高210像素，间隔17像素99算术均值滤波77算术均值滤波33算术均值滤波 几何均值滤波器几何均值滤波器mnStsxytsgyxf1,),(每一个被复原像素由子图像窗口中像素点乘积的 1/mn 次幂给出。几何均值滤波器所达到的平滑度可以与算术均值滤波器相比，但在滤波过程中会丢失更少的图像细节相对锐化。电路板的X射线图像高斯噪声污染图像(

23、=0, 2=400)33几何均值滤波后的图像33算术均值滤波后的图像图像更清晰图像变模糊 谐波均值滤波器谐波均值滤波器xyStstsgmnyxf),(),(1),(效果：适合处理高斯和均匀等随机噪声； 谐波均值滤波器对于“盐”噪声效果好，但不适用于“胡椒”噪声。原图：白条7像素宽，210像素高，白条间隔17像素99谐波均值滤波77谐波均值滤波33谐波均值滤波?白条带的宽度和高度白条带的宽度和高度xyxyStsQStsQtsgtsgyxf),(),(1),(),(),( 逆谐波均值滤波器逆谐波均值滤波器其中 Q 称为滤波器的阶数。当 Q 是正数时，滤波器用于消除 “椒” 噪声；当 Q 是负数时，

24、滤波器用于消除 “盐” 噪声。但它不能同时消除这两种噪声。当 Q = 0，逆谐波均值滤波器转变为算术均值滤波器；当 Q = - 1，逆谐波均值滤波器转变为谐波均值滤波器。“椒”噪声干扰的图像33 逆谐波均值滤波的结果 (Q=1.5)“盐”噪声干扰的图像33 逆谐波均值滤波的结果 (Q=-1.5)33 逆谐波均值滤波的结果 (Q=-1.5)33 逆谐波均值滤波的结果 (Q=1.5)使用逆谐波均值滤波器时，选择不当的使用逆谐波均值滤波器时，选择不当的 Q 值会带来严重的问题值会带来严重的问题!算术均值和几何均值滤波器适合于处理高斯或均匀等随机噪声；谐波均值滤波器适合于处理脉冲噪声。均值滤波器总结：

25、缺点：必须事先知道噪声是暗噪声还是亮噪声，以便于选择合适 的 Q 符号。4.3.2 统计排序滤波器统计排序滤波器顺序统计滤波器是空间域滤波器，它们的响应基于滤波器包围的图像区域中像素点的排序。滤波器在任何点的响应由排序结果决定。 中值滤波器 用该像素相邻像素的灰度中值来代替该像素的值),(),(),(tsgmedianyxfxySts特点：1）在相同尺寸下，比起均值滤波器引起的模糊少； 2) 对单极或双极脉冲 (椒盐) 噪声非常有效。椒盐噪声干扰图像1遍 33 中值滤波 2 遍 33 中值滤波 3 遍 33 中值滤波概率密度为 Pa = Pb =0.1 最大值滤波器：),(max),(),(t

26、sgyxfxySts特点：1) 对图像中的亮点敏感； 2) 最大值滤波器对于 “椒” 噪声具有良好消除效果。白条宽7像素，高210像素，间隔17像素99最大值滤波图像77最大值滤波图像33最大值滤波图像),(min),(),(tsgyxfxySts1) 对图像中的暗点敏感；2) 最小值滤波器对于 “盐” 噪声具有良好消除效果。 最小值滤波器： “椒” 噪声干扰图像“盐”噪声干扰图像33最小值滤波图像33最大值滤波图像概率密度为概率密度为Pa = Pb =0.1在范围内计算最大值和最小值之间的中点：),(min),(max),(),(),(tsgtsgyxfxyxyStsSts21特点：这种滤波

27、器结合了顺序统计和求平均 (均匀)，对于高斯和均匀随机分布噪声有最好的效果。 中点滤波器白条宽7像素，高210像素，间隔17像素99最大值滤波图像77最大值滤波图像33最大值滤波图像 修正后的Alpha均值滤波器xyStsrtsgdmnyxf),(),(1),(在 Sxy邻域内去掉 g(s,t) 最高灰度值的d/2和最低灰度值的d/2，用 gr(s,t) 代表剩余 mn-d 个像素，由这些剩余后的像素点的平均值形成的滤波器称为修正后的Alpha 均值滤波器：当 d = 0，蜕变为算术均值滤波；当 d = (mn-1)/2, 退变为中值滤波；当 d 取其它值时，适用于包括多种噪声的情况，例如高斯

28、噪声和椒盐噪声混合的情况。加性均匀噪声污染的图像加性椒盐噪声污染的图像由于脉冲噪声的存在，算术均值和几何均值滤波器没有起到好的作用.55几何均值滤波器55算术均值滤波器55中值均值滤波器D=5的修正后的Alpha均值滤波器（55）4.3.3 自适应滤波器自适应滤波器行为变化基于由mn矩形窗口Sxy定义的区域内图像的统计特性。与前述滤波器相比，性能更优；但也增加了算法复杂性。 自适应局部噪声消除滤波器 自适应中值滤波器 自适应局部噪声消除滤波器均值：给出了区域中灰度平均值的度量；方差：给出了区域的平均对比度的度量。 作用于局部区域 Sxy 的滤波器响应基于以下4个量： (a) g(x,y) -

29、噪声图像在点(x,y)上的值 (b) - 干扰 f(x,y)以形成 g(x,y) 的噪声方差 (c) mL - 在 Sxy 上像素点的局部均值 (d) - 在 Sxy 上像素点的局部方差22L1，如果 ，滤波器返回 g(x,y) 的值。 - 因为在 g(x,y) 下零噪声的情况等同于 f(x,y)。 滤波器的预期性能3，如果局部方差 与 是高相关的，滤波器返回一个g(x,y)的近似值。2，如果 ，滤波器返回区域 Sxy 上像素的算术均值。 - 这种情况发生在局部面积与全部图像有相同特性的条 件下，并且局部噪声简单地用求平均来降低。0222L22L唯一需要知道或估计的未知量是噪声方差；其它参数可

30、以从 Sxy 中的像素计算出来。 基于上述假定的自适应表达式L22myxgyxgyxfL),(),(),(自适应滤波器举例高斯噪声污染图像(均值为0，方差为1000)77的算术均值滤波77的几何均值滤波77的自适应滤波更加锐化 自适应中值滤波器zmin = Sxy 中的灰度级最小值zmax = Sxy 中的灰度级最大值zmed = Sxy 中的灰度级中值zxy = 坐标(x, y)处的灰度级Smax = Sxy 允许的最大尺寸主要目的： 除去“椒盐”噪声（冲激噪声） 平滑其它非冲激噪声 减少物体边界细化或粗化等失真 自适应中值滤波器算法工作在两个层次，定义为A层和B层 A层层：A1 = zme

31、d zminA2 = zmed zmax如果 A1 0 且 A2 0 且 B2 zmed zmin 说明 zmed 不是脉冲zmax zxy zmin 说明 zxy 不是脉冲zxy = zmax 或 zxy = zmin椒盐噪声干扰图像 (概率 Pa = Pb =0.25)77中值滤波器自适应中值滤波（Smax = 7）噪声被有效去除,但细节损失严重保持了点的尖锐性,细节清楚4.4 频域滤波削减周期噪声频域滤波削减周期噪声周期噪声的模型是二维正弦波，通过带阻、带通和陷波滤波器可以被有效去除。有关内容在上一章 “图像增强”里面已经详述，在此不再说明。 退化模型：g(x, y) = H f (x,

32、 y) + n(x, y)4.5 线性、位置不变的退化模型线性、位置不变的退化模型一副连续图像 f (x, y) 可以表示为：ddyxfyxf),(),(),( 式中 函数表示空间上点脉冲的冲激函数。),(),(ddyxfH ),(),(yxfHyxgddyxfHyxg),(),(),( 当退化函数 H f (x, y) 是线性算子时，可把 “加性特性” 扩展为积分：假设 n(x, y) = 0，则 g(x, y)= H f(x, y)退化输出图像函数：由于 独立于 x 和 y，使用 “均匀性” ，可得出),(fddyxfHyxg),(),(),( ddyxHf),(),( ),(),(yxH

33、yxh)(),(),(1ddyxhf 其中称为退化系统 H 的冲击响应函数(或点扩展函数PSF)。),(yxh(1) 式表明，如果退化系统 H 的冲击响应函数 h 已知，则对任意输入图像 的响应都可以求得 - 线性系统理论的核心。),(f当退化函数 H f (x, y) 是位置不变的系统时，),(),(yxHyxh)(),(),(),(2ddyxhfyxg 在有加性噪声存在的情况下，线性、位置不变的退化模型可表示为：),(),(),(),(yxnddyxhfyxg ),(),(*),(),(yxnyxfyxhyxg对上式取傅里叶变换，基于卷积定理，在频域内可得：),(),(),(),(vuNv

34、uFvuHvuG式(2)是卷积积分，输出 g 是冲击响应 h 和 输入图像 f 的简单卷积：),(*),(),(yxfyxhyxg许多退化类型可近似表示为线性的位置不变过程。这一方法的优点在于广泛线性系统理论对于解决图像复原问题很实用。由于退化模型为卷积的结果，且图像复原时需要滤波器，所以用于复原处理的滤波器称为 “去卷积滤波器”。4.6 离散图像退化的数学模型离散图像退化的数学模型 一维图像退化过程： 10eee1 1 0MmMxmxhmfxg,)()()(1010)()(eMxAAxxfxf1010)()(eMxBBxxhxh将输入数字图像和冲击响应函数扩展为：输出的降低质量的数字图像为：

35、输入的数字图像 f (x) ： 大小为 A冲激响应函数 h(x)： 大小为 B为避免卷积周期重叠误差，采用添零延拓方法将它们分别扩展成具有 M 个元素的周期函数，且 M A + B 1。 用矩阵形式表示)()()()()()()()()()()()()()()(110 021201110110eeeeeeeeeeeeeeeMfffhMhMhMhhhMhhhMgggHfg)()()()()()()()()(021201110eeeeeeeeehMhMhhhhhMhhH记忆方法： 将0 1 2 3 4 M-2 M-1画成一个圆圈，依次循环作为he(x)根据周期性： he(x) = he(x+M)循

36、环矩阵H= he(x) 1101010eNyBorMxAByandAxyxfyxf),(),(1101010eNyDorMxCDyandCxyxhyxh),(),(1 1 01 1 0 1010eeeNyMxnymxhnmfyxgMmNn,),(),(),( 二维退化过程 不考虑噪声，输出的降低质量的数字图像为：设输入的数字图像f (x, y) 的大小为AB，点扩展函数h(x,y)被均匀采样的大小为CD。采用添零延拓方法将其分别扩展成 M=A+C-1 和 N=B+D-1个元素的周期函数： 用矩阵形式表示：fHg H 是 MNMN 维矩阵， 即由MM个大小为NN的子矩阵组成， 032130122

37、1011210HHHHHHHHHHHHHHHHHMMMMMM子矩阵 Hj ( j=0, 1, 2, , M-1) 为分块循环矩阵，大小为NN：),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(0321301221011210eeeeeeeeeeeeeeeejhNjhNjhNjhjhjhjhjhjhNjhjhjhjhNjhNjhjhjH 考虑噪声1 1 01 1 0 1010eeeNyMxyxnnymxhnmfyxgMmNne,),(),(),(),(解解： H 是 MNMN 维矩阵，M=N=512， MNMN =512 512 512 512 =

38、236 = 68719476736例：对于512 512的图像，用离散退化模型的矩阵表示。这意味着要解出 f (x, y) 需要求解 236 (约687亿) 个联立方程组，上式表明，给定了退化图像 g(x,y)、退化扩展函数 h(x,y) 和噪声分布n(x,y)，就可以得到原始图像 f (x,y)的估计 。 ),(yxf可以证明在频率域中有：1 1 01 1 0NyMxvuNvuFvuHvuG,),(),(),(),(实际计算工作量巨大 (天文数字)，需要存储空间巨大 ! 两类解决办法两类解决办法(1) 通过对角化简化分块循环矩阵，再利用FFT快速算法可以大大降低计算量并能极大地节省存储空间。

39、(2) 分析退化的具体原因，找出退化函数 H 的具体简化形式 - 估计退化方法。4.7 估计退化函数估计退化函数 三种主要的估计退化函数的方法： (1) 图像观察估计法 (2) 试验估计法 (3) 模型估计法使用以某种方式估计的退化函数复原一幅图像的过程有时称为盲目去卷积，因为真正的退化函数很少能完全知晓。 假设空间不变的，由 推导出完全函数 构造一个估计图像，它与观察的子图像有相同大小和特性 表示观察子图像， 表示构造的子图像 和 为对应的傅立叶变换4.7.1 图像观察估计法图像观察估计法 给定一幅退化图像，但没有退化函数 H 的知识，那么估计该函数的方法之一就是收集图像自身的信息： 寻找简

40、单结构的子图像 寻找受噪声影响小的子图像),(),(),(vuFvuGvuHssS),(yxgs),(yxfs),(vuGs),(vuFs),(vuHs),(vuH4.7.2 试验估计法试验估计法 使用与被退化图像设备相似的装置，并得到一个脉冲的冲激响应，可以进行较准确的退化估计：AvuGvuH),(),(一个脉冲点一个脉冲点成像系统成像系统H 此处A是一个冲激的傅立叶变换，表示冲击强度，为一常数。 右图为一个放大的亮脉冲以及退化的冲激。),(yxg退化图像退化图像 4.7.3 模型估计法模型估计法 建立退化模型，考虑引起退化的环境因素22 5/6()( , )k uvH u ve 例如：Hu

41、fnagel 等 Stanley 的退化模型就是基于大气湍流的物理特性而提出来的，其中k为常数，与湍流特性相关。(除了指数5/6，该公式与高斯低通滤波形式相同.)大气湍流模型模拟退化模糊一幅图像：剧烈湍流(k=0.0025)中等湍流(k=0.001)轻微湍流(k=0.00025)可忽略的湍流22 5/6()( , )k uvH u ve 可以从基本原理开始推导出退化模型Tttyytxxfyxg000)d)(),(),(T: 采集 (曝光) 时间x方向运动分量 y 方向运动分量 当成像传感器与被摄景物之间存在足够快的相对运动时，所摄取的图像就会出现“运动模糊”。匀速直线运动造成的模糊，可用数学推

42、导出其退化函数。假定 f(x,y) 表示无运动模糊的清晰图象，相对运动用 x0(t) 和 y0(t) 表示，运动模糊的图像可表示为：令： yxvyuxyxgvuGd)d(j2)exp( )(,TttvytuxvuH000)d()(j2exp(),(所以，若 x0(t) 和 y0(t) 已知，传递函数H(u,v) 可以直接得到。进行傅立叶变换得到：)(exp),(),(0000j2vyuxvuFyyxxfTttvytuxvuFvuG000)d()(j2exp()( ,),( Ttyxvyuxtyytxxf000 d dd )(j2)exp(),(- )()()(vuFvuHvuG,则：平移性原始

43、图像运动引起的图像模糊 (a=b=0.1, T=1)假设当前图像只在 x 方向以给定的速度 x0(t) = at / T 做均匀直线运动。当 t = T 时，图像由总距离a取代。令 y0(t) = 0，则有：TtuxtvuH0)(j2de0),(若允许 y 分量也变化，y0(t) = bt / T 运动，则退化函数变为：)()(sin)(),(vbajevbavbaTvuHajTTuateaaTt)sin(/0j2de李抱朴, http:/ 如果已知系统的传递函数 ，则根据4.8 逆滤波复原法逆滤波复原法 vuFvuHvuG,vuH,vuHvuGvuF,可得复原图像的谱，经傅氏逆变换即可得到复

44、原图像在忽略噪声的影响，退化模型的傅氏变换为实际应用时存在病态的问题，即在 H(u,v) 等于零或非常小的数值点上， 将变成无穷大或非常大的数。 ),(vuF- 这就是逆滤波复原法 vuNvuFvuHvuG,vuHvuNvuFvuHvuNvuHvuGvuF,),(,系统中存在噪声时退化模型的傅立叶变换为：写成逆滤波复原的方式：1）即使知道退化函数，也不能准确复原图像，因为噪声函数 N(u,v) 是一个随机函数，其傅里叶变换未知。2）如果退化是零或非常小的值，噪声即使数值很小，但 N(u,v)/H(u,v) 之比 (上式第二项) 可能非常大，很容易决定的估计值。),(vuFdudvevuHvuN

45、yxfyxfvyuxj )(),(),(),(),(21解决退化是零或非常小的值的途径：限制滤波的频率，使其接近原点值。 在离频率平面离原点较远的地方，H(u,v)数值较小或为零，因此图像复原在原点周围的有限区域内进行，即将退化图像的傅立叶谱限制在没出现零点而且数值又不是太小的有限范围内，即通过将频率限制为接近原点分析，减少了遇到零值的几率。 剧烈湍流(k=0.0025)大气湍流模型模拟退化模糊一幅图像(前面已讲过内容)：可忽略的湍流652222/)/()/(),(NvMukevuH对退化函数 H(u,v) 进行精确取反，并进行逆滤波，结果如下图。全频直接逆滤波复原半径为40时截止H半径为70

46、时截止H半径为85时截止H结果表明：噪声明显影响了图像复原结果，一般直接逆滤波效果较差。剧烈湍流图(k=0.0025)退化图复原图像：目标： 寻找一个滤波器，使得复原后图像 与原始图像 的均方误差最小。4.9 最小均方误差复原法最小均方误差复原法 -Wiener滤波复原滤波复原逆滤波没有清楚说明如何处理噪声！),(yxf)(22ffEe误差度量：现讨论一种滤波复原法 - Wiener滤波复原： 综合考虑退化函数和噪声统计特征。E是宗量的期望值。),(yxf因此维纳滤波复原又称为最小均方误差复原。min),(),(2yxfyxfE),(),(/ ),(),(),(),(),(vuGvuSvuSv

47、uHvuHvuHvuFf221f误差函数的最小值在频域里可以通过近似图像 的傅里叶变换来计算：叶变换为复原近似图像的傅里换为退化图像的傅里叶变换为退化函数的傅里叶变其中：),(),(),(vuFvuGvuH为未退化图像的功率谱为噪声的功率谱22),(),(),(),(vuFvuSvuNvuSf),(),(),(*vuHvuHvuH2维纳滤波器(2) 未退化图像的功率谱难以知道，可用下式近似表示：KvuHvuHvuHvuHw221),(),(),(),(1) 如果噪声为 0，其功率谱消失，维纳滤波就退化为逆滤波。讨论：讨论：式中 K 是根据信噪比的某种先验知识确定的常数。),(/ ),(),(),(),(),(vuSvuSvuHvuHvuHvuHfw221维纳滤波复原：维纳滤波需要假定下述条件成立(或近似成立)： 系统为线性、空间不变；(1) 退化图像、原始图像和噪声都是均匀随机场，噪声的均值为零，且与图像不相关。 维纳滤波复原与逆滤波复原的比较全频逆滤波半径受限逆滤波维纳滤波复原 (交互选择K)逆滤波 维纳滤波运动模糊及加性噪声污染图像均值为0，方差为650的加性高斯噪声污染噪声幅度的方差比第一图小1个数量级噪声幅度的方差比第一图小5个数量级 维纳滤波的缺点： 未退化图像和噪声的功率谱必须是已知的； 功率比(信噪比)常数K 的估计一般还是没有合适的解。4.