張 瑞，湯心溢，朱雯青

1. 中國科學(xué)院大學(xué)，北京 100049 2. 中國科學(xué)院上海技術(shù)物理研究所，上海 200083 3. 中國科學(xué)院紅外探測與成像技術(shù)重點實驗室，北京 100049

引 言

生物光學(xué)窗口是我們根據(jù)生物組織對不同波段的光穿透能力的不同而定義的有利于生物組織進行光學(xué)成像的光學(xué)窗口。 根據(jù)光在生物組織中的衰減長度與波長之間的關(guān)系，一般地，我們定義第一生物光學(xué)窗口為700～900 nm的光波段，定義第二生物光學(xué)窗口為1 000～1 400 nm的光波段，定義第三生物光學(xué)窗口為1 500～1 700 nm的光波段。 在這三個光學(xué)窗口中，光在生物組織中的衰減長度都比較大，穿透生物組織的能力也比較強。

InGaAs短波紅外探測器響應(yīng)范圍很好的覆蓋了從900～1 700 nm的波段，精準的包含了第二和第三生物光學(xué)窗口。 在生物組織成像應(yīng)用中短波紅外具有以下的成像優(yōu)勢： 生物組織光學(xué)損傷小，成像深度大，成像信噪比高，空間和時間成像分辨率高。 根據(jù)這些特點配合適當(dāng)?shù)墓鈱W(xué)系統(tǒng)，可以實現(xiàn)針對小鼠靜脈的顯微短波紅外光譜觀測。

由于顯微成像的對焦區(qū)域集中的特點，使得小鼠靜脈熒光光譜反饋在單張圖像信息較少，因此我們設(shè)計了一種新穎的融合算法。 融合算法在圖像領(lǐng)域有了很好的發(fā)展。 主要分為兩大類，一類是基于深度學(xué)習(xí)的方法，另一類也是本工作引入的傳統(tǒng)算法。 基于深度學(xué)習(xí)的方法，需要大量的數(shù)據(jù)集，這一點在文獻11中有了很好的表述，而基于傳統(tǒng)算法的融合則可以盲設(shè)計。 在傳統(tǒng)算法中GF為基于梯度濾波的融合算法，其魯棒性差，且對細節(jié)保持較弱。 NSCT[1]算法則會在平滑區(qū)域引入噪聲，MFMGRW[2]、PCA[3]以及基于小波變換的融合算法fusedDctVar[4]、fusedDctVarCv[5]除了會在平滑區(qū)域引入噪聲之外，其在對焦區(qū)域檢測也缺失魯棒性。

針對短波紅外生物觀測的特點，設(shè)計了一種基于InGaAs探測器的紅外小鼠靜脈成像實驗，設(shè)計了一種高時空分辨率的短波紅外光譜生物探測系統(tǒng)，同時針對顯微成像的特點，提出了一種多尺度梯度域引導(dǎo)濾波的多焦距融合算法，很好的實現(xiàn)了對小鼠靜脈的觀測。

1 短波紅外光譜成像

設(shè)計的InGaAs短波紅外探測系統(tǒng)樣機如圖1(a)所示，其結(jié)構(gòu)原理如圖1(b)所示。 我們設(shè)計的短波相機最大分辨率達到640×512，其積分時間最大可以達到5 000 ms，可以最大程度接收弱光成像，以減少觀測中不必要的強激光對生物組織造成損傷。 配合不同的光學(xué)系統(tǒng)，可以實現(xiàn)對生物組織的宏觀和顯微成像。

圖1 InGaAs短波紅外小鼠靜脈生物成像(a)： 相機樣機; (b): 相機原理; (c): 現(xiàn)場熒光實驗圖; (d): 觀察所用的顯微鏡頭Fig.1 Vein imaging of mice with InGaAs camera(a)： The demo of camera; (b): Design framework of camera; (c): Fluorescence experiment; (d): Lens

實驗采用了顯微光學(xué)成像系統(tǒng)如圖1(d)所示，可以對小鼠的靜脈精準成像。 我們實現(xiàn)了1 100，1 250和1 350 nm等多個波段，不同積分時間下的小鼠靜脈數(shù)據(jù)采集，實驗現(xiàn)場如圖1(c)所示。 注射熒光顯影液后，激光照射在小鼠表面。 在顯微光學(xué)鏡頭下，由于對小鼠靜脈對焦區(qū)域較小，我們設(shè)計了針對靜脈圖像特點的多焦距融合算法。

2 單光譜多焦距融合算法

小鼠靜脈交織重疊，構(gòu)成了一張靜脈網(wǎng)絡(luò)。 單張圖像獲得的靜脈成像信息少，多張圖像又缺乏信息的連貫性。 有鑒于此，我們設(shè)計了一種基于梯度域引導(dǎo)濾波的對焦像素檢測的多焦距融合算法。 通過梯度域引導(dǎo)濾波最大程度的獲取原始圖像中的高頻靜脈成分，針對交織網(wǎng)絡(luò)場景的特點設(shè)計融合規(guī)則，來獲取信息量更加豐富的靜脈圖像。

