齊海超，宋延嵩,*，張 博，梁宗林，閆綱琦，薛佳音，張軼群，任 斌

（1.長(zhǎng)春理工大學(xué) 光電工程學(xué)院 空間光電技術(shù)研究所, 長(zhǎng)春 130022；2.鵬城實(shí)驗(yàn)室, 深圳 518052；3.中國(guó)空間技術(shù)研究院西安分院, 西安 710212）

1 引 言

多光譜圖像比RGB 圖像擁有更完整的空間和光譜信息，目前已廣泛應(yīng)用于醫(yī)學(xué)成像[1-3]、食品檢測(cè)[4]和人臉識(shí)別[5-6]等領(lǐng)域。多光譜圖像由多光譜濾光片陣列(Multispectral Filter Array,MSFA)成像獲得。該技術(shù)通過(guò)在圖像傳感器前端放置具有馬賽克模式的濾光片陣列[7]，經(jīng)一次曝光就可得到物體的全部光譜信息。由濾光片陣列直接得到的圖像稱(chēng)為Raw 圖像，在每個(gè)采樣點(diǎn)處估計(jì)其余波段丟失信息的過(guò)程稱(chēng)為多光譜去馬賽克過(guò)程[8]。與傳統(tǒng)彩色濾光片陣列(Color Filter Array, CFA)相比，由于MSFA 增加了波段數(shù)，每個(gè)波段Raw 圖像的空間采樣率更低，在去馬賽克過(guò)程中會(huì)丟失大量細(xì)節(jié)信息，且傳統(tǒng)CFA 去馬賽克方法無(wú)法直接應(yīng)用到MSFA 去馬賽克問(wèn)題中，因此多光譜去馬賽克問(wèn)題是近年來(lái)的一個(gè)難點(diǎn)。

與深度學(xué)習(xí)方法[9-11]不同，插值類(lèi)方法不需要大量的多光譜圖像來(lái)訓(xùn)練和驗(yàn)證模型，因此更適合應(yīng)用于某些實(shí)時(shí)場(chǎng)景中。目前，大多數(shù)去馬賽克工作中的濾光片陣列都是基于Miao[12]等人提出的二叉樹(shù)通用MSFA 模式而設(shè)計(jì)的。對(duì)于給定的波段數(shù)及其出現(xiàn)概率(Probability of Appearance, POA)，該方法通過(guò)對(duì)棋盤(pán)格進(jìn)行二元分解來(lái)得到所需陣列。針對(duì)該通用模式，Miao 等人利用邊緣信息[13]逐步估計(jì)每個(gè)波段的缺失像素值，當(dāng)陣列波段數(shù)較多時(shí)，其重建圖像會(huì)出現(xiàn)嚴(yán)重的模糊。Gupta[14]等人為不同波段數(shù)的MSFA 分別設(shè)計(jì)了特定的低通濾波器，并通過(guò)卷積操作來(lái)重建稀疏圖像，該方法速度快，但重建圖像質(zhì)量較差。隨后，Gupta 等人又提出了一種自適應(yīng)漸進(jìn)多光譜去馬賽克方法(Adaptive and Progressive Multispectral Image Demosaicking, APMID)[15]。該方法首先使用自適應(yīng)漸進(jìn)插值估計(jì)中間波段(POA 最高)，然后利用插值后的中間波段，采用漸進(jìn)譜差法對(duì)其它波段進(jìn)行插值，完成去馬賽克過(guò)程。Sun[16]等人設(shè)計(jì)了一種具有等空間概率比的八譜段MSFA 成像系統(tǒng)，作者利用鄰域采樣點(diǎn)的像素值和梯度值，結(jié)合泰勒多項(xiàng)式估計(jì)未采樣點(diǎn)并完成圖像的重建。RATHI[17]等人使用改進(jìn)的方向加權(quán)插值來(lái)估計(jì)缺失的像素值，該方法是通用的，作者在五到十六波段的MSFA 上驗(yàn)證了其性能。Monno[18]等人受CFA 的啟發(fā)，設(shè)計(jì)了一種特定模式的五波段MSFA，該陣列中G 波段的POA 為1/2。作者基于引導(dǎo)濾波方法(Guided Filter, GF)，首先使用自適應(yīng)高斯上采樣插值G 波段原始圖像，并將其作為引導(dǎo)圖像插值其余波段，該方法的重建圖像精度嚴(yán)重依賴(lài)于引導(dǎo)圖的質(zhì)量，且固定的懲罰因子對(duì)圖像的不同區(qū)域不能自適應(yīng)調(diào)整。

