王璟瑞，高 銳，邱 焓

（重慶真測(cè)科技股份有限公司，重慶 400044）

0 引 言

近年來(lái)，隨著計(jì)算機(jī)斷層成像技術(shù)的成熟，三維圖像的數(shù)字化處理逐漸成為研究熱點(diǎn)?，F(xiàn)有的三維圖像處理方法大多源自二維圖像領(lǐng)域，如形態(tài)學(xué)濾波算法、分水嶺算法、漫水填充算法、k均值聚類(lèi)算法、水平集算法以及Mean Shift算法等，在三維圖像處理上均有較好的運(yùn)用。Mean Shift算法是一種有效的迭代統(tǒng)計(jì)算法。Fukunaga[1]于1975年發(fā)表了一篇描述概率密度梯度函數(shù)估計(jì)相關(guān)理論的文章，提出了Mean Shift這一概念。在該概念的誕生之初，各國(guó)學(xué)者用Mean Shift來(lái)特指均值向量。從2000年開(kāi)始，Comaniciu[2-3]等人把Mean Shift成功運(yùn)用到特征空間分析中，在目標(biāo)跟蹤、圖像平滑和圖像分割方面取得了良好的效果。

Mean Shift分割方法具有良好的性能和無(wú)參化優(yōu)點(diǎn)，但其本身的計(jì)算復(fù)雜度高，屬于一種計(jì)算密集型算法。Mean Shift算法時(shí)間復(fù)雜度為O(kn2)[4]，k為每個(gè)數(shù)據(jù)點(diǎn)的迭代次數(shù)，n為數(shù)據(jù)集中數(shù)據(jù)點(diǎn)個(gè)數(shù)。計(jì)算機(jī)斷層成像技術(shù)重建而成的三維立體圖像，由若干層二維斷層圖像疊加而成，稱(chēng)為三維體數(shù)據(jù)。三維體數(shù)據(jù)的斷層往往多達(dá)成百上千層，數(shù)據(jù)量非常大，具有海量性[5]。Mean Shift算法的速度過(guò)慢，導(dǎo)致在三維體數(shù)據(jù)場(chǎng)景中難以得到實(shí)際應(yīng)用。

GPU強(qiáng)大的并行計(jì)算能力為復(fù)雜科學(xué)計(jì)算工作提供了有力支持，業(yè)界逐漸推出了多種并行計(jì)算開(kāi)發(fā)平臺(tái)。其中，開(kāi)放式計(jì)算語(yǔ)言（Open Computing Language，OpenCL）作為一個(gè)面向異構(gòu)計(jì)算平臺(tái)的通用編程框架[6]，為實(shí)現(xiàn)GPU并行計(jì)算程序的跨平臺(tái)移植提供了解決方案。為此，本文提出了一種基于OpenCL的Mean Shift三維體數(shù)據(jù)分割方法。

1 Mean Shift算法原理

對(duì)于給定的d維空間Rd的n個(gè)樣本點(diǎn)xi，i=1,2,…,n，在x點(diǎn)的Mean Shift向量的形式定義為：

其中，Sh指的是一個(gè)半徑為h的高維球區(qū)域，k表示在這n個(gè)樣本點(diǎn)xi中有k個(gè)點(diǎn)落入Sh區(qū)域中。Sh的定義為：

可以看出，式（1）中(xi-x)是樣本點(diǎn)xi相對(duì)于點(diǎn)x的偏移。那么，Mean Shift向量Mh(x)就是區(qū)域Sh中所有點(diǎn)相對(duì)于圓心的偏移向量求和之后的平均值，指向樣本集概率密度的梯度方向。

Mean Shift算法是一個(gè)迭代的步驟，即先算出當(dāng)前點(diǎn)的偏移均值，將該點(diǎn)移動(dòng)到此偏移均值，然后以此為新的起始點(diǎn)繼續(xù)移動(dòng)，直到滿足最終的條件。該算法的本質(zhì)是以自適應(yīng)的梯度上升搜索極值，通過(guò)反復(fù)迭代搜索樣本空間中的最密集區(qū)域，并計(jì)算收斂值，最后根據(jù)收斂值對(duì)樣本點(diǎn)進(jìn)行分類(lèi)[7]。

2 三維體數(shù)據(jù)的Mean Shift串行算法

在運(yùn)用并行加速技術(shù)前，首先設(shè)計(jì)基于CPU的串行計(jì)算框架中三維體數(shù)據(jù)的Mean Shift算法流程。

