趙明超 陳智斌 文有為

摘要：研究基于總變分（TV）的圖像去噪問題，針對中央處理器（CPU）計算速度較慢的問題，提出了在圖像處理器（GPU）上并行計算的方法?？紤]總變分最小問題的對偶模型，建立原始變量與對偶變量的關(guān)系，采用梯度投影算法求解對偶變量。數(shù)值實驗分別在GPU與CPU上進(jìn)行。實驗結(jié)果表明，總變分去噪模型對偶算法在GPU設(shè)備上執(zhí)行的效率高于在CPU上執(zhí)行的效率，并且隨著圖像尺寸的增大，GPU并行計算的優(yōu)勢更加突出。

關(guān)鍵詞：并行計算；總變分；圖像去噪； 圖像處理器

中圖分類號：TN911.73 文獻(xiàn)標(biāo)志碼：A

Abstract： The problem of Total Variation （TV）based image denoising was considered. Since the traditional serial computation speed based on Central Processing Unit （CPU） was low， a parallel computation based on Graphics Processing Unit （GPU） was proposed. The dual model of the total variationbased image denoising was derived and the relationship between the primal variable and the dual variable was considered. The projected gradient method was applied to solve the dual model. Numerical results obtained by CPU and GPU show that the algorithm implemented by GPU is more efficient than that by CPU， and with the increasing of image size， the advantage of GPU parallel computing is more outstanding.

Key words：parallel computation； Total Variation （TV）； denoising； Graphics Processing Unit （GPU）

0 引言

在物理成像系統(tǒng)以及傳輸過程中，圖像常常被噪聲污染而退化，從而造成視覺損傷，故圖像去噪是圖像處理過程中的一項重要操作。圖像去噪經(jīng)過幾十年的研究，已取得重大成果，但它仍然是充滿活力的領(lǐng)域[1]。隨著科技的發(fā)展，圖像尺寸持續(xù)增加，要求新算法提出的同時也需要計算硬件設(shè)備的更新。由過去單純提高單核時鐘頻率來提升設(shè)備計算性能的技術(shù)達(dá)到瓶頸后，利用多核并行計算提升設(shè)備計算性能已成為研究的熱點。

圖像處理器（Graphics Processing Unit， GPU）是早期為實現(xiàn)圖形實時渲染著色而開發(fā)的圖形處理設(shè)備，它擁有眾多計算核心和高帶寬，因此具有很高的計算吞吐量[2]。實驗表明，它是大規(guī)模并行計算得以實現(xiàn)的合適硬件設(shè)備。經(jīng)過數(shù)十年的開發(fā)研究，一種基于統(tǒng)一計算設(shè)備單元（Compute Unified Device Architecture， CUDA）新架構(gòu)構(gòu)建的GPU可以方便地使用C語言的擴展語言實現(xiàn)并行化計算[3]。CUDA C一經(jīng)面世，它便成為人們主要關(guān)注的對象。許多重要領(lǐng)域都積極研究開發(fā)基于GPU的應(yīng)用程序，例如，模式識別、基因（DeoxyriboNucleic Acid，DNA）序列校對、計算流體力學(xué)、量子力學(xué)和環(huán)境科學(xué)等。隨著CUDA技術(shù)的逐漸成熟，圖像處理領(lǐng)域的研究者也正在積極地將其引入到該領(lǐng)域中，特別是醫(yī)學(xué)圖像TechniScan的基于CUDA架構(gòu)的超聲波成像系統(tǒng)，使得醫(yī)生可在20min內(nèi)獲得患者高清三維圖像。

在圖像處理過程中，一幅灰度圖像可以對應(yīng)一個二維矩陣或一個列向量，列向量是由二維矩陣轉(zhuǎn)化來的。由Rudin，Osher和Fatemi提出的總變分去噪模型（RudinOsherFatemi， ROF）[4]，可以高效去除圖像噪聲，保留圖像的邊緣信息。ROF模型是根據(jù)式（1）得出：

近年來，研究者提出了許多關(guān)于總變分模型的算法，例如原對偶算法[6]、對偶算法[7-8]和Chambolle 算法[9]。這些算法包含相當(dāng)?shù)挠嬎懔?，隨著圖片像素的不斷提高以及圖片尺寸的增大，給實時圖像處理帶來巨大的挑戰(zhàn)。盡管GPU的新架構(gòu)CUDA可以實現(xiàn)并行計算，但并不是所有的算法可以不加更改地進(jìn)行并行計算。一個好的并行算法，可以將費時的計算劃分為一系列獨立的統(tǒng)一操作的計算，同時需要考慮數(shù)據(jù)之間的切換，因為GPU的存儲空間有限，數(shù)據(jù)間的交換同樣會花費大量的時間，那么在設(shè)計并行計算時，需要盡可能地將數(shù)據(jù)傳輸降到最低。

