肖漢，孫陸鵬，李彩林，周清雷

1.鄭州師范學(xué)院 信息科學(xué)與技術(shù)學(xué)院，鄭州450044

2.山東理工大學(xué) 建筑工程學(xué)院，山東 淄博255000

3.鄭州大學(xué) 計(jì)算機(jī)與人工智能學(xué)院，鄭州450001

隨著計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)字圖像處理系統(tǒng)應(yīng)用在安全監(jiān)控、智能識(shí)別和航天遙感等眾多領(lǐng)域。由于圖像采集容易受到光照以及設(shè)備質(zhì)量等因素的影響，有些圖像會(huì)存在例如，圖像對(duì)比度較低、模糊等質(zhì)量問題。這些質(zhì)量較差的圖像灰度動(dòng)態(tài)范圍窄，不能很好展現(xiàn)圖像的細(xì)節(jié)。圖像增強(qiáng)技術(shù)可以通過對(duì)圖像信息分布規(guī)律的再調(diào)整以改善圖像質(zhì)量。傳統(tǒng)的直方圖均衡方法是一種典型的基于統(tǒng)計(jì)學(xué)的圖像增強(qiáng)技術(shù)。該方法數(shù)學(xué)原理簡(jiǎn)單，計(jì)算方便，理論基礎(chǔ)清晰。但是它存在噪聲大、細(xì)節(jié)損失嚴(yán)重等不足，在處理過程中還易出現(xiàn)飽和失真等問題。在充分考慮信息分布與人眼一致性的前提下，各種改進(jìn)的直方圖統(tǒng)計(jì)局部增強(qiáng)方法可以對(duì)圖像信息實(shí)現(xiàn)較好的合理分配。圖像局部增強(qiáng)可以達(dá)到提高圖像細(xì)節(jié)的清晰度，突出圖像中需要的信息，改善圖像的視覺效果，有利于對(duì)圖像做進(jìn)一步的處理。然而，隨著圖像規(guī)模的增加，在不改變算法復(fù)雜度的情況下，算法的時(shí)間成本將會(huì)成倍增長(zhǎng)，已無法滿足大數(shù)據(jù)時(shí)代下快速處理的需求。

近些年來，CPU雖然有了很大的發(fā)展，但還是顯示出了其處理性能在功耗墻、存儲(chǔ)墻、頻率墻和過低的指令級(jí)并行等方面的限制。圖形處理器（graphics processing units，GPU）在數(shù)據(jù)并行上具有強(qiáng)大的浮點(diǎn)運(yùn)算能力和存 儲(chǔ)帶寬。2020 年5 月NVIDIA發(fā)布的基于最新一代GPU 架構(gòu)“Ampere”的Tesla A100，單精度和雙精度浮點(diǎn)計(jì)算能力分別可以達(dá)到19.5 TFLOPS和9.7 TFLOPS，帶寬1 555 GB/s。GPU 作為一種新型的高性能科學(xué)計(jì)算設(shè)備開始進(jìn)入人們的視野。與CPU 相比，GPU 具有突出的性能優(yōu)勢(shì)。

因此，本文將GPU 并行計(jì)算和直方圖統(tǒng)計(jì)算法相結(jié)合，著重研究在CPU+GPU 協(xié)同計(jì)算模型下，局部增強(qiáng)算法如何利用GPU 并行計(jì)算資源完成大像幅圖像的實(shí)時(shí)處理，從而驗(yàn)證在該協(xié)同計(jì)算結(jié)構(gòu)下圖像處理算法對(duì)于大規(guī)模圖像的運(yùn)行效率。該文的主要貢獻(xiàn)是通過分析算法的執(zhí)行特征，根據(jù)圖像數(shù)據(jù)的鄰接結(jié)構(gòu)，設(shè)計(jì)了一種基于統(tǒng)一計(jì)算設(shè)備架構(gòu)（compute unified device architecture，CUDA）的直方圖統(tǒng)計(jì)圖像增強(qiáng)并行算法，并實(shí)現(xiàn)了直方圖統(tǒng)計(jì)算法在CUDA 計(jì)算平臺(tái)上的naive 版本?；贕PU 的直方圖統(tǒng)計(jì)并行算法相比單線程CPU 直方圖統(tǒng)計(jì)算法獲得的加速比達(dá)到了261.35倍。

1 相關(guān)基礎(chǔ)與算法分析

1.1 研究現(xiàn)狀