三維體數(shù)據(jù)分割時(shí)，首先將三維圖像體素轉(zhuǎn)換成特征空間的樣本點(diǎn)，然后對(duì)樣本點(diǎn)進(jìn)行Mean Shift聚類(lèi)，其中特征空間中的聚類(lèi)對(duì)應(yīng)于圖像空間的分割?？紤]到圖像的空間信息和灰度信息，體素可轉(zhuǎn)換為特征空間中的一個(gè)4維樣本點(diǎn)。4個(gè)分量分別為體素的3個(gè)方向坐標(biāo)值和灰度值。

轉(zhuǎn)換樣本點(diǎn)后，選定初始樣本點(diǎn)x。根據(jù)式（1）計(jì)算樣本點(diǎn)x的迭代Mean Shift向量Mh(x)，每次迭代將樣本點(diǎn)x向量替換為Mh(x)。在這個(gè)過(guò)程中用到了兩個(gè)窗口參數(shù)hd和hg。hd是空間窗口，在分割過(guò)程中用于決定Sh半徑的大小。hd越大，Sh的范圍越大，落入Sh的樣本點(diǎn)越多。hg是灰度窗口，用于決定Sh范圍內(nèi)的樣本點(diǎn)是否參與計(jì)算。設(shè)置合適的hg，若落入Sh的點(diǎn)與當(dāng)前點(diǎn)的灰度值差別大于hg，則直接判斷不為同一類(lèi)，該點(diǎn)不參與計(jì)算。需要注意，hd和hg這兩個(gè)參數(shù)直接影響Mean Shift算法的執(zhí)行時(shí)間。

每次迭代計(jì)算后判斷是否收斂。若未收斂，則繼續(xù)執(zhí)行上述Mh(x)計(jì)算。收斂準(zhǔn)則由用戶根據(jù)實(shí)際情況自行決定，可以為指定的迭代次數(shù)，也可以判斷||Mh(x)-x||是否小于設(shè)定的ε值。本文將兩種方法結(jié)合使用，設(shè)定迭代次數(shù)為5次，ε為4，滿足任一條件則停止迭代。

當(dāng)前樣本點(diǎn)Mh(x)已收斂后，選定下一個(gè)樣本點(diǎn)，同樣執(zhí)行Mean Shift算法，直到所有體素處理完成為止。

以上是本文中Mean Shift串行算法的具體實(shí)現(xiàn)步驟?？梢钥闯觯畲笮阅芷款i在于體素的串行化處理。三維圖像的數(shù)據(jù)量十分龐大，僅常見(jiàn)的256×256×256圖像而言，體素?cái)?shù)量已經(jīng)達(dá)到千萬(wàn)級(jí)別。Mean Shift算法運(yùn)算量過(guò)大，性能問(wèn)題亟待解決。優(yōu)化方式可以從并行化入手，使圖像中的所有數(shù)據(jù)同時(shí)進(jìn)行計(jì)算。

3 OpenCL介紹

OpenCL是一個(gè)由非營(yíng)利性技術(shù)聯(lián)盟Khronos Group進(jìn)行管理的異構(gòu)編程框架，是一種開(kāi)放的工業(yè)標(biāo)準(zhǔn)。它可以在CPU、GPU和其他處理器上進(jìn)行通用目的的并行編程，使軟件開(kāi)發(fā)者可以更加方便高效地使用異構(gòu)處理平臺(tái)。OpenCL平臺(tái)模型如圖1所示，由一個(gè)主機(jī)（Host）和一個(gè)或更多的計(jì)算設(shè)備（Compute Device）組成。每個(gè)計(jì)算設(shè)備由一個(gè)或多個(gè)計(jì)算單元（Compute Unit）組成，每個(gè)計(jì)算單元進(jìn)一步分成一個(gè)或多個(gè)處理單元（Processing Element）[8]。

使用OpenCL技術(shù)發(fā)揮圖形處理器GPU內(nèi)大量處理單元的運(yùn)算能力，可有效實(shí)現(xiàn)Mean Shift算法的并行加速功能。OpenCL由兩部分組成，一部分是用于編寫(xiě)內(nèi)核程序（在OpenCL設(shè)備上運(yùn)行的代碼）的語(yǔ)言，另一部分是定義并控制平臺(tái)的API。使用OpenCL技術(shù)進(jìn)行并行化處理的邏輯過(guò)程如圖2所示。

圖1 OpenCL平臺(tái)模型

圖2 OpenCL邏輯過(guò)程

4 基于OpenCL的Mean Shift并行化實(shí)現(xiàn)

為實(shí)現(xiàn)Mean Shift算法的基于OpenCL技術(shù)的并行加速，本文設(shè)計(jì)的系統(tǒng)裝置如圖3所示。

4.1 創(chuàng)建GPU存儲(chǔ)空間

