角膜神經顯微圖像的自適應漸暈校正

李天宇，李光旭，2*，張 琛，李方烴，李德衡

（1.天津工業大學 電子與信息工程學院，天津 300387；2.天津市光電檢測技術與系統重點實驗室，天津 300387；3.天津醫科大學眼科醫院 眼視光學院，天津 300384；4.國家眼耳鼻喉疾病臨床醫學研究中心天津市分中心，天津 300384；5.天津市視網膜功能與疾病重點實驗室，天津 300384；6.北京大學人民醫院眼科，北京 100044；7.瑞達昇醫療科技有限公司，北京 101100）

1 引言

角膜中存在豐富的感覺神經末梢，這些角膜神經可調控淚液分泌，為上皮細胞和基質細胞提供營養支持［1］。同時，角膜神經功能障礙可引起結構變化的出現［2］。角膜也是人體唯一能夠直接觀察到神經末梢的組織。角膜神經的形態及密度變化可為神經類疾病，如糖尿病周圍神經病變［3］、退行性帕金森癥［4］等對神經系統的影響提供直接證據。此外，圓錐角膜等疾病也會引發角膜局部形態病變，導致角膜神經曲折度明顯增加，交織成網狀，臨床表現為不規則散光等［5］。

近年來，激光共聚焦顯微鏡（Laser Scanning Confocal Microscopy，LSCM）的廣泛應用［6］使對角膜神經形態特征的研究得以實現［7］。由于單張角膜共聚焦顯微圖像成像范圍較小［8］，為獲得更大范圍內角膜神經的整體形態分布，臨床上利用圖像拼接技術在保持相同分辨率的前提下擴大角膜圖像的可視區域［9］。然而，隨著離軸距離的增加，光線經過光學系統的有效孔徑減小，導致了從中心到邊緣非線性的緩慢的光衰減，圖像呈現出中心亮四周暗的特征，稱為光學漸暈效應［10-11］。該效應會引入拼接偽影，增大角膜神經的重構難度。

圖像漸暈校正算法可分為兩大類［12］：基于像素分析的校正方法和基于漸暈數學模型的校正方法。

基于像素分析的方法無須建立漸暈模型，直接在空間域或頻域對漸暈圖像進行校正替換操作。Zheng等人［13］提出真實圖像在像素梯度中服從稀疏概率分布，其特點是在梯度分布中有一個高峰度和兩個重尾。利用圖像的分段恒常性理論調整醫學圖像中的亮度不均衡現象，通過計算圖像梯度直方圖后擬合最大似然系數，使用迭代加權最小二乘法算法（IRLS）進行求解。Chernavskaia等人［14］提出了一種基于傅里葉變換的頻域濾波方法校正序列圖像漸暈。序列圖像空間域中的周期性漸暈導致頻域中出現周期性結構。通過改變頻域中的這種結構消除圖像空間域的漸暈現象。但顯微拼接圖像在空間域重疊的情況下不存在周期性，故不適用于角膜神經圖像。

基于漸暈數學模型的校正方法是先對漸暈效應建模，再將漸暈模型函數的逆函數作為增益函數，反向補償原圖像達到校正的作用。D'Angelo等人［15］參考cos4θ模型函數，建立了高次線形徑向漸暈模型，并用于相機光學系統漸暈效果校正。Goldman等 人［11］根 據 相 機 的 光 圈 大 小、焦距、輻射亮度和快門持續時間等外部固定參數來建立漸暈模型，但對于不同的相機拍攝圖像需進行不同的設置。Smith等人［16］提出了一種使用正則化能量最小化（CIDRE）的方法來校正圖像光照強度分布。首先構造一個線性光照強度模型，在求解模型時使用穩健的正則化能量函數，并使用Limited-memory Broyden-fletcher-goldfarbshan（L-BFGS）數值優化算法求解優化。Peng等人［17］提出了一種基于低秩和稀疏分解的圖像漸暈校正方法（BaSiC）。首先構造一個測量矩陣。然后分解成一個低秩矩陣和稀疏殘差矩陣，使用線性化的拉格朗日增廣方法迭代方式求解優化，最終獲得漸暈陰影模型。但需要多幅連續的圖像且每個要處理的圖像的前景與其他每個圖像的前景不相關，否則該方法會將視野中心始終較高的圖像強度視為局部亮度的增加，從而消除了真實的圖像亮度變化。