王化喆等采用基于直方圖統(tǒng)計(jì)的局部圖像增強(qiáng)算法，取得了較好的圖像處理質(zhì)量。Lei等通過對(duì)邊緣可靠度的直方圖統(tǒng)計(jì)，提出了一種快速可靠的二維相位圖展開算法，可以達(dá)到準(zhǔn)實(shí)時(shí)性能。Yang等提出了一種多屬性統(tǒng)計(jì)直方圖用于圖像的自動(dòng)配準(zhǔn)，對(duì)噪聲和變化的網(wǎng)格分辨率具有魯棒性。Chu等利用直方圖統(tǒng)計(jì)方法構(gòu)造了一個(gè)完整的紋理描述符，用于合成孔徑雷達(dá)圖像的分類，在斑點(diǎn)噪聲和極低的信噪比方面更具魯棒性。王勇等通過將亮度控制和直方圖分區(qū)域均衡進(jìn)行結(jié)合，實(shí)現(xiàn)了改進(jìn)的直方圖統(tǒng)計(jì)算法，具有良好的視覺效果。Ye等通過在全局方式的對(duì)比度差異的約束下進(jìn)行局部增強(qiáng)，提出了一種用于在線魚類行為監(jiān)測(cè)的自動(dòng)圖像增強(qiáng)方法，避免了假陽(yáng)性檢測(cè)。Zhu等運(yùn)用連通域的直方圖統(tǒng)計(jì)量來識(shí)別和統(tǒng)計(jì)儲(chǔ)糧中的昆蟲，實(shí)現(xiàn)了昆蟲識(shí)別計(jì)數(shù)方法。周啟雙等利用直方圖原理，實(shí)現(xiàn)了基于標(biāo)準(zhǔn)差改進(jìn)和局部均值的自適應(yīng)直方圖統(tǒng)計(jì)算法，影像局部對(duì)比度增強(qiáng)。He等采用texton 聚類和直方圖統(tǒng)計(jì)等策略，提出了一種新的用于合成孔徑雷達(dá)圖像分類的統(tǒng)計(jì)分布紋理特征。呂新正等提出了一種改進(jìn)的脈沖分選算法，算法采用直方圖統(tǒng)計(jì)和PRI（pulse repetition interval）變換相結(jié)合的方法，采用多DSP（digital signal processor）并行處理的方法提高了運(yùn)算速度。Yang等提出了一種結(jié)合強(qiáng)度梯度和強(qiáng)度統(tǒng)計(jì)直方圖的標(biāo)記邊緣檢測(cè)器，提高了道路標(biāo)記檢測(cè)的平均召回率和準(zhǔn)確率。胡英帥通過將極坐標(biāo)直方圖統(tǒng)計(jì)在CUDA 平臺(tái)上進(jìn)行并行化，Shape Context算法的速度得到很大提高。Karbowiak通過使用圖像直方圖實(shí)現(xiàn)了圖像銳化和更改顏色強(qiáng)度的變換，在GPU 上實(shí)現(xiàn)了自動(dòng)圖像標(biāo)簽算法并提高了效率。Gocho等提出了一種基于GPU 的高效SAR 圖像變化檢測(cè)分析器，用于對(duì)子圖像進(jìn)行直方圖分析以實(shí)現(xiàn)SAR 圖像的變化檢測(cè)，獲得了57.5 倍加速比。裴浩等提出了基于GPU 的直方圖并行算法，算法用于實(shí)現(xiàn)三維空間數(shù)據(jù)距離的計(jì)算，獲得了18 倍的性能提升。陳坤等利用現(xiàn)場(chǎng)可編程門陣列（field programmable gate array，F(xiàn)PGA）設(shè)計(jì)了圖像直方圖統(tǒng)計(jì)算法，提高了數(shù)據(jù)處理能力和處理速度。

目前國(guó)內(nèi)外大多數(shù)研究工作集中在將傳統(tǒng)的或改進(jìn)的直方圖統(tǒng)計(jì)算法應(yīng)用到各種專業(yè)領(lǐng)域中，從而改善了圖像增強(qiáng)效果，增強(qiáng)了圖像對(duì)比度。有部分研究成果利用GPU實(shí)現(xiàn)了直方圖統(tǒng)計(jì)算法，在應(yīng)用系統(tǒng)中取得了幾十倍的性能提高。還有部分學(xué)者基于FPGA平臺(tái)提出了直方圖統(tǒng)計(jì)并行算法，提高了算法處理速度?？傮w來看，在提高直方圖統(tǒng)計(jì)局部增強(qiáng)算法運(yùn)算效率方面的研究成果很少，速度的提升也不高。

1.2 CUDA 并行計(jì)算平臺(tái)

CUDA 的應(yīng)用系統(tǒng)是GPU 和CPU 的混合代碼系統(tǒng)。在執(zhí)行CUDA 系統(tǒng)時(shí)，主機(jī)端執(zhí)行的二進(jìn)制代碼在調(diào)用核函數(shù)時(shí)需要將設(shè)備端代碼通過CUDA API傳給設(shè)備端。GPU 傳給CUDA API 的設(shè)備端代碼不一定是二進(jìn)制代碼CUBIN，也可能是運(yùn)行于JIT 動(dòng)態(tài)編譯器上的匯編形式的PTX（parallel thread execution）代碼。最后傳到設(shè)備上的是適合具體GPU 的二進(jìn)制代碼，其中的信息多于PTX 或者CUBIN，這是因?yàn)镃UBIN 或者PTX 只包含了線程塊一級(jí)的信息，而不包括整個(gè)網(wǎng)格的信息。目前，在GPU 上可以運(yùn)行的指令長(zhǎng)度仍然有限制，不能超過兩百萬條PTX指令。GPU 端二進(jìn)制代碼主要包括網(wǎng)格的維度和線程塊的維度，每個(gè)線程塊使用的資源數(shù)量，要運(yùn)行的指令以及常數(shù)存儲(chǔ)器中的數(shù)據(jù)。

執(zhí)行內(nèi)核的線程可以訪問與線程的生存期相同的高性能線程寄存器。每個(gè)線程塊都有一個(gè)共享存儲(chǔ)器（對(duì)該塊的所有線程可見），具有該塊的生存期。只要同一塊的線程沒有沖突訪問，共享存儲(chǔ)器就與寄存器一樣快。所有線程都可以隨機(jī)訪問具有應(yīng)用程序生命周期的全局存儲(chǔ)器，即數(shù)據(jù)駐留在全局存儲(chǔ)器中，用于多次內(nèi)核啟動(dòng)。通常，全局存儲(chǔ)器比訪問寄存器或共享存儲(chǔ)器慢得多。當(dāng)相同halfwarp 的線程同時(shí)訪問全局存儲(chǔ)器時(shí)，如果被訪問的存儲(chǔ)器位于相同的全局存儲(chǔ)器段中，則可以將訪問合并為單個(gè)存儲(chǔ)器事務(wù)。如果不是，則同時(shí)訪問導(dǎo)致多個(gè)順序存儲(chǔ)器傳輸。

1.3 直方圖統(tǒng)計(jì)關(guān)鍵技術(shù)原理

設(shè)為在區(qū)間[0,-1]內(nèi)代表灰度值的一個(gè)離散隨機(jī)變量，并設(shè)(r)表示對(duì)應(yīng)于r值的歸一化直方圖分量。關(guān)于圖像均值的階矩定義為：

其中，是的全局均值，即圖像中像素點(diǎn)的平均灰度。

