利用拉氏算子鉴别散焦模糊图像点扩散函数

　CN4321258/TP ISSN10072130X 计算机工程与科学

COMPU TER EN GIN EERIN G&SCIENCE

2005年第27卷第9期　

　Vol127,No19,2005　

文章编号:10072130X(2005)0920040204

利用拉氏算子鉴别散焦模糊图像点扩散函数3

Identificatio n of t he Point Sp read

Functio n f ro m t he Defocus Image U sing

t he Laplacian Operators

陈前荣,陆启生,成礼智,刘泽金,舒柏宏,黎　全,王红霞

CHEN Q ian2rong,L U Q i2sheng,CHENGLi2zhi,L IU Z e2jin,Shu Bo2hong,Li Q u an,W ang H ong2xia

(国防科技大学理学院,湖南长沙410073)

(School of Science,N ational U niversity of Defense T echnology,Changsha410073,China)

摘　要:对于散焦模糊图像,通常把其散焦点扩散函数近似为圆柱形函数。若能有效地从模糊图像中鉴别出该圆柱形的半径,则等价于求出了散焦模糊点扩散函数。利用拉氏算子对散焦模糊图像进行无方向性的二阶微分,然后求微分图像的自相关,则三维显示的自相关图像上会出现一个环形槽。环形槽由数值为负的一系列相关峰连成,槽底位置连成鉴别圆。该鉴别圆以零频尖峰为圆心,直径等于作为散焦点扩散函数的圆柱形直径的两倍。通过插值求和,得到一条鉴别曲线,可以直接鉴别出其直径。这一新的鉴别方法计算量小,鉴别精度高,抗噪声能力较强。

Abstract:The point spread function(PSF)of the defocus image is approximately columniform.It is reasonable to equate identifying the PSF with identifying the radius of the column.A new method of identifying the radius is given.There are three important computing steps.The first step is to derivate the defocus image with a Laplacian operator,which is an iso2 tropic second difference arithmetic operator.The second step is to calculate the autocorrelation of the derivative image. There must be an annular groove in the3D show of the autocorrelation image.The bottom of the annular groove is linked by the autocorrelation peaks whose values are minus.The positions of the peaks link a cycle,whose radius is double that of the PSF column’s.The third step is to transform the X2Y coordinates into polar coordinates for the autocorrelation by the inter2 polation.An identifyication curve will be gained by integral calculus of the angle.The radius of the PSF column can be iden2 tified directly f rom the curve with high precision.This new method has the power to resist noise and does not cost much computing time.

关键词:图像恢复;散焦模糊图像;点扩散函数;鉴别;拉普拉斯算子

K ey w ords:image recovery;defocus image;point spread f unction;identification;Laplacian operator

中图分类号:TP391.4;TB133文献标识码:A

1　引言

散焦模糊图像恢复中,通常把散焦点扩散函数近似为圆柱形函数。若能有效地从模糊图像中鉴别出该圆柱形的半径,则等价于求出了散焦模糊点扩散函数,从而可以利用维纳滤波等复原算法很容易地恢复出清晰图像。

散焦模糊图像的频域图受一阶第一类Bessel函数调制,函数第一过零点的轨迹形成一个圆,利用这一特点可以鉴别散焦模糊点扩散函数,但该方法只实用于信噪比较高的情况。邹谋炎提出了“误差2参数分析法”[1],解决了抗噪声能力较弱问题。但是,从其所给误差鉴别曲线来看,真值附近那段曲线较平坦,不利于准确鉴别。

利用拉氏算子对散焦模糊图像进行无方向性的二阶微分,然后求微分图像的自相关,则三维显示的自相关图像上会出现一个环形槽。环形槽由一系列相关峰连成(连成环形槽的相关峰指向下,数值为负),槽底位置连成鉴别圆。该鉴别圆以零频尖峰为圆心,直径等于作为散焦点扩散函

3收稿日期:2004202209;修订日期:2004203208

基金项目:国防预研基金项目

作者简介:陈前荣(1972),男,湖北黄石人,博士,研究方向为图像恢复。

通讯地址:410073湖南省长沙市国防科技大学理学院203教研室;Tel:137********;E2mail:chenqianrong_cn@https://www.doczj.com/doc/cf14548732.html, Address:Section203,School of Science,National University of Defense Technology,Changsha,Hunan410073,P.R.China

数的圆柱形直径的两倍。通过插值,将微分图像的自相关

转化到极坐标系(ρ,θ

)(以零频尖峰为0点)表示,对角度θ求和,可以得到一条鉴别曲线,通过鉴别曲线即可鉴别出其直径。这一新的鉴别方法算法简单,计算量小,鉴别精度高,抗噪声能力较强。