圖4是OpenCL內(nèi)存模型，其中CPU有自己的主機(jī)內(nèi)存（Host Memory），而在Compute Device上則比較復(fù)雜。它分為3種類(lèi)型：私有內(nèi)存（Private Memory），由每個(gè)工作項(xiàng)獨(dú)享；局部?jī)?nèi)存（Local Memory），在工作組內(nèi)共享；全局/常量?jī)?nèi)存（Global/Constant Memory），對(duì)所有工作組可見(jiàn)[9]。內(nèi)存管理是顯式的，搬移數(shù)據(jù)必須從Host→Global→Local→Private，然后返回。

圖3 系統(tǒng)裝置

圖4 OpenCL內(nèi)存模型

進(jìn)行Mean Shift計(jì)算，需在GPU全局內(nèi)存中創(chuàng)建存儲(chǔ)空間，將待處理的圖像原始數(shù)據(jù)由主機(jī)內(nèi)存復(fù)制到GPU全局內(nèi)存。要?jiǎng)?chuàng)建的存儲(chǔ)空間分為源數(shù)據(jù)和目標(biāo)數(shù)據(jù)兩部分，使用函數(shù)clCreateBuffer()創(chuàng)建。源數(shù)據(jù)性質(zhì)為CL_MEM_COPY_HOST_PTR（創(chuàng)建GPU存儲(chǔ)空間的同時(shí)將數(shù)據(jù)從主機(jī)內(nèi)存復(fù)制到此處），目標(biāo)數(shù)據(jù)創(chuàng)建的存儲(chǔ)空間大小與源數(shù)據(jù)相同，性質(zhì)為CL_MEM_WRITE_ONLY（只可寫(xiě)）。

4.2 執(zhí)行Mean Shift算法的GPU kernel內(nèi)核程序

在創(chuàng)建內(nèi)核時(shí)并不需要確定目標(biāo)設(shè)備，這個(gè)內(nèi)核可以被發(fā)送到上下文中的任何設(shè)備上執(zhí)行，所以將內(nèi)核發(fā)送到設(shè)備上執(zhí)行前需要?jiǎng)?chuàng)建命令隊(duì)列，確定目標(biāo)設(shè)備，保障主機(jī)和設(shè)備之間的通信。確定目標(biāo)設(shè)備后，將內(nèi)核部署到命令隊(duì)列中時(shí)，所有的內(nèi)核會(huì)被發(fā)送到與這個(gè)隊(duì)列相關(guān)聯(lián)的設(shè)備上執(zhí)行。內(nèi)核執(zhí)行命令是唯一可以被發(fā)送到命令隊(duì)列中的命令。一條命令就是一條從主機(jī)發(fā)送來(lái)的消息，通知設(shè)備執(zhí)行某項(xiàng)操作。圖5是主機(jī)向3個(gè)設(shè)備發(fā)送命令的過(guò)程，其中每個(gè)設(shè)備都有自己的命令隊(duì)列。數(shù)據(jù)傳輸操作可能涉及讀、寫(xiě)設(shè)備上的數(shù)據(jù)，但是命令隊(duì)列中的命令卻只有一個(gè)方向，即從主機(jī)到設(shè)備。設(shè)備并不能向主機(jī)發(fā)送命令。

圖5 將命令發(fā)送到設(shè)備上

一個(gè)執(zhí)行內(nèi)核的命令，必須排隊(duì)到命令隊(duì)列，確定GPU為目標(biāo)設(shè)備后，使用函數(shù)clEnqueueNDRange Kernel()發(fā)送命令，執(zhí)行kernel內(nèi)核程序。該函數(shù)采用數(shù)據(jù)并行執(zhí)行模型，描述內(nèi)核執(zhí)行的索引空間，在設(shè)備上進(jìn)行數(shù)據(jù)劃分，需要NDRange維度和工作組大小的信息。主機(jī)程序通過(guò)一個(gè)叫做NDRange的索引空間調(diào)用內(nèi)核，NDRange的意思是“N維范圍”，可以是一維、二維或三維空間，由函數(shù)clEnqueueNDRangeKernel()的第3個(gè)參數(shù)work_dims指定。為了kerel程序中GPU處理單元索引方便，對(duì)于三維圖像一般將該參數(shù)設(shè)為3。劃分?jǐn)?shù)據(jù)定義N維計(jì)算域，在一個(gè)N維域中的每個(gè)獨(dú)立執(zhí)行單元被稱(chēng)為工作項(xiàng)。N維域定義的并行執(zhí)行的工作項(xiàng)總數(shù)等于全局工作大小。多個(gè)工作項(xiàng)可以被組織在一起稱(chēng)為工作組，組中的工作項(xiàng)可以相互通信，可以在組內(nèi)工作項(xiàng)之間同步執(zhí)行，以協(xié)調(diào)內(nèi)存訪問(wèn)。并行執(zhí)行多個(gè)工作組，映射全局工作大小（Global Work Size）到工作組，可以是顯式的或隱式的。