令(,) 表示給定圖像中任意像素的坐標(biāo)，S表示規(guī)定大小的以(,)為中心的鄰域。該鄰域中像素的均值為：

鄰域像素的方差為：

令(,)表示在圖像分辨率為×的像素點(diǎn)(,)處的灰度值，(,)表示增強(qiáng)后的圖像灰度值。對(duì)于=0,1,…,-1和=0,1,…,-1，有：其中，、、是小于1的正常數(shù)，是灰度放大系數(shù)。

使用模板參與計(jì)算，這在圖像處理的算法中有大量應(yīng)用。這些算法有一個(gè)相同的運(yùn)算特征：在對(duì)某個(gè)像素點(diǎn)進(jìn)行處理時(shí)，需要其模板覆蓋范圍下的鄰域像素點(diǎn)的信息。因此，算法在執(zhí)行時(shí)，需要依據(jù)模板大小遍歷模板范圍內(nèi)的圖像像素。在圖像的邊界區(qū)域運(yùn)用模板遍歷圖像時(shí)，判斷的情況尤為復(fù)雜。此時(shí)，模板覆蓋的圖像區(qū)域可能已經(jīng)越過了圖像邊界。

核函數(shù)的邏輯簡(jiǎn)潔，避免使用邏輯分支控制語句和循環(huán)體較短的循環(huán)語句，這是并行計(jì)算系統(tǒng)性能優(yōu)化的一個(gè)準(zhǔn)則。因此，對(duì)邊界的處理采用傳統(tǒng)的判斷邊界處理法將嚴(yán)重影響系統(tǒng)性能。根據(jù)模板的信息，并行算法采用顯式擴(kuò)充邊界法來處理圖像。圖像在完成邊界擴(kuò)充后，模板只需在原圖像大小范圍內(nèi)遍歷。此時(shí)，在擴(kuò)充圖像中不會(huì)發(fā)生遍歷越界問題。核函數(shù)的處理邏輯可實(shí)現(xiàn)歸一化，無需判斷邊界。直方圖統(tǒng)計(jì)并行算法采用偶拓展的顯式擴(kuò)邊方式進(jìn)行圖像擴(kuò)充，即第0列被復(fù)制到左拓展的第1 列，第1 列被復(fù)制到向左拓展的第2 列，以此類推進(jìn)行四條邊界上的拓展。采用邊緣像素點(diǎn)附近的點(diǎn)對(duì)原圖像邊界進(jìn)行擴(kuò)充，由于像素點(diǎn)鄰近，其相似度很高。并且使用圖像自身的像素點(diǎn)，既很好地解決了邊緣像素點(diǎn)的計(jì)算問題，也沒有違背圖像增強(qiáng)算法的理念。

圖1 圖像拓展示意圖Fig.1 Schematic diagram of image expansion

1.4 算法并行分析

對(duì)于總灰度級(jí)數(shù)為的圖像進(jìn)行的直方圖統(tǒng)計(jì)算法主要包括以下模塊：（1）初始化直方圖統(tǒng)計(jì)圖像。該操作需要對(duì)整體圖像中所有像素點(diǎn)都執(zhí)行一次，其時(shí)間復(fù)雜度為(×)。（2）灰度級(jí)概率估計(jì)計(jì)算。要將圖像中所有像素點(diǎn)的灰度值分配到某一個(gè)灰度級(jí)中進(jìn)行累加，需要執(zhí)行×次，其時(shí)間復(fù)雜度為(×)。然后求出每個(gè)灰度級(jí)的概率估計(jì)，其時(shí)間復(fù)雜度為()。因此，本步驟的時(shí)間復(fù)雜度為(×+)。（3）圖像均值求解。要求計(jì)算整體圖像的灰度均值和濾波窗口覆蓋下的子圖像塊的灰度均值，其時(shí)間復(fù)雜度為()。（4）圖像方差求解。要求計(jì)算整體圖像的灰度方差值和濾波窗口覆蓋下的子圖像塊的灰度方差值，其時(shí)間復(fù)雜度為()。（5）圖像擴(kuò)充。由于需要分別進(jìn)行原圖像左右拓展、上下拓展和保留原圖像區(qū)域，該步驟的時(shí)間復(fù)雜度為(×++×/2+×/2)。（6）局部直方圖統(tǒng)計(jì)計(jì)算。需要對(duì)×個(gè)像素點(diǎn)逐個(gè)判斷是否進(jìn)行灰度增強(qiáng)，其時(shí)間復(fù)雜度為(×)。同時(shí)計(jì)算某個(gè)像素點(diǎn)的局部均值和方差，其時(shí)間復(fù)雜度為(××(/4+2))。由此可得，直方圖統(tǒng)計(jì)算法的時(shí)間復(fù)雜度為(××)，步驟（6）是該算法的核心環(huán)節(jié)。表1 為圖像直方圖統(tǒng)計(jì)算法主要步驟在相同條件下的運(yùn)算時(shí)間。

表1 直方圖統(tǒng)計(jì)算法主要功能運(yùn)行時(shí)間占比Table 1 Proportion of running time of main functions of histogram statistical algorithm

由表1 可見，局部直方圖統(tǒng)計(jì)計(jì)算功能是算法中復(fù)雜度最高、最耗時(shí)的部分，該模塊占用了系統(tǒng)97%以上的資源。顯著降低直方圖統(tǒng)計(jì)算法復(fù)雜度的關(guān)鍵是采用高效的快速局部直方圖統(tǒng)計(jì)方法。因此，對(duì)局部直方圖統(tǒng)計(jì)計(jì)算功能結(jié)構(gòu)的優(yōu)化設(shè)計(jì)對(duì)直方圖統(tǒng)計(jì)并行算法的實(shí)現(xiàn)有著舉足輕重的意義。

1.5 直方圖統(tǒng)計(jì)算法的并行提取

