劉進(jìn)鋒

寧夏大學(xué) 數(shù)學(xué)計(jì)算機(jī)學(xué)院，銀川 750021

幾種CUDA加速高斯濾波算法的比較

劉進(jìn)鋒

寧夏大學(xué) 數(shù)學(xué)計(jì)算機(jī)學(xué)院，銀川 750021

1 引言

圖像濾波是圖像處理時(shí)常用的方法。圖像濾波總體上講包括空域?yàn)V波和頻域?yàn)V波。頻域?yàn)V波需要先進(jìn)行傅里葉變換至頻域處理，然后再反變換回空域還原圖像?？沼?yàn)V波是一種鄰域運(yùn)算，即輸出圖像中任何像素的值都是通過(guò)采用一定的算法，根據(jù)輸入圖像中一定鄰域內(nèi)像素的值得來(lái)的。如果輸出像素是輸入像素鄰域像素的線性組合則稱為線性濾波（例如最常見(jiàn)的均值濾波和高斯濾波），否則為非線性濾波（中值濾波、邊緣保持濾波等）。線性濾波器使用連續(xù)窗函數(shù)內(nèi)像素加權(quán)和來(lái)實(shí)現(xiàn)濾波，特別典型的是，同一模式的權(quán)重因子可以作用在每一個(gè)窗口內(nèi)，也就意味著線性濾波器是空間不變的，這樣就可以使用卷積模板來(lái)實(shí)現(xiàn)濾波。高斯濾波函數(shù)[1-2]為其中x，y是坐標(biāo)值，σ是正態(tài)分布的標(biāo)準(zhǔn)偏差。這個(gè)公式生成的曲面的等高線是從中心開(kāi)始呈正態(tài)分布的同心圓。濾波時(shí)中心像素的值有最大的高斯分布值，所以有最大的權(quán)重，周圍像素隨著距離中心像素越來(lái)越遠(yuǎn)，其權(quán)重也越來(lái)越小。這樣進(jìn)行處理比其他的均衡濾波器更好地保留了邊緣效果。

近十年來(lái)，計(jì)算機(jī)圖形處理器（Graphics Processing Unit，GPU）由原本只是處理計(jì)算機(jī)圖形的專用設(shè)備發(fā)展成為高并行度、多線程、多核的處理器。目前主流GPU的運(yùn)算能力已超過(guò)主流通用CPU，從發(fā)展趨勢(shì)上來(lái)看將來(lái)差距會(huì)越拉越大。GPU卓越的性能對(duì)開(kāi)發(fā)GPGPU（使用GPU進(jìn)行通用計(jì)算）非常具有吸引力。統(tǒng)一計(jì)算設(shè)備架構(gòu)（Compute Unified Device Architecture，CUDA）是NVIDIA公司推出的一種通用的GPU編程模型[3]，它使得通用計(jì)算過(guò)程不必再像過(guò)去的GPGPU那樣必須將計(jì)算映射到圖形API中，也就不需要學(xué)習(xí)艱澀的顯示芯片指令，因此開(kāi)發(fā)門檻大大降低。目前，基于GPU的通用計(jì)算成為高性能計(jì)算領(lǐng)域的熱點(diǎn)研究課題[4]。某些圖像處理對(duì)實(shí)時(shí)性要求很高，如何將GPU所提供的強(qiáng)大并行計(jì)算能力運(yùn)用到實(shí)時(shí)圖像處理中，成為一個(gè)很有意義的研究課題。

目前，基于CUDA的GPU加速圖像濾波的研究有不少成果。

盒型濾波器（boxFilter）將輸入圖像每個(gè)像素點(diǎn)用周圍像素平均值代替得到輸出圖像。文獻(xiàn)[5]給出使用滑窗技術(shù)實(shí)現(xiàn)快速盒型濾波器的方法。中值濾波是平滑與去噪應(yīng)用中常用的技術(shù)，文獻(xiàn)[6]利用CUDA實(shí)現(xiàn)了無(wú)分支向量化中值（BVM）濾波器。實(shí)驗(yàn)表明，使用該方法消除環(huán)狀偽影時(shí)處理速度是優(yōu)化了的基于CPU的方法的3.7倍。文獻(xiàn)[7]介紹了兩種基于CUDA加速的圖像濾波方法，其中直觀的實(shí)現(xiàn)方法是一個(gè)線程將相應(yīng)位置處與濾波核大小相同的圖像元素拿出來(lái)和濾波核做乘加計(jì)算，產(chǎn)生一個(gè)輸出像素，多個(gè)線程并行處理并輸出該圖像塊的濾波結(jié)果。另外一種是分離濾波器方法，將二維濾波分離為兩次一維濾波，即分兩趟處理。Nvidia提供了CUDA加速的CUFFT庫(kù)[8]，利用該庫(kù)可方便高效地實(shí)現(xiàn)頻域?yàn)V波[9]。此外，Nvidia的SDK提供了CUDA實(shí)現(xiàn)遞歸高斯濾波器的方法[10]。

