嚴發(fā)寶，蘇艷蕊，趙占鋒，左顥睿，柳建新

?

基于GPU的小尺寸FFT在實時圖像復(fù)原中的優(yōu)化

嚴發(fā)寶1，蘇艷蕊1，趙占鋒2，左顥睿3，柳建新4

(1. 山東大學(xué)(威海) 機電與信息工程學(xué)院，山東威海，264209；2. 哈爾濱工業(yè)大學(xué)(威海) 信息工程研究所，山東威海，264209；3. 中國科學(xué)院光電技術(shù)研究所，四川成都，610209；4. 中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南長沙，410083)

為滿足跟蹤識別系統(tǒng)對圖像復(fù)原的實時性需求，在圖形處理器(GPU)上進行高效實現(xiàn)小尺寸二維FFT的優(yōu)化策略研究。首先對二維FFT算法進行分析，根據(jù)圖形處理器的特點，提出基于圖形處理器的并行執(zhí)行模型?；谠撃Ｐ?，從算法的復(fù)雜度、跳轉(zhuǎn)指令的數(shù)量、共享存儲器的訪問沖突以及共享存儲器的訪問延遲及圖形處理器的利用效率這4個方面進行優(yōu)化策略的研究，提出相應(yīng)的優(yōu)化方法。在圖像復(fù)原的實驗中，先對基于GPU的小尺寸FFT優(yōu)化方法與基于CPU的MATLAB傳統(tǒng)算法進行計算精度對比，然后基于4種不同尺寸的圖像在相同的GPU平臺上再與NVIDIA公司提供CUFFT函數(shù)庫復(fù)原算法進行計算效率對比。研究結(jié)果表明：該優(yōu)化方法提供的圖像復(fù)原算法復(fù)原效果好，與MATLAB效果圖比較人眼觀察不出差異；在計算速率上，提出的優(yōu)化方法能夠在19.6 ms內(nèi)復(fù)原1幀128×128灰度模糊圖像，計算速度與直接采用CUFFT函數(shù)庫算法相比提高約1.8倍。

圖形處理器；小尺寸FFT；圖像復(fù)原；并行優(yōu)化；實時處理

圖像復(fù)原在目標跟蹤、目標識別系統(tǒng)中有著重要的應(yīng)用，但圖像復(fù)原算法特別是盲圖像復(fù)原算法對退化圖像均需要大量的二維快速傅里葉變換(二維FFT)迭代計算[1?4]。二維FFT的數(shù)值計算量大、數(shù)據(jù)依賴性強，同時具有計算密集型和存儲密集型的特點，在常規(guī)處理平臺如FPGA，DSP和通用CPU上難以滿足高速二維信號處理的需求, 因此，采用更高處理能力的圖形處理器(graphic process unit，GPU)成為當前研究的趨勢[5?8]。GPU能夠直接進行單精度和雙精度浮點計算，在圖像處理、雷達信號處理、物理建模等領(lǐng)域中有著廣泛的應(yīng)用[9?13]。ZHANG等[1?3]提出了噪雜退化圖像的盲解卷積及基于乘性迭代算法的多幀盲圖像復(fù)原，分別實現(xiàn)了高分辨率的圖像復(fù)原與多幀圖像的高分辨率復(fù)原；作為圖像復(fù)原情況的一種，ZHANG等[1?3]提出的基于GPU像素級優(yōu)化去霧準則實時單幀圖像去霧方法，以80幀/s的速率去處理百萬像素的大霧照片，效果優(yōu)良，同時該方法受非均勻照明的影響比傳統(tǒng)方法小；XIE等[5]提出的基于GPU的中頻盲圖像復(fù)原并行算法有效地提升了操作速率和復(fù)原的實時性；常方正等[8]提出的基于CUDA的遙感影像CVA變化檢測方法與傳統(tǒng)效果相比效率提高10倍左右；張全等[4]采用相位差算法對尺寸為256×256的圖像進行50次迭代，50次迭代單GPU與雙GPU分別耗時53 ms與45 ms，為目前搜索到的最高計算效率。實際中圖像復(fù)原計算度復(fù)雜，迭代次數(shù)要遠超50次，計算數(shù)據(jù)量大，而且53 ms遠遠低于所需實時性，實時圖像復(fù)原通常在原圖上選取特定區(qū)域?qū)嵤ǔ２淮笥?56×256，即使縮小了處理圖像的尺寸，在GPU上采用目前已有的計算方法，仍不能滿足工程中對圖像復(fù)原的實時需求[14]。為此，本文作者探討如何在GPU上高效實現(xiàn)小尺寸二維FFT的方法，并應(yīng)用到實時圖像復(fù)原處理中。

