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

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

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

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

1 算法概述

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

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

2 算法細節(jié)

2.1 圖像配準

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

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

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

2.2 聚焦度量

對于多焦圖像融合，需要先對圖像進行聚焦度量，再對聚焦區(qū)域進行分割和圖像融合。因此，選擇合適的聚焦度量方法對后續(xù)過程至關重要。筆者分析了幾種典型的度量方法：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外，其他方法的融合均出現(xiàn)了不同程度的白色偽影，并且采用SML方法能保留更多的細節(jié)。因此，筆者采用SML作為聚焦度量方法，并使用積分圖像進行優(yōu)化。

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

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

(10)

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

(11)

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

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

(12)

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

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

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

(13)

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

(14)

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

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

(15)

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

(16)

(17)

(18)

3 實驗結果與分析

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

3.1 實驗數(shù)據和誤差度量方法

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

圖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 張序列圖像，且焦面之間區(qū)分明顯，因此大部分區(qū)域融合質量接近，其他方法在不同區(qū)域的細節(jié)處均存在一些問題，尤其是DCT方法，在樣本1與樣本2的融合結果里均出現(xiàn)了大量錯誤的偽影，單獨以矩形標出，如圖8與圖9(a)所示。對于樣本3，由于生物顯微鏡下細胞切片可以被透視，很多區(qū)域沒有一個固定的焦面，即多聚焦圖像序列中某些區(qū)域沒有哪一處是非常清晰的，因此需要從多張圖像的失焦區(qū)域融合成相對清晰的結果，此類場景對于單純使用一張mask圖像來指導融合的方法(如DCT，GFF與ECNN等)是做不到的，所以針對樣本3，筆者方法可以得到更多的細節(jié)，與其他方法對比明顯，結果如圖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++作為實現(xiàn)語言，由于ECNN方法作者提供的源碼在600 W 分辨下測試便運行了數(shù)小時，因此不作運行效率參考。

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

4 結 論

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