這些CUDA加速的圖像濾波方法有的描述不清楚，也沒(méi)有人對(duì)這些算法的性能進(jìn)行詳盡的比較，這樣就給理解及應(yīng)用中如何選擇這些算法帶來(lái)了困難。這些算法應(yīng)用場(chǎng)合不完全相同，其中很多算法可用于高斯濾波。研究CUDA加速的高斯圖像濾波有很好的代表性。本文描述了幾種CUDA加速的圖像高斯濾波算法，包括直觀的實(shí)現(xiàn)方式、使用共享內(nèi)存的分離濾波器方法、使用紋理內(nèi)存的分離濾波器方法、基于CUFFT的卷積濾波以及遞歸高斯濾波器。強(qiáng)調(diào)了這些算法的核心思想，比較了它們的時(shí)間復(fù)雜度，通過(guò)實(shí)驗(yàn)對(duì)它們的性能進(jìn)行了分析。本研究成果對(duì)于理解及選擇合適的CUDA加速的圖像高斯濾波算法有較高的參考價(jià)值。

2 幾種基于CUDA的圖像高斯濾波方法

2.1 直觀的實(shí)現(xiàn)方法

總體上，該方法使用高斯函數(shù)生成高斯核，然后使用該核對(duì)圖像做卷積。離散的高斯卷積核 H的計(jì)算方法其中 σ為方差，k為核矩陣的半徑。比如：取k=1，σ=1，即可得到3×3的高斯卷積核如下：

進(jìn)行加權(quán)濾波，權(quán)系數(shù)之和必須為1，對(duì)上面所求出的高斯濾波核函數(shù)進(jìn)行歸一化后：

不同要求的其他高斯卷積核可以同法得到。

直觀的CUDA加速高斯濾波的具體實(shí)現(xiàn)方法是將高斯濾波核存儲(chǔ)在常量?jī)?nèi)存，將圖像劃分為多個(gè)塊[7]，一個(gè)線程塊加載圖像的一塊到共享內(nèi)存，對(duì)該塊濾波。該線程塊中，一個(gè)線程將共享內(nèi)存相應(yīng)位置處與濾波核大小相同的圖像元素拿出來(lái)和濾波核做乘加計(jì)算，產(chǎn)生一個(gè)輸出像素；多個(gè)線程并行處理并輸出該圖像塊的濾波結(jié)果。調(diào)入共享內(nèi)存的圖像塊邊緣像素在與濾波核做計(jì)算時(shí)，會(huì)依賴不在共享內(nèi)存中的像素。因此，在要做濾波計(jì)算的圖像塊周圍，需要一圈濾波器半徑寬度的像素“圍裙”。這樣每個(gè)線程塊既要加載需要濾波的像素，又要額外加載圍裙像素到共享內(nèi)存。如果每個(gè)像素加載到共享內(nèi)存都需要用一個(gè)線程，那么加載圍裙像素的線程在濾波器計(jì)算期間會(huì)處于空閑狀態(tài)。如果濾波器的半徑較大，空閑線程的比例就較大，這浪費(fèi)了GPU大量的并行計(jì)算能力，顯然，這種方式性能不會(huì)好。

2.2 基于CUDA的分離濾波器

通常二維卷積濾波器要為每個(gè)輸出像素做W×H次乘法；W和H是濾波器核的寬度和高度。分離濾波器是將二維濾波器表達(dá)為兩個(gè)一維濾波器，一個(gè)作用于圖像的行，另一個(gè)作用于圖像的列。分離濾波器只需為每個(gè)輸出像素做W+N次乘法。比如對(duì)圖像做二維濾波，與先做一維濾波，再做[- 1 0 1]一維濾波的效果一樣。

分離濾波器在執(zhí)行時(shí)可提供更多的靈活性，另外也減少了計(jì)算復(fù)雜度和計(jì)算每個(gè)像素點(diǎn)需要的帶寬。只有一部分二維濾波器可分離為兩個(gè)一維濾波器，高斯濾波器核是可分離的。

