大面陣航攝影像重建高效像素插值算法研究

胡海彥，余旭初，楊韞瀾，方 勇，江振治

(1.信息工程大學(xué),河南 鄭州 450001；2.西安測(cè)繪研究所，陜西 西安 710054)

復(fù)合型大面陣或超大面陣航攝相機(jī)提供給作業(yè)用戶的最終影像，其每一個(gè)像素都是從原始子影像重采樣而來，就計(jì)算耗時(shí)而言，根據(jù)經(jīng)驗(yàn)圖像重采樣所耗費(fèi)的時(shí)間至少要占到整個(gè)大面陣航攝影像預(yù)處理的1/3以上，所以對(duì)此環(huán)節(jié)進(jìn)行耗時(shí)優(yōu)化，必然會(huì)極大提升大幅面影像生成效率[1-3]。圖像內(nèi)插重建是基于一組等間距離散采樣點(diǎn)恢復(fù)連續(xù)空間圖像函數(shù)的過程，其是圖像縮放、子像素平移(非整像素)、圖像旋轉(zhuǎn)、圖像形變等基本圖像操作的技術(shù)基礎(chǔ)。當(dāng)這些圖像操作所在位置缺少既有圖像灰度值，即需要圖像插值進(jìn)行灰度填補(bǔ)。而且對(duì)于大幅面數(shù)字航攝相機(jī)而言，在影像復(fù)合過程中可通過高精度像素重采樣一并解決鏡頭畸變引起的圖像幾何失真、色彩配準(zhǔn)等圖像重建問題。由采樣理論可知，sinc函數(shù)是最為理想的圖像插值內(nèi)核函數(shù)(見圖1)。然而，實(shí)際應(yīng)用中并不能直接使用該函數(shù)進(jìn)行插值[4]，需要有限截?cái)嗖?shí)現(xiàn)離散化計(jì)算，同時(shí)要在計(jì)算速度與數(shù)學(xué)精度之間取得平衡，即在盡可能接近sinc函數(shù)理論精度的前提下，設(shè)計(jì)出一個(gè)內(nèi)核函數(shù)對(duì)sinc進(jìn)行函數(shù)逼近。通?？梢酝ㄟ^n階對(duì)稱分段多項(xiàng)式進(jìn)行內(nèi)核函數(shù)逼近，如常用的雙線性卷積核函數(shù)等即為其中一種特例，文獻(xiàn)[5]對(duì)此類內(nèi)核函數(shù)進(jìn)行了定量評(píng)估，相關(guān)研究也比較豐富[6～9]，但往往給出的插值算法精度接近像素級(jí)，也沒有挖掘此類算法的潛在效率，僅適用于圖像融合等某些常規(guī)應(yīng)用。本文從一維信號(hào)重構(gòu)原理著手，推導(dǎo)出n階對(duì)稱分段多項(xiàng)式內(nèi)核函數(shù)的數(shù)學(xué)形式，并給出多項(xiàng)式定義域變換后的加速算法實(shí)現(xiàn)，最后對(duì)所提算法實(shí)作性能進(jìn)行試驗(yàn)分析。算法在盡可能不引入額外像素插值誤差，從而避免損失或破壞影像自身固有幾何量測(cè)精度前提下，大幅度降低數(shù)據(jù)處理計(jì)算復(fù)雜度，滿足大幅面數(shù)字航攝影像重建對(duì)效率的要求。

圖1 空間域中的sinc函數(shù)

1 像素插值算法

1.1 內(nèi)核函數(shù)的多項(xiàng)式逼近

對(duì)于一維空間，近似sinc內(nèi)核函數(shù)的n階對(duì)稱分段多項(xiàng)式在區(qū)間[-(n-1),(n-1)]上的定義為

(1)

其中：i=0,1,…,m-1，m表示整數(shù)截?cái)喾秶?，滿足n=2m-1。這樣總共會(huì)有(n+1)m個(gè)aij系數(shù)。由于內(nèi)核函數(shù)具有對(duì)稱性及維度的可分離性，即二維內(nèi)核函數(shù)為兩個(gè)一維內(nèi)核函數(shù)的乘積，即h(x,y)=h(x)h(y)，所以為了描述的簡(jiǎn)潔性，本文均在一維情況下予以討論。

