基于三维高斯模型的参数盲解卷积算法

基于三维高斯模型的参数盲解卷积算法33

刘　莹,何小海33,陶青川,赵　佳

(四川大学电子信息学院,四川成都610064)

摘要:提出了一种新的基于三维高斯点扩展函数(PSF)模型的参数盲解卷积(PBD)算法,并将此算法用于显微光学切片的图像复原。由于PBD算法需要在估计样本函数的同时估计PSF的参数,一般采用的PSF的模型较为复杂,计算量大,收敛慢;而基于三维高斯PSF模型的算法只需要估计2个参数,因此计算量大大降低了。经过实验验证,该算法能够较好地复原光学切片的图像,并且估计出PSF的参数。关键词:计算光学切片显微术(COMS);三维高斯点扩展函数(PSF);参数盲解卷积(PBD)

中图分类号:TP391.4 文献标识码:A 文章编号:100520086(2006)0420493205

Parametric B lind Deconvolution B ased on32D G aussian Model

L IU Y ing,H E Xiao2hai33,TAO Qing2chuan,ZHAO Jia

(College of Electronic Information,Sichuan University,Sichuan Chengdu610064,China)

Abstract:A parametric blind deconvolution(PBD)algorithm was pre sented based on a new model of the32D Gaussian point spread function(PSF)to restore the micro scopic op2 tical slices.The PBD algorithm needed to simultaneously estimate the specimen function and the parameters of the PSF,while the PSF model wa s complicated,needed a large number of computation and converged slowly.Based on the Gaussian PSF model presen2 ted,this algorithm only needed to estimate two p arameters,so the amount of the computa2 tion wa s greatly reduced.The experiments proved that this algorithm can preferably re2 store the optical slices and estimate the parameters of the PSF.

K ey w ords:computational optical sectioning microscopy;32D Gaussian point spread function(PSF);para2 metric blind deconvolution(PBD)

1　引　言

显微镜是观察生物样本的主要工具[1,2]。在三维显微学里,通过将显微镜聚焦在生物样本的不同平面来获取二维序列图像。每个聚焦层面的图像,不仅包含本层面的信息,还包含其它层面的离焦模糊信息,从而产生降质。通常使用2种方法来复原图像:激光扫描共焦显微术(L SCM)[325];计算机光学切片显微术(COSM)[6]。前者可以直接得到物体的三维图像,但是不能得到高信噪比的图像,并很容易对生物样本产生漂白现象。后者要求分离图像序列中焦平面上和焦平面外的光信息,从而得到样本正确的三维形态,同前者相比,它具有信噪比高,并对样本不产生漂白现象,但运算速度慢。因此,本文提出了基于三维高斯点扩展函数(PSF)模型的参数盲解卷积(PBD)[7,8]算法,以减少运算量。实验表明,该算法能较快较好地复原图像,并估计出PSF的参数。

2　PBD算法的基本原理及步骤

三维光学切片的成像模型可表示为:通过显微镜观察到光学切片的图像g(X i),X i=(x i,y i,z i)为三维像空间的一点;恢复样本函数s(X o),X o=(x o,y o,

光电子・激光

第17卷第4期　2006年4月 J ournal of Optoelect ronics・L aser Vol.17No.4　Apr.2006 3收稿日期:2005208231　修订日期:2006201210

　3　基金项目:国家自然科学基金资助项目(60372079)

　33E2m ail:cias@

z o )为三维样本空间的一点。由成像公式有

g (x i )=h (x o ) s (x o )

(1)其中,h (x o )为显微镜的PSF ,并且由一组已知或未知

的参数决定。PBD 算法的目的是同时估计样本函数和PSF 的未知参数。这里使用最大似然函数(ML )算法来进行运算,而ML 是基于泊松分布概率的,即通过显微镜观察到的图像g (x i )是符合泊松分布[9],图像的均值为u (x i )=h (x o ) s (x o ),由此可得[10] P[g (x i )|s (x o )]=

∏

n

u (x i )g (x i )

e

-u (x i )

g (x i )!

(2)

在PBD 算法中,假设PSF 是由参数Θ={θ1,θ2,

θ3,…

}决定的,故s (x o )关于式(2)的对数ML 为[11]

L [s (x 0),Θ]=

-∫I

∫

O

h (x i -x o ;Θ)s (x o )d x o -g (x i )log

∫

o