針對(duì)Monno 等人方法的不足，本文提出了一種基于改進(jìn)引導(dǎo)濾波器的多光譜圖像去馬賽克方法。首先使用二維信號(hào)的自回歸模型對(duì)G 波段進(jìn)行插值，提高引導(dǎo)圖的精度；然后使用具有邊緣感知權(quán)重的引導(dǎo)濾波器重建其余稀疏帶，使本方法能夠自適應(yīng)處理圖像的邊緣和紋理區(qū)域，且重建圖像的細(xì)節(jié)信息更接近于真實(shí)圖像。

2 MSFA 模式設(shè)計(jì)與方法基本原理

2.1 基于二叉樹(shù)的5 波段MSFA 模式設(shè)計(jì)

MSFA 模式的選擇對(duì)后續(xù)去馬賽克方法重建效果至關(guān)重要，基于二叉樹(shù)法的通用MSFA 模式設(shè)計(jì)應(yīng)遵循以下準(zhǔn)則：(1)應(yīng)滿(mǎn)足光譜一致性要求，保證每個(gè)譜段采樣點(diǎn)鄰域內(nèi)的像素始終具有相同的排列方式，以避免光學(xué)串?dāng)_效應(yīng)的影響；(2)應(yīng)滿(mǎn)足空間均勻性要求，每個(gè)光譜帶的采樣點(diǎn)在整幅圖像內(nèi)應(yīng)盡可能均勻分布；(3)應(yīng)具有周期性以確保圖像重建的效率。

基于以上準(zhǔn)則，本文所設(shè)計(jì)的5 波段MSFA包括一個(gè)密集采樣波段G 和4 個(gè)稀疏采樣波段B、C、R、Y，其二叉樹(shù)分裂過(guò)程如圖1(a)（彩圖見(jiàn)期刊電子版）所示。二叉樹(shù)根節(jié)點(diǎn)的POA 為1，包含所有譜段的信息，從根節(jié)點(diǎn)出發(fā)開(kāi)始二元分裂，在每次分裂過(guò)程中，父通道平均分為兩個(gè)子通道，子通道的采樣率為父通道的一半，即位于第n級(jí)節(jié)點(diǎn)上的譜段空間采樣率為 2-n。在最后得到的MSFA 模式中，G 波段在第一級(jí)節(jié)點(diǎn)上，采樣率為1/2，其余4 個(gè)波段在第三級(jí)節(jié)點(diǎn)上，采樣率為1/8。MSFA 排列模式如圖1(b)（彩圖見(jiàn)期刊電子版）所示。

圖1 5 波段MSFA(a)二叉樹(shù)分裂過(guò)程及(b)排列模式Fig.1 (a) Binary tree splitting process and (b) arrangement of Five-band MSFA

2.2 G 波段去馬賽克

多光譜圖像鄰域內(nèi)像素點(diǎn)間通常存在著較大的相關(guān)性，自回歸模型是一種對(duì)該相關(guān)性建模的通用模型之一[19]。二維圖像信號(hào)Ym,n可以利用自回歸模型建模為：

其中，T是以像素點(diǎn) (m,n)為中心的空間鄰域窗口，Xm+p,n+q是T內(nèi)異于中心的任意像素點(diǎn)， αp,q為其自回歸模型參數(shù)， γm,n為誤差項(xiàng)?？紤]到存在變量不足、測(cè)量噪聲等不確定因素，該誤差會(huì)導(dǎo)致構(gòu)建的自回歸模型具有不完全準(zhǔn)確性，使Ym,n和真實(shí)值之間產(chǎn)生一定的差異。本方法引入軟決策自適應(yīng)插值(Soft decision Adaptive Interpolation,SAI)[20]思想，完成G 波段馬賽克圖像的重建，下面詳細(xì)闡述具體流程。

