角膜神經(jīng)顯微圖像的自適應(yīng)漸暈校正

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

（1.天津工業(yè)大學(xué) 電子與信息工程學(xué)院，天津 300387；2.天津市光電檢測技術(shù)與系統(tǒng)重點(diǎn)實(shí)驗(yàn)室，天津 300387；3.天津醫(yī)科大學(xué)眼科醫(yī)院 眼視光學(xué)院，天津 300384；4.國家眼耳鼻喉疾病臨床醫(yī)學(xué)研究中心天津市分中心，天津 300384；5.天津市視網(wǎng)膜功能與疾病重點(diǎn)實(shí)驗(yàn)室，天津 300384；6.北京大學(xué)人民醫(yī)院眼科，北京 100044；7.瑞達(dá)昇醫(yī)療科技有限公司，北京 101100）

1 引言

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

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

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

基于像素分析的方法無須建立漸暈?zāi)Ｐ?，直接在空間域或頻域?qū)u暈圖像進(jìn)行校正替換操作。Zheng等人［13］提出真實(shí)圖像在像素梯度中服從稀疏概率分布，其特點(diǎn)是在梯度分布中有一個高峰度和兩個重尾。利用圖像的分段恒常性理論調(diào)整醫(yī)學(xué)圖像中的亮度不均衡現(xiàn)象，通過計算圖像梯度直方圖后擬合最大似然系數(shù)，使用迭代加權(quán)最小二乘法算法（IRLS）進(jìn)行求解。Chernavskaia等人［14］提出了一種基于傅里葉變換的頻域?yàn)V波方法校正序列圖像漸暈。序列圖像空間域中的周期性漸暈導(dǎo)致頻域中出現(xiàn)周期性結(jié)構(gòu)。通過改變頻域中的這種結(jié)構(gòu)消除圖像空間域的漸暈現(xiàn)象。但顯微拼接圖像在空間域重疊的情況下不存在周期性，故不適用于角膜神經(jīng)圖像。

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

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

2 方法原理

2.1 漸暈?zāi)Ｐ偷臉?gòu)建與補(bǔ)償

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

圖1 圖像中的漸暈效應(yīng)Fig.1 Vignetting effect in the image

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

圖2 漸暈效應(yīng)衰減模型Fig.2 Vignetting effect decay model

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

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

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

2.2 漸暈參數(shù)存在條件

漸暈效果由中心點(diǎn)向四周遞增，像素值向四周單調(diào)遞減，所以增益函數(shù)必須單調(diào)遞增：

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

利用求根公式可得：

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

2.3 漸暈校正評價函數(shù)

當(dāng)距離r變化時，上述限定條件可以保證漸暈參數(shù)取值的合理性。然而，由于漸暈?zāi)Ｐ筒捎昧A多項式很難直接求得最優(yōu)解。為此，本研究通過評估圖像信息熵來控制參數(shù)優(yōu)化過程，以最小化圖像信息熵的方法優(yōu)化漸暈校正結(jié)果。

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

文獻(xiàn)［22］證明了對數(shù)信息熵在參數(shù)優(yōu)化過程中相較于信息熵的優(yōu)越性。圖3為信息熵與分布移動的關(guān)系比較。黑色曲線代表圖像中未受到漸暈效應(yīng)影響的區(qū)域（圖像中心部分）的熵值分布；灰色曲線代表圖像含有漸暈效應(yīng)的區(qū)域（圖像邊界部分）熵的分布。通過漸暈迭代校正，圖像亮度差別逐漸消失，即灰色曲線向黑色曲線靠近直至重疊。圖3的第一行表示直接計算的信息熵。在漸暈校正的過程中，會因校正漸暈引起圖像亮度增加而造成總信息熵增加，直到兩個分布重疊時總信息熵開始減小。圖3第二行表示對數(shù)信息熵計算結(jié)果?？梢钥闯鲈谝苿舆^程中對數(shù)信息熵基本保持不變，直到兩個分布開始重疊，對數(shù)信息熵才開始減小。因此，對數(shù)信息熵在優(yōu)化過程中具有單調(diào)性。

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

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

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

假設(shè)像素的對數(shù)信息熵趨于高斯分布，為減小量化誤差，使用標(biāo)準(zhǔn)高斯核Gσ平滑直方圖。用p?k表示每個直方圖區(qū)間的離散概率。

則對數(shù)信息熵H為：

2.4 模型優(yōu)化

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

輸入：漸暈圖像I。

輸出：校正圖像C。

Step1：設(shè)置漸暈參數(shù)初始值a1，a2，a3=（0，0，0），優(yōu)化步長δ=1/2，對數(shù)強(qiáng)度熵最小值Hmin=H（I），C=I，停止閾值T=1/256。Step2：當(dāng)δ＞T，計算6組漸暈參數(shù)Vi，

否則，結(jié)束優(yōu)化。

Step3：如果Vi滿足公式（6）中的條件之一，計算對數(shù)強(qiáng)度熵H=H（Vi）；若均不滿足條件，轉(zhuǎn)Step2。

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

否則δ=0.5*δ，轉(zhuǎn)至Step2。

實(shí)驗(yàn)過程中，遞進(jìn)步長δ的初始值為1，每次迭代減為之前的一半，即δi=0.5*δi-1。圖4顯示了一組迭代優(yōu)化過程中對數(shù)強(qiáng)度熵的變化。前期實(shí)驗(yàn)發(fā)現(xiàn)，當(dāng)?shù)螖?shù)大于40時對數(shù)強(qiáng)度熵趨于平穩(wěn)。

圖4 對數(shù)強(qiáng)度熵迭代過程Fig.4 Logarithmic intensity entropy iterative process

3 實(shí)驗(yàn)與結(jié)果分析

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

3.1 實(shí)驗(yàn)材料

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

3.2 主觀性評價

圖5為一組漸暈校正前后角膜顯微圖像對比圖。將單張原始圖像先進(jìn)行漸暈校正處理，再進(jìn)行手動拼接?？梢灾庇^看出，校正后圖像亮度分布均勻，且圖像拼接處無偽影，銜接自然。

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

3.3 客觀性評價

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

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

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

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

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

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

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

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

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

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

圖9展示了2組漸暈校正前后手動拼接的角膜神經(jīng)顯微圖像。第一行圖像是由未進(jìn)行漸暈校正的角膜神經(jīng)顯微圖像拼接成的圖像。第二行是由本文方法漸暈校正處理后進(jìn)行的拼接圖像?？梢钥闯霰疚姆椒軌蚋纳破唇訄D像銜接處的偽影問題。然而在實(shí)踐中我們也發(fā)現(xiàn)校正前后待拼接圖像間整體灰度值差別較大，影響后 續(xù)神經(jīng)提取。為此在圖像拼接前，我們對漸暈校正后圖像的像素均值進(jìn)行了圖像直方圖均衡化處理.

圖9 角膜神經(jīng)顯微拼接圖像Fig.9 Corneal nerve microscope stitching images

4 結(jié)論

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