引言

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

图像去模糊是图像处理方面基本问题之一，根据点扩散函数是否已知，将其分为盲去模糊与非盲去模糊。在非盲去模糊方面，Tiwari等人在假设没有高斯噪声影响的前提下，对比分析了不同的运动模糊参数估计方法[。同样在盲去模糊方面，Krishnan等人运用了一种归一化的稀疏度量方法来解决图像模糊问题[；Mai和Liu结合多种去模糊方法将它们得到的模糊核融合并提出了一种数据驱动方法，该方法能够有效地从训练集中学习得到融合后的模糊核[；Oliveira等人基于模糊图像的频谱和弱假设检验提出了对线性运动和散焦造成模糊参数的估计方法[，由于该方法不需要迭代，复原运算速度得到显著提升。对于OCT成像，目前已经出现了很多应用于图像的降噪和轮廓信息的提取等方面的技术

[。而在OCT图像去模糊方面，Kulkarni等人首次提出了用反卷积的方法去除图像模糊[；Liu等人同时针对纵向、横向相互独立的模糊核，利用维纳滤波和Lucy-Richardson滤波分别进行OCT图像去模糊，并对两种方法的效果进行了对比[；Liu等人采用了基于信息熵的方法对因失焦产生模糊OCT图像中的点扩散函数进行自动估计[。尽管如此，OCT图像去模糊的研究特别是结合OCT图像的特征属性进行针对性去模糊的处理依然比较少见。

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

1 基础理论

1.1 OCT技术原理及模糊原因

在眼科临床诊断中，OCT利用低相干光干涉成像原理，以光学组织切片的形式对视网膜成像。目前OCT系统都是基于低干涉光束，比如Michelson干涉仪。由OCT技术原理图(见[。

导致眼底OCT图像模糊的原因主要是运动模糊和散焦模糊。运动模糊通常是由于眼球运动造成，由于新型的OCT设备中带有自动追踪眼球运动的模块，因运动造成的模糊OCT图像比较少见。因此，本文主要关注由于散焦造成模糊的眼底OCT图像，如

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

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

1.2 传统图像退化模型

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

$

b\left( {x,y} \right) = I\left( {x,y} \right) * k\left( {x,y} \right) + n\left( {x,y} \right)

$

(1)

式中“*”为卷积运算。

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

[。线性运动的点扩散函数公式为

$

k\left( {x,y} \right) = \left\{ \begin{array}{l}

1/L\;\;\;\;\sqrt {{x^2} + {y^2}} \le \frac{L}{2}且\frac{x}{y} = - \tan \theta \\

0\;\;\;\;\;\;\;\;\;其他

\end{array} \right.

$

(2)

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

$

k\left( {x,y} \right) = \left\{ \begin{array}{l}

1/{\rm{ \mathsf{ π} }}{r^2}\;\;\;\;{x^2} + {y^2} \le {r^2}\\

0\;\;\;\;\;\;\;\;\;\;\;其他

\end{array} \right.

$

(3)

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

$

k\left( {x,y} \right) = \left\{ \begin{array}{l}

k{{\rm{e}}^{ - a\left( {\sqrt {{x^2} + {y^2}} } \right)}}\;\;\;\;\left( {x,y} \right) \in C\\

0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;其他

\end{array} \right.

$

(4)

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

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

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

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

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

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

$

b = {\rm{c}}\left( {k * I} \right) + n

$

(5)

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

$

{I_{{\rm{MAP}}}} = \arg \mathop {\max }\limits_I P\left( {I\left| {k,b} \right.} \right)

$

(6)

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

$

{I_{{\rm{MAP}}}} = \arg \mathop {\max }\limits_I \sum\limits_{m \in M} {P\left( {b,m\left| {k,I} \right.} \right)P\left( I \right)}

$

即

$

{I_{{\rm{MAP}}}} = \arg \mathop {\max }\limits_I \sum\limits_{m \in M} {P\left( {b\left| {m,k,I} \right.} \right)P\left( {m\left| {k,I} \right.} \right)P\left( I \right)}

$

(7)

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

(1) E过程

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

② 定义

$

Q\left( {I,{I^0}} \right) = E\left[ {\log P\left( {b,m\left| {k,I} \right.} \right)} \right] = E\left[ {\log P\left( {b\left| {m,k,I} \right.} \right) + \log P\left( {m\left| {k,I} \right.} \right)} \right]

$

(8)

式中E表示期望。

③ 根据文献[

$

P\left( {{b_x}\left| {m,k,I} \right.} \right) = \left\{ \begin{array}{l}

N\left( {{b_x}\left| {{f_x},\sigma } \right.} \right)\;\;\;\;\;\;{m_x} = 1\\

C\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{m_x} = 0

\end{array} \right.

$

(9)

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

$

P\left( {m\left| {k,I} \right.} \right){\rm{ = }}\prod\limits_x {P\left( {{m_x}\left| {k,I} \right.} \right)} = \prod\limits_x {P\left( {{m_x}\left| {{f_x}} \right.} \right)}

$

(10)

$

P\left( {{m_x} = 1\left| {{f_x}} \right.} \right) = \left\{ \begin{array}{l}

{P_{{\rm{in}}}}\;\;\;\;{f_x} \in \left[ {0,1} \right]\\

0\;\;\;\;\;\;其他

\end{array} \right.

$

(11)

因此有

$

Q\left( {I,{I^0}} \right) = E\left[ {\sum\limits_x {{m_x}\log N\left( {{b_x}\left| {{f_x},\sigma } \right.} \right)} } \right] = - \sum\limits_x {\frac{{E\left[ {{m_x}} \right]}}{{2{\sigma ^2}}}{{\left| {{b_x} - {f_x}} \right|}^2}}

$

(12)

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

④ 根据贝叶斯公式

$

P\left( {{m_x}\left| {b,k,{I^0}} \right.} \right) = \frac{{P\left( {{b_x}\left| {{m_x},k,{I^0}} \right.} \right)P\left( {{m_x}\left| {k,{I^0}} \right.} \right)}}{{\sum\limits_{{m_x} = 0}^1 {P\left( {{b_x}\left| {{m_x},k,{I^0}} \right.} \right)P\left( {{m_x}\left| {k,{I^0}} \right.} \right)} }}

$

(13)

因此可以进一步求得

$

E\left[ {{m_x}} \right] = \left\{ \begin{array}{l}

\frac{{N\left( {{b_x}\left| {f_x^0,\sigma } \right.} \right){P_{{\rm{in}}}}}}{{N\left( {{b_x}\left| {f_x^0,\sigma } \right.} \right){P_{{\rm{in}}}} + C{P_{{\rm{out}}}}}}\;\;\;\;\;\;f_x^0 \in \left[ {0,1} \right]\\

0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;其他

\end{array} \right.

$

(14)

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

(2) M过程

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

$

{I^n} = \arg \mathop {\max }\limits_I \left\{ {Q\left( {I,{I^0}} \right) + \log P\left( I \right)} \right\}

$

(15)

等价于求得

$

\sum\limits_x {\omega _x^m{{\left| {{b_x} - {{\left( {k * I} \right)}_x}} \right|}^2}} + \lambda \phi \left( I \right)

$

(16)

的最小值，从而近似为$\sum\limits_x {\omega _x^m{{\left| {{b_x}-{{\left( {k*I} \right)}_x}} \right|}^2} + \lambda \psi \left( I \right)} $。其中$\omega _x^m = E\left[{{m_x}} \right]/2{\sigma ^2}, \omega \left( I \right) = \sum {\left\{ {{{\left| {{\nabla ^h}I} \right|}^{\partial -2}} \cdot {{\left| {{{\left( {{\nabla ^h}I} \right)}_x}} \right|}^2} + {{\left| {{{\left( {{\nabla ^v}I} \right)}_x}} \right|}^{\partial -2}}{{\left| {{{\left( {{\nabla ^v}I} \right)}_x}} \right|}^2}} \right\}} $。

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。

http://yuzhikov.com/articles/BlurredImagesRestoration1.htm)、Pan等人[的算法以及经典的Lucy-Richardson滤波[对

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

为了客观分析基于EM反卷积算法在复原眼底OCT图像中的优势，[；而SSIM则是从结构、亮度和对比度衡量图像的相似性，值越大表明OCT图像复原后失真越小，与标准图像相似度更高。从

表 1(Tab. 1)

表 1 眼底OCT图像去模糊效果客观指标

Tab. 1 Deblurring indexes of OCT retinal image

指标

基于EM反卷积算法

SmartDeblur算法

文献[

Lucy-Richardson算法

PSNR

34.913 6

28.153 6

17.934 1

32.534 2

SSIM

0.938 1

0.872 8

0.598 8

0.935 5

表 1 眼底OCT图像去模糊效果客观指标

Tab. 1 Deblurring indexes of OCT retinal image

4 结束语

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