算法的可并行性大小與算法自身存在的任務(wù)依賴性有關(guān)。任務(wù)之間的關(guān)聯(lián)性越低，并行性效果越好，反之則效果越差。通過上面的算法分析，可以得出直方圖統(tǒng)計(jì)算法具有良好的并行性。具體來說，該算法包含三級(jí)并行性：

（1）像素點(diǎn)級(jí)并行：為了對(duì)圖像中包含的隱含特征進(jìn)行局部增強(qiáng)，需要將圖像進(jìn)行直方圖統(tǒng)計(jì)，對(duì)每一幅圖像的每一個(gè)像素點(diǎn)的處理是相互獨(dú)立的，可以并行執(zhí)行。

（2）均值和方差值級(jí)并行：每個(gè)像素點(diǎn)的局部均值和局部方差的計(jì)算是相互獨(dú)立的，可以并行執(zhí)行。

（3）窗口級(jí)并行：每個(gè)局部統(tǒng)計(jì)窗口的處理是相互獨(dú)立的，可以并行處理。

直方圖統(tǒng)計(jì)算法的并行性如圖2 所示。圖2（1）中灰色區(qū)域?yàn)樵即幚韴D像區(qū)域，進(jìn)行局部處理時(shí)的鄰域大小為3×3，白色區(qū)域?yàn)閳D像擴(kuò)充區(qū)域。每個(gè)方塊代表一個(gè)像素點(diǎn)。在局部直方圖統(tǒng)計(jì)計(jì)算步驟中，先對(duì)像素點(diǎn)進(jìn)行處理，即求出子圖像塊1 中的均值和方差，并進(jìn)行是否增強(qiáng)暗區(qū)的判斷等一系列操作。然后依次對(duì),,…,,像素點(diǎn)進(jìn)行處理。子圖像塊1…子圖像塊一系列子圖像塊之間沒有相互依賴性，因此，局部直方圖統(tǒng)計(jì)計(jì)算功能可進(jìn)行并行化處理。

圖2 局部直方圖統(tǒng)計(jì)的并行性示意圖Fig.2 Parallelism diagram of local histogram statistics

由此可見，直方圖統(tǒng)計(jì)算法的計(jì)算量非常大并且具備良好的并行性。因此，算法特性非常適合GPU 計(jì)算平臺(tái)大規(guī)模并行的架構(gòu)特點(diǎn)。同時(shí)，根據(jù)1.1 節(jié)對(duì)有關(guān)研究成果的分析，本文試圖采用CUDA平臺(tái)來解決在異構(gòu)計(jì)算設(shè)備上利用直方圖統(tǒng)計(jì)圖像增強(qiáng)算法對(duì)海量圖像實(shí)時(shí)處理的問題。實(shí)驗(yàn)數(shù)據(jù)表明，基于GPU 的直方圖統(tǒng)計(jì)圖像增強(qiáng)并行算法在保持算法精度的同時(shí)，取得了兩個(gè)數(shù)量級(jí)的性能提高。

2 基于GPU 的直方圖統(tǒng)計(jì)并行映射模型

2.1 并行算法描述

根據(jù)CUDA 架構(gòu)的特點(diǎn)，設(shè)在單指令多線程（single instruction multiple thread，SIMT）模型中啟動(dòng)×個(gè)線程，每一個(gè)線程只負(fù)責(zé)處理一個(gè)像素的局部直方圖統(tǒng)計(jì)計(jì)算，然后把灰度增強(qiáng)值存儲(chǔ)到對(duì)應(yīng)的像素點(diǎn)位置。具體直方圖統(tǒng)計(jì)并行算法形式描述如下。

SIMT 模型上的直方圖統(tǒng)計(jì)并行算法

假設(shè)圖像像幅為×，采用GPU 的多線程對(duì)局部直方圖統(tǒng)計(jì)計(jì)算部分進(jìn)行并行計(jì)算，一個(gè)線程負(fù)責(zé)處理一個(gè)像素點(diǎn)，算法的計(jì)算時(shí)間復(fù)雜度降為()。若在一次內(nèi)核函數(shù)中沒有處理完所有像素點(diǎn)，那么每個(gè)線程需要執(zhí)行內(nèi)核函數(shù)至少(×)/次，其中是活動(dòng)線程的數(shù)量。此時(shí)，直方圖統(tǒng)計(jì)圖像增強(qiáng)并行算法時(shí)間復(fù)雜度為((×)/×)。需要注意的是，一般GPU 中可保持的活動(dòng)線程量很大，即是一個(gè)很大的值。因此，存在直方圖統(tǒng)計(jì)圖像增強(qiáng)并行算法時(shí)間復(fù)雜度((×)/×)?(××)。

2.2 計(jì)算過程

直方圖統(tǒng)計(jì)并行算法的執(zhí)行流程如圖3 所示。

圖3 直方圖統(tǒng)計(jì)并行算法執(zhí)行模式Fig.3 Execution mode of histogram statistical parallel algorithm

直方圖統(tǒng)計(jì)并行運(yùn)算主要實(shí)現(xiàn)步驟如下：

（1）獲取設(shè)備信息和初始化設(shè)備。在當(dāng)前機(jī)器上獲得支持CUDA 的設(shè)備數(shù)量，枚舉設(shè)備信息等，并設(shè)置工作設(shè)備。

（2）圖像預(yù)處理。將圖像從文件讀入主機(jī)內(nèi)存。將文件頭結(jié)構(gòu)、信息頭和位圖信息讀入特定的數(shù)據(jù)結(jié)構(gòu)體中，從而獲得圖像的高度、寬度和大小等數(shù)據(jù)。

（3）計(jì)算級(jí)灰度級(jí)的概率密度函數(shù)，計(jì)算圖像整體統(tǒng)計(jì)指標(biāo)全局均值和方差。