1 DFT與二維FFT計算流程

1.1 DFT算法

一維點向量對應(yīng)的函數(shù)()的正反離散傅里葉變換(DFT)公式為[15]：

1.2 二維FFT計算流程

FFT算法是傅里葉算法的快速計算方法，詳細計算方法見文獻[16?17]。圖像的數(shù)字化一般都是用二維矩陣進行描述，這里只介紹二維FFT計算的步驟，對于輸入矩陣(，)，具體的計算流程如圖1所示。

圖像數(shù)據(jù)的二維矩陣(，)輸入到處理器，首先對矩陣中的行一維數(shù)據(jù)進行FFT計算，即

得到一維計算結(jié)果(，)后，需要進行FFT計算結(jié)果倒位序計算，得到1(，)，然后進行轉(zhuǎn)置 計算：

再對矩陣(，)進行如式(3)一樣的一維計算，得到2(，)，然后對該矩陣進行倒位序計算后，轉(zhuǎn)置得到最后的二維FFT計算結(jié)果(，)。

圖1 二維FFT計算流程圖

2 GPU上的二維FFT計算

2.1 算法并行特征分析

二維矩陣(，)的FFT的計算主要由一維FFT構(gòu)成，主要分析一維點(=2)FFT的可并行性。

計算分為層，第層計算需要?1層的計算結(jié)果(1＜≤)，層與層之間的計算只能串行；

每層中參與計算的數(shù)據(jù)均為個，每層的計算次數(shù)一樣，每一層內(nèi)的計算可以并行執(zhí)行；

旋轉(zhuǎn)因子W根據(jù)輸入點序號進行計算，任意點之間沒有相關(guān)性，具有良好的可并行性；

二維FFT計算由多個一維FFT計算構(gòu)成，一維FFT之間的計算是相互獨立的，具有良好的可并行性。

2.2 算法并行計算模型

根據(jù)2.1節(jié)的分析，二維FFT計算具有較好的可并行性，能夠較好地利用GPU的并行計算能力，圖2所示為在GPU上對1個×的二維矩陣執(zhí)行二維FFT的并行化模型。

1) 并行任務(wù)劃分。對于×的二維矩陣，要進行+次一維FFT計算，以次點一維FFT計算為例，將行劃分為個GPU任務(wù)分塊，每個塊中包含/2個線程，每個線程操作2個數(shù)據(jù)，基于該任務(wù)劃分，個塊可以并行執(zhí)行。

2) 一維FFT計算。以點一維FFT為例，計算共分為(=log2)個串行層次，每個層次內(nèi)部由/2個線程并行計算得到結(jié)果，每層計算完成后必須對/2個線程進行同步操作，才能保證下一層所需輸入數(shù)據(jù)的正確性。

3 基于GPU的并行優(yōu)化

在實際計算的時候，主要考慮從算法本身的優(yōu)化、指令在GPU中執(zhí)行的優(yōu)化、數(shù)據(jù)存儲的優(yōu)化以及GPU中計算部署的優(yōu)化等4個方面進行討論。

3.1 算法優(yōu)化

在GPU中進行并行計算，想要提升實時性，首先需要從算法本身入手。雖然現(xiàn)有用于GPU的算法也有一定的優(yōu)化效果，但是還可以從以下2個方面進行進一步優(yōu)化。

