基于圖形處理器的形態(tài)學重建系統(tǒng)

何希 吳炎桃 邸臻煒 陳佳

摘 要：形態(tài)學重建是醫(yī)學圖像處理中非?；A和重要的操作。它根據(jù)掩膜圖像的特征對標記圖像反復進行膨脹操作，直到標記圖像中的像素值不再變化為止。對于傳統(tǒng)基于中央處理器（CPU）的形態(tài)學重建系統(tǒng)計算效率不高的問題，提出了使用圖形處理器（GPU）來加速形態(tài)學重建。首先，設計了適合GPU處理的數(shù)據(jù)結構：并行堆集群;然后，基于并行堆集群，設計和實現(xiàn)了一套基于GPU的形態(tài)學重建系統(tǒng)。實驗結果表明，相比傳統(tǒng)基于CPU的形態(tài)學重建系統(tǒng)，基于GPU的形態(tài)學重建系統(tǒng)可以獲取超過20倍的加速比?；贕PU的形態(tài)學重建系統(tǒng)展示了如何把基于復雜數(shù)據(jù)結構的軟件系統(tǒng)高效地移植到GPU上。

關鍵詞：圖形處理器;形態(tài)學重建;并行計算;并行堆;并行數(shù)據(jù)結構

Abstract： Morphological reconstruction is a fundamental and critical operation in medical image processing， in which dilation operations are repeatedly carried out on the marker image based on the characteristics of mask image， until no change occurs on the pixels of the marker image. Concerning the problem that traditional CPU-based morphological reconstruction system has low computational efficiency， using Graphics Processing Unit （GPU） to quicken the morphological reconstruction was proposed. Firstly， a GPU-friendly data structure： parallel heap cluster was proposed. Then， based on the parallel heap cluster， a GPU-based morphological reconstruction system was designed and implemented. The experimental results show that compared with traditional CPU-based morphological reconstruction system， the proposed GPU-based morphological reconstruction system can achieve speedup ratio over 20 times. The proposed system demonstrates how to efficiently port complex data structure-based software system onto GPU.

Key words： Graphics Processing Unit （GPU）; morphological reconstruction; parallel computing; parallel heap; parallel data structure

0 引言

形態(tài)學重建（morphological reconstruction）是醫(yī)學圖像處理中非常基礎和重要的操作。它根據(jù)掩膜圖像的特征對標記圖像反復進行膨脹操作，直到標記圖像中的像素值不再變化為止。其中的標記圖像（Marker Image）和掩膜圖像（Mask Image）是兩張尺寸相同的高像素醫(yī)學圖像。在標記圖像的膨脹操作中，符合擴散條件的像素點會不斷地向鄰接的像素點傳播其像素值。假設I（e）、I（f）是標記圖像中相鄰的兩個像素點，J（e）是掩膜圖像中與I（e）位置相對應的像素點，那么I（f）對I（e）發(fā)生擴散的條件如下：

如圖1（a）所示，標記圖像中像素點O的像素值（32）比其鄰接像素點C的像素值（17）大，而且C在掩膜圖像中位置相對應的像素點C1的像素值（45）也比C的數(shù)值大，于是O的像素值在膨脹操作中擴散到了它的鄰接像素點C（如圖1（c））。接著，C的像素值又會通過其鄰接像素點進一步地擴散。這種像素值的擴散會不斷迭代地進行下去，直到標記圖像中的像素值不再變化為止。形象地理解形態(tài)學重建，標記圖像中較大的像素值會像“水”一樣四處擴散。如果沒有掩膜圖像，那么當形態(tài)學重建結束時標記圖像里所有的像素值都會變成標記圖像中某個最大的像素值，而掩膜圖像就像“高山”一樣阻擋了標記圖像中像素值任意的傳播，使得像素值擴散只在若干個隔絕的區(qū)域中發(fā)生。由于標記圖像，掩膜圖像的超高分辨率（有些圖像分辨率可以達到51200×51200，102400×102400），并且隨著傳感器技術和掃描設備的提高，圖像的分辨率會越來越高，形態(tài)學重建是一個計算量很大的問題?？紤]到形態(tài)學重建里像素值擴散的迭代特性，傳統(tǒng)的算法使用了隊列來跟蹤像素點的擴散過程，然而隊列是按照像素值的擴散順序而不是像素值的大小來決定下一次迭代中像素值擴散的先后順序，使用隊列會導致頻繁出現(xiàn)對于同一個像素點的多次像素值擴散。形態(tài)學重建的計算量會因此大幅度提高。這個問題促使了在形態(tài)學重建中使用以像素值大小為優(yōu)先級的優(yōu)先隊列，由優(yōu)先隊列先選出具有較大像素值的像素點讓其像素值先擴散以避免較小像素值的擴散，從而可以大幅度減少形態(tài)學重建的計算量。