首先，對(duì)于每個(gè)待估計(jì)像素，以其為中心選定一個(gè)大小為 7×7 的 局部窗口T，為便于表達(dá)，對(duì)T內(nèi)像素點(diǎn)進(jìn)行編號(hào)，令xi表示已知像素點(diǎn)，yj表示未知像素點(diǎn)，已知像素點(diǎn)集合為X={x1,.x2,···,x24} ， 未知像素點(diǎn)集合為Y={y1,y2,···,y25}。如圖2 所示，窗口內(nèi)實(shí)心圓點(diǎn)表示已知像素，空心圓點(diǎn)表示未知像素。

圖2 鄰域T 內(nèi)像素排列Fig.2 Pixel arrangement in neighborhood T

然后，利用自回歸模型對(duì)窗口T內(nèi)像素點(diǎn)間的相關(guān)性在水平-垂直方向上進(jìn)行建模，如圖3 所

圖3 水平-垂直方向上的自回歸模型Fig.3 Autoregressive model in the horizontal-vertical direction

式中 γj， γi分別為yj和xi的估計(jì)誤差，在假設(shè)自回歸模型參數(shù)bj,t，bi,t給定的前提下，通過(guò)最小化窗口T內(nèi)像素點(diǎn)的估計(jì)誤差總和，可以得到未知像素點(diǎn)yj的最優(yōu)估計(jì)值，該過(guò)程等價(jià)于最小二乘問(wèn)題：

求解上述最小二乘問(wèn)題，得到局部窗口T內(nèi)一組待估計(jì)像素的最優(yōu)解后，將局部窗口在圖像空間內(nèi)進(jìn)行遍歷，完成整幅圖像的插值過(guò)程。

考慮到SAI 在每個(gè)局部窗口內(nèi)采用固定的模型參數(shù)，當(dāng)邊緣或紋理尺寸小于窗口時(shí)，估計(jì)的模型參數(shù)會(huì)產(chǎn)生較大誤差。因此，本文采用漸進(jìn)方法獨(dú)立估計(jì)每一個(gè)像素點(diǎn)的模型參數(shù)b，使自回歸模型對(duì)圖像的紋理和細(xì)節(jié)具有更好的適應(yīng)能力。基于像素之間的空間連續(xù)性，首先使用已知像素點(diǎn)強(qiáng)度之間的高斯梯度函數(shù)計(jì)算待估計(jì)像素參數(shù)bj,t：

式中， σj是待估計(jì)像素yj四鄰域內(nèi)已知像素點(diǎn)間的標(biāo)準(zhǔn)差。如圖4(a)所示，xj,t和xj,s為yj在水平-垂直方向上最近鄰的已知采樣像素點(diǎn)，計(jì)算出bj,t后，利用式(2)可以得到未知像素的初步估計(jì)值：

圖4 模型參數(shù)估計(jì)Fig.4 Estimation of model parameter

如圖4(b)所示，對(duì)于已知像素xi的模型參數(shù)bi,t，可利用其與四鄰域內(nèi)未知像素的初步估計(jì)值間的相似性得到：

其中，y和x分別是窗口T內(nèi)的未知像素點(diǎn)和已知像素點(diǎn)集合，I13×13為13 階單位矩陣，其余子矩陣定義如下：

則公式(8)的解析解可以表示為：

本方法輸出式(13)中與窗口T內(nèi)中心像素初步估計(jì)值差異最小的y值作為最后的插值結(jié)果。

2.3 其余四波段去馬賽克