2) 矩陣轉(zhuǎn)置優(yōu)化。標準的二維FFT計算需要進行2次轉(zhuǎn)置，對于在頻域計算后還將變換回時域計算的情況，可不進行第2次轉(zhuǎn)置直接后續(xù)計算，而轉(zhuǎn)置計算的存儲復(fù)雜度較高，每減少1次轉(zhuǎn)置，內(nèi)存訪問次數(shù)降低2××。

從理論上分析，進行了這2種優(yōu)化后，在計算速度上將會有所提高。

3.2 指令優(yōu)化

指令執(zhí)行是并行計算關(guān)鍵的一步，因此，指令運行導(dǎo)致的流水線問題以及一些復(fù)雜的指令運行都是優(yōu)化的目標。

1) 解循環(huán)優(yōu)化。在一維FFT計算時，通過循環(huán)對層計算進行流程控制，循環(huán)能夠有效降低代碼占用的存儲空間，但是執(zhí)行循環(huán)時，需要多次執(zhí)行判斷和跳轉(zhuǎn)指令，破壞了GPU的執(zhí)行流水線，降低了GPU的效率。對于實際工程，輸入數(shù)據(jù)的長度是已知的，在計算前就能夠計算出一共需要多少層計算，因此可以將整個循環(huán)展開，去除一維FFT計算中流程控制中的判斷和跳轉(zhuǎn)。

2) 復(fù)雜指令優(yōu)化。

①取模計算優(yōu)化：在一維FFT計算、矩陣轉(zhuǎn)置中采用了大量取模計算進行下標計算[16]，取模計算需要進行多次除法計算，計算復(fù)雜。在實際工程中，輸入數(shù)據(jù)的長度是已知的，而且通常是2的指數(shù)，這就可以用執(zhí)行效率高的位運算代替取模計算。當=2時，(MOD)可以用&(?1)代替，可以實現(xiàn)1個時鐘周期輸出計算結(jié)果，無需常規(guī)取模計算所需的數(shù)百個時鐘周期。

圖2 二維FFT并行執(zhí)行模型

②三角函數(shù)計算優(yōu)化：在計算旋轉(zhuǎn)因子W時，需要計算三角函數(shù)cos和sin，直接采用GPU提供的快速硬件指令取代函數(shù)實現(xiàn)，可以實現(xiàn)1個時鐘周期輸出計算結(jié)果，無需常規(guī)直接計算所需的數(shù)百個時鐘周期，能夠有效提高計算效率。

3.3 存儲優(yōu)化

1) 共享存儲器優(yōu)化。在GPU共享存儲器的訪問時，必須盡量避免訪問沖突，當共享存儲器的訪問模式為32位時，訪問一個復(fù)數(shù)需要2次操作，這導(dǎo)致了較多的共享存儲器訪問沖突；而采用64bit共享存儲器訪問方式，能減少沖突訪問的次數(shù)，提高共享存儲器的訪問效率。

2) 用查表法取代復(fù)雜計算。倒位序操作涉及到多次循環(huán)和條件判斷[16?17]，若每個點直接計算對應(yīng)位序，則GPU難以達到較高的處理效率。在實際工程中，可根據(jù)實際計算的點數(shù)，事先計算好倒位序表，在倒位序計算時直接查表，對數(shù)據(jù)進行倒位序操作。

3.4 部署優(yōu)化

GPU由于部署不同，在執(zhí)行效率上會產(chǎn)生較大的差異。一般而言，必須使每個塊中有盡可能多的線程，同時也要保證足夠的分塊數(shù)量，才能充分利用GPU的計算能力。在圖2所示的并行模型中，每個塊中的線程數(shù)量為數(shù)據(jù)長度的一半，每個線程僅操作2個數(shù)據(jù)，無法有效地屏蔽存儲器訪問延遲。因此，在保證分塊數(shù)足夠的前提下，應(yīng)使每個線程操作更多的數(shù)據(jù)，如4個或8個數(shù)據(jù)。在實際工程中，可通過試差法，選擇最優(yōu)的部署。