2　拉氏算子

拉氏算子是一个无方向性的二阶微分算子,在极坐标

系中表示为: 2=929ρ2+1ρ99ρ+1ρ292

9ρ2,在X 2Y 坐标系中表

示为: 2=

92

9x 2+92

9y 2

。其八邻域算子为:1111-811

1

1

。

3　散焦模糊点扩散函数鉴别计算流程

设M ×M 模糊图像为g (i ,j ),从模糊图像出发,利用

拉氏算子鉴别散焦模糊点扩散函数的圆柱形半径,计算流程如下:

(1)利用拉氏算子对散焦模糊图像g (i ,j )进行微分。

(2)利用快速二维傅立叶变换的方法计算出微分图像的自相

关S 。

(3)通过插值,将自相关S 的分布由直角坐标系转化到极坐标系(ρ,θ)(以零频尖峰为0点),然后对角度θ求和,得到一条鉴别曲线S add 。

(4)画图和鉴别。通过鉴别曲线S add 即可鉴别出散焦模糊点扩散函数的圆柱形直径2r 。当信噪比很高时,还可以三维显示出微分图像的自相关S ,从中可以非常直观地看出图中含有一个由一系列负的尖峰所连成的环形槽。它以零频尖峰(自相关图像中心处最高最强的正尖峰)为圆心,其直径是散焦模糊点扩散函数的圆柱形直径的两倍,即4r ,鉴别出4r 即知散焦模糊点扩散函数。

4　散焦模糊点扩散函数鉴别结果

图1是一幅散焦模糊图像,它对应的散焦模糊点扩散函数的圆柱形半径真实值为r 0=7pixels 。将该散焦模糊图像与拉氏算子进行二维卷积运算,得到其二阶差分图像,利用二维快速傅立叶变换方法求得二阶差分图像的自相关函数,并将其归一化,得到归一化自相关函数S 。图2是S 的三维显示(只显示靠近中心的部分

),可以看出有一个环形槽,其直径是散焦模糊点扩散函数的圆柱形直径的两倍,即4r 。

图1　散焦模糊图像 图2　微分图像归一化自相关函数S

图3是过零频尖峰的S 的切面图,从中可以看出,图中

两个负的尖峰比其它值都要小许多,很容易鉴别出来。鉴别出两负尖峰的位置依次为x 1=169pixels ,x 2=197pixels ,故鉴别出来的散焦模糊点扩散函数的圆柱形半径为:

r =(x 2-x 1)/4=7(pixels )(1)

利用鉴别结果,即可得到散焦点扩散函数。利用Wiener 滤波对图1进行恢复,得到清晰的恢复图像如图

4所示。恢复计算中应用循环边界法处理边界[1],以防止G ibbs 效应。

图3　微分图像归一化自相关函数S 图4　恢复图像

虽然对模糊图像进行微分会一定程度上放大噪声,但由于计算自相关及求和过程会大大抑制噪声,这一新的鉴别方法具有较强的抗噪声能力。

将图2用加性高斯白噪声污染,其峰值信噪比PS N R =33dB 。此时,由于噪声干扰作用,零频尖峰及其近邻几点值非常大,S 的三维显示图中环形槽不明显,要画出鉴别曲线S add 加以鉴别(ρ<3时S add 值很大,为了体现负尖峰处的波动,一般不画ρ<3的点)。如图5所示,鉴别得2r =13.6pixels ,即r =6.8pixels ,误差为-0.2pixels 。当模糊图像的峰值信噪比PS N

R <30dB 时(噪声污染严重的情况),该鉴别方法随噪声增大而逐渐失效。

图5　加噪散焦模糊图像鉴别曲线S add

(PS N R =33dB )

5　理论分析

图像模糊的物理模型通常写为:

g (x ,y )=f (x ,y )3h (x ,y )+n (x ,y )

(2)其中,g (x ,y )为模糊图像,f (x ,y )为原清晰图像,h (x ,y )为散焦模糊点扩散函数,n (x ,y )为加性噪声,3表示二维卷积。先不考虑噪声,则有:

g (x ,y )=f (x ,y )3h (x ,y )

(3)对g (x ,y )进行无方向性的二阶微分,有:

2g (x ,y )= 2

∫∫f (α,β

)h (α-x ,β-y )d αd β=∫∫

f (α,β) 2

h (α-x ,β-y )d αd

β=f (x ,y )3 2h (x ,y )

(4)

故有:

S = 2g (x,y ) 2g (x,y )= ( 2g (-x,y ))33 2g (x,y )=