2.2.1 使用共享內(nèi)存實(shí)現(xiàn)分離濾波

基于CUDA的分離濾波器方法是將二維濾波分離為兩次一維濾波[7]，即分兩趟處理，每趟一維濾波將需要處理的圖像塊首先加載到共享內(nèi)存后再做濾波計(jì)算。由于一維濾波計(jì)算邊緣像素時(shí)需要的“圍裙”像素較少，因此可以減少不必要的數(shù)據(jù)加載量。

實(shí)際上，程序在處理時(shí)，更多的是受到線程塊大小的限制，而不是共享內(nèi)存的大小。為了達(dá)到更高的效率，每個(gè)線程必須處理多個(gè)輸出像素。分離濾波的方法能更靈活地增加單個(gè)線程塊所處理圖像塊的寬度，這會(huì)帶來(lái)顯著的性能改進(jìn)。

設(shè)備內(nèi)存的帶寬理論上遠(yuǎn)遠(yuǎn)高于主機(jī)內(nèi)存。為了實(shí)現(xiàn)高內(nèi)存吞吐量，GPU會(huì)嘗試將多個(gè)線程訪問(wèn)合并為單個(gè)內(nèi)存事務(wù)，這稱為內(nèi)存合并優(yōu)化（Coalescence）。如果一個(gè)warp中所有的線程（32個(gè)線程）同時(shí)讀取連續(xù)的字，GPU可以用優(yōu)化方式讀取內(nèi)存，加塊內(nèi)存訪問(wèn)速度。如果讀取32個(gè)隨機(jī)內(nèi)存地址，讀取性能會(huì)降低很多，只能達(dá)到設(shè)備內(nèi)存總帶寬的一小部分。為了使0號(hào)線程總是從正確對(duì)齊的地址讀取從而滿足所有warp內(nèi)存對(duì)齊的約束，行濾波器中使用的方法是在要處理的一條數(shù)據(jù)的邊緣前端增加線程。這看起來(lái)像是在浪費(fèi)線程，但這沒(méi)太大關(guān)系，因?yàn)楫?dāng)數(shù)據(jù)塊足夠大時(shí)，圍裙像素相比輸出像素是很少的。在列濾波時(shí)，圍裙并不影響合并的對(duì)齊要求，只要圖像條的寬度是16的倍數(shù)。

分離的行濾波和列濾波對(duì)應(yīng)的CUDA內(nèi)核都分為兩個(gè)階段。第一階段將數(shù)據(jù)從全局內(nèi)存加載到共享內(nèi)存中，第二階段執(zhí)行濾波，并將結(jié)果寫回到全局內(nèi)存。

行濾波過(guò)程：在加載階段，每個(gè)活動(dòng)的線程加載一個(gè)像素。在計(jì)算階段，每個(gè)線程的循環(huán)次數(shù)是濾波器半徑寬度的兩倍加1，每個(gè)像素乘以存儲(chǔ)在常量?jī)?nèi)存中的相應(yīng)的濾波器核系數(shù)。在半warp中的每個(gè)線程訪問(wèn)相同的常量?jī)?nèi)存地址，因此，沒(méi)有常量?jī)?nèi)存bank沖突開(kāi)銷。此外，連續(xù)的線程始終訪問(wèn)連續(xù)的共享內(nèi)存地址，所以也沒(méi)有共享內(nèi)存bank沖突。

列濾波過(guò)程：列濾波操作很像行濾波。主要的區(qū)別在于線程號(hào)的增加跨過(guò)了濾波區(qū)域。和行濾波一樣，單個(gè)半warp中的線程始終訪問(wèn)不同的bank共享內(nèi)存，但下一個(gè)地址根據(jù)列寬遞增，而不是加1。列濾波加載階段，沒(méi)有行濾波那樣為了合并對(duì)齊而加上的處于非活動(dòng)狀態(tài)的線程，因?yàn)榧俣袑挾仁呛喜⒆x取大小的倍數(shù)。為了減少圍裙像素與輸出像素的比率，希望圖像列盡可能長(zhǎng)，為了盡量利用共享內(nèi)存，就要使圖像列盡可能窄，綜合考慮，使用16列。

在行和列濾波的內(nèi)部循環(huán)中，計(jì)算較少，循環(huán)分支的開(kāi)銷很大，為了提高性能，可以將濾波計(jì)算的循環(huán)展開(kāi)，這樣會(huì)獲得一定的性能提升。

