引言
光學相幹層析成像技術(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 [email protected] 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圖像中點擴散函數這兩個問題，還需要後續研究[。