（4）圖像數(shù)據(jù)擴(kuò)展和傳輸。直方圖統(tǒng)計(jì)算法在遍歷計(jì)算到圖像邊界時(shí)，會(huì)發(fā)生圖像數(shù)據(jù)數(shù)組越界尋址的問題，程序中可能會(huì)出現(xiàn)空指針異常，從而導(dǎo)致意想不到的錯(cuò)誤。因此，采用圖像數(shù)據(jù)擴(kuò)展方式避免異常發(fā)生。

系統(tǒng)將存儲(chǔ)在主存中的圖像數(shù)據(jù)傳輸?shù)紾PU 全局存儲(chǔ)器中。由于帶寬的限制，制約系統(tǒng)整體性能提升的最大瓶頸是主存與顯存之間的數(shù)據(jù)通信。因此，并行算法中采用了整體數(shù)據(jù)傳輸，而非多次分塊傳送。根據(jù)圖像像幅大小在設(shè)備上申請(qǐng)分配存儲(chǔ)空間，用于存放由GPU 計(jì)算出來的更新灰度值數(shù)據(jù)。

（5）系統(tǒng)調(diào)度內(nèi)核函數(shù)執(zhí)行直方圖統(tǒng)計(jì)計(jì)算。主要執(zhí)行下面幾個(gè)步驟：

①確定內(nèi)核函數(shù)的執(zhí)行配置。根據(jù)圖像規(guī)模大小，設(shè)置網(wǎng)格維度和線程塊維度。

②分解圖像數(shù)據(jù)。每個(gè)線程塊完成所包含的線程需要計(jì)算的相關(guān)像素點(diǎn)對(duì)應(yīng)點(diǎn)子圖像塊的有關(guān)計(jì)算。

③發(fā)起kernel，進(jìn)行線程間的并行計(jì)算。每個(gè)線程完成與之相應(yīng)的像素點(diǎn)的有關(guān)計(jì)算。

（6）更新數(shù)據(jù)讀出并顯示圖像增強(qiáng)結(jié)果。GPU運(yùn)算結(jié)束后，還需將更新后的像素灰度值從GPU 傳回CPU。

2.3 算法的加速策略

設(shè)計(jì)內(nèi)核函數(shù)kernel<<>>，主要實(shí)現(xiàn)鄰域內(nèi)像素的灰度級(jí)概率密度函數(shù)、局部均值、局部方差和像素局部增強(qiáng)判定。詳細(xì)設(shè)計(jì)如下：

（1）內(nèi)核函數(shù)參數(shù)block 的維度為(T,T)，grid 的維度為((+T-1)/T,(+T-1)/T)。、分別是圖像的高和寬，每個(gè)線程塊處理一個(gè)子圖像塊，塊內(nèi)每個(gè)線程處理一個(gè)像素及其鄰域。系統(tǒng)開啟線程總數(shù)滿足((+T-1)/T×(+T-1)/T×T×T)＞(×)。

（2）每個(gè)線程依據(jù)索引地址計(jì)算相應(yīng)像素鄰域的灰度級(jí)概率密度函數(shù)，根據(jù)鄰域的灰度級(jí)概率密度函數(shù)計(jì)算出像素的局部均值和局部方差。根據(jù)像素點(diǎn)鄰域的局部統(tǒng)計(jì)特性數(shù)據(jù)判斷其是否滿足灰度值局部增強(qiáng)條件，并將相應(yīng)更新的灰度值數(shù)據(jù)寫入主存空間的像素值，如圖4 所示。為了實(shí)現(xiàn)每個(gè)線程對(duì)像素點(diǎn)鄰域的局部統(tǒng)計(jì)特性的計(jì)算，在GPU 設(shè)備函數(shù)中設(shè)計(jì)了局部均值計(jì)算函數(shù)和局部方差計(jì)算函數(shù)。__device__函數(shù)是一種設(shè)備端的子函數(shù)，是由GPU 中一個(gè)線程調(diào)用的函數(shù)。

圖4 直方圖統(tǒng)計(jì)劃分Fig.4 Histogram statistical segmentation

在CUDA 的執(zhí)行模型中，需要滿足以下三方面的限制：（1）一個(gè)流多處理器（streaming multiprocessor，SM）所能并行運(yùn)行的線程總數(shù)的限制；（2）一個(gè)SM 所能并行運(yùn)行的線程塊總數(shù)的限制；（3）一個(gè)線程塊所能并行運(yùn)行的線程總數(shù)的限制。一個(gè)SM 所能并行運(yùn)行的線程總數(shù)是確定的，為了獲得圖像直方圖統(tǒng)計(jì)算法最大的并行計(jì)算效率，合理地設(shè)計(jì)線程塊中所包含的線程數(shù)目是關(guān)鍵。在對(duì)圖像直方圖統(tǒng)計(jì)算法進(jìn)行GPU 實(shí)現(xiàn)的過程中，GPU 中每個(gè)SM最多可以激活的線程數(shù)目為2 048。由于一個(gè)warp一次調(diào)度32 個(gè)線程，為了提高并行計(jì)算效率，一個(gè)線程塊所包含的線程數(shù)應(yīng)該是32 的倍數(shù)。同時(shí)考慮線程塊中寄存器等資源的限制，設(shè)計(jì)線程塊的尺寸為16×16=256，小于線程塊中的最大線程數(shù)1 024 的限制。每個(gè)SM 運(yùn)行2 048/256=8 個(gè)線程塊，小于每個(gè)SM 上最多可以激活的線程塊數(shù)32 的限制。因此，GPU 中的所有線程塊和線程都可以同時(shí)處于并行計(jì)算的狀態(tài)。

在圖像增強(qiáng)中進(jìn)行局部直方圖處理時(shí)，由于設(shè)置的線程數(shù)量大于圖像像素的數(shù)量，每個(gè)線程和處于、平面中的像素成一一對(duì)應(yīng)關(guān)系。以式（6）計(jì)算得到的線程索引作為相應(yīng)像素的、維坐標(biāo)，就可實(shí)現(xiàn)線程坐標(biāo)向像素坐標(biāo)的轉(zhuǎn)換。即在GPU 上進(jìn)行圖像空間的像素坐標(biāo)變換計(jì)算如式（6）所示。