2.2.2 使用CUDA紋理內(nèi)存實(shí)現(xiàn)分離濾波

CUDA紋理內(nèi)存是CUDA內(nèi)存的重要組成部分，在很多圖像處理中都能用到，使用紋理內(nèi)存是加速CUDA程序的有效途徑之一。

在執(zhí)行CUDA程序前，都要把數(shù)據(jù)先從主機(jī)內(nèi)存復(fù)制到設(shè)備內(nèi)存中；一般情況下，主要使用設(shè)備的全局內(nèi)存進(jìn)行存取。但在某些應(yīng)用場(chǎng)合，使用紋理內(nèi)存代替全局內(nèi)存進(jìn)行數(shù)據(jù)存取能有效地提高性能。

相對(duì)于全局內(nèi)存和常量?jī)?nèi)存來(lái)說(shuō)，紋理內(nèi)存主要有以下優(yōu)點(diǎn)：

（1）與常量?jī)?nèi)存類似，紋理內(nèi)存也具有數(shù)據(jù)緩存能力，但容量比常量?jī)?nèi)存的64 kB大很多。

（2）對(duì)于隨機(jī)存取數(shù)據(jù)的效能可能更好。

（3）可以套用filter，使用內(nèi)建的功能來(lái)做內(nèi)插取值，并且可以設(shè)定要用clamping或repeat模式來(lái)處理邊界數(shù)據(jù)（即如果訪問(wèn)數(shù)據(jù)的索引地址在邊界以外，可以選擇取0或重復(fù)邊界數(shù)據(jù)，該特性對(duì)于卷積等操作，編程上方便了很多）。

在CUDA中，使用紋理的大流程分為三步：Host端的紋理建立部分；Device端使用紋理；Host端刪除紋理。

使用CUDA紋理內(nèi)存實(shí)現(xiàn)分離濾波器的具體實(shí)現(xiàn)過(guò)程不是像2.2.1節(jié)介紹的將全局內(nèi)存中的數(shù)據(jù)分塊使用共享內(nèi)存來(lái)加速訪問(wèn)，而是使用CUDA紋理內(nèi)存的方式。首先將存儲(chǔ)在主機(jī)內(nèi)存中的數(shù)據(jù)傳到紋理內(nèi)存，然后執(zhí)行行濾波，將行濾波的結(jié)果再傳到紋理內(nèi)存，最后執(zhí)行列濾波。

由于紋理內(nèi)存的特點(diǎn)，所以在實(shí)現(xiàn)上不需要使用共享內(nèi)存時(shí)要考慮的邊緣像素，也不需要特別考慮全局內(nèi)存訪問(wèn)對(duì)齊，因而實(shí)現(xiàn)代碼相對(duì)簡(jiǎn)單。但從實(shí)驗(yàn)結(jié)果看，性能比使用共享內(nèi)存的分離濾波器要差一些。

2.3 基于CUFFT的卷積濾波

假設(shè)圖像大小為 dataW×dataH，卷積核的大小為kW×kH。傅里葉變換二維卷積濾波分為以下五個(gè)步驟：

（1）要對(duì)圖像和卷積核分別作擴(kuò)展[9，11]，擴(kuò)展后的大小都是 fftW×fftH，要求 fftW≥dataW+kW-1，fftH≥dataH+ kH-1。如果不擴(kuò)展就直接進(jìn)行傅里葉變換和乘積，會(huì)產(chǎn)生折疊誤差（卷繞），不能得到準(zhǔn)確結(jié)果。

卷積核擴(kuò)展：卷積核擴(kuò)展方式如圖1所示，左邊是卷積核，圖中標(biāo)號(hào)為12的位置是空域中卷積操作時(shí)卷積核與圖像數(shù)據(jù)點(diǎn)對(duì)準(zhǔn)的位置，kx和ky是相對(duì)于標(biāo)號(hào)12的值。右圖是卷積核擴(kuò)展后的狀況，中間數(shù)據(jù)補(bǔ)0。

圖1 卷積核擴(kuò)展方式實(shí)例

圖像數(shù)據(jù)擴(kuò)展：卷積計(jì)算會(huì)到圖像邊界以外，圖像邊界以外數(shù)據(jù)或者以0擴(kuò)展或者是將邊界數(shù)據(jù)復(fù)制向外擴(kuò)展。圖像數(shù)據(jù)擴(kuò)展如圖2所示，圖示的是邊界數(shù)據(jù)復(fù)制向外擴(kuò)展的情況，如果是0擴(kuò)展，就將源圖像以外的數(shù)據(jù)填0。