與密集波段G 相比，其余4 個(gè)稀疏帶的采樣率很低，僅利用空間相關(guān)性進(jìn)行重建無(wú)法得到滿(mǎn)意的結(jié)果。基于多光譜圖像不同波段間邊緣信息高度相似這一假設(shè)，利用改進(jìn)引導(dǎo)濾波器來(lái)對(duì)稀疏帶圖像去馬賽克。以稀疏帶R 為例，圖5(彩圖見(jiàn)期刊電子版)說(shuō)明了其去馬賽克流程。引導(dǎo)濾波器[21]由局部線(xiàn)性模型導(dǎo)出，其通過(guò)考慮引導(dǎo)圖像的內(nèi)容來(lái)計(jì)算濾波器的輸出，能夠?qū)⒁龑?dǎo)圖像的結(jié)構(gòu)信息轉(zhuǎn)移到輸出圖像中，因此具有良好的邊緣保持特性。

圖5 基于加權(quán)引導(dǎo)濾波的去馬賽克流程Fig.5 Demosaicing process based on weight-guided filtering

傳統(tǒng)引導(dǎo)濾波器對(duì)整幅圖像使用統(tǒng)一的線(xiàn)性模型和固定的正則化參數(shù)，無(wú)法適應(yīng)圖像自身不同區(qū)域的紋理特性[22]，因此將加窗固有變分系數(shù)(Windowed Inherent Variation, WIV)[23]作為邊緣感知權(quán)重引入到引導(dǎo)濾波器中，WIV 系數(shù)能夠突出圖像邊緣和紋理上的像素，捕捉窗口內(nèi)的空間變化。對(duì)于一幅輸入圖像I，Lx(p)和Ly(p)分別為其上像素點(diǎn)p在x和y方向上的WIV 系數(shù)，表達(dá)式為：

式中，Rr(p) 是以像素點(diǎn)p為中心，以r為半徑的局部鄰域窗口，q是窗口Rr(p) 內(nèi) 異于p的任意像素點(diǎn)， ?x和 ?y分別表示x和y方 向上的偏導(dǎo)數(shù)，gp,q是根據(jù)空間相似性定義的高斯權(quán)重函數(shù)：

其中 (xp,yp) 和 (xq,yq) 是窗口內(nèi)像素點(diǎn)p和q的空間坐標(biāo)， σ是控制窗口的空間尺度。最終p點(diǎn)的WIV系數(shù)定義為：

邊緣處像素點(diǎn)的L值通常大于紋理和平滑區(qū)域內(nèi)像素點(diǎn)的L值，這是因?yàn)榧y理區(qū)域內(nèi)的?I值因方向不同相互抵消，而平滑區(qū)域內(nèi) ?I值的模長(zhǎng)很小，所以WIV 系數(shù)能夠有效地檢測(cè)出邊緣信息。

本文將重建后的G 波段作為引導(dǎo)圖像P，并使用加權(quán)引導(dǎo)濾波來(lái)重建其余波段的圖像。引導(dǎo)濾波的關(guān)鍵假設(shè)是引導(dǎo)圖像P和 輸出q在一個(gè)二維局部窗口內(nèi)滿(mǎn)足線(xiàn)性關(guān)系，如下：

式中， ωk是以像素k為中心的窗口，i是窗口內(nèi)的任意像素點(diǎn)。 (ak,bk) 是 ωk內(nèi)一組假定為常數(shù)的線(xiàn)性系數(shù)。對(duì)等式兩邊同時(shí)求梯度可以得到?q=a?P， 說(shuō)明輸出q和 引導(dǎo)圖P有相似的梯度，即引導(dǎo)濾波器具有邊緣保持特性。

為了確定線(xiàn)性系數(shù) (ak,bk)，需要確保濾波器輸出q和原始圖像I之間的差異最小，通過(guò)最小化下述代價(jià)函數(shù)來(lái)實(shí)現(xiàn)該目的：

式中，Mi是二元掩模，在原始圖像I內(nèi)已采樣點(diǎn)處的值為1，其余位置值為0，第二項(xiàng)是正則化項(xiàng)，用來(lái)對(duì)系數(shù)ak進(jìn)行懲罰， ε 是 懲罰因子， ψ(k)是基于WIV 系數(shù)的邊緣感知權(quán)重，該權(quán)重能夠自適應(yīng)地調(diào)整懲罰因子 ε，對(duì)于紋理或邊緣區(qū)域，系數(shù)ak的值較大[24]，需要較小的懲罰系數(shù)；對(duì)于平滑區(qū)域，系數(shù)ak值較小，需要較大的懲罰系數(shù)。 ψ(k)定義為：

