基于CUDA的遙感圖像快速超分辨率重建算法研究

陳 昕 胡 杰 李 映

(西北工業(yè)大學(xué) 西安 710129)

0 引言

圖像超分辨率重建技術(shù)的基本思想是采用信號(hào)處理的方法，在改善圖像質(zhì)量的同時(shí)，重建成像系統(tǒng)截止頻率之外的信息，從而在不改變硬件設(shè)備的前提下，獲取高于成像系統(tǒng)空間分辨率的圖像。如圖1所示，由于成像系統(tǒng)中存在不同的降質(zhì)過(guò)程［1］，理想的高分辨率圖像x在成像過(guò)程中會(huì)受到形變、模糊(包括光學(xué)、傳感器和幀內(nèi)運(yùn)動(dòng)模糊)、抽樣和噪聲等因素的影響，造成圖像質(zhì)量的下降。實(shí)際應(yīng)用中，對(duì)同一場(chǎng)景的多次成像會(huì)得到k幀低分辨率圖像:

式中，M(k)為第k幀低分辨率圖像成像過(guò)程中由運(yùn)動(dòng)引起的形變矩陣;B(k)為點(diǎn)擴(kuò)散函數(shù)或衍射極限引起的模糊矩陣;D為下采樣矩陣，反映了成像系統(tǒng)對(duì)原始高分辨率場(chǎng)景的下采樣過(guò)程;n(k)為噪聲，通常被認(rèn)為是加性高斯白噪聲。

從圖1 可以看出，圖像的超分辨率重建過(guò)程就是根據(jù)獲取的一幀或多幀低分辨率圖像y(k{)} 重建理想的高分辨率圖像x的過(guò)程，即在對(duì)成像系統(tǒng)通帶內(nèi)獲取的圖像低頻信息進(jìn)行復(fù)原的基礎(chǔ)上，恢復(fù)其截止頻率以上的高頻信息的過(guò)程。

圖1 圖像獲取過(guò)程示意圖

由于在頻域中很難有效表達(dá)圖像中目標(biāo)的運(yùn)動(dòng)信息，近些年大部分研究工作都集中在空域法。其中主要包括反向迭代投影法(IBP，Iterative Back Projection)［2］，凸集投影法(POCS，Projection Onto Convex Sets)［3］，最大后驗(yàn)概率估計(jì)法(MAP，Maximum A Posteriori)［4］和一些濾波方法，相比傳統(tǒng)的插值方法，能夠更好的重建圖像的高頻信息(如邊緣，結(jié)構(gòu)信息等)，獲得更好的視覺(jué)效果。在稀疏表示框架下基于字典學(xué)習(xí)的圖像超分辨率重建算法［5-6］提出了一種更加新穎的思路，從信號(hào)、圖像底層的表示模型出發(fā)，討論高、低分辨率圖像之間的關(guān)系，進(jìn)而確定高、低分辨率圖像之間的映射關(guān)系以及相應(yīng)的約束，該類方法獲得了高質(zhì)量的重建結(jié)果。然而由于字典學(xué)習(xí)、稀疏表示等過(guò)程具有較高的計(jì)算復(fù)雜度，很難直接應(yīng)用在實(shí)際中。在保證高質(zhì)量的圖像超分辨率重建效果前提下，獲得接近或達(dá)到工業(yè)級(jí)(25 幀每秒)的處理速度，一直是圖像/視頻超分辨率重建領(lǐng)域的一大難題，但同時(shí)也是一個(gè)具有重大研究和實(shí)際意義的研究方向。2007年6月NVIDIA 公司推出了統(tǒng)一設(shè)備計(jì)算架構(gòu)(CUDA，Compute Unified Device Architecture)。CUDA 是一種將GPU 作為數(shù)據(jù)并行計(jì)算設(shè)備的軟硬件體系，被廣泛應(yīng)用于天文計(jì)算、圖像處理、音視頻編解碼等領(lǐng)域，獲得了幾倍、幾十倍，乃至上百倍的加速比。CUDA 為加速圖像/視頻處理提供了一個(gè)很好的平臺(tái)，使得圖像實(shí)時(shí)/快速超分辨率重建成為可能。文章重新設(shè)計(jì)并實(shí)現(xiàn)了Non-local IBP［7］算法，實(shí)現(xiàn)了基于CUDA 的遙感圖像快速超分辨率重建。

