實時多聚焦顯微圖像融合算法

(1.浙江工業大學 計算機科學與技術學院，浙江 杭州 310023；2.浙江在線 新聞中心，浙江 杭州310039；3.浙江中控技術股份有限公司 研發中心，浙江 杭州 310052)

顯微鏡在生物醫療、工業界等不同領域有著廣泛的應用，因為顯微鏡平臺高倍率條件下景深很小，難以觀看樣品全貌，所以需要研究多聚焦顯微鏡圖像融合問題。目前，針對圖像融合問題已經有很多研究[1]，然而多聚焦顯微鏡圖像需要細節表現力強的算法，以保證融合圖像的質量，目前多聚焦顯微鏡圖像融合仍是一個挑戰性難題。

針對圖像融合的研究廣義上主要可以分為兩類：基于變換域的方法和基于空間域的方法。基于變換域的方法主要使用金字塔分解、小波變換和多尺度幾何分析。這些方法包括拉普拉斯金字塔[2]、梯度金字塔[3]、離散小波變換[4]、對偶樹復小波變換[5]、曲波變換[6]和ND濾波[7]等。基于變換域的多焦點圖像融合方法的基本框架是：首先將源圖像分解為一個多尺度域，在多個尺度上進行合成，然后通過反變換對融合后的圖像進行重構。變換域的選擇很重要，高頻低頻域的融合規則也很重要。基于空間域的方法通過融合聚焦度量方法直接估算圖像的銳度。一般來說，其可進一步分為基于像素[8-10]的方法和基于塊[11-12]的方法。這兩類方法各有優缺點，基于像素的方法雖然可以顯著提高融合后圖像的視覺效果，但是它對噪聲極其敏感，因為噪聲點處也有高頻信號，導致合成圖像空間上不連續，文獻[13]使用塊的一致性校驗，文獻[14]使用去除面積小于1%的塊以得到更連續的融合結果，不過在不同類型圖像中需要調整參數才能得到較好的結果。基于塊的方法通常將圖像分為M×N塊，通過分塊計算塊的銳度，可有效避免少數噪聲點帶來的影響，但是此類方法容易產生塊狀效應，且塊的大小對融合結果的影響很大。因為空間域的方法是直接估算圖像銳度，所以聚焦度量方法在此處很重要。通常通過計算圖像的梯度能量(EOG)、拉普拉斯能量(EOL)、改進拉普拉斯和(SML)、空間頻率(SF)和方差等來度量圖像的聚焦程度。近些年由于深度學習的發展，又出現了各類通過訓練神經網絡進行融合的方法[15-18]，其思路與傳統方法一致，只不過使用神經網絡決策得到聚焦度量，一般有將整幅圖像作為輸入、將圖像切分為多個塊進行輸入以及將邊緣檢測后的圖像進行輸入等，輸出一張決策圖或者某一點的聚焦結果。基于神經網絡的方法雖然在多數圖像上可以獲得很好的效果，但是由于存在大量的參數，計算復雜度很高，難以用于實時應用場景。針對上述問題，筆者提出了一種快速多聚焦圖像融合算法。該算法采用拉普拉斯金字塔分解的同時，利用一種聚類算法來引導金字塔圖像的融合。實驗結果表明：算法可以有效提高融合圖像的連續性，減小噪聲帶來的影響，并很好地保留圖像的細節特征。此外，算法多處采用了積分圖像的思想，大大提高了計算性能，在高分辨率下也可滿足實時融合的應用場景。

1 算法概述

在實時系統中，算法的輸入為相機數據流，當前輸入的圖像與上一張相似度很高時應不參與融合，由于在體視顯微鏡下序列圖像存在較大程度的偏移，所以融合前還需要進行圖像配準，配準完成后再進行圖像融合。算法流程圖如圖1所示。

圖1 算法流程圖Fig.1 Flow chart of algorithm

2 算法細節

2.1 圖像配準

針對圖像配準有基于灰度模板和特征點的算法。基于灰度模板的算法主要通過灰度匹配，使用誤差函數來判斷匹配程度；基于特征點的方法有SIFT，SURF和ORB等，此類方法基本思想是：首先找出圖像中穩定的關鍵點，根據特征點給它一個向量用于描述，然后通過這個描述進行特征點的匹配，估計出圖像的變換矩陣。然而，在顯微圖像融合的應用場景中有如下兩個特殊之處：

