【作 者】王吉喆，嚴(yán)壯志,，溫軍玲

1 上海大學(xué)通信與信息工程學(xué)院，上海市，200444

2 上海大學(xué)上海生物醫(yī)學(xué)工程研究所，上海市，200444

隨著世界人口老齡化的加劇，阿爾茨海默癥（Alzheimer's Disease，AD）患者逐年增加。AD作為一種起病隱匿的進(jìn)行性發(fā)展神經(jīng)系統(tǒng)退行性疾病，目前尚無(wú)特效治療方法。所以只能通過(guò)早診斷，早治療延緩病癥。通過(guò)磁共振圖像（Magnetic Resonance Imaging，MRI）測(cè)量大腦中海馬體的萎縮程度，是診斷AD患者或者輕度認(rèn)知功能障礙患者的最有效的方法之一[1]。所以針對(duì)海馬體的形態(tài)學(xué)分析在臨床診斷和監(jiān)控疾病上具有十分重要的意義。海馬體的輪廓提取作為形態(tài)學(xué)分析的重要步驟，在臨床上以人工逐層分割為主，完成一組海馬圖片的手動(dòng)勾勒的時(shí)間在2個(gè)小時(shí)左右[2]。人工分割具有時(shí)間效率低，可重復(fù)性差等缺點(diǎn)，而且逐層分割失去了數(shù)據(jù)的空間連續(xù)性等三維空間的性質(zhì)。直接進(jìn)行三維分割，又由于海馬體復(fù)雜的邊界形狀，其計(jì)算代價(jià)較高。因此，研究適用于海馬體三維分割的高效算法是十分必要的。

近幾年，學(xué)者們提出很多針對(duì)海馬體的分割算法。一種是基于圖譜的方法，Tong等[3]使用多通道大形變微分同胚度量圖在配準(zhǔn)模板和待分割圖像之間尋找微分同胚對(duì)應(yīng)關(guān)系。Kichang等[4]使用單圖譜的方法對(duì)海馬體進(jìn)行粗分割，并將分割結(jié)果作為下一步使用圖像切割的方法進(jìn)行細(xì)分割的初始化和先驗(yàn)信息。Wu等[5-6]通過(guò)融合不同尺度的圖譜圖像塊，來(lái)增加圖譜與待分割圖像的相似度。Suppa等[7]已經(jīng)開始嘗試將基于圖譜的全自動(dòng)分割、測(cè)量方法用于臨床輔助診斷，并證明該方法在預(yù)測(cè)輕度認(rèn)知功能障礙轉(zhuǎn)化為AD患者的過(guò)程中表現(xiàn)良好。另一種方法是基于偏微分方程的形變模型，但由于其僅依賴于從圖片中提取的信息，存在著一些不足。而先驗(yàn)知識(shí)的引入可以彌補(bǔ)這些不足。最近研究學(xué)者們提出了名為邊界梯度分配（Gradient Distribution on Boundary，GDB）的方法。通過(guò)機(jī)器學(xué)習(xí)建立起一個(gè)區(qū)域權(quán)重圖。這張權(quán)重圖成為GDB[8]，GDB被用作在活動(dòng)輪廓模型演化時(shí)候的先驗(yàn)信息。并且通過(guò)文獻(xiàn)[9]的擴(kuò)展，GDB模型成為了基于活動(dòng)輪廓模型的邊界演化的自適應(yīng)邊界梯度分配（Adaptive Gradient Distribution on Boundary，AGDB）。文獻(xiàn)[10]構(gòu)建了三相GDB模型來(lái)解決迭代次數(shù)與細(xì)節(jié)分割的矛盾。

基于形變模型的方法具有良好的空間擴(kuò)展性，可以將成熟的模型直接擴(kuò)展至三維空間。但求解形變模型中的偏微分方程，是一個(gè)十分耗時(shí)的過(guò)程。格子玻爾茲曼（Lattice Boltzmann，LB）方法作為一種快速求解偏微分方程的方法，其模型經(jīng)過(guò)離散化處理，求解過(guò)程簡(jiǎn)單，易于進(jìn)行編程實(shí)現(xiàn)；本身具有三維模型，可以直接與三維形變模型相結(jié)合，加速計(jì)算求解；計(jì)算求解只與其周圍像素有關(guān)，可以通過(guò)并行化，進(jìn)一步加速計(jì)算，所以便有學(xué)者將其應(yīng)用于圖像處理領(lǐng)域[11-12]。