為了保證節(jié)點(diǎn)處函數(shù)擬合的連續(xù)性與光滑性，除保證多項(xiàng)式與sinc函數(shù)在節(jié)點(diǎn)上的取值相同外，還需要對(duì)式(1)附加導(dǎo)數(shù)連續(xù)性約束條件。以下將3階對(duì)稱分段多項(xiàng)式內(nèi)核函數(shù)作為樣例，給出推導(dǎo)過程，更高階多項(xiàng)式類同。在區(qū)間[-2,2]上sinc函數(shù)的近似3階對(duì)稱多項(xiàng)式內(nèi)核函數(shù)為

hC(x)=

(2)

由約束條件(2)總共可列出7個(gè)方程，其中含8個(gè)未知數(shù)系數(shù)。進(jìn)一步令a13=α作為調(diào)整參數(shù)，則未知數(shù)個(gè)數(shù)減1，使得方程組可解。有

a03=(α+2)，a02=-(α+3)，a01=0，a00=1

a13=α,a12=-5α，a11=8α，a10=-4α.

表1 平坦性約束下的3～9階多項(xiàng)式α參數(shù)

可以看出不同的α值，對(duì)應(yīng)不同的插值計(jì)算，常用的3次卷積對(duì)應(yīng)α=-1/2，也稱之為基斯(Keys)插值內(nèi)核函數(shù)[10](見圖2)。

hkeys=

(3)

圖2 基斯3次卷積內(nèi)核函數(shù)

1.2 多項(xiàng)式定義域轉(zhuǎn)換

內(nèi)插權(quán)重計(jì)算即是利用分段多項(xiàng)式求得采樣點(diǎn)集上對(duì)應(yīng)的多項(xiàng)式函數(shù)值，以此求得卷積核模板，與采樣點(diǎn)進(jìn)行卷積運(yùn)算便完成插值點(diǎn)函數(shù)插值。所以需要先行計(jì)算插值點(diǎn)到采樣點(diǎn)的相應(yīng)距離ξ。在3次卷積核的情況下，這些距離值分別為ξ，1-ξ，1+ξ及2-ξ，其中ξ∈[0,1]?？紤]計(jì)算效率與便捷性，可將所有多項(xiàng)式的定義域區(qū)間變換為區(qū)間[0,1]，利用多項(xiàng)式的分段對(duì)稱性(圖2(a)已給出圖示)，可對(duì)式(1)中的多項(xiàng)式p0(x)在ξ=x,ξ=1-x及p1(x)在ξ=x+1,ξ=2-x上進(jìn)行分段處理。對(duì)于更高的n階多項(xiàng)式內(nèi)核函數(shù)，實(shí)現(xiàn)的策略是相同的，這樣可以得到n+1個(gè)相同定義域的分段多項(xiàng)式，對(duì)于偶數(shù)i,即(0,2,4,6…)。令k=i/2，即(0,1,2,3…)。則有

(4)

其中，

(5)

對(duì)圖2進(jìn)行3次多項(xiàng)式定義域變換后對(duì)應(yīng)的函數(shù)圖形如圖3所示。

圖3 定義域變換為[0,1]的Keys 3次卷積核

為了節(jié)約計(jì)算時(shí)間，給出關(guān)于變換后多項(xiàng)式系數(shù)屬性的5條規(guī)則。圖4為圖形化說明(以5階插值內(nèi)核函數(shù)系數(shù)表為例，α=3/64；表2為3階插值內(nèi)核函數(shù)系數(shù)，α=-1/2)。

圖4 5階插值多項(xiàng)式內(nèi)核函數(shù)系數(shù)圖示

規(guī)則1：僅第一個(gè)多項(xiàng)式f0(x)的最低階項(xiàng)系數(shù)為1，其余多項(xiàng)式fi(x)最低階項(xiàng)系數(shù)ai0均為0；

規(guī)則2：多項(xiàng)式fi(x)和fi+1(x)，最高階項(xiàng)系數(shù)ain成對(duì)相等、符號(hào)相反；

規(guī)則3：所有j≤-2階項(xiàng)系數(shù)成對(duì)出現(xiàn)，且奇數(shù)階項(xiàng)大小相等符號(hào)相反、偶數(shù)階項(xiàng)大小相等且符號(hào)相同；

規(guī)則4：多項(xiàng)式fn(x)的所有j≤n-2階數(shù)項(xiàng)的系數(shù)anj為0；