現(xiàn)代的NVIDIA圖形處理器（Graphics Processing Unit， GPU）是一個可編程的基于眾核架構的處理器。每一個圖形處理器包含若干個流處理器（Streaming Multiprocessors）。每個流處理器上有32個或者更多的計算單元以及存取速度很快但數(shù)量有限的寄存器和共享內(nèi)存（Shared Memory）。同時，每個圖形處理器上還配有所有流處理器都可以訪問的全局內(nèi)存（Global Memory）。運行在圖形處理器上的CUDA（Compute Unified Device Architecture）程序可以啟動成百上千線程同時執(zhí)行任務。這些線程是以線程組的形式組織起來的，每個線程組會由圖形處理器的硬件調(diào)度系統(tǒng)分配到最空閑的流處理器上執(zhí)行。圖形處理器強大的計算能力使得很多傳統(tǒng)問題的計算時間縮短了十幾倍甚至是幾十倍至原來的十幾分之一甚至幾十分之一此處表述不當，應該為縮短至原來的十幾分之一甚至幾十分之一。形態(tài)學重建問題涉及的是高像素圖像以及由此產(chǎn)生的很多的細粒度任務，它非常適合使用圖形處理器來加速其計算過程，而如何設計基于圖形處理器并且適合形態(tài)學重建的并行優(yōu)先隊列決定了能否高效地把形態(tài)學重建以及別的類似應用移植到圖形處理器上進行加速處理。

1 相關工作

1.1 形態(tài)學重建

文獻[1]詳細描述了傳統(tǒng)上形態(tài)學重建的幾種算法，本文簡單地介紹三種。

第一種算法稱為順序重建法（Sequential Reconstruction），這種算法依賴不斷地正向掃描和反向掃描標記圖像，直到在一輪掃描以后再也找不到符合擴散條件的像素點為止。正向掃描指的是從圖像中最左上的像素點開始，到最右下的像素點為止，以從上到下、從左到右的順序掃描圖像中的每一個像素點，而反向掃描恰恰和正向掃描相反，它是從圖像中最右下的像素點開始，到最左上的像素點為止，以從下到上、從右到左的順序掃描圖像中的每一個像素點。在掃描過程中，一旦發(fā)現(xiàn)某個像素點符合擴散的條件，就會立即執(zhí)行像素值擴散操作，將其像素值賦給鄰接像素點。

第二種算法稱為快速混合重建法（Fast Hybrid Reconstruction）。這種算法首先會運行一次正向掃描和一次反向掃描，對標記圖像執(zhí)行初步的膨脹操作，并把有擴散可能的像素點收集起來放進一個先進先出隊列中，然后，對于隊列里的像素點，快速混合重建法會逐個取出來，嘗試把像素點上的像素值向其鄰接像素點擴散，并且把被擴散的鄰接像素點添加到隊列里。上述操作會持續(xù)進行下去，直到隊列為空為止。對比順序重建法，快速混合重建法更加有效率，因為它只需要處理有擴散可能的像素點，而不需要像順序重建法一樣反復進行全局的掃描，但是快速混合重建法也有可以提高的地方，例如其中很多像素值的擴散操作是對同一個像素點執(zhí)行的，因而是重復的計算。

第三種算法稱為下坡過濾器（Downhill Filter）[2]，與第二種算法不一樣的是，它使用了一個隊列，把有擴散可能的像素點按像素值大小進行排序，先處理像素值較大的像素點，避免了像素值較小的像素點的擴散，從而減少了整個形態(tài)學重建的計算量。

1.2 基于圖形處理器的并行數(shù)據(jù)結構

基于中央處理器的數(shù)據(jù)結構已經(jīng)被研究了很多年，但是基于圖形處理器的并行數(shù)據(jù)結構研究還是非常有限?？紤]到本文涉及了圖形處理器上并行數(shù)據(jù)結構，首先對這個領域作一個簡述。