但在處理三維數(shù)據(jù)時(shí)，若僅僅使用中央處理器（Central Processing Unit，CPU）進(jìn)行串行計(jì)算，不能充分發(fā)揮LB方法的并行性，所以仍存在計(jì)算耗時(shí)過(guò)長(zhǎng)的問(wèn)題。作為專門為了處理圖像顯示而設(shè)計(jì)的圖形處理器（Graphics Processing Unit，GPU），由于其設(shè)計(jì)的初衷就是能夠并行處理大量數(shù)據(jù)，利用GPU進(jìn)行通用并行計(jì)算便應(yīng)運(yùn)而生。CUDA（Compute Unified Device Architecture）作為一種GPU并行計(jì)算開發(fā)環(huán)境，以C語(yǔ)言為基礎(chǔ)，對(duì)C語(yǔ)言進(jìn)行擴(kuò)展，使得開發(fā)人員不需要去學(xué)習(xí)特定的顯示芯片指令或者特殊的結(jié)構(gòu)就能完成GPU并行程序的編寫。

CUDA一經(jīng)推出便受到了很多研究LB學(xué)者的青睞，Balla-Arabé等[13]和Jones等[14]將LB方法與GPU并行計(jì)算相結(jié)合，加速了圖像處理的速度。這充分地顯示出LB方法與GPU并行計(jì)算的匹配性。所以結(jié)合LB方法的算法優(yōu)勢(shì)和GPU的硬件平臺(tái)優(yōu)勢(shì)，可以大大加快處理三維數(shù)據(jù)的速度。本文利用LB方法結(jié)合現(xiàn)有的形變模型與并行計(jì)算平臺(tái)，實(shí)現(xiàn)對(duì)海馬體的三維快速分割。

1 基于先驗(yàn)信息的三維LB模型

1.1 三維LB分割模型

LB方法是近二十年來(lái)發(fā)展起來(lái)的一種系統(tǒng)建模和模擬的方法，由于其對(duì)偏微分方程的快速求解和圖像處理領(lǐng)域相關(guān)方法的發(fā)展，已在圖像去噪、分割、修復(fù)、增強(qiáng)等方面取得了顯著的成就。本文利用三維LB模型，結(jié)合腦部MRI數(shù)據(jù)，完成對(duì)海馬體的快速分割。

圖1 三維格子玻爾茲曼的空間模型Fig.1 Three dimensions lattice Boltzmann model

三維LB的離散空間模型如圖1所示。模型中的每個(gè)節(jié)點(diǎn)都可以看作是三維數(shù)據(jù)中的體素，每個(gè)節(jié)點(diǎn)中的粒子都可沿著矢量方向遷移至相鄰節(jié)點(diǎn)，即LB模型中的演化過(guò)程；隨后對(duì)各個(gè)節(jié)點(diǎn)的粒子進(jìn)行密度重分配，即LB模型中的碰撞過(guò)程。

在有源項(xiàng)（即外力項(xiàng)）存在的情況下，LB模型的一般方程為：

其中fa(r+ea△t, t+△t)表示在t+△t時(shí)刻在節(jié)點(diǎn)r+ea.△t上的粒子分布函數(shù)，ea是粒子運(yùn)動(dòng)的方向，△t是時(shí)間步長(zhǎng)，為松弛因子，(r, t)是平衡態(tài)分布函數(shù)，F(xiàn)a(r, t) 是有源項(xiàng)。對(duì)于三維LB來(lái)說(shuō)，本文設(shè)定的各個(gè)適量方向上的速度ea為

通過(guò)對(duì)等式(1)左邊使用泰勒展開，等式右邊使用Chapmann-Enskog展開，可寫作下面的形式：

其中，F(xiàn)為外力項(xiàng)，D為擴(kuò)散張量