上述算法雖然可以解決單個或組合序列圖像的漸暈問題，但需針對特定圖像進行參數設置。如設置不當，會出現圖像欠校正或過度校正等問題。對于角膜神經圖像而言，醫生采樣角膜神經進行鏡頭移動時，采樣移動方向不固定，導致角膜神經圖像序列不規律。此外，因顯微鏡采樣速度較快，圖像間相關性較高等原因，導致采樣圖像通常無法滿足以上算法條件。因此，本文提出一種角膜神經顯微圖像自適應漸暈校正算法，針對共聚焦顯微鏡的漸暈效應，提出使用六階多項式對漸暈模型進行擬合，更好地保留原圖像未受漸暈影響區域的特征并校正漸暈區域。在校正過程中使用對數信息熵對建模效果進行評估反饋，增強了漸暈校正效果的穩定性。

2 方法原理

2.1 漸暈模型的構建與補償

圖1（a）是含有漸暈效應的角膜神經共聚焦顯微圖像，直觀上可以看出圖像邊界區域比中心區域亮度低。漸暈效果以圖像光學中心為原點，向外漸暈強度逐漸增加。漸暈呈對稱性分布，且漸暈強度隨徑向距離的增加而增大。圖1（b）為拼接后角膜神經顯微圖像。圖中漸暈失真會在圖像拼接處產生偽影，造成角膜神經結構不清晰。圖1（c）表示圖1（a）圖中對角線處（黃色虛線）像素的亮度變化趨勢。藍色實線表示對角線像素值曲線，紅色虛線表示像素變化擬合曲線。圖像中心部分區域為非漸暈區域，漸暈校正方法的目的為調整漸暈區域像素值，使其接近非漸暈區域像素值。

圖1 圖像中的漸暈效應Fig.1 Vignetting effect in the image

共聚焦顯微鏡的光源光路與成像光路相同，其成像亮度值呈現中心高并向四周逐漸減弱的特點［18］，如圖2所示。通常，受光路影響顯微鏡組的光學中心（漸暈中心）與獲取的圖像中心有坐標偏差，其偏移量最大可達25個像素［19］。然而，共聚焦顯微鏡采用平移掃描的方式獲取圖像。且單次成像范圍極小，成像位置偏差可以忽略。因此，可將光學中心與圖像中心視為同軸。

圖2 漸暈效應衰減模型Fig.2 Vignetting effect decay model

基于漸暈數學模型的校正方法是將漸暈函數的逆函數作為增益函數M，反向補償原圖像像素灰度值。若用Iorg（x，y）表示包含光學漸暈效應圖像。漸暈校正后圖像為Ires（x，y），則：

本研究利用相機響應模型EMoR［20］來擬合漸暈效果，該模型為含有偶次項的六階多項式函數，可以用非常少的參數來模擬范圍廣泛的響應函數，并且在準確性方面優于其他多項式和非參數模型。進而求得其增益補償函數為：

其中：r是圖像像素點到圖像中心點的歸一化距離，令圖像中心點處r=0，圖像四個角點處的r=1；系數ai（i=1，2，3）稱為漸暈參數，通過迭代優化求解。

2.2 漸暈參數存在條件

漸暈效果由中心點向四周遞增，像素值向四周單調遞減，所以增益函數必須單調遞增：

若令q=r2，則式（3）可化簡為：

利用求根公式可得：

根據解的分布，得出當漸暈參數滿足如下條件之一時，不等式（4）成立。

2.3 漸暈校正評價函數

當距離r變化時，上述限定條件可以保證漸暈參數取值的合理性。然而，由于漸暈模型采用六階多項式很難直接求得最優解。為此，本研究通過評估圖像信息熵來控制參數優化過程，以最小化圖像信息熵的方法優化漸暈校正結果。

Likar等［21］提出通過最小化圖像的信息熵對圖像像素強度一致性進行評估，并將其用于漸暈偽影的校正模型的評估。通常，漸暈圖像周圍的暗角會使圖像產生額外的信息熵；當漸暈被校正同時信息熵應該逐漸減小。然而，在此過程中過度校正會引入部分使圖像增亮的信息熵，導致校正算法無法實現全局最優。