當圖形處理器開始應用于解決一般問題時，研究人員嘗試把各種傳統(tǒng)數(shù)據(jù)結構轉(zhuǎn)變成適合眾核架構的并行數(shù)據(jù)結構。文獻[3]按照廣度優(yōu)先的順序構造KD（請補充KD的英文全稱?；貜停哼@個從原來的英文名字翻譯過來的。原來的名字叫kd-trees或者k-d tree. 國內(nèi)的文獻一般叫KD樹或者kd-樹（在中國知網(wǎng)上搜索"KD樹"可以看到相關文章）. 我覺得就叫KD樹吧。）-樹，并且對大節(jié)點采用了新穎的構造算法以充分利用圖形處理器的計算能力。文獻[4]針對前者在構造KD-樹過程中消耗了過多圖形處理器內(nèi)存的問題提出了按照部分廣度優(yōu)先的順序構造KD-樹，犧牲一部分并發(fā)性以換取內(nèi)存的大量節(jié)省。文獻[5]中提出了一種在圖形處理器上表面積啟發(fā)式構建KD樹的并行方法。其他基于樹的數(shù)據(jù)結構，例如八叉樹（Octree）[6-12]、決策樹[13-14]、層次包圍盒（Bounding Volume Hierarchy）[15]也已經(jīng)在文獻中被討論了如何移植到圖形處理器上。文獻[16]通過把包含指針的操作轉(zhuǎn)化為數(shù)組操作，設計出了一種高效率、適合圖形處理器的跳表結構（Skiplist）。在文獻[17]里，研究者們在完美哈希表（Perfect Hash）[18]和布谷鳥哈希表（Cuckoo Hash）[19]的基礎上設計出了一個適合圖形處理器的并行哈希表。文獻[20]中提出了在圖形處理器中實現(xiàn)布隆過濾器（Bloom Filter）的方案。對于集合，文獻[21]中提出了圖形處理器中新穎的集合數(shù)據(jù)組織方式以及實現(xiàn)集合相交的算法。在圖方面，文獻[22]使用了多層隊列的方式在圖形處理器中加速了廣度優(yōu)先圖的遍歷。文獻[23]解決了如何并行計算最小生成樹的問題。

1.3 并行優(yōu)先隊列

單線程優(yōu)先隊列有不同的實現(xiàn)方式，包括二元堆（Binary Heap）[24]、二項堆（Binomial Heap）[25]、斐波那契堆（Fibonacci Heap）[26]等。一般來說有三種方式來并行化優(yōu)先隊列。

第一種方式是使用多個計算單元來并行化單個元素的入隊列或者出隊列操作[27]，縮短了單個操作的時間。

第二種方式是允許優(yōu)先隊列里同時有多個單線程入隊列和出隊列操作[28]，減少了單個操作的平均時間。

第三種方式則是既允許入隊列操作和出隊列操作同時進行，又使用多個計算單元來并行化入隊列和出隊列操作[29]。

本文的主要工作有兩點：

1）在圖形處理器上設計和實現(xiàn)了并行堆集群。并行堆是并行優(yōu)先隊列在圖形處理器上的實現(xiàn)，而并行堆集群則是并行堆的集合，出于性能考慮而設計，可以理解為并線優(yōu)先隊列在圖形處理器上的近似實現(xiàn)。

2）基于并行堆集群，在圖形處理器上設計和實現(xiàn)了形態(tài)學重建系統(tǒng)。

2 基于圖形處理器的形態(tài)學重建系統(tǒng)

2.1 系統(tǒng)概述

為了充分利用圖形處理器的并行計算能力以縮短形態(tài)學重建的計算時間，本文設計開發(fā)了一套基于圖形處理器的形態(tài)學重建系統(tǒng)，取名為MR_GPU（Morphological Reconstruction_GPU）。MR_GPU設計的一個重要原則是把計算量大、適合使用并行計算的操作放在圖形處理器端處理，而中央處理器端則負責系統(tǒng)的初始化、輸入輸出以及與圖形處理器的協(xié)調(diào)工作。圖2展示的是MR_GPU的流程。r是并行堆節(jié)點可容納的最大像素點個數(shù)。系統(tǒng)的輸入是一張標記圖像和一張掩膜圖像，輸出是一張膨脹后的標記圖像。系統(tǒng)內(nèi)的處理主要分為三個階段：準備階段、建堆階段和迭代膨脹階段。在準備階段，系統(tǒng)會完成初始化和數(shù)據(jù)準備的工作。初始化工作主要包括配置系統(tǒng)參數(shù)，初始化系統(tǒng)中引用的各種外部類庫和在中央處理器和圖形處理器兩端計算并分配存儲空間。數(shù)據(jù)準備工作則包括：1）從標記圖像和掩膜圖像文件里把圖像信息數(shù)據(jù)讀取出來，并存放在二維數(shù)組里;2）對圖像信息進行必要的類型轉(zhuǎn)換和數(shù)據(jù)清洗;3）把圖像信息傳輸?shù)綀D像處理器端的全局內(nèi)存;4）掃描圖像信息，找出符合擴散條件的種子像素點，并且把這些種子像素點均勻地分組。在建堆階段，上一階段留下的每一個像素點組，會被按照像素值進行排序以構造一個最大并行堆。在迭代膨脹階段，每一個最大并行堆都會獨立地、迭代地執(zhí)行膨脹操作，其步驟如算法1所示。細節(jié)會在介紹并行堆時一起討論。