ωα為粒子速度函數(shù)，稱為平衡態(tài)加權(quán)系數(shù)，當(dāng)α=0時(shí)ω=，當(dāng)α=1, 2, ...,6時(shí)ω=。αα

腦部MRI中，海馬體由于其體積較小，且與周圍的腦灰質(zhì)部分邊界不是很明顯，在分割的時(shí)候經(jīng)常出現(xiàn)過(guò)分割的情況。所以本文將邊緣信息、局部信息和先驗(yàn)形狀信息引入至演化方程的外力項(xiàng)，來(lái)彌補(bǔ)分割精度不足的問(wèn)題。

1.2 基于邊緣信息和局部信息的外力項(xiàng)

根據(jù)文獻(xiàn)[15]，可得D=g(│△I│)，則方程可改寫為

F1=αg為分割提供邊緣信息，其中g(shù)(│△I│)是邊緣停止函數(shù)，定義為：

其中，│△Gσ*I │是經(jīng)過(guò)高斯平滑的圖像梯度，Gσ是高斯核函數(shù)，σ是高斯函數(shù)的方差。

F2=β(λ1e1+λ2e2)為分割提供了區(qū)域信息，λ1和λ2是演化曲線的內(nèi)部和外部區(qū)域信息權(quán)重。其中e1，e2定義為：

其中，H(.)是Heaviside函數(shù)，本文采用平滑函數(shù)Hε近似為：

1.3 基于先驗(yàn)形狀信息的外力項(xiàng)

在使用先驗(yàn)形狀分割模型時(shí)應(yīng)考慮到個(gè)體差異性的存在，需要針對(duì)先驗(yàn)信息進(jìn)行配準(zhǔn)和融合，獲得一個(gè)具有普適性的先驗(yàn)形狀信息，隨后將先驗(yàn)信息配準(zhǔn)至待分割圖像，作為分割所使用的外力項(xiàng)之一。本文采用基于最大互信息的配準(zhǔn)方法，利用Powell算法進(jìn)行尋優(yōu)。將所得20組配準(zhǔn)結(jié)果進(jìn)行加和平均，求得融合先驗(yàn)信息，結(jié)果如圖2所示。

先驗(yàn)形狀的描述方法有很多種，其中，距離符號(hào)函數(shù)能夠?qū)⑿螤畹妮喞c主動(dòng)輪廓模型進(jìn)行融合[16]，在變分法的水平集中有廣泛的應(yīng)用。本文采用歐式距離作為相似度衡量指標(biāo)，構(gòu)建的形狀先驗(yàn)外力項(xiàng)，即：

圖2 配準(zhǔn)前后位置的對(duì)比Fig.2 Comparison of the position before and after registration

綜上，本文使用的LB模型的方程為：

其中α，β，μ是有源項(xiàng)的權(quán)重系數(shù)，可根據(jù)不同的圖片進(jìn)行調(diào)整。

2 三維LB分割算法的CUDA實(shí)現(xiàn)

基于CUDA的程序設(shè)計(jì)流程為：首先在CPU端和GPU端進(jìn)行內(nèi)存的分配，初始化CPU端相關(guān)數(shù)據(jù)，并將初始化數(shù)據(jù)復(fù)制到GPU端中；設(shè)置grid和block的大小，并根據(jù)核函數(shù)，在GPU端完成并行計(jì)算過(guò)程；最后，在當(dāng)計(jì)算滿足條件的時(shí)候，計(jì)算結(jié)束，并將數(shù)據(jù)由GPU端復(fù)制回CPU端，完成并行計(jì)算過(guò)程。

由于LB方法的計(jì)算只與其周圍網(wǎng)格有關(guān)，每個(gè)格點(diǎn)的計(jì)算相對(duì)簡(jiǎn)單，因此影響程序性能的主要原因就是數(shù)據(jù)的復(fù)制。所以，應(yīng)盡可能地將計(jì)算所需數(shù)據(jù)經(jīng)由一次復(fù)制過(guò)程復(fù)制到GPU端。其次便是資源的分配。由于并行線程的數(shù)量是不可能無(wú)限增加的，所以，如何根據(jù)具體情況進(jìn)行線程的分配，也是并行計(jì)算的重中之重。

2.1 算法流程