1 Non－local IBP 算法

Non-local IBP［7］算法的主要思想是利用圖像含有非局部冗余信息的特點(diǎn)來(lái)改善超分辨率重建圖像的質(zhì)量，并且通過(guò)結(jié)合IBP［2］來(lái)減小重建圖像的誤差。如圖2所示，大多數(shù)自然圖像含有許多重復(fù)的結(jié)構(gòu)。這些重復(fù)的結(jié)構(gòu)，即非局部冗余信息對(duì)于求解圖像處理問(wèn)題方面的逆問(wèn)題都有很大的輔助作用。算法流程如圖3所示。

其中:

圖2 圖像中重復(fù)的結(jié)構(gòu)

圖3 Non-local IBP 算法流程

其中‖·‖1是L1范數(shù)，sum_pixel 是方形鄰域中像素點(diǎn)的個(gè)數(shù)。如果圖像塊之間的距離d(i，j)比預(yù)設(shè)的閾值小，說(shuō)明兩個(gè)圖像塊是相似的，并把像素y(i，j)歸于集合S(i0，j0)。S(i0，j0)中包含以及和其相似的像素。在文章的實(shí)驗(yàn)部分中，只選擇保存7 個(gè)最小的MAD 值所對(duì)應(yīng)的像素點(diǎn)y(i，j)。同時(shí)，記錄下和y(i，j)坐標(biāo)之差，即vector_x= i0-i，vector_y= j0-j。盡管塊匹配處理可以在整個(gè)圖像上進(jìn)行，但是計(jì)算量過(guò)大。因此，在實(shí)驗(yàn)中，我們把搜索限定在一個(gè)相對(duì)較大(21×21)的窗口中。

b.圖像重建誤差:

Il是低分辨率圖像;Ih是高分辨率圖像;D是下采樣矩陣;* 是卷積算子;G是點(diǎn)擴(kuò)散方程(PSF)［3］，通常是一個(gè)高斯核;是圖像重建誤差。

c.利用運(yùn)動(dòng)矢量進(jìn)行圖像融合:第tth次迭代時(shí)計(jì)算圖像的重建誤差，用重建誤差來(lái)更新:

其中F(·)是一個(gè)bicubic 插值算子，p是一個(gè)向后投影核來(lái)控制收斂的速度，對(duì)于高分辨率像素的更新，用運(yùn)動(dòng)矢量估計(jì)得到的7 個(gè)低分辨率像素進(jìn)行更新。

2 Non －local IBP 算法在CUDA 上的實(shí)現(xiàn)

CUDA 程序的執(zhí)行包含兩部分，一部分是在CPU 上串行執(zhí)行的主機(jī)(Host)代碼，另一部分是在GPU 上并行執(zhí)行的設(shè)備(Device)代碼。運(yùn)行在GPU 上的函數(shù)稱為核函數(shù)(Kernel)，而運(yùn)行在CPU上的串行代碼主要完成Kernel 啟動(dòng)之前的數(shù)據(jù)準(zhǔn)備和設(shè)備初始化工作，以及在Kernel 之間進(jìn)行的一些串行運(yùn)算。CUDA 程序一般包括以下過(guò)程:

在CPU 上初始化數(shù)據(jù);

將數(shù)據(jù)傳輸給GPU;

調(diào)用Kernel 進(jìn)行并行計(jì)算;

將結(jié)果數(shù)據(jù)傳回CPU。