其中，和為某個(gè)線程塊中的線程分別在和方向的索引，.和.為某個(gè)線程塊在網(wǎng)格中的和方向的索引，grid在和方向的維度分別為.=(+T-1)/T，.=(+T-1)/T。

2.4 優(yōu)化設(shè)計(jì)

在CUDA 上實(shí)現(xiàn)直方圖統(tǒng)計(jì)并行算法的關(guān)鍵是減少數(shù)據(jù)在主機(jī)和設(shè)備之間的數(shù)據(jù)傳輸。由于主機(jī)與設(shè)備間的數(shù)據(jù)傳輸帶寬遠(yuǎn)低于設(shè)備和設(shè)備之間的數(shù)據(jù)傳輸顯存帶寬，如果主機(jī)和設(shè)備之間數(shù)據(jù)傳輸次數(shù)過多，就會(huì)進(jìn)入傳輸瓶頸而不能充分發(fā)揮GPU的并行運(yùn)算能力，造成計(jì)算效率下降。為了減少數(shù)據(jù)的傳輸過程，充分利用GPU 的計(jì)算能力，把局部直方圖統(tǒng)計(jì)計(jì)算功能中的全部計(jì)算過程全部映射到GPU 上，減少了中間數(shù)據(jù)在主機(jī)和設(shè)備之間的換入換出操作，并把圖像增強(qiáng)前后的數(shù)據(jù)一次性在主機(jī)和設(shè)備之間進(jìn)行傳輸。這樣就大大減少了傳輸次數(shù)，充分利用了GPU 的計(jì)算能力，提高了計(jì)算密集度。

在CUDA 計(jì)算模型中，系統(tǒng)為了最大限度利用SM 的計(jì)算資源，將會(huì)定期從一個(gè)warp 切換到另一個(gè)warp 進(jìn)行調(diào)度線程。由于同一個(gè)warp 中的32 個(gè)線程是被綁定在一起執(zhí)行同一指令，應(yīng)以32 的整數(shù)倍設(shè)置每個(gè)線程塊中線程數(shù)量，并根據(jù)任務(wù)的具體情況確認(rèn)線程塊每個(gè)維度上的大小。表2 中顯示了圖像像幅大小為4 357×4 872 時(shí)，在線程塊中設(shè)置不同數(shù)量線程時(shí)的系統(tǒng)運(yùn)行時(shí)間。從中可以看出，當(dāng)線程塊維度為16×16 時(shí)，系統(tǒng)性能最優(yōu)。

表2 線程塊維度對(duì)運(yùn)算速度的影響Table 2 Influence of thread block dimension on computing speed

3 性能評(píng)測(cè)和分析

3.1 實(shí)驗(yàn)條件

硬件平臺(tái)配置如表3 所示。

表3 實(shí)驗(yàn)環(huán)境配置Table 3 Experimental environment configuration

軟件平臺(tái)中操作系統(tǒng)為Microsoft Windows 10 64 bit；MATLAB R2008b；集成開發(fā)環(huán)境IDE 為Microsoft Visual Studio 2013；多核處理器支持環(huán)境OpenMP 3.0；環(huán)境支持為CUDA Toolkit 10。

3.2 實(shí)驗(yàn)數(shù)據(jù)與實(shí)驗(yàn)描述

串行算法使用的是C 語言實(shí)現(xiàn)的圖像直方圖統(tǒng)計(jì)圖像增強(qiáng)系統(tǒng)。直方圖統(tǒng)計(jì)并行算法是采用CPU多核多線程和GPU 眾核分別實(shí)現(xiàn)的OpenMP 并行系統(tǒng)和CUDA 并行系統(tǒng)。每種算法共執(zhí)行7 組圖像，像幅分別為547×754、1 246×1 652、2 341×2 684、3 241×3 685、4 357×4 872、6 471×6 832 和8 856×8 476。圖5分別給出了不同像幅大小的直方圖統(tǒng)計(jì)算法串行執(zhí)行與并行執(zhí)行的圖像增強(qiáng)結(jié)果。

圖5 直方圖統(tǒng)計(jì)圖像增強(qiáng)效果Fig.5 Histogram statistics image enhancement effect

表4 給出了直方圖統(tǒng)計(jì)算法CPU 串行和并行執(zhí)行時(shí)間的統(tǒng)計(jì)數(shù)據(jù)。GPU 并行執(zhí)行時(shí)間分為三部分：初始化、計(jì)算和結(jié)束。因此，系統(tǒng)總運(yùn)行時(shí)間可以用這三部分時(shí)間之和來表示，即=++。其中，是指系統(tǒng)初始化、讀取圖像數(shù)據(jù)到內(nèi)存、為數(shù)據(jù)分配顯存空間和將數(shù)據(jù)從主存加載到GPU 等的執(zhí)行時(shí)間，是指GPU 計(jì)算以及CPUGPU 間數(shù)據(jù)傳輸時(shí)間，是指把計(jì)算得到的結(jié)果寫回到指定位置的時(shí)間。

加速比是指CPU 執(zhí)行串行算法的總運(yùn)算時(shí)間與并行系統(tǒng)執(zhí)行并行算法總運(yùn)算時(shí)間的比值。加速比反映了相應(yīng)并行計(jì)算架構(gòu)下并行算法相比CPU 串行算法整體提升效率情況，可用于對(duì)實(shí)際系統(tǒng)速度方面的客觀評(píng)價(jià)，如表5 所示。根據(jù)表4 和加速比的定義，由表4 中第二列數(shù)據(jù)與第三列相應(yīng)數(shù)據(jù)的比值可得表5 中基于OpenMP 的直方圖統(tǒng)計(jì)并行算法的加速比。由表4 中第二列數(shù)據(jù)與第四列相應(yīng)數(shù)據(jù)的比值可得表5 中基于CUDA 的直方圖統(tǒng)計(jì)并行算法的加速比。