2.2 并行堆

并行堆是專門為圖形處理器設計的優(yōu)先隊列，可以看成是二元堆的升級版本。與二元堆類似，并行堆實際上是一棵完全二叉樹，并且同樣維護堆的屬性，即每個節(jié)點里元素的值都比孩子節(jié)點元素的值要大（最大堆）或者?。ㄗ钚《眩?。在本文中，所討論的并行堆都為最大堆。與二元堆不同的是，并行堆里的節(jié)點可以包含多個元素。實際上，在本文的形態(tài)學重建系統(tǒng)里，并行堆節(jié)點的元素數(shù)是幾十甚至上百，這樣可以方便分配多個線程同時對并行堆進行操作，也符合圖形處理器中通過多線程調(diào)度來解決內(nèi)存存取延遲問題的硬件特性。

與二元堆類似，并行堆也有入堆和出堆操作。出堆操作從緩沖區(qū)內(nèi)取回r個元素并放到根節(jié)點上，然后執(zhí)行出堆調(diào)整操作，使并行堆保持堆的屬性。具體的調(diào)整方案是把根節(jié)點上的r個元素和它的兩個孩子節(jié)點n1、n2內(nèi)的2r元素合并。假設n1里的最小元素比n2里的最小元素要小，那么合并后最大的r個元素存放在根節(jié)點，最小的r個元素放在n2節(jié)點，其余的r個元素放在n1節(jié)點。可以證明只有n2節(jié)點及它的子節(jié)點需要繼續(xù)調(diào)整[28]。出堆調(diào)整操作會繼續(xù)迭代地調(diào)整n2節(jié)點及其子節(jié)點，直到到達了葉子節(jié)點。

入堆操作是要把新的r元素添加到并行堆現(xiàn)有元素之后，并且調(diào)整并行堆使其保持堆的屬性。具體的調(diào)整方案如下：計算一條從根節(jié)點到待插入節(jié)點之間由上往下的入堆路線。從根節(jié)點開始，新的r個元素與根節(jié)點的r個元素合并。較大的r個元素留在根節(jié)點，較小的r個元素繼續(xù)流向入堆路線上的下一個節(jié)點里，然后在下一個節(jié)點上繼續(xù)合并，保留較大的r個元素，讓較小的r個元素繼續(xù)流向再下一個節(jié)點。上述步驟持續(xù)進行，直到較小的元素到達待插入節(jié)點。

傳統(tǒng)的二元堆入堆調(diào)整由下往上進行，而出堆調(diào)整是由上往下進行的。由于在多線程環(huán)境下會存在多個出堆調(diào)整和入堆調(diào)整線程，這種出堆、入堆調(diào)整方向不一致的情況會導致死鎖問題的出現(xiàn)。實際上，并行堆根據(jù)圖形處理器同步的特點，采用了流水線（Pipeline）的并行策略，堆中每一個層次同時都有一組入堆調(diào)整線程和一組出堆調(diào)整線程在運行，因此，在并行堆中，把出堆、入堆調(diào)整的方向都設計為由上到下以避免死鎖的出現(xiàn)。

在形態(tài)學重建系統(tǒng)里需要一個緩沖區(qū)，緩沖區(qū)內(nèi)存放的是最近產(chǎn)生的待擴散的像素點。每次迭代中，并行堆根節(jié)點的像素點會被取出來，與緩沖區(qū)內(nèi)的像素點合并成有序序列，然后最大的r個像素點會被取出進行擴散處理，次大的r個像素點會被放在并行堆的根節(jié)點，然后進行由上到下的出堆調(diào)整操作。另外的r像素點則會沿著計算好的路徑由上而下，最終把經(jīng)過調(diào)整后較小的像素點添加到并行堆的待插入節(jié)點中。