4.3 GPU kernel內(nèi)核程序的實(shí)現(xiàn)

在GPU kernel內(nèi)核程序中實(shí)現(xiàn)Mean Shit的并行運(yùn)算。這一部分的工作由GPU執(zhí)行，是整個(gè)Mean Shift分割算法系統(tǒng)加速的核心。kernel內(nèi)核使用工作項(xiàng)函數(shù)get_global_id()獲取該處理單元在各維度上的全局工作項(xiàng)ID，處理單元根據(jù)工作項(xiàng)ID同步訪問(wèn)GPU全局內(nèi)存中圖像源數(shù)據(jù)的相應(yīng)位置上的體素點(diǎn)，將該體素點(diǎn)及Sh范圍內(nèi)的數(shù)據(jù)復(fù)制到工作項(xiàng)私有內(nèi)存，并轉(zhuǎn)換為4維樣本點(diǎn)數(shù)據(jù)對(duì)該體素點(diǎn)作相應(yīng)的Mean Shift運(yùn)算。同樣地，運(yùn)算結(jié)果也根據(jù)工作項(xiàng)ID存放在目標(biāo)數(shù)據(jù)相應(yīng)位置。對(duì)于尺寸為M×N×H的三維體數(shù)據(jù)，若全局工作項(xiàng)ID的3個(gè)維度的值是i、j、k，那么對(duì)應(yīng)的源數(shù)據(jù)體素點(diǎn)的地址偏移量計(jì)算公式為：

hd窗口內(nèi)的數(shù)據(jù)地址偏移量計(jì)算公式為：

u、v、w的定義域?yàn)閇-hd,hd]，目標(biāo)數(shù)據(jù)對(duì)應(yīng)體素點(diǎn)的地址偏移量與源數(shù)據(jù)相同。OpenCL中可以使用函數(shù)distance()計(jì)算||Mh(x)-x||。本部分的內(nèi)容由各處理單元同時(shí)運(yùn)行完成，算法不需要作原理上的改變，只是由各體素點(diǎn)依次串行處理順序變?yōu)椴⑿羞M(jìn)行。

內(nèi)核程序執(zhí)行完成后，CPU使用函數(shù)clEnqueue ReadBuffer()將上述運(yùn)算結(jié)果從GPU全局內(nèi)存讀回主機(jī)內(nèi)存。

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

以256×256×105的三維體數(shù)據(jù)分別進(jìn)行OpenCL加速前后的測(cè)試，參數(shù)hd和hg分別為10、18，硬件平臺(tái)為CPU i5-4590（3.3 GHz，4核）、內(nèi)存8G、GPU GTX 970（1 664個(gè)流處理器，顯存4G）。經(jīng)測(cè)試，加速前所耗時(shí)間為246.09 s，加速后所耗時(shí)間為6.754 s，加速比為36.44。

分割前后的效果如圖6、圖7所示，三維體數(shù)據(jù)通常用4個(gè)窗口顯示，右下窗口顯示三維立體圖，另外3個(gè)窗口分別顯示3個(gè)方向的平面切片圖。可以看出，經(jīng)Mean Shift分割后的目標(biāo)更加明顯，并能濾除部分噪聲。

另外，本文對(duì)該體數(shù)據(jù)在參數(shù)hd和hg不同取值時(shí)的運(yùn)行效率進(jìn)行測(cè)試。參數(shù)(hd,hg)分別取(5,10)、(8,10)、(8,18)、(10,18)、(10,25)時(shí)， 加 速 前和加速后的執(zhí)行時(shí)間如圖8所示（縱坐標(biāo)取對(duì)數(shù)，單位：s）。通過(guò)測(cè)試數(shù)據(jù)可以發(fā)現(xiàn)：Mean Shift算法的執(zhí)行時(shí)間（加速前和加速后）隨著參數(shù)hd、hg的增加而增加，其中hd的影響尤為顯著；不同參數(shù)條件下，加速比均能保持在30～40的水平，滿足三維圖像處理對(duì)實(shí)時(shí)性的要求。

圖6 分割前圖像

圖7 分割后圖像

圖8 執(zhí)行時(shí)間對(duì)比

6 結(jié) 語(yǔ)

本文基于Mean Shift算法的原理，充分發(fā)揮GPU內(nèi)大量處理單元并行運(yùn)算的特性，極大地提高了三維圖像處理的速度，實(shí)現(xiàn)了三維體數(shù)據(jù)的快速分割算法，且基于OpenCL技術(shù)，開(kāi)發(fā)簡(jiǎn)捷，用戶核心程序改動(dòng)小，具有平臺(tái)無(wú)關(guān)性。