圖2 圖像數(shù)據(jù)的擴(kuò)展

（2）對(duì)擴(kuò)充后的圖像和卷積核分別進(jìn)行傅里葉變換。

（3）傅里葉變換后的結(jié)果相乘。

（4）相乘的結(jié)果進(jìn)行傅里葉逆變換，并取實(shí)部。

（5）將左上部的矩形修剪為原始大小。

CUDA架構(gòu)下基于FFT卷積的方法在實(shí)現(xiàn)時(shí)，在對(duì)圖像和卷積核作擴(kuò)展時(shí)，采用了GPU并行計(jì)算方法，加快了速度。

FFT變換自身具有可“分治”實(shí)現(xiàn)的特點(diǎn)，因此能高效地在GPU平臺(tái)上實(shí)現(xiàn)。CUDA提供了封裝好的FFT實(shí)現(xiàn)庫(kù)，稱為CUFFT[8]，能夠讓使用者輕易地挖掘GPU的強(qiáng)大浮點(diǎn)處理能力，又不用自己去實(shí)現(xiàn)專門的FFT內(nèi)核函數(shù)。CUFFT提供了與CPU上廣泛使用的FFTW庫(kù)相似的接口[12]，CUFFT與FFTW的性能比較請(qǐng)參考文獻(xiàn)[13]。

CUFFT使用步驟：

（1）使用cufftPlan類型的函數(shù)設(shè)置狀態(tài)，可以分別設(shè)置1維、2維、3維和多維狀態(tài)。

（2）使用cufftExec類型的函數(shù)執(zhí)行FFT，可以實(shí)現(xiàn)實(shí)數(shù)、復(fù)數(shù)的FFT變換，支持單精度和雙精度。

（3）使用cufftDestroy函數(shù)釋放GPU資源。

2.4 基于CUDA的遞歸高斯濾波器

高斯窗常用于對(duì)圖像進(jìn)行濾波，但是隨著高斯半徑的增加，時(shí)間消耗會(huì)逐級(jí)增加。如高斯半徑為N時(shí)，計(jì)算每個(gè)輸出采樣點(diǎn)需要的乘法次數(shù)為(2N+1)×(2N+1)，這種情況下，當(dāng)N=100甚至更大時(shí)，計(jì)算量是非常大的。使用CPU計(jì)算時(shí)，即使進(jìn)行SIMD指令集優(yōu)化，很多情況下仍然不能滿足應(yīng)用要求，比如N=100時(shí)，優(yōu)化后的匯編代碼的執(zhí)行時(shí)間也通常在幾百毫秒以上，遠(yuǎn)不能達(dá)到實(shí)時(shí)處理的要求。

高斯窗方法是使用高斯窗口對(duì)準(zhǔn)原理實(shí)現(xiàn)的，屬于FIR（有限長(zhǎng)單位沖激響應(yīng)）型濾波，對(duì)于半徑大于N的像素點(diǎn)，其權(quán)重取為0，即對(duì)當(dāng)前點(diǎn)無(wú)貢獻(xiàn)。然而在實(shí)際情況下，即使在3倍標(biāo)準(zhǔn)差外的像素也應(yīng)該對(duì)中心點(diǎn)有貢獻(xiàn)的，雖然很小。

由于高斯濾波器應(yīng)用普遍，對(duì)它的性能優(yōu)化非常必要，這樣就出現(xiàn)了IIR（有限長(zhǎng)單位沖激響應(yīng)）型的高斯濾波器，也稱為遞歸高斯濾波器（Recursive Gaussian Filter）。

遞歸高斯濾波器的工作過(guò)程[14-15]：

考慮xn是輸入信號(hào)，yn是輸出信號(hào)。高斯濾波器的二階遞歸實(shí)現(xiàn)可以使用如下的因果序列和非因果序列相加來(lái)實(shí)現(xiàn)，其中：

具體處理方式是圖像水平方向處理一遍，垂直方向處理一遍，每遍都有前后兩趟。當(dāng)然，第一遍和第二遍的處理次序也可以顛倒。

第一遍，按水平方向計(jì)算，分為兩趟。

第一趟，從左到右，按照公式（1）展開(kāi)：y1(n)=a0× x(n)+a1×x(n-1)-b1×y1(n-1)-b2×y1(n-2)。

第二趟，從右到左，按照公式（2）展開(kāi)：y2(n)=a2× x(n)+a3×x(n+1)-b1×y2(n+1)-b2×y2(n+2)。