2.1 梯度域引導(dǎo)濾波金字塔對焦像素檢測

我們設(shè)計的融合算法框架如圖2所示，通過多尺度梯度域引導(dǎo)濾波提取圖像對焦像素區(qū)域，然后根據(jù)獲得的對焦區(qū)域求解融合決策函數(shù)，最后進一步改進了融合函數(shù)，可以獲得更好的融合魯棒性。 下面具體介紹算法的內(nèi)容。

圖2 融合算法總體框架流程Fig.2 Framework of algorithm

2.1.1 梯度域引導(dǎo)濾波

梯度域引導(dǎo)濾波GDGF(gradient domain guided image filter)由引導(dǎo)濾波[6]改進而來。 通過對邊緣感知加權(quán)函數(shù)的改進，可以更大程度的提取圖像的邊緣細節(jié)。 原理如下：

GDGF的損失函數(shù)如式(1)所示

(1)

式(1)中，γk為邊緣感知濾波因子，ΓI, k為邊緣感知加權(quán)濾波因子，根據(jù)引導(dǎo)圖像I計算像素k的重要性。γk與ΓI, k分別定義如式(2)和式(3)

(2)

(3)

其中，N是引導(dǎo)圖像的總像素數(shù)，ε=(0.001×L)2，L是輸入圖像的動態(tài)范圍。Sk為引導(dǎo)圖像I像素k周圍兩個窗口像素的方差的積，即窗口因子Ωr(k)的方差，Ωr(k)表示以像素k為中心半徑為r的正方形窗口，S表示所有的Si(i=1，…，N)，μS表示Si的均值。

Sk=σI, 1(k)σI, r(k)

(4)

Sk同時計算了兩個窗口函數(shù)的方差，這兩個窗口函數(shù)分別是3×3和(2r+1)×(2r+1)。

最優(yōu)化值ak和bk分別由式(5)和式(6)得到

(5)

bk=μp, r(k)-akμI, r(k)

(6)

最后輸出值為

qi=μa(i)Ii+μb(i)

(7)

其中μa(i)和μb(i)分別表示ak與bk在窗口處的均值，按如式(8)和式(9)計算

(8)

(9)

梯度域引導(dǎo)濾波相對引導(dǎo)濾波可以更好的提取圖像的細節(jié)，針對窗口大小的控制，可以獲取不同尺度的圖像頻率信息，基于此特點，我們設(shè)計了多尺度梯度域金字塔對焦像素提取算法，同時可以利用梯度域引導(dǎo)濾波對邊緣的感知特點，針對小鼠靜脈的邊緣信息可以更好的保存。

2.1.2 單光譜對焦像素檢測

在一般的光電成像系統(tǒng)中，系統(tǒng)函數(shù)是一個低通濾波器。 模糊區(qū)域的系統(tǒng)函數(shù)比清晰區(qū)域的高頻分量更低[7-8]，因此在對焦區(qū)域有更多的高頻信息。 在多焦距融合算法中首要的就是對焦像素區(qū)域檢測，我們可以通過對高頻分量的檢測實現(xiàn)對焦區(qū)像素區(qū)域的檢測[9]。

針對靜脈圖像條狀曲線紋理的特點，設(shè)計的基于GDGF金字塔對焦像素檢測[10-12]，主要分為以下兩個步驟。 首先是多尺度高頻分量提取，基于GDGF保持高頻細節(jié)的特點，通過多尺度的形式，增加算法對圖像細節(jié)提取的魯棒性，其原理如圖2黃色與藍色框所示。

其中I是輸入圖像，初始的濾波半徑是rinitial=1，eps為濾波器的參數(shù)，設(shè)定為0.05，于是得到圖像細節(jié)的分解函數(shù)如式(10)所示

(10)

式(10)中，ri作為迭代因子，其迭代方式為ri=2×r(i-1)，于是根據(jù)多尺度分解函數(shù)，得到圖像細節(jié)金字塔Di，其中(i=0, 1, 2, …，n)。 最后通過式(11)獲取n的數(shù)值

(11)

于是得到IA,IB兩張待融合圖像的對焦像素細節(jié)分量(focus pixel decision map，F(xiàn)PDM)為

(12)

2.2 融合規(guī)則

在多焦距融合中，清晰的區(qū)域一定含有最多的對焦像素，而模糊的區(qū)域則含有最少的對焦像素。 FPDMA與FPDMB分別表示待融合圖像IA,IB的對焦像素區(qū)域，初步的粗融合規(guī)則CF(coarse-fusion)定義如式(13)

(13)

由于生物組織中靜脈的交織重疊，圖像中所反映的單一焦距下靜脈成像會交錯出現(xiàn)在不同的焦距的圖像中。 這些圖像中難免會出現(xiàn)小區(qū)域孔洞，而圍成這些孔洞的靜脈圖像是不可以被拋棄的，如圖3所示。 再者，如果直接使用這個CF粗融合規(guī)則作為最后的融合規(guī)則的話，最后得到的融合圖像又會不可避免的引入融合噪聲，所以對這個CF融合規(guī)則做一個精細化處理，如式(14)—式(16)所示