h (x i

-x o

;Θ

)s (x o

)d x o

d x i (3)

因此,基于ML 的PBD 算法即是:求出使式(3)最大

的s (x o )和Θ。

为了获得ML 的最大值,使用一组协调的迭代算法进行计算。即,在每次迭代中,当PSF 的参数不变时,更新样本函数;当样本函数不变时,更新PSF 的参数。 首先使用EM 算法[12]更新s (x o )。对式(3)求导,得

5L

5s

=-∫

I

∫O h (-x o )d x o -g (x i )

h (x o ) s (x o )

h (-x o )d x i =0

(4)

由此可推出 ^s (x )

(k +1)

=^s k

(x )

H o ×h (-x ) g (x )^g k

(x )

(5)

这里:^g k (x

)=^s k (x ) ^h k

(x );H o

=

∫

O

h (x o

)d x

o

。

其次,为了更新PSF 的参数,用共轭梯度算法[13]

来优化运算,有

^Θk +1=^Θk +α

τk (6)这里,τk =γk +β

k -1τk -1

,是迭代的方向。其中 γk = ΘL k

(7)

是式(3)关于PSF 参数的导数。有

5L

5Θ

=-∫

I

h Θ s -g ×

h Θ s

h s

d x i

(8)

其中:h Θ=5h/5Θ,是PSF 对参数的导数。βk -1=〈

γk

-γk -1,γk 〉/〈γk -1,γk -1

〉,〈・〉表示内积。使用线性

搜索法决定式(6)中的步长α,即求使L (^Θk +α

τk

)-L (^Θk )最大的α。如果

‖(L k -L k -1)/L k ‖<ε(9)

满足,则此次迭代的参数更新结束,进入下次迭代运算。反之,则共轭梯度算法要不断的重复,直到满足式(9)。其中,ε为任意小的正值。

3　基于简化PSF 模型的PB D 算法

典型的非共焦显微镜的PSF 的形式为[14]

h (X o ;Θ

)=∫

10

J 0[(2π/λ)a ρr o /z d ]P (ρ,z o ,Θ)ρd

ρ2

(10)

其中:P (ρ,z o ,Θ)=A (ρ

)exp [j W (ρ,z o )]为物镜的瞳孔函数,A (ρ

)是瞳孔函数的幅度,ρ为物镜后焦平面的径向坐标,W (ρ,z o )=

2πλ

O PD (ρ,z o )是相位偏差,λ

是波长,O PD (ρ,z o )是高斯参考平面和实际的波阵面的光程差;J o 为阶数为0的第一类贝塞尔函数;a =

z d N A /

M 2

-N A 2

≈z d N A /M ,N A 和M 分别是物

镜的数值孔径和放大率;z d 是物镜到探测器的距离;

r o =

x 2

o +y 2

o 。

在应用中,该PSF 存在不足:首先,它需要大量的参数,但实际的光路比较复杂,因此这些参数不易获得;其次,由于它形式较复杂,因此运算量大,速率低。因此,需要提出另一种形式的PSF 来进行序列图像复原。研究表明,三维PSF 的能量集中在其中心的沙漏部分,因此提出了沙漏形状简单的PSF 模型。考虑到高斯函数和贝塞尔函数的相似性,且它的计算简单,因此使用高斯函数作为系统近似的PSF ,不同层面可以用不同方差的高斯函数来模拟。 三维高斯函数PSF 的模型为

h (x ,y ,z )=1

2πσ2z

e

-

x 2+y 2

2σ2z

(11)

这里

σ2

z =σz 0+α

|z -z 0|(12)

即,方差同此平面和中心位置的距离差成正比。在进

行共轭梯度运算时,式(7)即为γk =[γk σz 0,γk

α],而相应的式(8)为5L/5Θ=5L/5σz 0,5L/5

α,其中h Θ=5h/5σz 0,5h/5α=[h σz 0,h α]。由式(11),PSF 对σz 0

和α分别求导,得

^h σz 0=-1

2πσ4

z

exp [-

x 2+y

2

2σ2z

]+

x 2+y

2

4πσ6

z

exp [-

x 2+y

2

2σ2

z

]

(13)

^h α=-

1

2πσ4

z

exp [-

x 2+y

2

2σ2

z

]|z -z 0|+

x 2+y

2

4πσ6

z

exp

・

494・ 光电子

・激光　

2006年　第17卷