4 實驗結(jié)果及分析

高速小尺寸二維FFT優(yōu)化方法的提出是為了提高圖像復(fù)原算法的計算速度。因此，在實驗中用乘性迭代算法(MIA)作為基礎(chǔ)算法來對文中的小尺寸二維FFT算法進行精度及運算效率評測[2]。

實驗測試平臺如下：GPU 為NVIDIA公司的GeForce Titan Black(顯存6GB)；CPU是i7 4770K，主頻3.5G(4核)；內(nèi)存16 G；測試程序在VC2010和CUDA6.5環(huán)境下編譯完成。首先對算法的計算精度進行實驗。

實驗1的結(jié)果如圖3所示。其中，圖3(a)所示為原始圖像，尺寸為128×128；圖3(b)所示為人工卷積模糊后的圖像；圖3(c)所示為在CPU上采用Matlab執(zhí)行MIA算法復(fù)原后的圖像；圖3(d)所示為GPU上執(zhí)行MIA算法進行200次迭代復(fù)原后的圖像。

實驗2的結(jié)果如圖4所示。其中，圖4(a)所示為原始圖像，尺寸為256×256；圖4(b)所示為在CPU上采用Matlab執(zhí)行MIA算法復(fù)原后的圖像；圖4(c)所示為GPU上執(zhí)行MIA算法進行200次迭代復(fù)原后的圖像。

從復(fù)原效果來看，Matlab復(fù)原圖像基本與GPU復(fù)原圖像一致，肉眼評價無法觀測出二者之間的差異。表1所示為圖3(c)與圖3(d)、圖3(b)與圖3(c)的差異。從表1可知：GPU復(fù)原圖像的精度非常好，其計算精度和Matlab結(jié)果相差不超過0.2%。

(a) 原始圖像；(b) 降質(zhì)圖像；(c) Matlab復(fù)原圖像；(d) GPU復(fù)原圖像

(a) 原始圖像；(b) Matlab復(fù)原圖像；(c) GPU復(fù)原圖像

表1 GPU復(fù)原圖像與Matlab復(fù)原圖像比較

然后對復(fù)原算法的實時性進行實驗。

為了與現(xiàn)有應(yīng)用于GPU的復(fù)原CUFFT庫圖像復(fù)原算法進行對比，采用32×32，64×64，128×128和256×256這4種尺寸的圖像進行模糊后，分別在相同的GPU平臺上利用CUFFT庫算法和MIA算法進行圖像復(fù)原計算，其中迭代次數(shù)為200次，同時，在目前主流的CPU上對同樣的算法處理速度進行對比實驗。

表2所示為在GPU上使用CUFFT庫算法和本文介紹方法對不同尺寸圖像執(zhí)行MIA算法的時間。算法需要執(zhí)行1 200次二維FFT，執(zhí)行時間為純計算時間，沒有包括從CPU到GPU的數(shù)據(jù)傳輸時間。CUFFT為采用標準CUFFT庫函數(shù)的實驗結(jié)果，CPU_S為采用標準C代碼后的實驗結(jié)果，GPU_F為采用本文優(yōu)化策略后的實驗結(jié)果。

從表2可知：當數(shù)據(jù)尺寸越小，采用本文優(yōu)化策略相比標準CUFFT庫函數(shù)取得的加速比越大；當數(shù)據(jù)規(guī)模為256×256時，本文優(yōu)化策略仍能取得1.32倍的加速比，相比標準CPU至少可以取得1.81倍的加速比，與目前已知基于GPU的最快計算效率相比，仍然優(yōu)越得多。在處理實際使用的128×128圖像時，處理時間控制在19.6 ms以內(nèi)，即50幀/s的處理速率，達到了實時圖像復(fù)原處理速度的要求，并已在實際光電跟蹤系統(tǒng)中得到了應(yīng)用。

表2 耗時測試結(jié)果

5 結(jié)論