Kernel 以線程網(wǎng)格(Grid)的形式組織，每個(gè)線程網(wǎng)格包含若干個(gè)線程塊(Block)，而每個(gè)線程塊包含若干個(gè)線程(Thread)，如圖4所示。實(shí)質(zhì)上，Kernel 是以Block 為單位執(zhí)行的，CUDA 引入Grid 只是用來(lái)表示一系列可以被并行執(zhí)行的Block 的集合。不同Block 之間可以并行執(zhí)行［8］，但無(wú)法進(jìn)行通信，也沒(méi)有執(zhí)行順序;同一Block 內(nèi)的線程不僅能夠并行執(zhí)行，而且能夠通過(guò)共享存儲(chǔ)器(shared memory)和_syncthreads()函數(shù)相互通信。在Kernel 函數(shù)啟動(dòng)之前需要確定Grid 和Block 的維度。CUDA 提供了dim3 類型(基于uint3 定義的矢量類型，由3 個(gè)unsigned int 組成的結(jié)構(gòu)體)的內(nèi)建變量threadIDx 和blockIDx。這樣，就可以使用一維、二維或三維的索引來(lái)標(biāo)識(shí)線程，構(gòu)成一維、二維或三維的線程塊。

如圖3所示，Non-local IBP 算法主要由bicubic 插值，運(yùn)動(dòng)矢量估計(jì)，高斯處理，圖像差值和圖像融合這5 個(gè)步驟構(gòu)成，對(duì)應(yīng)到CUDA 架構(gòu)下相應(yīng)的核函數(shù)依次為bicubic Kernel，運(yùn)動(dòng)矢量估計(jì)Kernel［9］，高斯模糊Kernel，圖像差值Kernel 和圖像融合Kernel 這5 個(gè)Kernel 函數(shù)。在核函數(shù)被調(diào)用前，必須配置核函數(shù)的執(zhí)行參數(shù)。例如，對(duì)于352×288的低分辨率遙感圖像進(jìn)行放大因子為2 的超分辨率重建，bicubic Kernel 核函數(shù)用來(lái)對(duì)遙感圖像并行地進(jìn)行雙三次插值。線程塊Block 的大小blockSize 可取(blockSize.x= 16，blockSize.y= 16)，則線程網(wǎng)格Grid 的大小gridSize 為(gridSize.x= 352/block-Size.x，gridSize.y= 288/blockSize.y)。bicubic Kernel 的配置為bicubic Kernel＜＜＜gridSize，blockSize＞＞＞，一個(gè)線程處理圖像的一個(gè)像素點(diǎn)。由于GPU實(shí)際上是以線程束(Warp)為單位來(lái)調(diào)度線程執(zhí)行計(jì)算的，因此通常設(shè)定M為線程束大小(即32)的倍數(shù)，以更加方便、有效地對(duì)線程進(jìn)行管理。CUDA程序執(zhí)行如圖5所示，首先將圖像數(shù)據(jù)讀入主機(jī)端內(nèi)存，再將數(shù)據(jù)從內(nèi)存?zhèn)魉偷皆O(shè)備端顯存，然后設(shè)備端的核函數(shù)，即bicubic Kernel，運(yùn)動(dòng)矢量估計(jì)Kernel，高斯模糊Kernel，圖像差值Kernel 和圖像融合Kernel 依次對(duì)顯存數(shù)據(jù)做并行處理，最后將得到的超分辨率重建的圖像數(shù)據(jù)從設(shè)備端顯存?zhèn)魉突刂鳈C(jī)端內(nèi)存。

圖4 線程層次結(jié)構(gòu)

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

本實(shí)驗(yàn)實(shí)現(xiàn)了基于CUDA 的Non-local IBP 算法，對(duì)真實(shí)遙感圖像進(jìn)行了尺度大小為2 的超分辨率重建，并與用傳統(tǒng)bicubic 插值算法得到的結(jié)果進(jìn)行了對(duì)比。如圖6、圖7 和圖8 三種狀態(tài)所示，本文提出的算法重建得到的超分辨率圖像的細(xì)節(jié)更豐富，對(duì)比度更高，邊緣信息更加清晰，圖像的結(jié)構(gòu)特征基本保持完整。

圖5 CUDA 程序執(zhí)行

表1 是Non-local IBP 算法在CPU 和GPU 兩種環(huán)境下的仿真實(shí)現(xiàn)的速度對(duì)比。如表1所示，相較于CPU 的遙感圖像超分辨率重建速度，GPU 的處理速度有了質(zhì)的提升，實(shí)現(xiàn)了對(duì)遙感圖像的快速超分辨率重建。

