基于EM算法的眼底OCT图像反卷积去模糊技术*

闫　芳　宋　双　连　剑　任衍具　尹义龙　郑元杰

(1.山东师范大学信息科学与工程学院，济南，250014； 2.山东师范大学心理学院，济南，250014； 3.山东大学计算机科学与技术学院，济南，250100)

引　　言

光学相干层析成像技术(Optical coherence tomography, OCT)作为一种发展迅速的视网膜检查手段，具有分辨率高、灵敏性强、成像速度快、检查方便以及对眼睛无伤害等特点[1]。然而在获取OCT视网膜图像时，由于眼球运动或者散焦作用会导致图像产生模糊现象。对于此类模糊图像，退化信息隐藏在图像当中，这要求我们找出其中的退化信息来恢复图像。与自然图像相似，在不考虑噪声存在的情况下，OCT模糊图像可以看作是清晰图像与点扩散函数卷积的结果。

图像去模糊是图像处理方面基本问题之一，根据点扩散函数是否已知，将其分为盲去模糊与非盲去模糊。在非盲去模糊方面，Tiwari等人在假设没有高斯噪声影响的前提下，对比分析了不同的运动模糊参数估计方法[2]。同样在盲去模糊方面，Krishnan等人运用了一种归一化的稀疏度量方法来解决图像模糊问题[3]；Mai和Liu结合多种去模糊方法将它们得到的模糊核融合并提出了一种数据驱动方法，该方法能够有效地从训练集中学习得到融合后的模糊核[4]；Oliveira等人基于模糊图像的频谱和弱假设检验提出了对线性运动和散焦造成模糊参数的估计方法[5]，由于该方法不需要迭代，复原运算速度得到显著提升。对于OCT成像，目前已经出现了很多应用于图像的降噪和轮廓信息的提取等方面的技术[1,6]。而在OCT图像去模糊方面，Kulkarni等人首次提出了用反卷积的方法去除图像模糊[7]；Liu等人同时针对纵向、横向相互独立的模糊核，利用维纳滤波和Lucy-Richardson滤波分别进行OCT图像去模糊，并对两种方法的效果进行了对比[8]；Liu等人采用了基于信息熵的方法对因失焦产生模糊OCT图像中的点扩散函数进行自动估计[9]。尽管如此，OCT图像去模糊的研究特别是结合OCT图像的特征属性进行针对性去模糊的处理依然比较少见。

由于视网膜具有结构特殊、各组织层次联系密切和局部边缘不明显等特点，导致对眼底OCT图像成像的精准度要求更高。此外，因为OCT模糊图像中含有一定数量的异常值干扰，例如非高斯噪声等，均为OCT模糊图像的复原带来了挑战。为了在不破坏图像关键信息的前提下，尽量保持眼底OCT图像轮廓特征，本文基于已知点扩散函数信息，提出了一种新的基于最大期望(Expectation-maximization, EM)方法的反卷积去除OCT图像模糊算法。该算法的思想是：将图中像素分为符合传统模型和使传统模型失效两种类型，在此基础上建立一个鲁棒的非盲去模糊反卷积方法。实验结果表明该方法可以抑制由异常值干扰引起的振铃效应，从而显著提高眼底OCT图像的清晰度。

1　基础理论

1.1OCT技术原理及模糊原因

图1　OCT技术原理图Fig.1　Schematic diagram of OCT

图2　散焦造成的眼底模糊图像Fig.2　　　　OCT image with blur caused by 　　　out-of-focus

在眼科临床诊断中，OCT利用低相干光干涉成像原理，以光学组织切片的形式对视网膜成像。目前OCT系统都是基于低干涉光束，比如Michelson干涉仪。由OCT技术原理图(见图1)可以看出，干涉仪内的光纤耦合器将入射光分为参考光和测定光，前者通过平面反射镜被反射回来，而后者经过眼底各层反射，产生具有时间差的反射光，最终经光纤耦合器合为一束光进入探测器。当参考光和测定光再次重合时会产生干涉进而能够测定物体距离。通过测定反射光的强度和时间偏差，再换算成空间位置关系就可以获得眼底断层图像[10]。