其中，k為中心像素，i取遍圖像的所有像素，N為像素總數(shù), λ 是 取值為 (0.001×D)2的常數(shù)，D為圖像I的動(dòng)態(tài)范圍(8 位圖像的動(dòng)態(tài)范圍為0～255)。式(19)的解為：

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

為了評(píng)估本方法的性能，本文使用以下兩個(gè)公開(kāi)的、著名的多光譜圖像數(shù)據(jù)集進(jìn)行驗(yàn)證：

(1) CAVE[25]數(shù)據(jù)集：CAVE 數(shù)據(jù)集由哥倫比亞大學(xué)提出，該數(shù)據(jù)集由32 類(lèi)場(chǎng)景圖片組成，每種場(chǎng)景共包含31 個(gè)波段，波長(zhǎng)范圍在400 ～700 nm 之間，光譜分辨率為10 nm，空間分辨率為512×512。該數(shù)據(jù)集是由VariSpec 液晶可調(diào)濾波器和冷CCD 拍攝的。

(2) TokyoTech[26]數(shù)據(jù)集：該數(shù)據(jù)集包含30種場(chǎng)景圖像，這些圖像具有豐富的色彩和紋理。每種多光譜圖像包含31 個(gè)波段，波長(zhǎng)范圍為420～720 nm，光譜間隔為10 nm，空間分辨率為500×500。

首先在每種數(shù)據(jù)集中選擇10 種不同場(chǎng)景，然后通過(guò)31 波段原始圖像模擬得到5 波段真實(shí)圖像，并根據(jù)指定的MSFA 對(duì)真實(shí)圖像進(jìn)行下采樣，生成Raw 圖像，最后采用GF 方法、APMID 方法作為對(duì)比方法來(lái)重建馬賽克圖像。實(shí)驗(yàn)使用線(xiàn)性模型估計(jì)光譜反射率[27]，并基于文獻(xiàn)[26]中的方法將5 波段重建圖像變換到sRGB 顏色空間中。為避免卷積操作跨越圖像邊界帶來(lái)的影響，本方法從原始圖像和去馬賽克圖像中排除了5 個(gè)邊界像素。

3.1 去馬賽克圖像質(zhì)量客觀評(píng)價(jià)

表1 和表2 分別為3 種方法在CAVE 數(shù)據(jù)集和TokyoTech 數(shù)據(jù)集上的重建圖像精度對(duì)比。在表1 中，本方法的sRGB PSNR、sRGB SSIM 分別較GF 方法、APMID 方法提高了1.36 dB、0.82 dB，0.0129、0.0029；在表2 中，本方法的sRGB PSNR、sRGB SSIM 比 GF 方法、APMID 方法高2.25 dB、1.27 dB，0.0218、0.0137；CIEDE 2000 衡量不同圖像之間的色差，數(shù)據(jù)值越小，色差越小，本方法重建圖像的色差均小于另外兩種方法?？梢钥闯霰痉椒ㄔ趦煞N數(shù)據(jù)集上的去馬賽克圖像質(zhì)量大多優(yōu)于另外兩種方法，這是因?yàn)椋?/p>

表1 CAVE 數(shù)據(jù)集上3 種方法的客觀評(píng)價(jià)指標(biāo)Tab.1 Objective evaluation metrics of the three methods on the CAVE dataset

表2 TokyoTech 數(shù)據(jù)集上3 種方法的客觀評(píng)價(jià)指標(biāo)Tab.2 Objective evaluation metrics of the three methods on the TokyoTech dataset

(1) 本方法提出的自回歸模型插值基于漸進(jìn)的思想逐點(diǎn)估計(jì)模型參數(shù)，根據(jù)圖像梯度信息自適應(yīng)計(jì)算插值權(quán)重。并且引入了誤差反饋機(jī)制，該機(jī)制通過(guò)最小化未知像素估計(jì)值和已知像素值之間的誤差使得估計(jì)值更接近于真實(shí)值。