(f (-x,-y )3 2h (-x,-y ))33(f (x,y )3 2h (x,y ))= (f 3(-x,-y )3f (x,y ))3(( 2h (-x,-y ))33 2h (x,y ))= (f (x,y ) f (x,y ))3( 2h (x,y ) 2h (x,y ))= S f 3S 2h

(5)

S f =f (x,y ) f (x,y ),S 2h = 2h (x,y ) 2h (x,y )

这里, 表示二维相关,S 是微分模糊图像的自相关,S f

是原始清晰图像f (x ,y )的自相关,S 2h 是微分点扩散函数的自相关。公式(5)表明,S 等于S f 与S 2h 的二维卷积。

图6是r =7时的散焦点扩散函数,对它用拉氏算子进行微分,然后求微分结果的自相关,得到S 2h 。如图7所示,可以观察到S 2h 存在一个由负的尖峰连成的环形槽,它以零频尖峰为圆心,半径为2r 。图8是S 2h 的过零频尖峰的剖面图。可以看出,负尖峰处的值比所有其它位置的值小许多,这保证了负尖峰的位置易于通过数值鉴别出来。计算过程中还表明,S 2

h 的二维积分等于0。

由于原始清晰图像f (x ,y )通常可以看作为无方向性的一阶马尔科夫过程,它的自相关三维显示通常表现为:中心有一个尖锐的圆锥状相关峰,远离中心部分起伏较小,如图

9,它是清晰图像图4的自相关S f 。

S 等于S f 与S 2h 的二维卷积。在求解过程中,当S 2h

远离S f 的零频尖峰时,由于S f 的起伏较小,且S 2h 的二维积分等于0,故此时对应的卷积结果在0值附近,波动也不会很大;当S 2h 靠近S f 的零频尖峰时,可以将S f 的零频尖峰及其邻域S ′f 看作点扩散函数,微分点扩散函数的自相关

S 2h 看作图像,则此时求S 即可看作S ′f 对S 2h 的模糊。

由图8知,S ′f 可近似认为尖锐圆锥状,故S ′f 对S 2h 的模糊不会改变S 2h 的基本轮廓,即三维显示时,在零频尖峰邻域,S 与S 2h 的基本轮廓相同。故在S 中必然也像S 2h 一样存在一个由负的尖峰连成的环形槽,它同样以零频尖峰为圆心,半径为2r,且半径值易于通过数值鉴别出来,正如图2、图3所示的那样。

6　对S 2

h 的进一步分析

这一节将解释为什么微分散焦点扩散函数的自相关

S 2h 中必定含有一个半径为2r 、由负尖峰连成的环形槽。

在X 2Y 坐标系中,散焦点扩散函数可写为:

h (x ,y )=

1,x 2+y 2≤r 20,其它

(6)

而在极坐标系(ρ∈(-∞,+∞),θ∈(-π2,π

2

))中,散焦点扩散函数可写为:

h (

ρ,θ)=1,-r ≤

ρ≤r 0,其它

(7)

上述散焦点扩散函数表达式都未作归一化,但不影响下面相对强度量的分析。

在极坐标系中讨论微分:

2h (

ρ,θ)=929

ρ2+1ρ99ρ+1ρ292

9θ2

h (ρ,θ)=99

ρ(δ(ρ-r )-δ(ρ+r ))+1ρ(δ(ρ-r )-δ(ρ+r ))=

99ρδ(ρ-r )+1ρδ(ρ-r )

-99

ρδ(ρ+r )+1ρδ(ρ+r )

记:

a (

ρ)=99

ρδ(ρ)+1ρδ(ρ)

(8)

显然有:

a (ρ)=0,当ρ|(-ε,+ε),ε?0+时a (ρ)≥0,当ρ∈(-ε,ε1],ε?0+,0<ε1<ε时a (ρ)<0,当ρ∈(ε1,ε),ε?0+,0<ε1<ε时(9)

则:

2h (

ρ,θ)=a (ρ-r )-a (ρ+r )=a (

ρ)3δ(ρ-r )-a (ρ)3δ(ρ+r )(10) 由式(10)知, 2h (ρ,θ

)的取值与θ无关,即 2h (ρ,θ)的取值具有圆对称性。式(10)中,E 1=a (ρ

)3δ(ρ-r )构成一个右半圆环,它只在r -ε<ρ

(ρ)3δ(ρ+r )构成一个左半圆环,它只在-r -ε<ρ<-r +ε

范围内取非0值。两项合起来,刚好组成一个圆环。圆环的

取值内正外负,即r -ε<|ρ|≤r +ε1时, 2

h (ρ,θ)≥0;r +ε1

<|ρ|

)<0。数值计算出来的 2h 如图10所示,由 2h 的三维显示及下面的等高线可以看出,计算