导致眼底OCT图像模糊的原因主要是运动模糊和散焦模糊。运动模糊通常是由于眼球运动造成，由于新型的OCT设备中带有自动追踪眼球运动的模块，因运动造成的模糊OCT图像比较少见。因此，本文主要关注由于散焦造成模糊的眼底OCT图像，如图2所示。从图2可以看出，由于拍摄聚焦在视网膜色素上皮层中心区域，使得聚焦点外本应清晰分层的结构变得模糊，如图方框中区域所示。

1.2　传统图像退化模型

通常图像模糊可以看作原清晰图像I(x,y)与造成模糊的点扩散函数k(x,y)经过卷积运算后，加入外来噪声n(x,y)，退化成模糊图像b(x,y)的过程，表达式为

b(x,y)=I(x,y)*k(x,y)+n(x,y)

(1)

式中“*”为卷积运算。

由此可见，在不考虑噪声的影响下，图像去模糊可以看作是图像模糊的逆过程，即求解出点扩散函数，然后经过反卷积运算，从而预测求得原图像I(x,y)。常用的点扩散函数k(x,y)主要有线性运动退化函数、散焦模糊退化函数以及光学系统衍射、像差等引起的高斯退化函数[11]。线性运动的点扩散函数公式为

(2)

式中：L为运动模糊长度；θ为运动模糊方向。散焦点扩散函数是一个均匀分布的圆盘函数，表示为

(3)

式中r为散焦半径，也就是模糊半径。高斯退化函数表达式为

(4)

式中：k为常数；a为正常数；C为k(x,y)的圆形支持域。

1.3　常用图像去模糊技术分类

图像去模糊技术通常根据模糊的区域差别、信息源的分类、恢复手段不同以及点扩散函数是否已知分为4种类型[12]。非盲去模糊是指已知点扩散函数信息对图像清晰化处理的过程，它是盲去模糊的基础，只要估计出造成图像模糊的点扩散函数，所有的非盲去模糊方法都可以应用在盲去模糊中。相较于非盲去模糊，盲去模糊因缺乏点扩散函数信息而变得复杂：一方面，在盲去模糊的过程中，不同的点扩散函数会产生不同的清晰图像，使得一张模糊图像对应许多清晰图像，因此很难预测合适的点扩散函数；另一方面，图像盲去模糊算法运行时间一般较长，而非盲去算法效率较高。鉴于盲去模糊过程中产生的这些难题，本文采用了一种非盲去模糊的算法处理眼底OCT图像，取得了较好的效果。

2　基于EM算法的OCT图像反卷积技术

OCT图像在成像过程中会出现异常值，一方面，由于OCT眼底图像在成像过程中会引入噪声，而已有的很多方法都是针对特定的噪声去模糊，如高斯噪声，这对于散焦引起的眼底OCT图像去模糊并不适合；另一方面，在OCT成像过程中，由于某些原因OCT图像会产生过饱和现象，而相机具有精准度范围限制，会对OCT图像进行裁剪，在这两种异常值存在的情况下采用反卷积算法去模糊容易产生振铃效应。因为振铃的存在，会对OCT各组织层次造成干扰，而眼底OCT图像在层次结构分明、轮廓清晰等方面有着较高要求，使得在去模糊同时有效抑制振铃效应，防止层次边缘失真成为研究的难题。

本文针对以上情况，采用了一种基于EM算法的OCT图像反卷积技术，不仅能够去模糊，而且可以最大程度上抑制振铃效应。本文技术根据Cho等人[13]的思想引入了一种新的表示图像退化的模型，它能够更加准确地描述包括异常情况在内的图像模糊过程。根据这一模型，可以将模糊图像中的像素分为满足式(1)的和不满足式(1)的两种情况。为此再引入一个二值映射关系m:当mx=1时表示该像素满足式(1)，而当mx=0时表示该像素不满足式(1)，其中x为像素的索引。假设相机传感器捕获到没有噪声的模糊图像，过饱和的像素值经裁剪限定在相机的允许范围之内，模糊图片中出现的噪声等异常值便被加入到被裁剪的模糊图像中。针对此情况，本文算法提出了一种新的模型，此过程可以表示为

b= c (k*I) +n

(5)