基于上述的模型，本文提出了三維LB分割算法。為了充分發(fā)揮LB方法并行化的特性，使用了GPU平臺(tái)對(duì)運(yùn)算過(guò)程進(jìn)行加速。算法流程圖見(jiàn)圖3。

步驟1：設(shè)置初始化距離函數(shù)φ；

步驟2：設(shè)置參數(shù)α，β，μ，λ1，λ2，根據(jù)式(6)計(jì)算邊緣信息；

步驟3：將邊緣信息、距離函數(shù)、待分割圖像及參數(shù)從CPU端復(fù)制到GPU端；

圖3 三維格子玻爾茲曼分割算法單GPU平臺(tái)實(shí)現(xiàn)流程圖Fig.3 Flow chart of three dimensions lattice Boltzmann model on GPU platform

步驟4：根據(jù)式(7)計(jì)算區(qū)域信息，根據(jù)式(9)計(jì)算先驗(yàn)信息；

步驟5：計(jì)算外力項(xiàng)等其他必要參量；

步驟6：根據(jù)式(10)執(zhí)行碰撞遷移過(guò)程，并且更新距離函數(shù)的值；

步驟7：由CPU判斷是否達(dá)到迭代步數(shù)，若沒(méi)有達(dá)到，則重新啟動(dòng)核函數(shù)，返回步驟4；若已經(jīng)達(dá)到，則將計(jì)算得到的新的距離函數(shù)從GPU端復(fù)制回CPU端，進(jìn)行邊界處理，完成分割過(guò)程。

GPU特化了計(jì)算能力的代價(jià)就是做分支判斷的能力很低，所以一般情況下，分支判斷的部分都交由CPU端運(yùn)行。由于邊界處理涉及到很多判斷過(guò)程，所以在計(jì)算的時(shí)候，將邊界處理放在了完成迭代之后，交由CPU一并完成。

雙GPU平臺(tái)能夠解決因內(nèi)存固化而導(dǎo)致的顯存不足的問(wèn)題，但是又出現(xiàn)了一個(gè)新的問(wèn)題：數(shù)據(jù)集的拆分。三維LB算法存在著層與層之間的數(shù)據(jù)交換，如果只是簡(jiǎn)單的進(jìn)行上下拆分，那么在計(jì)算的過(guò)程中，位于邊界上的兩層將無(wú)法獲得來(lái)自上面一層（下面一層）擴(kuò)散來(lái)的粒子，造成計(jì)算錯(cuò)誤。而對(duì)于迭代求解的方法，錯(cuò)誤是會(huì)累積的，所以需要將GPU1中向下擴(kuò)散的部分搬運(yùn)至GPU2中，完成GPU1中的最后一層向GPU2中第一層擴(kuò)散的過(guò)程。但是CUDA無(wú)法支持GPU對(duì)GPU通信，所以需要CPU作為中介，協(xié)助兩組GPU完成擴(kuò)散過(guò)程。雙GPU平臺(tái)算法流程見(jiàn)圖4。

相較于單GPU，雙GPU在執(zhí)行值步驟6時(shí)，會(huì)分別取GPU1中的最后一層的數(shù)據(jù)向下擴(kuò)散的緩存，和GPU2中的第一層數(shù)據(jù)上擴(kuò)散的數(shù)據(jù)緩存進(jìn)行交換，完成遷移過(guò)程。

3 驗(yàn)證實(shí)驗(yàn)結(jié)果與討論

3.1 驗(yàn)證實(shí)驗(yàn)設(shè)計(jì)

圖4 三維格子玻爾茲曼分割算法雙GPU平臺(tái)實(shí)現(xiàn)流程圖Fig.4 Flow chart of three dimensions lattice Boltzmann model on dual GPU platform