文獻［22］證明了對數信息熵在參數優化過程中相較于信息熵的優越性。圖3為信息熵與分布移動的關系比較。黑色曲線代表圖像中未受到漸暈效應影響的區域（圖像中心部分）的熵值分布；灰色曲線代表圖像含有漸暈效應的區域（圖像邊界部分）熵的分布。通過漸暈迭代校正，圖像亮度差別逐漸消失，即灰色曲線向黑色曲線靠近直至重疊。圖3的第一行表示直接計算的信息熵。在漸暈校正的過程中，會因校正漸暈引起圖像亮度增加而造成總信息熵增加，直到兩個分布重疊時總信息熵開始減小。圖3第二行表示對數信息熵計算結果。可以看出在移動過程中對數信息熵基本保持不變，直到兩個分布開始重疊，對數信息熵才開始減小。因此，對數信息熵在優化過程中具有單調性。

圖3 信息熵與分布移動的關系［22］Fig.3 Relationship between information entropy and distribution movement［22］

信息熵與對數信息熵計算公式雖然相同，在計算對數信息熵時，首先要將圖像的像素值L進行對數映射，如公式7所示。

其中：N為映射后的像素階數，i為對數映射后的像素值。在統計其直方圖時，nk表示每一階像素直方圖的統計值。運算符■■與■■分別表示向上取整和向下取整。

假設像素的對數信息熵趨于高斯分布，為減小量化誤差，使用標準高斯核Gσ平滑直方圖。用p?k表示每個直方圖區間的離散概率。

則對數信息熵H為：

2.4 模型優化

采用Levenberg-Marquardt非線性優化方法獲取參數的最優解。獲得所需漸暈校正模型需要設定漸暈參數ai的初始值及遞進步長δ。每個參數都會獨立地增減相同的δ。每輪優化會產生6組計算參數，記作Vi（i=1，2…，6）。選取6組參數中最小對數信息熵對應的漸暈參數，直至H收斂。具體校正算法如下：

輸入：漸暈圖像I。

輸出：校正圖像C。

Step1：設置漸暈參數初始值a1，a2，a3=（0，0，0），優化步長δ=1/2，對數強度熵最小值Hmin=H（I），C=I，停止閾值T=1/256。Step2：當δ＞T，計算6組漸暈參數Vi，

否則，結束優化。

Step3：如果Vi滿足公式（6）中的條件之一，計算對數強度熵H=H（Vi）；若均不滿足條件，轉Step2。

Step4：如果H＜Hmin，Hmin=H；C=I·MVi；

否則δ=0.5*δ，轉至Step2。

實驗過程中，遞進步長δ的初始值為1，每次迭代減為之前的一半，即δi=0.5*δi-1。圖4顯示了一組迭代優化過程中對數強度熵的變化。前期實驗發現，當迭代次數大于40時對數強度熵趨于平穩。

圖4 對數強度熵迭代過程Fig.4 Logarithmic intensity entropy iterative process

3 實驗與結果分析

本文以角膜共聚焦顯微鏡圖像為研究對象，完成單幅圖像的漸暈校正及全景圖像拼接。為了更好地驗證算法的有效性與可靠性，將提出的算法與針對相同圖像特性（非自然圖像）的算法進行對比，包括基于醫學圖像漸暈校正的非建模算 法［13］，基 于 顯 微 圖 像 漸 暈 校 正 的CIDRE算法［16］和BaSiC算法［17］。

3.1 實驗材料

實驗選取5名患者雙眼的角膜共聚焦顯微圖像，每組包含圖像的數量為300～900張不等。圖像由海德堡視網膜激光斷層掃描系統（Ⅱ代）獲取。視場大小為400 μm×400 μm，最高分辨率為1 μm。

3.2 主觀性評價

圖5為一組漸暈校正前后角膜顯微圖像對比圖。將單張原始圖像先進行漸暈校正處理，再進行手動拼接。可以直觀看出，校正后圖像亮度分布均勻，且圖像拼接處無偽影，銜接自然。

圖5 漸暈校正前后圖像Fig.5 Images before and after vignetting correction

3.3 客觀性評價

為了客觀評估漸暈校正對于角膜神經圖像的有效性，我們對比了校正前后2幅圖像局部灰度值變化。如圖6（a）所示，以圖像中心為原點，求取固定像素寬度的環形區域內像素的平均值。圖6（b）為當環形區域寬度取10個像素值時，校正前后像素均值的變化曲線。由圖6（b）漸暈圖像曲線可以看出像素均值由內到外呈遞減趨勢。表明漸暈效果是從圖像中心向外亮度逐漸降低，且隨著距離增大，漸暈強度逐漸增加，符合漸暈效果的物理性質。BaSiC方法與CIDRE方法雖然可以較好地校正漸暈，但中心區域像素明顯降低，改變了原始圖像的亮度，會對后續醫生診斷產生影響。非建模方法雖然中心區域像素均值下降不明顯，但隨著圈層數增加，漸暈校正效果有所降低。本文提出的方法校正后圖像周圍區域像素無明顯增減趨勢，亮度均衡。可證明漸暈校正效果良好。圖7為不同漸暈校正方法的結果。