式中c为一个限幅函数。对于符合式(1)的像素默认增加的噪声为高斯噪声，而对于不符合式(1)的像素则认为增加的噪声与点扩散函数、原图像均相互独立。为了得到原图像I，在上述分析与假设的基础上，建立如下基于最大后验概率的公式

(6)

然后，根据贝叶斯公式对式(6)求解得到

(7)

式中M表示所有可能出现m的排列。由于P(b,m|k,I)求解困难，而且很难获得m的真实值，因此Cho等人利用EM方法估计P(b,m|k,I)的期望以此发现I的估计值。具体的EM解决过程如下：

(1)E过程

①先对原图像I进行估计，初始估计值为I0。

②定义

Q(I,I0)=E[logP(b,m|k,I)]=E[logP(b|m,k,I)+logP(m|k,I)]

(8)

式中E表示期望。

③根据文献[14]可以得知

(9)

式中：f=k*I；N表示高斯分布；σ表示标准方差；C为常数。此外，假设m空间独立，则

(10)

(11)

因此有

(12)

式中E[mx]=P(mx=1|b,k,I0)。

④根据贝叶斯公式

(13)

因此可以进一步求得

(14)

式中：f0=k*I0；Pout=1-Pin；Pin∈[0,1]。

(2)M过程

该过程用于对原图像I的估计In进一步修正。令

(15)

等价于求得

(16)

3　实验与分析

本文实验采用清晰眼底OCT图像，高斯模糊卷积核大小是5×5，标准差值为2。处理器Inter(R) Core(TM) i5-2450M CPU@2.50 Hz 2.49 GHz,内存容量4 GB，操作系统为Windows 7(64 bit)，处理平台为Matlab R2015a。