規(guī)則5：多項(xiàng)式f0(x)的所有j≤n-2(j為奇數(shù))階數(shù)項(xiàng)的系數(shù)a0j為0。

利用這5個(gè)規(guī)則，可以明顯提升計(jì)算效率。首先對(duì)于所有多項(xiàng)式,距離x具有相同的值，因此可以對(duì)x冪進(jìn)行預(yù)先計(jì)算。而且成對(duì)相同的系數(shù)及與x冪相乘的積僅需計(jì)算1次，同時(shí)可跳過具有零系數(shù)的計(jì)算項(xiàng)。

表2 多項(xiàng)式變換后的3次卷積核系數(shù)

所以，式(3)給出的Keys 3次卷積核，根據(jù)上述理論進(jìn)行多項(xiàng)式轉(zhuǎn)換，可等價(jià)得到[0,1]空間上的4個(gè)多項(xiàng)式為

(6)

2 算法性能分析及試驗(yàn)結(jié)果

2.1 性能分析

表3給出一維3階多項(xiàng)式卷積核插值權(quán)重的計(jì)算次數(shù)。對(duì)于一維3階卷積，這里提出的方法總共只需16個(gè)計(jì)算操作，若替代傳統(tǒng)方法的26個(gè)計(jì)算操作，可使計(jì)算加速1.625。表4為n階多項(xiàng)式卷積核運(yùn)算次數(shù)，隨著多項(xiàng)式階數(shù)的增加，加速倍率逐漸降低(例如9階多項(xiàng)式的加速倍率為1.52)，對(duì)于無窮階多項(xiàng)式，加速倍率穩(wěn)定為4/3。由于分段n階卷積是可分離的(卷積級(jí)聯(lián))，所以高維N-D卷積的計(jì)算耗時(shí)是N×1D卷積時(shí)間。

當(dāng)然，要對(duì)整個(gè)插值操作進(jìn)行評(píng)估時(shí)，除需要考慮插值多項(xiàng)式權(quán)重計(jì)算外，還有采樣點(diǎn)與這些權(quán)重的卷積計(jì)算。該卷積計(jì)算的時(shí)間復(fù)雜度由(n+1)N個(gè)加法、乘法和計(jì)算機(jī)內(nèi)存訪問耗時(shí)共同決定。這也意味著當(dāng)進(jìn)行高維插值(N>2)時(shí)，卷積計(jì)算量占主導(dǎo)地位。但通過將多項(xiàng)式定義域變換為[0，1]區(qū)間后，對(duì)任意維度下的插值都具加速效果，只是此方法對(duì)于一、二維低維度插值而言更為高效。

表3 3階多項(xiàng)式一維卷積核運(yùn)算次數(shù)

表4 n階多項(xiàng)式卷積核最大運(yùn)算次數(shù)估計(jì)

2.2 試驗(yàn)結(jié)果

用C++程序語言對(duì)上述高階多項(xiàng)式內(nèi)核卷積進(jìn)行算法實(shí)現(xiàn)。在配置為Intel Core i3-370M(2.4 GHz)的計(jì)算機(jī)上完成試驗(yàn)測(cè)試。為了進(jìn)一步加速高階多項(xiàng)式的計(jì)算效率，編程中同時(shí)采用查找表(LUT)數(shù)據(jù)結(jié)構(gòu)對(duì)內(nèi)核函數(shù)進(jìn)行預(yù)先計(jì)算，為保證插值精度，查找表中將單位長度設(shè)置為10 000個(gè)網(wǎng)格點(diǎn)(即插值像素點(diǎn)位可精確到0.000 1)。在3階多項(xiàng)式情況下，由于內(nèi)存訪問相對(duì)多項(xiàng)式直接計(jì)算耗時(shí)，所以采用查找表比直接計(jì)算反而略顯低效——構(gòu)建查找表耗時(shí)18 982 ms，整個(gè)插值過程耗時(shí)51 193 ms，但更高階內(nèi)核插值試驗(yàn)表明，采用查找表仍是必要的，所以在編碼實(shí)現(xiàn)時(shí)建議分3階以下/以上多項(xiàng)式兩種情況進(jìn)行程序設(shè)計(jì)。為了排除計(jì)算的偶然性，獲取插值時(shí)長的平均值，對(duì)指定影像總共反復(fù)執(zhí)行100次插值處理，二維插值的試驗(yàn)結(jié)果如表5所示。

表5 試驗(yàn)結(jié)果 ms