1) 根據(jù)跟蹤識別系統(tǒng)對圖像復(fù)原的實時性要求，介紹了在GPU上高效并行實現(xiàn)小尺寸二維FFT的方法，通過分析二維FFT的計算流程和算法特點，對二維FFT計算的可并行性進行了梳理；給出了二維FFT計算在GPU上的并行執(zhí)行模型；提出了算法優(yōu)化、指令優(yōu)化、存儲優(yōu)化以及部署優(yōu)化等基于GPU的并行優(yōu)化策略，充分挖掘了GPU的性能。

2) MIA圖像復(fù)原算法的運算實驗效果表明本文的小尺寸二維FFT計算方法的精度基本與Matlab的一致，在多次迭代后仍能達到較高的精度。文中小尺寸二維FFT算法在處理不同尺寸的圖像時，均達到了較優(yōu)越的運算性能。

3) 在滿足計算精度的基礎(chǔ)上，計算速度滿足實時圖像復(fù)原處理速度的要求，相關(guān)優(yōu)化策略已經(jīng)在多項光電跟蹤系統(tǒng)中得到應(yīng)用，顯示了本文方法的有效性和優(yōu)越性。

[1] ZHANG Jianlin, ZHANG Qiheng, HE Guangming. Blind deconvolution of a noisy degraded image[J]. Applied Optic, 2009, 48(12): 2350?2355.

[2] ZHANG Jun, HU Shiqiang. A GPU-accelerated real-time single image de-hazing method using pixel-level optimal de-hazing criterion[J]. Journal of Real-Time Image Processing, 2014, 9(4): 661?672.

[3] ZHANG Jianlin, ZHANG Qiheng, HE Guangming. Multiframe blind image restoration based on a multiplicative iterative algorithm[J]. Optical Engineering, 2009, 48(2): 027004.

[4] 張全, 鮑華, 饒長輝, 等. 相位差算法在多GPU平臺上的并行化實現(xiàn)[J]. 光電工程, 2016, 43(3): 66?72.ZHANG Quan, BAO Hua, RAO Changhui, et al. Parallel implementation of phase diversity algorithm on multi-GPU[J]. Opto-Electronic Engineering, 2016, 43(3): 66?72.

[5] XIE Lang, LUO Yihan, BAO Qiliang. GPU-based parallel algorithm for blind image restoration using midfrequency-based methods[C]//ISPDI 2013-Fifth International Symposium on Photoelectronic Detection and Imaging. Beijing, 2013: 89101R?89101R-10.

[6] MEYER-BAESE U, MEYER-BAESE A, GONZáLEZ D, et al. Code obfuscation using very long identifiers for FFT motion estimation models in embedded processors[J]. Journal of Real-Time Image Processing, 2016, 11(4): 817?827.

[7] PUCHA?A D, STOKFISZEWSKI K. Effectiveness of fast fourier transform implementations on GPU and CPU[J]. Przeglad Elektrotechniczny, 2016, 92(7): 69?71.

[8] 常方正, 趙銀娣, 劉善磊. 遙感影像CVA變化檢測的CUDA并行算法設(shè)計[J]. 遙感學(xué)報, 2016, 20(1): 114?128.CHANG Fangzheng, ZHAO Yindi, LIU Shanlei. CUDA parallel algorithm for CVA change detection of remote sensing imagery[J]. Journal of Remote Sensing, 2016, 20(1): 114?128.

[9] XU Ming, CHEN Feiguo, LIU Xinhua, et al. Discrete particle simulation of gas-solid two-phase flows with multi-scale CPU-GPU hybrid computation[J]. Chemical Engineering Journal, 2012, 207/208(10): 746?757.

[10] YANG Canqun, WU Qiang, HU Huili, et al. Fast weighting method for plasma PIC simulation on GPU-accelerated heterogeneous systems[J]. Journal of Central South University, 2013, 20(6): 1527?1535.

[11] CHU H, LI T, WANG P. IP address lookup by using GPU[J]. IEEE Transactions on Emerging Topics in Computing, 2016, 4(2): 187?198.

