基于圖形處理器的形態學重建系統

何希 吳炎桃 邸臻煒 陳佳

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

關鍵詞：圖形處理器;形態學重建;并行計算;并行堆;并行數據結構

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 引言

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

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

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

1 相關工作

1.1 形態學重建

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

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

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

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

1.2 基于圖形處理器的并行數據結構

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

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

1.3 并行優先隊列

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

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

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

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

本文的主要工作有兩點：

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

2）基于并行堆集群，在圖形處理器上設計和實現了形態學重建系統。

2 基于圖形處理器的形態學重建系統

2.1 系統概述

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

2.2 并行堆

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

與二元堆類似，并行堆也有入堆和出堆操作。出堆操作從緩沖區內取回r個元素并放到根節點上，然后執行出堆調整操作，使并行堆保持堆的屬性。具體的調整方案是把根節點上的r個元素和它的兩個孩子節點n1、n2內的2r元素合并。假設n1里的最小元素比n2里的最小元素要小，那么合并后最大的r個元素存放在根節點，最小的r個元素放在n2節點，其余的r個元素放在n1節點?？梢宰C明只有n2節點及它的子節點需要繼續調整[28]。出堆調整操作會繼續迭代地調整n2節點及其子節點，直到到達了葉子節點。

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

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

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

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

2.3 并行策略

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

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

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

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