可以看出，新的多項(xiàng)式計(jì)算耗時(shí)相比常規(guī)計(jì)算提速比為1.73倍，接近理論值1.625。整個(gè)二維3次卷積計(jì)算與傳統(tǒng)方法相比提速比為1.15倍。這里所用影像數(shù)據(jù)為西安測(cè)繪研究所研制的DMZ大面陣航攝相機(jī)實(shí)攝航空影像，raw數(shù)據(jù)為12個(gè)子影像，內(nèi)插合成后的影像像幅為24 K×24 K(見圖5)。

圖5 DMZ大面陣航攝相機(jī)實(shí)攝航空影像內(nèi)插合成

3 結(jié) 論

本文提出分段多項(xiàng)式核函數(shù)高效插值算法，數(shù)學(xué)形式簡(jiǎn)潔統(tǒng)一。相比較于常規(guī)插值方法，新的3階多項(xiàng)式算法耗時(shí)提速比為1.625倍,對(duì)于任意高階多項(xiàng)式，耗時(shí)提速比上限為1.3倍。這說明新插值算法對(duì)常規(guī)插值算法具有一定的替代性。

對(duì)于大幅面數(shù)字航空攝像機(jī)而言，因?yàn)橛糜趫D像重建的插值計(jì)算量十分巨大，因此，本文方法對(duì)更為快速高效地進(jìn)行原始攝影數(shù)據(jù)預(yù)處理意義非凡，已對(duì)西安測(cè)繪研究所研制的國產(chǎn)DMZ大面陣航測(cè)相機(jī)進(jìn)行相關(guān)試驗(yàn)，由于航測(cè)相機(jī)要求的量測(cè)精度極高(最低子像元級(jí))，利用文中的高階多項(xiàng)式內(nèi)核函數(shù)，可保證所引入的插值誤差被忽略不計(jì)(實(shí)踐表明采用5階多項(xiàng)式已足夠有效)，這使得大幅面影像生成兼顧效率與精度，進(jìn)而保障相機(jī)的使用效能。還需說明的是，所提方法不限于大幅面數(shù)字航空相機(jī)的子影像復(fù)合拼接應(yīng)用，需要高精度、高效率的圖像插值應(yīng)用場(chǎng)合，該算法都可發(fā)揮一定的作用。當(dāng)然，可將此算法進(jìn)一步結(jié)合當(dāng)代GPU并行編程技術(shù)予以實(shí)現(xiàn)，使得實(shí)際作業(yè)效率提升更為顯著。

參考文獻(xiàn)：

[1] 郭海青.信息化航空攝影測(cè)量生產(chǎn)與管理系統(tǒng)建設(shè)[J].測(cè)繪通報(bào),2014(11):53-72.

[2] 劉翠.大面陣CCD成像系統(tǒng)關(guān)鍵技術(shù)[J].科學(xué)與財(cái)富,2015(2):127-128.

[3] 葛明,許永森,沈宏海,等.航空相機(jī)大面陣CCD多自由度拼接方法研究[J].紅外與激光工程,2015(3):923-928.

[5] MEIJERING E,NIESSEN W,VIERGEVER M.Quantitative evaluationofconvolution-basedmethodsformedicalimageinterpolation[J].Medical Image Analysis,2001,5:111-126.

[6] 鐘寶江,陸志芳,季家歡.圖像插值技術(shù)綜述[J].數(shù)據(jù)采集與處理,2016(6):1083-1096.

[7] 楊文波.航空?qǐng)D像超分辨率重構(gòu)技術(shù)研究[D].長春：中國科學(xué)院長春光學(xué)精密機(jī)械與物理研究所,2014.

[8] 喬少華,李潤鑫,劉輝,等.基于統(tǒng)計(jì)量的加權(quán)函數(shù)圖像重建方法[J].傳感器與微系統(tǒng),2017(9):53-56.

[9] 趙秀影,王洪玉,從福仲.亞像元幾何圖像成像與超分辨恢復(fù)算法[J].光學(xué)技術(shù),2010,36(5):795-798.

[10] KEYS R G.Cubic convolution interpolation for digital image processing[C].IEEE Transactions on Acoustics,Speech,and Signal Processing,1981,29(6):1153-1160.

[11] Weisstein E W.Horner’s rule.Eric Weisstein’s World of Mathematics.2003c,http://mathworld.wolfram.com/HornersRule.html