基于GPU的小尺寸FFT在實時圖像復原中的優化

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

?

基于GPU的小尺寸FFT在實時圖像復原中的優化

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

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

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

圖形處理器；小尺寸FFT；圖像復原；并行優化；實時處理

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

1 DFT與二維FFT計算流程

1.1 DFT算法

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

1.2 二維FFT計算流程

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

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

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

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

圖1 二維FFT計算流程圖

2 GPU上的二維FFT計算

2.1 算法并行特征分析

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

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

每層中參與計算的數據均為個，每層的計算次數一樣，每一層內的計算可以并行執行；

旋轉因子W根據輸入點序號進行計算，任意點之間沒有相關性，具有良好的可并行性；

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

2.2 算法并行計算模型

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

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

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

3 基于GPU的并行優化

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

3.1 算法優化

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

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

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

3.2 指令優化

指令執行是并行計算關鍵的一步，因此，指令運行導致的流水線問題以及一些復雜的指令運行都是優化的目標。

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

2) 復雜指令優化。

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

圖2 二維FFT并行執行模型

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

3.3 存儲優化

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

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

3.4 部署優化

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

4 實驗結果及分析

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

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

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

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

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

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

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

表1 GPU復原圖像與Matlab復原圖像比較

然后對復原算法的實時性進行實驗。

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

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

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

表2 耗時測試結果

5 結論

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

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

3) 在滿足計算精度的基礎上，計算速度滿足實時圖像復原處理速度的要求，相關優化策略已經在多項光電跟蹤系統中得到應用，顯示了本文方法的有效性和優越性。

[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平臺上的并行化實現[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并行算法設計[J]. 遙感學報, 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] 夏健明，魏德敏. 圖形處理器在大規模力學問題計算中的應用進展[J]. 力學進展, 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] 李仕, 王晶, 孫輝. 基于圖形處理器的實數FFT 在圖像處理中的應用[J]. 光學精密工程, 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

國家科技基礎性工作專項(2013FY110800)；中國博士后科學基金資助項目(2016M600538)；國家自然科學基金資助項目(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)

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