隨著總變分模型應(yīng)用于不同的圖像恢復(fù)任務(wù)，該模型也在不斷變換，以適應(yīng)新的環(huán)境，如文獻(xiàn)[8，10]。而在文獻(xiàn)[11]中，研究的是在GPU上實現(xiàn)TVL1正則模型的原始與對偶方法；文獻(xiàn)[12]探究了TV模型在醫(yī)學(xué)領(lǐng)域核磁共振圖像去噪的應(yīng)用，并且提出了解決TV正則化參數(shù)估計問題的方法。本文主要研究TVL2正則模型對偶算法的并行計算問題，將其在GPU的CUDA架構(gòu)上實現(xiàn)，并與TVL2正則化模型對偶算法在CPU上的執(zhí)行進(jìn)行對比。

2 GPU實現(xiàn)

對偶算法是在GPU的新架構(gòu)CUDA上實現(xiàn)的。在CUDA上運行的函數(shù)稱為Kernel（內(nèi)核函數(shù)）。Kernel以線程網(wǎng)格的形式組織，每個線程網(wǎng)格由若干個塊組成，每個線程塊又由若干線程組成。對于給定型號的GPU，一個塊中可開辟最大線程數(shù)量是固定的。CUDA將計算任務(wù)映射為大量可以并行執(zhí)行的程序，并由硬件動態(tài)調(diào)度和執(zhí)行這些線程。只要聲明了執(zhí)行參數(shù)（設(shè)備的限制內(nèi)），GPU設(shè)備會自動將數(shù)據(jù)分配到相應(yīng)的處理單元上。

圖像去噪空間差分算子和它的轉(zhuǎn)置以及離散偏差運算都是主要的耗時運算操作。本文將在GPU上開啟n個線程同時進(jìn)行運算，每一個線程返回一個單精度的結(jié)果。當(dāng)進(jìn)行空間操作時，如空間差分第（i， j）個位置進(jìn)行運算時需要用到它相鄰位置的數(shù)據(jù)，不同塊線程之間不能進(jìn)行數(shù)據(jù)共享，故此處應(yīng)用了GPU的緩存技術(shù)[2]。

當(dāng)算法操作執(zhí)行一個數(shù)據(jù)集的歸約運算[2]，如對偶投影的計算，需要謹(jǐn)慎處理這樣的操作，因為程序?qū)懙貌划?dāng)可能使運算時間加倍或者產(chǎn)生線程訪問數(shù)據(jù)沖突。在這樣的運算中，應(yīng)將數(shù)據(jù)先分成大小合適的塊，每個塊執(zhí)行部分?jǐn)?shù)據(jù)求和，然后再將塊中部分?jǐn)?shù)據(jù)和再求和，使得GPU處理器眾核得到高效利用。當(dāng)進(jìn)行分配塊時，需要考慮數(shù)據(jù)傳輸?shù)膸捪拗坪虶PU本身資源的限制如寄存器、內(nèi)存等。

本文算法在實現(xiàn)一個全局求和的過程中，每個塊有256個線程，每個線程讀取全局內(nèi)存中的1個數(shù)據(jù)，然后將它們放到一個大小16×16的共享內(nèi)存塊中。每個塊執(zhí)行相同加和操作7次，然后將塊中的256個線程中的數(shù)據(jù)再求和放到事先開辟的部分和全局向量（n2/2048）內(nèi)存中，這樣每個塊可以執(zhí)行2048個數(shù)據(jù)求和，最后將部分和全局向量內(nèi)存中的數(shù)據(jù)求和即可完成本次操作。

在執(zhí)行不依賴彼此的運行過程中，本文用到了GPU并行計算當(dāng)中的流處理[15]。每個流是CUDA當(dāng)中的一個同步單元。一個同步點可以定義多個由CUDA函數(shù)cudaThreadSynchronize（）控制的流。使用流，可以同時更新對偶變量，因為不同行的更新是獨立的。

3 數(shù)值結(jié)果

理論上，在GPU和CPU上執(zhí)行的代碼迭代次數(shù)（Iterations，Iter）應(yīng)該是相同的，但是當(dāng)達(dá)到一定精度后，就會出現(xiàn)迭代次數(shù)的差異，這是由于硬件的限制，算法在CPU上執(zhí)行的數(shù)值運算是雙精度，而在GPU上執(zhí)行的是單精度。圖1～2中也反映出，隨λ的取值變化，算法在GPU上與CPU上的執(zhí)行時間的變化不一致。因為在GPU上執(zhí)行時精度比較低，因此它對參數(shù)變化的靈敏度比較低，而在CPU上執(zhí)行時精度比較高，故它對參數(shù)變化的靈敏度較高，所以在CPU上執(zhí)行算法時時間隨參數(shù)的變化比較明顯。