图3(a)为视网膜黄斑图像，它是由OCT直线扫描人眼视网膜黄斑区获取。图3(b)则是输入的模糊图像。图3(c～f)分别显示了基于EM反卷积方法、利用Vladimir Yuzhikov的SmartDeblur软件(http://yuzhikov.com/articles/BlurredImagesRestoration1.htm)、Pan等人[15]的算法以及经典的Lucy-Richardson滤波[16]对图3(b)图像去模糊的效果图。从图3(c)可以看出EM反卷积算法复原的图像边缘轮廓相对清晰，几乎没有出现振铃效应；图3(d，f)虽然恢复出清晰的组织结构，但同时伴有明显的振铃效应造成干扰；从图3(e)可以看出黄斑区各层间轮廓突出，但具体细节仍需优化。

为了客观分析基于EM反卷积算法在复原眼底OCT图像中的优势，表1通过峰值信噪比(Peak signal to noise ratio, PSNR)和图像相似度(Structural similarity, SSIM)对各种方法复原模糊眼底OCT图像质量进行评估。这些评价指标广泛应用于图像处理和计算机视觉研究中，能够直接测算复原后结果与标准图像的视觉差异。其中，PSNR值越大表明去噪效果越明显，算法引入噪声越少[17]；而SSIM则是从结构、亮度和对比度衡量图像的相似性，值越大表明OCT图像复原后失真越小，与标准图像相似度更高。从表1所有算法处理后的PSNR和SSIM对比结果可以看出，基于EM反卷积算法的PSNR与SSIM值均较高，从客观上表明该算法可以很好地克服异常值对复原图像带来的干扰，从而能够恢复出清晰图像。

图3　不同算法去模糊效果图Fig.3　Comparison for reconstruction of retinal image using different deblurring algorithms

指标基于EM反卷积算法SmartDeblur算法文献[15]算法Lucy⁃Richardson算法PSNR34．913628．153617．934132．5342SSIM0．93810．87280．59880．9355

4　结束语

本文研究列举了传统图像退化模型，针对视网膜组织结构复杂对成像标准要求高这一特性，具体分析了眼底OCT设备工作原理以及造成眼底OCT图像模糊的原因，根据眼底OCT图像自身特点(不符合高斯噪声分布和过饱和情况)，详细介绍并分析了基于EM的非盲去模糊的反卷积算法，并结合3种方法作参照进行仿真实验。实验结果表明，相对于传统图像去模糊技术，该算法在抑制振铃效应及保留细节方面有着显著效果。然而生活中利用OCT获取的眼底图像一般不容易确定点扩散函数信息，且图像处理后的PSNR值较低，因此如何解决夹杂大量噪声的OCT模糊图像的恢复、准确测算实际拍摄眼底OCT图像中点扩散函数这两个问题，还需要后续研究[18]。

参考文献:

[1]刘涛.眼底OCT图像降噪及边缘检测算法研究[D].北京：清华大学，2010.

Liu Tao. Research on noise degradation and contour extraction algorithm for eye fundus OCT images[D]. Beijing: Tsinghua University, 2010.

[2]Tiwari S, Shukla V P, Biradar S R, et al. Review of motion blur estimation techniques[J]. Journal of Image and Graphics, 2013,1(4):176-184.

[3]Krishnan D, Tay T, Fergus R. Blind deconvolution using a normalized sparsity measure[C]//Computer Vision and Pattern Recognition (CVPR), 2011 IEEE Conference on. [S.l.]: IEEE, 2011:233-240.

[4]Mai L, Liu F. Kernel fusion for better image deblurring[C]//Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. [S.l.]: IEEE, 2015:371-380.

[5]Oliveira J P, Figueiredo M A T, Bioucas-Dias J M. Parametric blur estimation for blind restoration of natural images: Linear motion and out-of-focus[J]. Image Processing, IEEE Transactions on, 2014,23(1):466-477.

[6]刘新文，王惠南，钱志余.小波变换对OCT图像的降噪处理[J].光子学报,2006(6):935-939.

Liu Xinwen, Wang Huinan, Qian Zhiyu. Denoising of OCT image by wavelet transform[J]. Acta Photonica Sinica, 2006(6):935-939.

[7]Kulkarni M D, Izatt J A, Sivak M V. Image enhancement in optical coherence tomography using deconvolution: U.S. Patent 5,994,690[P]. 1999-11-30.

[8]Liu Y, Liang Y, Mu G, et al. Deconvolution methods for image deblurring in optical coherence tomography[J]. JOSA A, 2009,26(1):72-77.

[9]Liu G, Yousefi S, Zhi Z, et al. Automatic estimation of point-spread-function for deconvoluting out-of-focus optical coherence tomographic images using information entropy-based approach[J]. Optics Express, 2011,19(19):18135-18148.

[10] Huang D, Swanson E A, Lin C P, et al. Optical coherence tomography[J]. Science, 1991, 254(5035): 1178-1181.

[11] 刘扬阳，金伟其，苏秉华. 数字图像去模糊处理算法的对比研究[J]. 北京理工大学学报，2004(10):905-909.

Liu Yangyang,Jin Weiqi,Su Binghua.Comparison and study of image deconvolution algorithms[J]. Transactions of Beijing Institute of Technology,2004(10): 905-909.

[12] 田宇. 单幅运动模糊图像去模糊研究[D].长春：吉林大学，2011.

Tian Yu. Research on a single image motion deblurring[D]. Changchun: Jili University,2011.

[13] Cho S, Wang J, Lee S. Handling outliers in non-blind image deconvolution[C]//Computer Vision (ICCV), 2011 IEEE International Conference on. [S.l.]: IEEE, 2011: 495-502.

[14] Levin A, Fergus R, Durand F, et al. Image and depth from a conventional camera with a coded aperture[C]//ACM Transactions on Graphics (TOG). [S.l.]: ACM, 2007, 26(3): 70.

[15] Pan J, Liu R, Su Z, et al. Kernel estimation from salient structure for robust motion deblurring[J]. Signal Processing: Image Communication, 2013, 28(9): 1156-1170.

[16] Tai Y W, Tan P, Brown M S. Richardson-Lucy deblurring for scenes under a projective motion path[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2011, 33(8): 1603-1618.

[17] 王田芳,李浩,温四林,等. 基于Gauss滤波和Euler修复模型的SAR图像去噪[J]. 数据采集与处理,2016,31(3):562-569.

Wang Tianfang, Li Hao, Wen Silin, et al.Speckle removement of SAR image based on Euler′s elastica model and Gauss filtering[J]. Journal of Data Acquisition and Processing, 2016,31(3):562-569.

[18] Liu G, Yousefi S, Zhi Z, et al. Automatic estimation of point-spread-function for deconvoluting out-of-focus optical coherence tomographic images using information entropy-based approach[J]. Optics Express, 2011,19(19):18135-18148.