其中，x(n-1)、x(n)和 x(n+1)是第n-1個(gè)、第n個(gè)和n+1個(gè)輸入，y1(n-2)、y1(n-1)和 y1(n)是第一趟的第n-2個(gè)、n-1個(gè)和第n個(gè)輸出，y2(n)、y2(n+1)和 y2(n+2)是第二趟的第n個(gè)、n+1個(gè)和n+2個(gè)輸出。a0、a1、a2、a3、b1、b2為濾波系數(shù)，它們的值可以根據(jù)公式（1）和（2）及k值算出來(lái)。

將兩趟的輸出結(jié)果相加得到第一遍的處理結(jié)果，該結(jié)果作為第二遍的輸入。

第二遍，按垂直方向計(jì)算，分為兩趟。第一趟，從上到下計(jì)算，計(jì)算公式與水平方向的第一趟相同。

第二趟，從下到上計(jì)算。計(jì)算公式與水平方向的第二趟相同。將兩趟的輸出相加，得到最終結(jié)果。

二維無(wú)限濾波器在該遞歸實(shí)現(xiàn)時(shí)，不管高斯半徑為多少，每個(gè)輸出采樣點(diǎn)的計(jì)算量是相同的。每個(gè)輸出元素只需要16次乘法和14次加法，計(jì)算結(jié)果與無(wú)限情況誤差很小，這就是遞歸高斯濾波器的優(yōu)點(diǎn)。

CUDA實(shí)現(xiàn)遞歸高斯濾波器的方法與上述方法基本相同，為了提高CUDA實(shí)現(xiàn)的效率，在具體實(shí)現(xiàn)時(shí)，做了一些并行化工作[10]。首先，第一遍處理時(shí)，每個(gè)進(jìn)程計(jì)算一列數(shù)據(jù)，先從上到下，再?gòu)南碌缴?，這樣的處理可以滿足全局內(nèi)存合并的要求。第一遍處理完后，第二遍本來(lái)需要按行計(jì)算。為了提高程序效率，滿足全局內(nèi)存合并的要求，將第一遍的計(jì)算結(jié)果做矩陣轉(zhuǎn)置，轉(zhuǎn)置后，仍然按照與第一遍相同的程序做處理。第二遍得到的結(jié)果需要再轉(zhuǎn)置，得到最終的結(jié)果。

3 幾種高斯濾波實(shí)現(xiàn)的時(shí)間復(fù)雜度分析

假設(shè)圖像的大小為M×N，濾波核的半徑是r。

直觀實(shí)現(xiàn)方法的時(shí)間復(fù)雜度是O(M×N×r2)，也就是，對(duì)每個(gè)像素而言，復(fù)雜度是O(r2)，與濾波器大小成二次關(guān)系。如果增加濾波器大小，算法會(huì)明顯減慢。

分離濾波器的算法復(fù)雜度是O(M×N×r)，針對(duì)每個(gè)像素，濾波器的算法復(fù)雜度為O(r)。

基于FFT卷積的計(jì)算復(fù)雜性僅取決于填充后的圖像大小，而填充圖像的大小由圖像和卷積核的大小決定。如果填充后的圖像大小為(fftH，fftW)，每個(gè)元素的計(jì)算復(fù)雜度是：

增加卷積核的大小，仍然假設(shè)它比圖像小得多（使用中多半如此），可以看到基于 FFT卷積的計(jì)算復(fù)雜度大約保持不變。由于直觀實(shí)現(xiàn)方法的計(jì)算復(fù)雜度按卷積核中元素的數(shù)目同比增加，因此，對(duì)于“足夠大”的輸入的圖像，從某一卷積核大小開(kāi)始，基于FFT的卷積比直觀實(shí)現(xiàn)方法的性能越來(lái)越好。當(dāng)濾波器大小較小時(shí)，該算法不如使用可分離濾波，而當(dāng)濾波器大小較大(r＞lb(M×N))時(shí)，基于FFT卷積的速度會(huì)比較快。

對(duì)于遞歸高斯濾波器，每個(gè)輸出采樣點(diǎn)的計(jì)算與高斯半徑?jīng)]有關(guān)系，而6個(gè)濾波系數(shù)是高斯半徑的函數(shù)，只計(jì)算一次，這樣，對(duì)高斯半徑為50、100、300等的處理，每個(gè)輸出像素點(diǎn)的計(jì)算量是相同的，都是16次乘法，14次加法，計(jì)算量大幅下降，很多時(shí)候能滿足圖像處理的性能需求，并且質(zhì)量不會(huì)下降。

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