為了驗(yàn)證三維LB分割算法的有效性，本文將該算法與AGDB方法[11]進(jìn)行對(duì)比。AGDB方法將自適應(yīng)邊界分布圖和先驗(yàn)知識(shí)引入主動(dòng)輪廓模型，對(duì)海馬體進(jìn)行三維分割。本文的實(shí)驗(yàn)硬件平臺(tái)CPU為Intel?Core? i7-4790K @ 4.00 GHz，GPU為兩塊NVIDIA GeForce GTX 660，8 GB內(nèi)存。軟件平臺(tái)為Microsoft Visual Studio 2012，CUDA版本為6.5。實(shí)驗(yàn)圖像由 ADNI(The Alzheimer's Disease Neuroimaging Initiative, ADNI)數(shù)據(jù)庫(kù)獲得，EADC-ADNI 工作組有國(guó)際知名的研究海馬體勾畫的專家團(tuán)體，該組織定義了一個(gè)關(guān)于海馬體勾畫的標(biāo)準(zhǔn)協(xié)議[17]，為診斷研究、臨床試驗(yàn)和算法驗(yàn)證工作提供了統(tǒng)一的標(biāo)準(zhǔn)。為了保證實(shí)驗(yàn)具有統(tǒng)計(jì)學(xué)意義，本文選取了20組AD患者的腦部MRI圖像進(jìn)行分割，每組圖像尺寸均為166×256×30，且每組圖像都有相對(duì)應(yīng)的、由專家進(jìn)行勾畫的海馬體的金標(biāo)準(zhǔn)。

具體的實(shí)驗(yàn)步驟如下：

（1）取20組腦部MRI海馬體分割的金標(biāo)準(zhǔn)，進(jìn)行圖像配準(zhǔn)及融合，獲得一組166×256×30的海馬體的先驗(yàn)信息，如圖5所示。

圖5 融合先驗(yàn)信息的三維重建Fig.5 3-D reconstruction of fuse prior information

（2）使用上述先驗(yàn)信息，對(duì)20組腦部MRI圖像進(jìn)行分割，并按照不同情況，設(shè)置α，β，μ，λ1，λ2，進(jìn)行迭代分割。

（3）使用同樣的參數(shù)，分別在CPU平臺(tái)、單GPU平臺(tái)、雙GPU平臺(tái)上進(jìn)行計(jì)算。

（4）針對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行評(píng)估實(shí)驗(yàn)。

3.2 評(píng)估方法

算法的評(píng)估包括對(duì)分割準(zhǔn)確性的評(píng)估和時(shí)間效率的評(píng)估。為了定量比較算法分割的準(zhǔn)確性，本文使用目前在衡量分割精度之中最常采用的Dice相似系數(shù)(Dice Similarity Coefficient，DSC)[18]和豪斯多夫距離(Hausdorあ Distance，HD)進(jìn)行定量比較。時(shí)間效率的評(píng)估標(biāo)準(zhǔn)為算法達(dá)到穩(wěn)定后所需要的時(shí)間和迭代步數(shù)。

DSC可以衡量分割結(jié)果與金標(biāo)準(zhǔn)之間的重疊率，其計(jì)算公式為：

其中：A是分割算法獲得的結(jié)果，M是金標(biāo)準(zhǔn)，運(yùn)算符∩計(jì)算A與M的重疊的像素?cái)?shù)量，│A│+│M│計(jì)算A與M中，全部像素的數(shù)量。當(dāng)兩個(gè)區(qū)域完全重疊的時(shí)候，DSC=1；當(dāng)兩個(gè)區(qū)域完全不重疊的時(shí)候，DSC=0。

HD計(jì)算結(jié)果為兩個(gè)輪廓曲面之間最接近點(diǎn)的最大值，其計(jì)算公式為：

其中，A和B表示兩個(gè)曲面，h(A, B)表示曲面A上的所有點(diǎn)到曲面B上的所有點(diǎn)的距離的最小值，h(B, A)表示曲面B上的所有點(diǎn)到曲面A上的所有點(diǎn)的距離的最小值，這兩者的最大值，即為豪斯多夫距離。HD越小，說(shuō)明曲面A與曲面B越接近。

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

由于迭代次數(shù)會(huì)影響到二者的分割精度，在進(jìn)行三維LB算法分割質(zhì)量的評(píng)估之前，應(yīng)確定最佳迭代方案。所以本文根據(jù)DSC隨時(shí)間變化的情況，來(lái)確定分割的最佳迭代次數(shù)。