圖3 靜脈紋理細節(jié)的網(wǎng)孔說明Fig.3 Illustration of texture detail

ICF(x,y)=CF(x,y)IA(x,y)+

(1-CF(x,y))IB(x,y)

(14)

FFM(x,y)=GDGF(ICF(x,y), CF(x,y))

(15)

Ioutput=FFM(x,y)IA(x,y)+

(1-FFM(x,y))IB(x,y)

(16)

其中FFM(x,y)表示最終的決策融合映射函數(shù)(final fusion map，F(xiàn)FM)。

3 實驗與結(jié)果討論

3.1 成像

通過不同波段的短波紅外激光照射，實現(xiàn)1 100，1 250和1 350 nm等不同波段的小鼠靜脈觀測，InGaAs短波探測器成像結(jié)果如圖4所示。

圖4 采集到的不同波段的靜脈成像圖Fig.4 Vein images of different wavelengths

從圖4中我們可以清楚的看到小鼠的靜脈短波紅外光譜紋理(高亮細節(jié))，四周深色的須狀為小鼠的毛發(fā)，但是由于光學(xué)對焦在不同的區(qū)域，單張圖像中無法全面的顯示小鼠靜脈信息，并且不同的波段的紅外光在小鼠機體組織中穿透力也有不同。

3.2 單光譜多焦距圖像融合

將獲得的不同譜段InGaAs短波紅外小鼠靜脈圖像部署到多個算法平臺上包括GF，NSCT，MFMGRW，PCA，fusedDctVar以及fusedDCTVarCv等6個算法，并與我們的算法結(jié)果進行比較。 同時采用一些評價指標包括Qabf，SCD，MS_SSIM以及SSIM等4個指標[13]，來綜合評價我們的算法。

我們設(shè)計的基于梯度域引導(dǎo)濾波的金字塔圖像融合算法結(jié)果如圖5—圖7所示。 在視覺效果上本算法兩焦距融合的邊界平滑度更好，且針對原始圖起到了一定的降噪作用。 此外，從1 100，1 250和1 350 nm三個波段的圖像融合結(jié)果中可以看出，其他的幾種算法包括GF，NSCT，MFMGRW，PCA，fusedDctVar以及FusedDctVarCv等在不同程度上引入額外的噪聲，同時針對對焦區(qū)域檢測上也缺乏魯棒性。 在GF算法中，1 250 nm的的融合圖像中直接丟失了A圖像的對焦區(qū)域，其他幾種算法也有不同程度的細節(jié)丟失。 這些問題在我們設(shè)計的算法上得到了很好的解決。 采用了Qabf，SCD，MS_SSIM, SSIM等四個融合圖像評價指標，其結(jié)果如表1—表3所示，我們的算法綜合評價最高，其中在MS_SSIM這個指標上我們得到了全優(yōu)。 我們設(shè)計的算法成功的將兩個光譜圖像的對焦區(qū)域融合到了一張圖中，使得我們在單張光譜圖像中獲取更多的小鼠靜脈信息。

表1 1 100 nm 5 000 ms各個融合算法評價指標得分Table 1 Evaluation index scores of each fusion algorithm at 1 100 nm 5 000 ms

表2 1 250 nm 5 000 ms各個融合算法評價指標得分Table 2 Evaluation index scores of each fusion algorithm at 1 250 nm 5 000 ms

表3 1 350 nm 5 000 ms各個融合算法評價指標得分Table 3 Evaluation index scores of each fusion algorithm at 1 350 nm 5 000 ms

4 結(jié) 論

設(shè)計了InGaAs短波紅外小鼠靜脈熒光光譜觀測系統(tǒng)，通過短波紅外光譜照射，實現(xiàn)了針對小鼠靜脈短波紅外光譜熒光觀測成像。 小鼠靜脈交織重疊，為了在單張圖像中呈現(xiàn)更多小鼠靜脈光譜的細節(jié)，提出的一種基于梯度域引導(dǎo)濾波金字塔的圖像融合算法，先通過梯度域引導(dǎo)濾波金字塔算法檢測出圖像的對焦像素區(qū)域，在基于梯度域引導(dǎo)濾波設(shè)計新的融合規(guī)則。 最終在短波紅外光譜的照射下，實現(xiàn)了針對小鼠靜脈的不同波段紅外光譜下的廣域觀測。 由于多焦距成像是因為景深的因素，即短波紅外光譜對生物組織的穿透特性，未來可以開展針對生物組織靜脈的深度短波紅外光譜成像研究工作，進一步探索短波紅外光譜對生物組織觀測的作用。

圖6 1 250 nm波長5 000 ms積分時間融合結(jié)果Fig.6 Fusion result of 5 000 ms integration time at 1 250 nm wavelength

圖7 1 350 nm 5 000 ms各個融合算法結(jié)果Fig.7 Fusion result of 5 000 ms integration time at 1 350 nm wavelength