4.1 實(shí)驗(yàn)結(jié)果

實(shí)驗(yàn)是在NVIDIA GeForce GTX9800+上實(shí)現(xiàn)的，該GPU有128個(gè)頻率為1.836 GHz的流處理器，分為16個(gè)SM，896 MB顯存，裝配在CPU是Intel Core 2 Duo E7400，時(shí)鐘頻率為2.8 GHz，內(nèi)存為1 GB的計(jì)算機(jī)上。

基于CUDA的濾波實(shí)驗(yàn)分別做了2.2.1節(jié)的采用共享內(nèi)存的分離濾波器，2.2.2節(jié)采用紋理內(nèi)存的分離濾波器，2.3節(jié)CUFFT卷積以及2.4節(jié)遞歸高斯濾波。為了做對(duì)比，在CPU上也實(shí)現(xiàn)了分離濾波器濾波。為了比較濾波核大小對(duì)算法的影響，這些方法又分別采用了7×7和41×41兩種大小的濾波核，圖像的大小分別為1 024×1 024，2 048× 2 048，3 072×3 072，4 096×4 096。

圖3是濾波核大小為7×7的幾種高斯濾波方法性能比較，圖4是濾波核大小為41×41的幾種高斯濾波方法性能比較，圖5是CPU上使用分離濾波器方法濾波的性能。由于CPU上實(shí)現(xiàn)的濾波耗時(shí)太長(zhǎng)，如果和GPU上基于CUDA的幾種方法放在同一個(gè)圖中比較會(huì)比例失調(diào)，所以它單獨(dú)成圖。

圖3 濾波核大小為7×7的幾種高斯濾波方法性能比較

基于CUDA的幾種濾波方法所消耗的時(shí)間都沒(méi)有包括申請(qǐng)?jiān)O(shè)備內(nèi)存、將數(shù)據(jù)從主機(jī)內(nèi)存拷貝到設(shè)備內(nèi)存以及釋放設(shè)備內(nèi)存的時(shí)間。

4.2 結(jié)論與分析

（1）CPU與GPU上實(shí)現(xiàn)的高斯濾波性能差距巨大。兩種濾波核大小，四種圖像大小的組合情況下，CPU上分離濾波器比GPU上分離濾波器都慢了上百倍。這一方面是因?yàn)榇_實(shí)GPU加速效果明顯，另外也是由于基于CUDA的方法在耗時(shí)上沒(méi)有包括數(shù)據(jù)輔助處理時(shí)間。

（2）GPU上基于CUDA的四種方法中，采用共享內(nèi)存的分離濾波器是最快的，采用紋理內(nèi)存的分離濾波器次之，兩者性能比較接近。CUFFT卷積與循環(huán)高斯濾波方法性能差不多。

分離濾波器算法本身時(shí)間復(fù)雜度就低，再加上CUDA實(shí)現(xiàn)時(shí)滿足了共享內(nèi)存bank和設(shè)備內(nèi)存合并訪問(wèn)的要求，提高了吞吐量，因此性能優(yōu)越。

（3）由于基于CUDA的高斯濾波器可以快速地實(shí)現(xiàn)，完全能夠達(dá)到實(shí)時(shí)顯示的要求，因此，可以采用CUDA與OpenGL結(jié)合的技術(shù)[16]，隨時(shí)動(dòng)態(tài)改變高斯濾波器的參數(shù)，觀察濾波效果。

（4）基于CUDA的濾波實(shí)現(xiàn)方式中，直觀的實(shí)現(xiàn)方法和基于CUFFT卷積的方法適用性最廣，除了可用于高斯濾波器外，還可用于其他濾波器。分離濾波器只能用于二維濾波可分離為兩個(gè)一維濾波的情況，使用場(chǎng)合有限。

圖4 濾波核大小為41×41的幾種高斯濾波方法性能比較

圖5 CPU上采用分離濾波器濾波的性能

[1]Haddad R A，Akansu A N.A class of fast Gaussian binomial filtersforspeech and image processing[J].IEEE Transactions on Acoustics，Speech and Signal Processing，1991，39：723-727.

[2]Nixon M S，Aguado A S.Feature extraction and image processing[M].[S.l.]：Academic Press，2008：88-89.

[3]NVIDIA Corporation.NVIDIA CUDA programming guide version 3.2[EB/OL].[2011-03-27].http：//developer.nvidia.com/cuda.