圖6記錄了DSC隨迭代次數(shù)變化的曲線。從圖中可以得知，兩種方法的DSC均隨著迭代次數(shù)的增加而增加，并逐漸達(dá)到穩(wěn)定。但是三維LB分割算法在迭代30次左右，DSC就可以達(dá)到一個(gè)相對(duì)穩(wěn)定的值。相較于AGDB方法的70次迭代，三維LB分割算法能夠在迭代次數(shù)更少的前提下獲得一個(gè)較好的分割效果。

圖6 三維格子玻爾茲曼分割算法和AGDB方法與金標(biāo)準(zhǔn)的Dice相關(guān)系數(shù)隨迭代次數(shù)變化曲線Fig.6 The curve of DSC based on standard for three dimensions lattice Boltzmann method and AGDB

使用得到的最佳迭代次數(shù)后，使用本文提出的三維LB分割算法，在20名AD患者的腦部MRI圖像上進(jìn)行分割結(jié)果精度的驗(yàn)證。結(jié)果見(jiàn)圖7。圖中白色部分均為金標(biāo)準(zhǔn)。紅色部分分別為三維LB分割算法和AGDB方法的分割結(jié)果。從圖7可以很直觀地看出三維LB分割算法的分割精度要高于AGDB方法。而表1記錄的兩種方法的平均DSC和HD，也從定量的角度上印證了三維LB分割算法分割的準(zhǔn)確程度，要高于AGDB方法。

圖7 三維格子玻爾茲曼分割算法與AGDB方法分割結(jié)果對(duì)比Fig.7 Segmentation results between 3-D LB and AGDB

表1 三維格子玻爾茲曼分割算法和AGDB方法的分割結(jié)果與金標(biāo)準(zhǔn)的Dice相關(guān)系數(shù)與豪斯多夫距離Tab.1 The result of DSC and HD for three dimensions lattice Boltzmann method and AGDB with standard

三維LB分割算法分割效率的評(píng)估是建立在最優(yōu)資源分配的前提下，線程塊（block）大小對(duì)計(jì)算效率的影響。線程塊是GPU進(jìn)行計(jì)算時(shí)的獲得硬件資源的最小單位。所以以一組AD患者腦部MRI圖像為例，輸入實(shí)際分割時(shí)所使用的參數(shù)，進(jìn)行十次迭代，只計(jì)算GPU端核函數(shù)的計(jì)算時(shí)間，以觀察block大小對(duì)計(jì)算效率的影響，結(jié)果如表2所示。

表2 不同block大小對(duì)同規(guī)模并行計(jì)算所需時(shí)間的影響Tab.2 The eあect of block sizes in parallel computation

從表2中可以看出，在一定范圍內(nèi)，隨著block的變大，計(jì)算時(shí)間在減少，說(shuō)明硬件資源在逐步飽和，而當(dāng)block的大小達(dá)到256時(shí)，硬件資源達(dá)到飽和，計(jì)算時(shí)間不會(huì)進(jìn)一步縮短。這也與多數(shù)程序中推薦的block維度相同。

在上述實(shí)驗(yàn)的基礎(chǔ)上，使用已經(jīng)達(dá)到最優(yōu)化的block參數(shù)，針對(duì)同一組腦部MRI進(jìn)行分割，記錄達(dá)到最大DSC時(shí)，分割所需要的時(shí)間，結(jié)果如表3所示。

表3 不同平臺(tái)下三維格子玻爾茲曼方法和AGDB方法完成分割所需時(shí)間Tab.3 The time required to complete the segment between AGDB and three dimensions lattice Boltzmann method on diあerent platform

從表3中我們可以看出，在完成分割的時(shí)候三維LB分割算法的速度明顯優(yōu)于AGDB方法，而且在使用GPU并行加速后，計(jì)算速度又有了進(jìn)一步的提升。說(shuō)明三維LB分割算法在時(shí)間效率上要優(yōu)于AGDB方法，并且，在GPU并行計(jì)算加入后，計(jì)算效率提升了一個(gè)數(shù)量級(jí)，說(shuō)明了LB方法與GPU并行計(jì)算有著很好的契合性。