表4 直方圖統(tǒng)計(jì)算法運(yùn)行時(shí)間對(duì)比Table 4 Comparison of running time of histogram statistical algorithm

表5 直方圖統(tǒng)計(jì)算法性能對(duì)比Table 5 Performance comparison of histogram statistical algorithm

3.3 驗(yàn)證并行算法結(jié)果

圖像直方圖統(tǒng)計(jì)并行處理的目的是縮短圖像處理時(shí)間，獲得更高局部圖像增強(qiáng)速度。然而，如果以損失圖像質(zhì)量為前提，就沒有了并行化處理的作用。下面以三組實(shí)驗(yàn)圖像數(shù)據(jù)為例進(jìn)行局部圖像增強(qiáng)效果分析。

如圖5 可見，原始鎢絲圖像中部的鎢絲和支架都很清楚并且容易分析。圖像右側(cè)即暗區(qū)域，有部分比較細(xì)微的圖像幾乎看不見。原始圖像經(jīng)過串行/并行直方圖統(tǒng)計(jì)圖像增強(qiáng)處理后，在圖像右側(cè)的鎢絲，甚至鎢絲上的紋路都清晰可見。原始人體脊椎骨折圖像左部存在大片的暗淡區(qū)域，在經(jīng)過串行/并行直方圖統(tǒng)計(jì)圖像增強(qiáng)處理后，在圖像左部顯示出了更多的人骨表面細(xì)節(jié)。原始月球北極圖像右上部存在大片的模糊區(qū)域，在經(jīng)過串行/并行直方圖統(tǒng)計(jì)圖像增強(qiáng)處理后，在圖像右上部顯示出了更多的月球表面細(xì)節(jié)。

以上可以看出，原始圖像經(jīng)過串行/并行直方圖統(tǒng)計(jì)處理后增強(qiáng)了其暗色區(qū)域，亮灰度區(qū)域被完整地保留，串行/并行處理后的圖像展示效果完全一致。

從微觀層面上看，將串行系統(tǒng)和并行系統(tǒng)處理三組實(shí)驗(yàn)圖像數(shù)據(jù)得到的圖像灰度值用直方圖表示，如圖6 所示。經(jīng)過分析對(duì)比，各組串/并行處理后的圖像直方圖均相同。

圖6 串/并行算法直方圖對(duì)比Fig.6 Histogram comparison of series/parallel algorithm

由此可以看出，無論是從宏觀層面還是微觀層面來看，直方圖統(tǒng)計(jì)串行算法和并行算法雖然在方法設(shè)計(jì)和運(yùn)行時(shí)間上不同，但是在圖像處理的結(jié)果上保持一致，從而驗(yàn)證了算法的正確性和可行性。

3.4 討論實(shí)驗(yàn)結(jié)果

在直方圖統(tǒng)計(jì)并行算法的處理過程中，需要原始圖像數(shù)據(jù)的×次存儲(chǔ)器讀取，直方圖統(tǒng)計(jì)圖像增強(qiáng)數(shù)據(jù)的×次存儲(chǔ)器寫入操作。設(shè)圖像數(shù)據(jù)大小×=3 241×3 685，每個(gè)像素值分配存儲(chǔ)空間大小是2 Byte，因此，存儲(chǔ)器存取數(shù)據(jù)總量約為0.048 GB。除以kernel實(shí)際執(zhí)行的時(shí)間0.000 22 s，得到的帶寬數(shù)值是218.18 GB/s，這已經(jīng)接近GeForce GTX 1070 顯示存儲(chǔ)器的256 GB/s 帶寬了。因此，可以看出，基于CUDA 架構(gòu)的直方圖統(tǒng)計(jì)并行算法的效率受限于全局存儲(chǔ)器帶寬。

從表5 可以看出，GPU 加速的直方圖統(tǒng)計(jì)算法性能提升明顯，但GPU并行算法的加速比隨著圖像大小的增加呈現(xiàn)緩慢下降的趨勢(shì)。主要原因是在CUDA并行算法實(shí)現(xiàn)中，CPU 負(fù)責(zé)讀取和輸出圖像數(shù)據(jù)，而這一過程并沒有加速。隨著被處理圖像大小的增加，讀取和輸出圖像數(shù)據(jù)所花費(fèi)的時(shí)間也在增加。因此，CUDA 架構(gòu)下的直方圖統(tǒng)計(jì)并行算法的性能瓶頸是顯存帶寬和主存與顯存之間數(shù)據(jù)傳輸?shù)膸挕?/p>

從圖7 中可以看出，同樣的算法，當(dāng)計(jì)算規(guī)模較小時(shí)，CUDA 并行計(jì)算方式加速效果較為明顯，甚至出現(xiàn)260 倍的加速比效果。如圖像數(shù)據(jù)計(jì)算規(guī)模為3 241×3 685 時(shí)，串行算法計(jì)算時(shí)間為205 714.00 ms，OpenMP 算法計(jì)算時(shí)間為84 655.96 ms，CUDA 并行算法計(jì)算時(shí)間為787.11 ms，CUDA 并行方式計(jì)算時(shí)間遠(yuǎn)小于傳統(tǒng)串行方式和OpenMP 并行方式計(jì)算時(shí)間。當(dāng)圖像規(guī)模較大時(shí)，串行算法耗時(shí)呈現(xiàn)近直線上升趨勢(shì)，OpenMP 算法耗時(shí)也是一種緩慢上升勢(shì)頭，而CUDA 并行算法耗時(shí)只是出現(xiàn)了十分平緩的上升。雖然CUDA 加速效果有所放緩，但仍舊保持了兩個(gè)數(shù)量級(jí)的加速。

圖7 不同架構(gòu)下直方圖統(tǒng)計(jì)算法運(yùn)算時(shí)間對(duì)比Fig.7 Comparison of running time of histogram statistical algorithm under different architectures