[12] 夏健明，魏德敏. 圖形處理器在大規(guī)模力學(xué)問題計算中的應(yīng)用進展[J]. 力學(xué)進展, 2010, 40(1): 57?63.XIA Jianming, WEI Deming. Advances in graphices processing units’ application to the computation of large-scale mechanical problems[J]. Advances In Mechanics, 2010, 40(1): 57?63.

[13] XIA Chao, GUAN Qingxiao, ZHAO Xianfeng, et al. Highly accurate real-time image steganalysis based on GPU[J]. Journal of Real-Time Image Processing, 2016, 11(Special Issue): 1?14.

[14] 李仕, 王晶, 孫輝. 基于圖形處理器的實數(shù)FFT 在圖像處理中的應(yīng)用[J]. 光學(xué)精密工程, 2008, 16(12): 2414?2420. LI Shi, WANG Jing, SUN Hui. Real FFT based on graphic processing unit for image processing[J]. Optics and Precision Engineering, 2008, 16(12): 2414?2420.

[15] CAVICCHI T. DFT time-domain interpolation[J]. IEE Proceedings F: Radar & Signal Processing, 1992, 139(3): 207?211.

[16] COOLEY J, TUKEY J. An algorithm for the machine calculation of complex Fourier series[J]. Mathematics Computation, 1965, 19: 296?301.

[17] ZHANG Xinxue, WANG Guizeng. Parallel FFT architecture consisting of FFT chips[J]. Journal of Circuits and Systems, 2000, 5(2): 38?42.

(編輯 楊幼平)

Optimization on FFT of small size in real-time image restoration based on GPU

YAN Fabao1, SU Yarui1, ZHAO Zhanfeng2, ZUO Haorui3, LIU Jianxin4

(1. School of Mechanical, Electrical & Information Engineering, Shandong University, Weihai, Weihai 264209, China;2. Institute of Information Engineering, Harbin Institute of Technology, Weihai, Weihai 264209, China;3. Institute of Optics and Electronics, Chinese Academy of Sciences, Chengdu 610209, China;4. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China)

To meet the real-time demand of image restoration for recognition and tracking system, an optimization research on two-dimensional FFT of small size realized in graphics processor unit(GPU) efficiently was done. An analysis of two-dimensional FFT algorithm was analyzed first. And according to the characteristics of GPU, a parallel execution model based on graphics processor was proposed. Based on this model, the optimization research was done considering the aspects of algorithm complexity, the number of jump instructions, access conflict and access latency of the shared memory, and the utilization efficiency of GPU. And two-dimensional FFT computation of small size was realized in the GPU. In image restoration experiment, comparison on the calculation accuracy of two-dimensional FFT of small size optimization algorithm based on GPU and the traditional algorithm in MATLAB based on CPU was done. And a comparison on the computational efficiency of optimization algorithm proposed and the library function image restoration algorithm of CUFFT provided by NVIDIA Corp in four different sizes based on the same GPU platform was made. The results indicate that this optimization algorithm has excellent recovery performance, and human vision system could not distinguish the difference between the results and the MATLAB demonstrations. And the optimization algorithm can recover a frame of 128×128 gray fuzzy image within 19.6 ms, while the computing speed increases 1.8 times approximately compared with that using library function of CUFFT directly.

graphic processing unit (GPU); FFT of small size; image restoration; parallel optimization; real-time computation

10.11817/j.issn.1672?7207.2017.10.019

TP39

A

1672?7207(2017)10?2691?06

2016?10?13；

修回日期：2017?01?16

國家科技基礎(chǔ)性工作專項(2013FY110800)；中國博士后科學(xué)基金資助項目(2016M600538)；國家自然科學(xué)基金資助項目(41674080，41574123，21505028)(Project(2013FY110800) supported by the National Science and Technology Basic Work; Project (2016M600538) funded by China Postdoctoral Science Foundation; Projects(41674080, 41574123, 21505028) supported by the National Science Foundation of China)

蘇艷蕊，博士，講師，從事跟蹤控制平穩(wěn)性研究、嵌入式系統(tǒng)研究，E-mail：suyanrui@126.com