但是同時(shí)發(fā)現(xiàn)，在使用雙GPU計(jì)算的情況下，計(jì)算效率反而比單GPU計(jì)算時(shí)的效率要低。這是因?yàn)樵谑褂秒pGPU在計(jì)算的時(shí)候，涉及到層間信息交換的過(guò)程，這個(gè)過(guò)程是交由CPU來(lái)完成的。數(shù)據(jù)的GPU端與CPU端互相復(fù)制和CPU的串行執(zhí)行，是導(dǎo)致雙GPU計(jì)算效率低于單GPU計(jì)算效率的主要原因。雖然存在著這樣的弊端，但是雙GPU計(jì)算的優(yōu)勢(shì)在于更大的顯存和并行計(jì)算規(guī)模，所以在針對(duì)更大數(shù)據(jù)量的計(jì)算時(shí)，雙GPU可以完成單GPU所無(wú)法完成的工作。

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

本文提出了結(jié)合先驗(yàn)信息的LB三維算法，并且證明該模型在分割弱邊緣，低對(duì)比度的海馬體MRI圖像時(shí)，表現(xiàn)良好。并且通過(guò)GPU平臺(tái)的加入計(jì)算，速度得到了明顯的提升。原因分析如下：

（1）先驗(yàn)信息作為約束項(xiàng)的加入，使得模型在針對(duì)圖像邊界模糊，待分割區(qū)域與背景灰度相似等復(fù)雜背景圖像進(jìn)行分割時(shí)，有著較強(qiáng)的準(zhǔn)確性和魯棒性。特別是針對(duì)海馬體體積較小，邊緣模糊，組織對(duì)比度低，結(jié)構(gòu)復(fù)雜的器官，具有較強(qiáng)的分割能力。

（2）LB方法的求解過(guò)程只與該像素點(diǎn)周圍的區(qū)域有關(guān)，自身帶有可并行化的特點(diǎn)，通過(guò)GPU平臺(tái)的使用，大大減少了計(jì)算耗時(shí)。

（3）通過(guò)對(duì)block中線程數(shù)目的優(yōu)化，達(dá)到最佳的硬件資源分配；通過(guò)對(duì)算法流程分析，減少兩端復(fù)制過(guò)程，能夠讓GPU在不被打斷的情況下，完成高密集度的計(jì)算過(guò)程，進(jìn)一步提升了計(jì)算效率。

4 結(jié)論

在腦部MRI中能夠快速、準(zhǔn)確地分割海馬體的結(jié)構(gòu)，對(duì)于AD等疾病的臨床診斷具有重大意義。本文提出的三維LB分割算法，直接在三維空間中對(duì)海馬體進(jìn)行分割，避免了逐層分割帶來(lái)的層間信息損失和三維重建的過(guò)程。并針對(duì)腦部MRI中，海馬體成像邊界模糊、不連續(xù)等問(wèn)題，引入了先驗(yàn)形狀等信息作為三維LB分割算法的外力項(xiàng)，使得活動(dòng)曲面能夠快速、準(zhǔn)確地演化至海馬體的邊界，獲得海馬體的三維輪廓信息。在實(shí)驗(yàn)中，通過(guò)與AGDB方法的比較，三維LB分割算法在DSC、HD和計(jì)算速度等指標(biāo)上，均優(yōu)于AGDB方法。而針對(duì)三維分割帶來(lái)的耗時(shí)問(wèn)題，根據(jù)LB方法自身可高度并行化的特點(diǎn)，通過(guò)對(duì)算法流程的優(yōu)化，在單GPU平臺(tái)實(shí)現(xiàn)了算法的并行計(jì)算，進(jìn)一步提高了計(jì)算速度，解決了三維數(shù)據(jù)帶來(lái)的計(jì)算代價(jià)高的問(wèn)題。為了應(yīng)對(duì)可能出現(xiàn)的板載顯存不足的問(wèn)題，針對(duì)雙GPU平臺(tái)計(jì)算的特點(diǎn)，重新設(shè)計(jì)了算法的實(shí)現(xiàn)過(guò)程，使得LB方法能夠在并行平臺(tái)上進(jìn)行更大規(guī)模的計(jì)算。進(jìn)一步擴(kuò)展了LB方法在快速圖像處理領(lǐng)域中的應(yīng)用前景。