圖3（a）顯示的是一個簡化版本的并行堆例子。它共有11個節(jié)點，每個節(jié)點都有一個編號，最多可以包含2個像素點，緩沖區(qū)內(nèi)有上次迭代收集的4個像素點（62，45，38，35）。在當前迭代中，出堆操作會從并行堆取出根節(jié)點內(nèi)的像素點（57，59）放到緩沖區(qū)內(nèi)，然后對緩沖區(qū)內(nèi)的像素點的值進行排序，然后，如圖3（b）所示，對于像素值最大的2個像素點（62，59），系統(tǒng)會對它們進行擴散處理，并把擴散所涉及的鄰接像素點收集到緩沖區(qū)為下一輪迭代做準備。對于緩沖區(qū)像素值排第3第4的像素點（57，45），系統(tǒng)會把它們放回到并行堆的根節(jié)點，然后從根節(jié)點開始對并行堆進行出堆調(diào)整。至于剩下的像素點（38，35），則會執(zhí)行入堆操作。調(diào)整后具有較小像素值的兩個像素點（35，38）將會被放到節(jié)點12上去。因為待插入節(jié)點已知，系統(tǒng)由此可以計算入堆路線，并由上往下進行入調(diào)整。圖3（c）顯示是入堆、出堆調(diào)整后的并行堆。

2.3 并行策略

行之有效的并行策略是MR_GPU中重要的一環(huán)。在準備階段，分配了數(shù)量眾多的線程來并行掃描標記圖像和掩膜圖像。由于圖形處理器獨特的硬件線程調(diào)度實現(xiàn)，數(shù)量眾多的線程不但不會因為線程調(diào)度而影響性能，反而可以隱藏讀寫內(nèi)存帶來的延遲從而提高整體的性能。在建堆階段和迭代膨脹階段，適用的并行策略可以分為三個層次。

并行策略的第一層次是在圖形處理器集群上。圖形處理器集群上部署了多個圖形處理器，可以切割大圖像文件并分配到每個圖形處理器上，由這些圖形處理器獨立地、并行地進行形態(tài)學重建。同時，圖形處理器間有同步策略，可以在不同圖像文件分割塊間進行同步。

并行策略第二層次是在每個圖形處理器上運行的并行堆集群上，可以配置多個并行堆同時執(zhí)行形態(tài)學重建中的膨脹操作。在最開始的設計中只維護了一個并行堆。一個并行堆可以保證具有較大像素值的像素點首先可以得到擴散的機會，從而避免很多較小像素值不必要的擴散，但是問題在于如果一個圖形處理器中只維護一個并行堆的話，圖形處理器的計算能力遠遠得不到充分的利用，整個形態(tài)學重建的效率并沒有得到最大化?，F(xiàn)在MR_GPU中設計的并行堆集群中各個并行堆是相互獨立的，保證了效率不會因為并行堆之間的相互依賴關系而下降。每個并行堆內(nèi)部也是保證具有較大像素值的像素點首先可以得到擴散的機會從而減少了不必要的計算，而并行堆之間由于是相互獨立，沒有同步的機制，會導致一些重復的像素值擴散操作。然而，相對于采用并行堆集群而獲取的圖形處理器計算能力的充分利用，這些重復計算的代價是可以承受的。

并行策略的第三層次是在并行堆內(nèi)部。在建堆階段，需要對初始化階段找出的首先擴散的像素點進行并行排序從而構造并行堆。在圖形處理器上并行排序是一個已經(jīng)有很多研究者在研究的課題[30-33]?？梢灾苯硬捎靡呀?jīng)優(yōu)化的并行排序方案。在迭代膨脹階段，采用了流水線（pipeline）的并行機制，出堆和入堆調(diào)整線程從根節(jié)點開始，一層一層往下對并行堆進行調(diào)整。當一組出堆調(diào)整和一組入堆調(diào)整線程完成了對根節(jié)點的調(diào)整，開始對并行堆的第二層進行調(diào)整時，系統(tǒng)會啟動另外一組出堆調(diào)整和一組入堆調(diào)整線程執(zhí)行下一輪出堆入堆操作。一般來說，對于一個有n層的并行堆，會有n組出堆調(diào)整線程和n組入堆調(diào)整線程。每一個出堆或入堆調(diào)整線程組內(nèi)會有r個線程，r是并行堆節(jié)點里最大的元素數(shù)，而為節(jié)點里每一個元素分配一個線程無論從實現(xiàn)角度還是效率角度來看都是不錯的選擇。出堆和入堆調(diào)整線程還涉及了并行合并操作和并行收集像素點操作，對于前者，別的研究已經(jīng)提供了解決方案[30，34]。對于并行收集像素點到緩沖區(qū)問題，本文的策略是為每個產(chǎn)生待擴散的像素點的線程分配臨時空間存放待擴散像素點，然后統(tǒng)計像素點的個數(shù)。線程同步以后使用并行掃描的方式統(tǒng)計像素點總的個數(shù)以確定每個線程產(chǎn)生的像素點在緩沖區(qū)內(nèi)的位置，然后讓每個線程分別把各自的像素點從臨時空間拷貝到緩沖區(qū)內(nèi)。