1) 由于光源是穩定的，視角是固定的，不會有透視變換和旋轉，只會有位移以及光圈彌散，即視覺上的輕微放大。大量實驗表明光圈引起的放大不會對圖像有顯著影響。

2) 由于聚焦的變化，特征點檢測也會有較大的變化，假設有兩幅圖像，第1幅圖像左半邊比較清晰，右半邊比較模糊，第2幅圖像則相反，左半邊比較模糊，右半邊比較清晰，這時候在進行圖像配準時，用基于特征點的方法會發現第1幅圖像找出來的特征點集中在左邊，而第2幅圖像特征點集中在右邊，此時特征點匹配會比較困難。同時，由于基于特征點的方法可以應對相機各種類型的變換，如平移、縮放、旋轉和透視等，所有的這些參數組合進算法中，使得算法的復雜度提高，計算量巨大，不太適用于實時應用程序。因此，筆者采用了基于多級下采樣與最大相關性方法，由粗到精進行匹配，大大減少了計算復雜度。

2.2 聚焦度量

對于多焦圖像融合，需要先對圖像進行聚焦度量，再對聚焦區域進行分割和圖像融合。因此，選擇合適的聚焦度量方法對后續過程至關重要。筆者分析了幾種典型的度量方法：EOG，EOL，SF和SML。

EOG的計算式為

(1)

式中：I為源圖像；w為窗口大小；(i,j)為窗口中心位置。

EOL和SF的計算式分別為

(2)

(3)

式中：RF為行頻率；CF為列頻率。其計算式分別為

(4)

(5)

SML的計算式為

(6)

(7)

4 種聚焦度量方法只有SML在計算完每個像素點的銳度之后進一步考慮了鄰域內的銳度信息，因此使用SML可以獲得更好的融合結果。如圖2所示，除SML外，其他方法的融合均出現了不同程度的白色偽影，并且采用SML方法能保留更多的細節。因此，筆者采用SML作為聚焦度量方法，并使用積分圖像進行優化。

圖2 不同聚焦度量方法融合結果對比Fig.2 Comparison of the fusion results by different focus measurements

圖3 積分圖像計算示意圖Fig.3 Schematic diagram of integral image calculation

SML(x,y)=S(A)+S(D)-S(B)-S(C)

(8)

2.3 掩膜圖像

因為聚焦的圖像擁有更大的銳度，所以圖像SML(x,y)越大則它的銳度越大，初始的掩膜圖像M0計算式為

M0(x,y)=argnmax|SMLn(x,y)|-1n∈[1,2]

(9)

因為M0易受到噪聲的影響，導致某些位置的掩膜圖像數據錯誤，影響融合結果的連續性，所以通過對以R為半徑的窗口內的類別進行聚類以提高連續性，圖像對比如圖4所示。圖4(a)為聚類前的掩膜圖像，圖4(b)為經過聚類修正后的掩膜圖像，可以看到修正后的圖像中很多小的顆粒都消失了，表明使用聚類后的掩膜圖像可以得到更連續的融合結果。

圖4 聚類前后掩膜圖像對比Fig.4 Comparison of the mask images before and after clustering

(10)

因為M0中只包含0和1，所以可以先計算M0的積分圖像，積分圖像計算方法見2.2節，然后只需要通過3 個加法計算就可以得到s。此時mask(x,y)等于0或1則代表該點是一個種子點，進而對每個種子點執行運算，把周圍的點歸到簇中，運算公式為

(11)

式中：x,y為mask圖像每個種子點的坐標位置。將未歸類的點，即mask圖像中為-1的點歸到最近棋盤距離的簇中，距離度量計算式為

d=max(|x1-x2|,|y1-y2|)

(12)

至此，mask為一張只有0和1的二值圖像，每個點代表聚焦程度更高的圖像索引。

2.4 拉普拉斯金字塔圖像融合

首先進行高斯金字塔的分解，高斯金字塔中的圖像由源圖像通過高斯濾波并以去掉偶數行列的下采樣方式得到，故下一層圖像的大小永遠是上一層圖像的1/4，其中：Gl為高斯金字塔中的第l層圖像，G0為原始圖像。具體計算式為

(13)

式中：n為最大金字塔層數；Rl，Cl分別為第l層金字塔圖像的行和列數；w為高斯核，其計算式為

(14)

由于在采樣的過程中已經丟失了一些信息，為了在之后合成時減少這些信息的丟失，所以需要將這些信息保存下來，這些信息就形成了拉普拉斯金字塔。