圖6 漸暈校正圖像像素法評估Fig.6 Pixel method evaluation of vignetting corrected image

圖7 不同漸暈校正方法的結果Fig.7 Results of different vignetting correction methods

本文收集了5例不同患者的角膜神經圖像，每例圖像的光照和采集深度條件均不同。為了客觀評估漸暈校正效果的有效性和可靠性，我們以正常圖像為基準圖像，再對基準圖像人為添加漸暈效果，獲得模擬漸暈圖像供校正使用。將正常圖像與校正后圖像比較計算均方誤差（Mean Square Error，MSE）、峰值信噪比（Peak Signal to Noise Ratio，PSNR）、結構相似性（Structural Similarity，SSIM）［23］的平均值，對比評估不同漸暈校正方法結果的圖像質量。

本論文引入SSIM對校正前后圖片進行相似性計算，從亮度、對比度、結構三個方面評估漸暈校正算法的可靠性，分析算法對圖像質量的影響。公式（12）為圖像I1，I2間的結構相似性指數的計算公式。

其中：μ1，μ2，σ1，σ2分別為對應圖像的均值和方差；σ12為協方差。用均值、方差、協方差作為圖像的亮度、對比度、結構的估計。SSIM取值范圍為［0，1］。SSIM值越大表示兩幅圖像相似性越高，亮度與結構誤差越小。

本文同時采用MSE計算圖像像素差異的L2范數。MSE值越小代表圖像像素誤差越小，圖像質量越好。PSNR能夠反應圖像的失真程度，PSNR值越大，圖像失真越小。我們將圖像PSNR范圍限定在0～100 dB。

5組圖像的數量分別為587、943、371、593和754組，其中病例A和病例B為正常采集組，其余病例圖像采集特性分別為采集亮度低、采集亮度高、采集深度淺（部分顯微圖像出現噪聲）。五名患者校正效果評估平均結果如表1所示，本文方法處理的5名患者校正后圖像MSE、PSNR、SSIM平均值均高于其他方法，即校正后圖像與基準圖像最相似，具有最佳的校正效果。由此也可以表明本文方法漸暈校正后圖像并未出現較大程度失真，不會改變原始圖像中的亮度與結構等信息。圖8繪制了5組圖像校正后的評估指標箱型圖。從圖中可以看出本文方法相較于其他算法更加穩定，不會受其他外部因素影響，魯棒性較高。

表1 五名患者校正效果評估平均結果Tab.1 Mean results of correction effect evaluation in five patients

3.4 角膜共聚焦顯微圖像拼接

圖9展示了2組漸暈校正前后手動拼接的角膜神經顯微圖像。第一行圖像是由未進行漸暈校正的角膜神經顯微圖像拼接成的圖像。第二行是由本文方法漸暈校正處理后進行的拼接圖像。可以看出本文方法能夠改善拼接圖像銜接處的偽影問題。然而在實踐中我們也發現校正前后待拼接圖像間整體灰度值差別較大，影響后 續神經提取。為此在圖像拼接前，我們對漸暈校正后圖像的像素均值進行了圖像直方圖均衡化處理.

圖9 角膜神經顯微拼接圖像Fig.9 Corneal nerve microscope stitching images

4 結論

本文根據共聚焦顯微鏡成像的特點，提出一種角膜神經顯微圖像的漸暈校正方法，介紹了漸暈效應并使用六階偶次多項式建立漸暈效應的數學模型，采用Levenberg-Marquardt優化方法進行參數迭代優化，并設置限制條件進行約束。在優化過程中使用對數信息熵對漸暈數學模型進行評估。實驗證明本文所述算法能很好地消除角膜神經圖像的漸暈效果，避免角膜拼接圖像產生偽影。校正后圖像MSE、PSNR、SSIM評估指標平均值分別達到0.004 2、72.225 1 dB、0.960 0。漸暈校正后具有大視野范圍的角膜神經顯微拼接圖像滿足醫生臨床診斷需求。