[4]Hwu Wen-mei.GPU computing gems[M].New York：Morgan Kaufmann，2011.

[5]NVIDIA Corporation.Boxfilter[EB/OL].[2011-03-27].http：// developer.nvidia.com/gpu-computing-sdk/.

[6]Chen Wei.High performance median filtering using commodity graphics hardware[C]//Nuclear Science Symposium Conference Record（NSS/MIC），2009：4142-4147.

[7]Podlozhnyuk V.Image convolution with CUDA[EB/OL]. [2011-03-20].http：//developer.nvidia.com/gpu-computing-sdk/.

[8]NVIDIA Corporation.CUFFT Library documentation[EB/OL]. [2011-03-27].http：//developer.nvidia.com/gpu-computing-sdk/.

[9]Podlozhnyuk V.FFT-based 2D Convolution[EB/OL].[2011-03-20]. http：//developer.nvidia.com/gpu-computing-sdk/.

[10]NVIDIA Corporation.Recursive Gaussian[EB/OL].[2011-03-27]. http：//developer.nvidia.com/gpu-computing-sdk/.

[11]William K.Digital image processing[M].4th ed.Hoboken，New Jersey：Wiley-Interscience，2007：203-247.

[12]Frigo M，Johnson S.The Faster Fourier Transform in the West（FFTW）[EB/OL].[2012-01-21].http：//www.fftw.org/.

[13]Merz H.CUFFT vs FFTW comparison[EB/OL].[2012-01-25]. http：//www.sharcnet.ca/～merz/CUDA_benchFFT/.

[14]Alvarez L，Deriche R，Santana F.Recursivity and PDE’s in image processing[C]//Proceedings of the International ConferenceonPatternRecognition（ICPR 2000），Barcelona，Spain，2000：242-248.

[15]Bharti B.IIR Gaussian blur filter implementation using Intel? advanced vector extensions[EB/OL].[2010-09-17].http：//software.intel.com/sites/default/files/m/d/4/1/d/8/Gaussian_Filter.pdf.

[16]劉進(jìn)鋒，郭雷.CUDA和OpenGL互操作的實(shí)現(xiàn)及分析[J].微型機(jī)與應(yīng)用，2011（23）：40-43.

LIU Jinfeng

School of Mathematics and Computer,Ningxia University,Yinchuan 750021,China

There are some image filtering algorithms based on CUDA,but some of them are not clearly described,and no one to compare the performance of these algorithms,which brings difficulties for understanding and using these algorithms.This paper discusses five different Gaussian image filters based on CUDA,they are naive method,separable share memory method,separable texture memory method,FFT convolution filtering and recursive Gaussian filter.Core ideas are emphasized,time complexities are compared,and performances are analyzed through experiments.

Gaussian filter;separable filter;recursive Gaussian filter;Compute Unified Device Architecture（CUDA）;Graphics Processing Unit（GPU）

目前已有幾種CUDA加速的圖像高斯濾波算法，但這些算法有的描述不清楚，也沒(méi)有人對(duì)它們的性能進(jìn)行詳盡的比較，這給理解及應(yīng)用帶來(lái)了困難。描述了幾種CUDA加速的圖像高斯濾波算法，包括直觀的實(shí)現(xiàn)方式、使用共享內(nèi)存的分離濾波器方法、使用紋理內(nèi)存的分離濾波器方法、基于CUFFT的卷積濾波以及遞歸高斯濾波器。強(qiáng)調(diào)了這些算法的核心思想，比較了它們的時(shí)間復(fù)雜度，通過(guò)實(shí)驗(yàn)對(duì)它們的性能進(jìn)行了分析。

高斯濾波；可分離濾波器；遞歸高斯濾波器；統(tǒng)一計(jì)算設(shè)備架構(gòu)；圖形處理器

A

TP391

10.3778/j.issn.1002-8331.1306-0035

LIU Jinfeng.Comparation of several CUDA accelerated Gaussian filtering algorithms.Computer Engineering and Applications,2013,49（23）：14-18.

寧夏自然科學(xué)基金（No.NZ12163）。

劉進(jìn)鋒（1971—），男，在讀博士，副教授，研究領(lǐng)域?yàn)镚PU通用計(jì)算、圖形學(xué)、圖像處理。E-mail：ljf@sjtu.org

2013-06-07

2013-07-22

1002-8331（2013）23-0014-05

CNKI出版日期：2013-07-29 http://www.cnki.net/kcms/detail/11.2127.TP.20130729.1104.002.html