表2 是Non-local IBP 和bicubic 兩種算法得到的遙感圖像超分辨率重建質(zhì)量的定量評(píng)價(jià)。BRISQUE［10］(Blind/Referenceless Image Spatial Quality Evaluator，盲圖像空間質(zhì)量評(píng)估)的值越小(最小值為0，最大值為100)，說(shuō)明重建圖像的質(zhì)量越好。由圖6，圖7，圖8，表1，表2 可以看出，基于CUDA 實(shí)現(xiàn)的超分辨率重建算法不僅滿足了快速處理的要求，相比bicubic 插值得到的重建圖像，質(zhì)量要高出很多，增加了圖像的對(duì)比度，邊緣保持得很清晰，人工痕跡很少。

圖6 遙感圖像超分辨率重建對(duì)比

圖7 遙感圖像超分辨率重建對(duì)比

圖8 遙感圖像超分辨率重建對(duì)比

表1 GPU 和CPU 實(shí)現(xiàn)Non-local 算法的速度對(duì)比(單位:幀/秒)

表2 Non-local IBP 和bicubic 算法超分辨率重建圖像的BRISQUE

4 結(jié)束語(yǔ)

文章基于Non-local IBP 算法，在CUDA 架構(gòu)下實(shí)現(xiàn)了遙感圖像的快速超分辨率重建算法，合理利用了Non-local 特性在圖像處理中對(duì)描述圖像結(jié)構(gòu)特征的優(yōu)勢(shì)，并能同時(shí)包含全局和局部結(jié)構(gòu)特征信息，使得Non-local IBP 超分辨率重建算法得到的圖像細(xì)節(jié)更豐富，對(duì)比度更高，邊緣更加清晰;并且充分地利用GPU 超強(qiáng)的并行計(jì)算能力，使得超分辨率重建算法的速率得到了質(zhì)的提高。

［1］Park S C，Parkm K，Kangm G.Super-resolution image reconstruction:A technical review［J］.IEEE Signal Processing Magazine，2003，20(3):21-36.

［2］S.Dai，M.Han，Y.Wu，and Y.Gong.Bilateral back projection for single image super resolution［C］.In Proc.Int.Conf.on Multimedia and Expo，2007，1039-1042.

［3］H.Stark，P.Oskoui.High resolution image recovery from image-plane arrays，using convex projections［J］.J.Opt.Soc.Am.A，1989，6:1715-1726.

［4］M.Elad，A.Feuer.Restoration of a single superresolution image from several blurred，noisy and undersampled measured images［J］.IEEE Transaction on Image Processing，1997，6(12):1646-1658.

［5］Jianchao Yang，John Wright，Thomas Huang，and Yi Ma.Image super-resolution as sparse representation of raw image patches ［C］.In Proceedings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition，2008，1-8.

［6］Jianchao Yang，John Wright，Thomas Huang，and Yi Ma.Image super-resolution via sparse representation［J］.IEEE Transactions on Image Processing，2010，19(11):2861-2873.

［7］Dong W，Zhang D，Shi G，et al.Nonlocal back-projection for adaptive image enlargement［C］.Image Processing(ICIP)，2009 16thIEEE International Conference on.IEEE，2009，349-352.

［8］張舒，褚艷利，趙開(kāi)勇，張玉勃.GPU 高性能計(jì)算之CUDA ［M］.北京:中國(guó)水利出版社，2009.

［9］E.Monteiro，B.Vizzotto，C.Diniz，B.Zatt，and S.Bampi.Applying CUDA architecture to accelerate full search block matching algorithm for high performance motion estimation in video encoding［C］.2011 23rd International Symposium on Computer Architecture and High Performance Computing(SBAC-PAD)，2011，128-135.

［10］A.Mittal，A.K.Moorthy，and A.C.Bovik.No-reference image quality assessment in the spatial domain ［J］.IEEE Trans.Image Process.，2012，21(12):4695-4708.