(2) 得益于引入的邊緣感知權(quán)重WIV，改進(jìn)的引導(dǎo)濾波器能夠根據(jù)圖像不同結(jié)構(gòu)自適應(yīng)改變懲罰系數(shù)ε，更好地保留了圖像的高頻分量部分。

表3 為3 種方法分別在CAVE 數(shù)據(jù)集和TokyoTech 數(shù)據(jù)集上的平均運(yùn)行時(shí)間，用時(shí)從長(zhǎng)到短排序依次為：GF 方法>本文方法>APMID 方法。APMID 方法運(yùn)行速度最快，這是因?yàn)槠浜?jiǎn)單地使用了鄰域加權(quán)插值，而本方法引入誤差反饋機(jī)制，通過(guò)犧牲一定的實(shí)時(shí)性換取了重建圖像更高的峰值信噪比、結(jié)構(gòu)相似度以及更低的色差。

表3 不同方法在兩種數(shù)據(jù)集上的運(yùn)行時(shí)間Tab.3 Running times of different methods on the two datasets (s)

3.2 去馬賽克圖像視覺(jué)效果評(píng)價(jià)

分別選取CAVE 數(shù)據(jù)集的Balloons 場(chǎng)景和TokyoTech 數(shù)據(jù)集的CD、Party 場(chǎng)景作為原始圖像，對(duì)3 種方法的去馬賽克圖像進(jìn)行視覺(jué)效果評(píng)價(jià)。如圖6、7、8（彩圖見(jiàn)期刊電子版）所示?？梢钥闯觯篏F 方法的重建圖像中，邊緣附近存在著明顯的顏色偽影，且平滑區(qū)域內(nèi)出現(xiàn)了虛假紋理和顯著的噪聲，造成了圖像的失真，這是因?yàn)镚F 方法的懲罰系數(shù)固定，不能適應(yīng)邊緣信息；而APMID 方法由于僅使用鄰域內(nèi)像素均值來(lái)插值未知像素，因此當(dāng)未知像素點(diǎn)在邊緣附近時(shí)，插值操作會(huì)引入跨越邊緣的法向誤差，導(dǎo)致APMID 方法的重建圖像中邊緣附近有明顯的噪聲出現(xiàn)；相比于以上兩種方法，本方法的重建圖像中僅有少量的邊緣偽影，且圖像不同區(qū)域之間過(guò)渡更加平滑，更接近于真實(shí)圖像，說(shuō)明本方法能夠提高重建圖像的精度，對(duì)邊緣信息有更好的適應(yīng)能力，綜合性能優(yōu)于GF 方法和APMID 方法。

圖6 不同算法Balloons 場(chǎng)景重建圖像對(duì)比Fig.6 Comparison of Balloons images reconstructed by different algorithms

圖7 不同算法CD 場(chǎng)景重建圖像對(duì)比Fig.7 Comparison of CD images reconstructed by different algorithms

圖8 不同算法Party 場(chǎng)景重建圖像對(duì)比Fig.8 Comparison of party images reconstructed by different algorithms

4 結(jié) 論

本文提出了一種基于改進(jìn)引導(dǎo)濾波器的多光譜圖像去馬賽克方法。本方法首先采用逐像素點(diǎn)漸進(jìn)估計(jì)模型參數(shù)的自回歸模型插值密集采樣波段G，提高了G 波段Raw 圖像的重建質(zhì)量，并間接提高了引導(dǎo)圖像的精度；然后基于邊緣感知權(quán)重WIV，提出了一種邊緣自適應(yīng)的加權(quán)引導(dǎo)濾波器來(lái)插值其余波段。根據(jù)CAVE 數(shù)據(jù)集和TokyoTech 數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果表明，本方法的綜合性能優(yōu)于GF 方法和APMID 方法，其中重建圖像的峰值信噪比和結(jié)構(gòu)相似度分別提高了3.40%， 2.02%， 1.34%， 0.30%和 6.11%， 5.95%，2.28%，1.42%，且本方法有效地抑制了邊緣偽影現(xiàn)象，減少了噪聲的出現(xiàn)，重建圖像的色差更小，視覺(jué)效果更接近于真實(shí)圖像。