[1] Dubois B, Feldman H H, Jacova C, et al. Research criteria for the diagnosis of Alzheimer's disease: revising the NINCDS-ADRDA criteria[J]. Lancet Neurol, 2007, 6(8): 734-746.

[2] Dill V, Franco A R, Pinho M S. Automated methods for hippocampus segmentation: the evolution and a review of the state of the art[J]. Neuroinformatics, 2015, 13(2):133-150.

[3] Tong T, Wolz R, Coupé P, et al. Segmentation of MR images via discriminative dictionary learning and sparse coding: application to hippocampus labeling[J]. Neuroimage, 2013, 76(1):11-23.

[4] Kichang K, Uicheul Y, Dong-Kyun L, et al. Fully-automated approach to hippocampus segmentation using a graphcuts algorithm combined with atlas-based segmentation and morphological opening[J]. Magn Reson Imag, 2013, 31(7):1190-1196.

[5] Wu G, Shen D. Hierarchical label fusion with multiscale feature representation and label-specif i c patch partition[C]. Int Conf Med Imag Comput Computer-Assisted Intervent, 2014:299-306.

[6] Wu G, Kim M, Sanroma G, et al. Hierarchical multi-atlas label fusion with multi-scale feature representation and label-specific patch partition[J]. Neuroimage, 2015, 106: 34-46.

[7] Suppa P, Hampel H, Spies L, et al. Fully automated atlas-based hippocampus volumetry for clinical routine: validation in subjects with mild cognitive impairment from the ADNI cohort[J]. J Alzheimers Dis, 2015, 46(1):199-209.

[8] Zarpalas D, Zafeiropoulos A, Daras P, et al. Brain structures segmentation using optimum global and local weights on mixing active contours and neighboring constraints[C]. Int Symp Appl Sci Biomed Commun Tech, ACM, 2011: 127-132

[9] Zarpalas D, Gkontra P, Daras P, et al. Hippocampus segmentation through gradient based reliability maps for local blending of ACM energy terms[C]. Int Symp Biomed Imaging, IEEE, 2013: 53-56.

[10] Zarpalas D, Gkontra P, Daras P, et al. Gradient-Based Reliability Maps for ACM-Based Segmentation of Hippocampus[J]. IEEE Trans Biomed Eng, 2014, 61(4): 1015-1026.

[11] 林笑曼, 嚴(yán)壯志, 蔣皆恢,等. 基于三維格子玻爾茲曼模型的海馬結(jié)構(gòu)MRI快速分割[J]. 生物醫(yī)學(xué)工程學(xué)進(jìn)展, 2014, 35(1): 1-7.

[12] Wen J, Jiang J, Yan Z. A new lattice Boltzmann algorithm for assembling local statistical information with MR brain imaging segmentation applications[J]. Multidimen Syst Sign Proc, 2016:1-17.

[13] Balla-Arabé S, Gao X, Wang B. GPU accelerated edge-region based level set evolution constrained by 2D gray-scale histogram[J].IEEE Trans Imag Proc, 2013, 22(7): 2688-2698.

[14] Jones B D, Feng Y T. Effect of image scaling and segmentation in digital rock characterisation[J]. Comput Part Mech, 2016, 3(2):201-213.

[15] Wen J, Yan Z, Jiang J. Novel lattice Boltzmann method based on integrated edge and region information for medical image segmentation[J]. Biomed Mater Eng, 2014, 24(1):1247-1252.

[16] Bresson X, Vandergheynst P, Thiran J P. A Variational model for object segmentation using boundary information and shape prior driven by the Mumford-Shah functional[J]. Int J Comput Vis, 2006,68(2):145-162.

[17] ADNI 數(shù)據(jù)庫(kù)[R/OL]. http://www.hippocampal-protocol.net/SOPs/LINK_PAGE/Harmonized Protocol ACPC User Manual_biblio.pdf.

[18] Cabezas M, Oliver A, Lladó X, et al. A review of atlas-based segmentation for magnetic resonance brain images[J]. Comput Meth Program Biomed, 2011, 104(3): 158-177.