然后將上一層的圖像行列擴張為兩倍，新增的行列以0填充，再使用先前的卷積核乘以4對圖像作卷積，計算式為

(15)

新的圖像與下一層的原圖相減即為拉普拉斯金字塔保存的丟失信息，計算式為

(16)

(17)

(18)

3 實驗結果與分析

實驗中，分別采用LP[2]，GFF[8]，DCT[11]，ECNN[15]以及筆者方法對3 組圖像進行融合測試，除筆者方法外的其他實現均來自于公開的代碼，實驗運行環境為i5 8400的CPU+8 G內存，操作系統為Win 10，使用C++實現。

3.1 實驗數據和誤差度量方法

為了比較的一致性，樣本配準后對融合結果進行比較。總共有3 組圖像(包括2 組電路板與1 組植物切片)，電路板每組有5 張不同焦面且尺寸為5 472×3 648的24 位真彩色圖像，樣本數據如圖5，6所示。因為植物切片是一個斜面，所以序列圖像較多，共11 張，取序列中第1，3，5，7，9，11張圖像，如圖7所示。圖5，6采集于體視顯微鏡，可以直觀地看出立體的形態。圖5的焦面從頂部往下依次變化。圖6的焦面從底部往上依次變化。圖7采集于生物顯微鏡，景深很小，由于植物切片崎嶇不平，焦面分散在圖像的各個區域。

圖5 樣本1電路板Fig.5 PCB of sample 1

圖6 樣本2電路板Fig.6 PCB of sample 2

圖7 樣本3序列部分圖像Fig.7 Sequence of partial images of sample 3

采用方差和信息熵來度量融合質量，方差公式為

(19)

(20)

3.2 可視化比較和分析

對于樣本1與樣本2，由于只有5 張序列圖像，且焦面之間區分明顯，因此大部分區域融合質量接近，其他方法在不同區域的細節處均存在一些問題，尤其是DCT方法，在樣本1與樣本2的融合結果里均出現了大量錯誤的偽影，單獨以矩形標出，如圖8與圖9(a)所示。對于樣本3，由于生物顯微鏡下細胞切片可以被透視，很多區域沒有一個固定的焦面，即多聚焦圖像序列中某些區域沒有哪一處是非常清晰的，因此需要從多張圖像的失焦區域融合成相對清晰的結果，此類場景對于單純使用一張mask圖像來指導融合的方法(如DCT，GFF與ECNN等)是做不到的，所以針對樣本3，筆者方法可以得到更多的細節，與其他方法對比明顯，結果如圖10所示。

圖8 樣本1的融合結果Fig.8 Fusion results of sample 1

圖9 樣本2的融合結果Fig.9 Fusion results of sample 2

圖10 樣本3的融合結果Fig.10 Fusion results of sample 3

為了進一步定量分析融合后的結果，采用方差、信息熵對每組融合圖像進行客觀度量，計算結果如表1所示，由表1可以看出筆者方法均獲得了最佳度量結果。

表1 樣本1,2的電路板以及樣本3的植物切片融合結果對比

算法運行效率如表2所示，由表2可以看出：隨著圖像分辨率的增加，筆者算法運行時間雖然在增加，但是即使在1 600 W這種高分辨率情況下，運行時間也只有511 ms，而如果不采用積分圖像的方法，算法效率明顯低很多，特別是當分辨率達到1 600 W時，所需時間超過95 s，無法支持實時計算。DCT，LP，GFF和ECNN方法均采用作者提供的源碼進行實驗，編程語言為Python/Matlab，筆者方法采用C++作為實現語言，由于ECNN方法作者提供的源碼在600 W 分辨下測試便運行了數小時，因此不作運行效率參考。

表2 算法運行效率Table 2 Efficiency of algorithm

4 結 論

筆者提出了一種實時多聚焦圖像融合算法，實驗結果表明：一方面，該算法能有效抑制偽影，且在不同場景下都可以保留顯微圖像的局部細節特征；另一方面，筆者采用了積分圖像的思想，提高了計算速度，可滿足實時融合的應用場景。對于圖像融合，都需要先進行圖像配準，誤配準將會嚴重影響融合結果，所以在后續工作中將繼續研究相關配準算法以保證配準的準確性。此外，在顯微鏡的應用場景下，序列圖像之間存在著強烈的相關性，充分利用好這部分相關信息可以提高融合算法的效果，這也是未來可以研究的課題。