產(chǎn)生該現(xiàn)象的原因是，主機(jī)內(nèi)存和GPU 存儲(chǔ)器之間數(shù)據(jù)交互需要一定的時(shí)間開銷，圖像規(guī)模不大時(shí)這部分開銷對(duì)最終的計(jì)算時(shí)間影響較小，而當(dāng)圖像數(shù)據(jù)規(guī)模很大時(shí)，數(shù)據(jù)交互時(shí)間所占比例有所提高，GPU 計(jì)算時(shí)間不足以抵過傳輸延遲帶來的時(shí)間開銷，CUDA 加速的效果才出現(xiàn)減緩的趨勢(shì)。

從圖8 可以看出，當(dāng)圖像像幅大小在1 246×1 652以內(nèi)時(shí)，OpenMP 并行方式的加速比曲線斜率變化不大，而CUDA 并行方式的加速比曲線斜率變得較為陡峭；當(dāng)像幅大小從1 246×1 652 擴(kuò)展到3 241×3 685時(shí)，OpenMP并行方式曲線斜率變化仍不明顯，CUDA并行方式的加速比曲線斜率則出現(xiàn)一種緩慢上升的特點(diǎn)；當(dāng)像幅大小從3 241×3 685 擴(kuò)展到5 361×5 763時(shí)，OpenMP 并行方式曲線斜率仍保持著非常平穩(wěn)的上升態(tài)勢(shì)，而CUDA 并行方式的加速比曲線則呈現(xiàn)出一種緩慢下降的趨勢(shì)。因此，從圖8中可見，隨著像幅規(guī)模的增加，CUDA 曲線斜率的變化都大于OpenMP方式曲線斜率變化。

圖8 直方圖統(tǒng)計(jì)并行算法加速比趨勢(shì)Fig.8 Speedup trend of histogram statistical parallel algorithm

曲線斜率的大小，在一定程度上可反映出數(shù)據(jù)計(jì)算規(guī)模與計(jì)算時(shí)間的關(guān)系，即在相同的計(jì)算規(guī)模下，曲線斜率越大，說明該并行計(jì)算方式的時(shí)間消耗變化越劇烈。當(dāng)曲線正斜率較大時(shí)，數(shù)據(jù)計(jì)算規(guī)模稍微擴(kuò)大，僅能夠?qū)е掠?jì)算時(shí)間的平穩(wěn)增加，這時(shí)計(jì)算規(guī)模與消耗時(shí)間的性價(jià)比極高，形成了計(jì)算效率的高峰期，但擴(kuò)展性較差。

根據(jù)上述分析，可以看出CUDA 并行計(jì)算方式的擴(kuò)展性不如OpenMP 方式，CUDA 方式更容易形成計(jì)算瓶頸。但是，由于CUDA 并行計(jì)算方式帶來的巨大的加速比優(yōu)勢(shì)，在圖像像幅大小不斷擴(kuò)大的趨勢(shì)下，仍然具有OpenMP 傳統(tǒng)并行計(jì)算方式無法企及的加速效果，因此，CUDA 并行計(jì)算方式仍更有優(yōu)勢(shì)。

由圖8 還可以看出，使用在CUDA 架構(gòu)上進(jìn)行GPU 計(jì)算，速度遠(yuǎn)高于基于OpenMP 傳統(tǒng)模式的CPU 計(jì)算，且隨著圖像像幅的不斷增大，這個(gè)速度差距始終保持。這是由于CPU 核心數(shù)的限制，CPU 處于滿載狀態(tài)時(shí)，很難有較大的性能提升，同時(shí)還存在線程創(chuàng)建和調(diào)度的額外開銷。而在一定的計(jì)算量范圍之內(nèi)，GPU 中每個(gè)線程的計(jì)算時(shí)間大致相同，計(jì)算時(shí)間的增加僅僅由于更多的線程和線程塊與硬件之間交互造成的必要時(shí)間開銷。

4 結(jié)束語

本文提出了一種基于CPU+GPU 協(xié)同計(jì)算的直方圖統(tǒng)計(jì)的圖像增強(qiáng)并行算法。首先進(jìn)行核函數(shù)基本配置設(shè)計(jì)和任務(wù)劃分映射，使該并行算法在GPU上可執(zhí)行。然后通過優(yōu)化數(shù)據(jù)傳輸和優(yōu)化線程塊資源配置，進(jìn)一步將并行算法的執(zhí)行效率提高。實(shí)驗(yàn)結(jié)果表明，在NVIDIA GTX 1070和Intel Core i5-7400 CPU組成的實(shí)驗(yàn)平臺(tái)上，使用并行算法進(jìn)行了三組直方圖統(tǒng)計(jì)圖像增強(qiáng)實(shí)驗(yàn)，相比單線程CPU 算法最大加速比達(dá)到了261.35 倍，可以滿足大數(shù)據(jù)圖像的實(shí)時(shí)處理需求。此外，將經(jīng)過本文方法與傳統(tǒng)CPU 串行算法處理后的圖像采用直方圖計(jì)算比較，二者在宏觀和微觀層面結(jié)果均保持一致，驗(yàn)證了方法的正確性。文中采用的并行優(yōu)化方法，對(duì)其他圖像處理算法在GPU 計(jì)算架構(gòu)下的實(shí)現(xiàn)和優(yōu)化也具有很好的指導(dǎo)意義和參考價(jià)值。

本文所完成的工作還有進(jìn)一步優(yōu)化的空間，還有待做更深入的研究：數(shù)字圖像本質(zhì)上就是一個(gè)二維數(shù)組，每個(gè)像素的處理過程是相互獨(dú)立并且完全一樣的計(jì)算過程。因此，可以在GPU 集群上采用MPI 和CUDA 相結(jié)合的技術(shù)，由MPI 完成更大圖像塊之間的并行，由每個(gè)節(jié)點(diǎn)上的GPU 完成圖像塊內(nèi)的并行，通過GPU 集群可以使處理速度更快，爭(zhēng)取在更短的時(shí)間內(nèi)完成更大尺寸的圖像處理工作。