由于λ的取值也影響著算法迭代的時間，為了方便在不同設(shè)備上運行的時間進(jìn)行對比，那么λ的取值盡量使得在GPU與CPU上運行迭代次數(shù)相同或差異控制在5%以內(nèi)。表1實驗結(jié)果表明，上述算法在GPU設(shè)備上執(zhí)行的效率高于在CPU上，并且隨著計算規(guī)模的增加，GPU加速效果更加明顯。

4 結(jié)語

本文借助GPU計算設(shè)備，實現(xiàn)總變分圖像去噪對偶算法的并行計算，解決總變分模型在實時圖像處理時計算量大、耗時較長的問題。本文僅研究了灰度圖像的并行計算，可以更進(jìn)一步探究彩色圖像的并行計算。在應(yīng)用總變分圖像去噪時沒有考慮模型中參數(shù)選擇的問題，參數(shù)選擇在圖像處理中是一個難點和熱點問題。GPU并行計算并不僅限于圖像處理這一領(lǐng)域，它主要針對具有計算規(guī)模大，耗時且實時性要求高的問題，對于這樣的問題，都可以在GPU設(shè)備上進(jìn)行并行計算研究。

參考文獻(xiàn)：

[1]CHATTERJEE P， MILANFAR P. Is denoising dead？[J]. IEEE Transactions on Image Processing， 2010， 19（4）： 895-911.

[2]張舒， 褚艷利. GPU高性能運算之CUDA[M]. 北京：中國水利水電出版社，2009：5-13，141-189.（ZHANG S， CHU Y L. CUDA High Performance Computing GPU[M]. Beijing： China Water & Power Press， 2009： 5-13，141-189.）

[3]SANDERS J， KANDROT E. GPU高性能編程CUDA實戰(zhàn)[M]. 聶學(xué)軍，譯. 北京： 機械工業(yè)出版社， 2011： 3-8. （SANDERS J， KANDROT E. CUDA by Example： an Introduction to GeneralPurpose GPU Programming[M]. NIE X J， translated. Beijing： China Machine Press， 2011： 3-8.）

[4]RUDIN L， OSHER S， FATEMI E. Nonlinear total variation based noise removal algorithms[J]. Physica D： Nonlinear Phnomena， 1992， 60（1）： 259-268.

[5]CHAMBOLLE A， LIONS P L. Image recovery via total variation minimization and related problem[J]. Numerische Mathematik， 1997， 76（2）： 167-188.

[6]CHAN T F， GOLUB G H， MULET P. A nonlinear primaldual method for total variation based image restoration[J]. SIAM Journal on Scientific Computing， 1996， 20（6）： 1964-1997.

[7]CARTER J L. Dual methods for total variationbased image restoration[D]. Los Angeles： University of California， 2001： 8-41.

[8]ZHU M， WRIGHT S J， CHAN T F. Dualitybased algorithms for total variation image restoration[J].Computational Optimization and Applications， 2010， 47（3）： 377-400.

[9]CHAMBOLLE A. An algorithm for total variation minimization and applications[J]. Journal of Mathematical Imaging and Vision， 2004， 20（1/2）： 89-97.

[10]CHAN T F， ESEDOGLU S， PARK F， et al. Total Variation Image Restoration： Overview and Recent Developments[M]. New York： Springer US， 2006： 17-31.

[11]POCK T， UNGER M， CREMERS D， et al. Fast and exact solution of total variation models on the GPU[C]// Proceedings of the 2008 IEEE Computer Society Conference on Computer Vision and Pattern Recognition Workshops. Washington， DC： IEEE Computer Society， 2008： 1-8.

[12]LIU L R， SHI L， HUANG W H， et al. Generalized total variationbased MRI Rician denoising model with spatially adaptive regularization parameters[J]. Magnetic Resonance Imaging， 2014， 32（6）： 702 -720.

[13]BOYD S， VANDENBERGHE L.凸優(yōu)化[M].王書寧，許鋆， 黃曉霖，譯. 北京：清華大學(xué)出版社， 2013： 85-89.（BOYD S， VANDENBERGHE L. Convex Optimization[M].WANG S N， XU Y， HUANG X L， translated. Beijing： Tsinghua University Press， 2013： 85-89.）

[14]BERTSEKAS D P. Nonlinear Programming[M]. 2nd ed. Nashua： Athena Scientific， 1999： 9.

[15]KIRK D B， HWU W M.大規(guī)模并行處理器編程實戰(zhàn)[M].陳曙暉， 熊淑華，譯. 北京：清華大學(xué)出版社，2010：65-78.（KIRK D B， HWU W M. Programming Massively Parallel Processors： a Handson Approach[M]. CHEN S H， XIONG S H， translated. Beijing： Tsinghua University Press， 2010：65-78.）

[16]DERIN B S， RAFAEL M， KATSAGGELOS A K. Parameter estimation in TV image restoration using variational distribution approximation[J].IEEE Transactions on Image Processing，2008，17（3）：326-339.