结果与前面的解析描述基本一致。当然,由于计算中所用的拉氏算子只是对 2的数值近似,故图10的各峰值有一定的

起伏,而不是象解析描述那样完全相等(圆对称)。

呈图9数值分布样式的 2h 的自相关S 2h ,为何必定含有一个半径为2r 、由负尖峰连成的环形呢?求自相关的计算在X 2Y 坐标系,虽然可以转化到极坐标系中求解,但最终很难得出一个包含明确物理意义的表达式来。但是,有一点是可以明确的,即由 2h 于具有圆对称性,故 2h 的自相关

S 2h 必然也有圆对称性。所以,分析 2h 的值分布形式,只

须分析过S 2h 的零频尖峰的X 轴方向切面的值分布即可。分析这一点,可以借助信息光学中用于目标识别的联合变换相关器J TC [2～4]的思想原理。

在J TC 中,输入面有位于(x 1,y 1)的参考像和位于(x 2,

图10 2h的三维显示 图11　借J TC思想分析

2h的相关S 2h

y2)的参考像。当参考象和目标像相同时,输出面(其复光分布等于输入面自相关)将有一对共轭的相关峰,且其关于零频尖峰对称分布,位置分别为(x2-x1,y2-y1)、(x1-x2,y1 -y2)。

对于 2h,将其分成等份的小段微元p1、p2、p3、p4等,其中p1、p3在X轴上关于y轴对称分布,p2、p4关于原点对称分布,如图11所示。假设输入面只有p1、p3,则在表征输入面的自相关的输出面有一对共轭的相关峰,位置分别为(2 2r,0)、(2r,0)。又由式(10)知,p3=-p1,故相关峰为负值。实际上,若只考虑p1、p3,由式(10)有(此时y∈(-dy/2,dy/ 2),故只写一维形式):

2h(x)=a(x)3δ(x-r)-a(x)3δ(x+r)(11)从而:

2h(x) 2h(x)=2a(x) a(x)-a(x) a(x)3(δ(x -2r)-a(x) a(x)3δ(x+2r)(12)在式(12)中,第一项为正的零频尖峰,幅值为自相关峰的两倍,位置为(0,0);第二项U2=-a(x) a(x)3δ(x-2r)表示一个负的自相关峰,位置为(2r,0);第三项U3=-a(x) a(x)3δ(x+2r)表示一个负的自相关峰,位置为(22r,0)。

同理,p2、p4之间也会形成一对共轭的负相关峰,它们关于原点对称,与原点的距离为2r,…。由S 2h的圆对称特性知,这些负的相关峰连成一个环形槽,该圆半径为2r,以零频尖峰为圆心。

上述分析过程中没有考虑到p1和p2、p1和p4等之间的相关峰。实际上:(1)它们间的互相关峰位置不在半径为2r的圆上;(2)由于p1和p2间相对旋转了一个角度,而由J TC的性质知道,该相关器不具有旋转不变性,即p1和p2间的互相关峰的值(绝对值)很小。由于上述两点原因,p1和p2、p1和p4等之间的互相关峰不会对用于鉴别2r的环形的槽形成干扰。

所以,S 2h中必定含有一个半径为2r、由负尖峰连成的环形槽,且其半径易于数值鉴别。

7　实验结果

图12是利用窄视场长焦距光学系统拍摄的约600米远处的散焦模糊图像(像方焦距偏离1.3mm);图13是它的鉴别曲线,鉴别得2r=17.2pixels,其圆柱形散焦点扩散函数如图14所示;利用鉴别结果恢复图像,得到图15。与图12相比,恢复效果是很明显的,

这验证了鉴别方法的有效性。8　结束语

利用拉氏算子求微分,然后计算微分图像的自相关,得到一种新的散焦模糊图像点扩散函数鉴别方法。该鉴别方法算法简单,计算量小,鉴别精度高,抗噪声能力较强。文中借用信息光学中用于目标识别的光学联合变换相关器J TC 概念,给出了该鉴别方法的解析理论分析。

参考文献:

[1]　邹谋炎.反卷积和信号复原[M].北京.国防工业出版社,2001.

[2]　S J utamulia.Joint T ransform Correlators and Their Applications

[A].Proc SPIE.Vol1812[C].1992.2332243.

[3]　Zhang Y imo,Qin Yuwen,Huang Zhanhua,et al.Real2Time High2

E fficiency Joint T ransform C orrelator[J].Acta Optica S inica,

1995,15(7):8832887.

[4]　Zhai Hongchen,Chen Z ikuan,Zhang Tiequn,et al.A Joint T rans2

form Correlator Based on Phase Modulation of the Power Spec2 trum[J].Chinese Journal of Lasers,1998,A25(2):1392144.