一種基于格子玻爾茲曼前向模型的GPU并行加速熒光擴(kuò)散斷層成像的方法

吳煥迪，嚴(yán)壯志，岑星星

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

0 引言

熒光擴(kuò)散斷層成像技術(shù)（Fluorescent Diffuse Optical Tomography，F(xiàn)DOT）具有無創(chuàng)、定量、靈敏度高、成本低等優(yōu)點，已成為醫(yī)學(xué)計算中一種強(qiáng)力的分子成像技術(shù)，在癌癥診斷、藥物機(jī)理研究、腦功能成像方面具有巨大的潛力[1-2]。通常選擇基于輻射傳輸方程（Radiative Transfer Equation，RTE）或擴(kuò)散方程（Diffusion Equation，DE）模型求解FDOT。前者是一種復(fù)雜的積分微分方程，用于描述光在生物組織等介質(zhì)中傳播的精確模型；后者是前者的一階球諧簡化模型。盡管許多研究人員致力于縮短RTE和DE的計算時間，然而即使是處理小尺寸圖像，求解RTE或DE模型也仍然不能滿足動態(tài)成像的要求[3]。

FDOT的求解方法通常比較耗時。因此，迫切需要一種新穎有效的傳播模型。格子玻爾茲曼方法（Lattice Boltzmann Method，LBM）與FDOT的協(xié)同發(fā)展有可能解決快速成像問題。LBM是一種基于宏觀和微觀之間的介觀尺度模型的偏微分方程數(shù)值求解工具，具有計算穩(wěn)定、物理含義明確、并行結(jié)構(gòu)簡單、易實現(xiàn)的特點。

目前LBM已經(jīng)廣泛應(yīng)用于流體力學(xué)和熱輻射，但是在光學(xué)傳輸領(lǐng)域LBM還是一個新興、極具前景的方法[4-5]。例如PATIDAR等[6]提出了基于LBM的輻射傳輸方程瞬態(tài)解并且利用二維組織模型進(jìn)行了驗證；ZHANG等[7]提出了一種用基于格子玻爾茲曼方法求解前向模型的方法，與利用COMSOL方法求解相比，該方法的求解速度比用COMSOL快2.3倍。但是目前為了獲得準(zhǔn)確的目標(biāo)位置和豐富的成像信息需要提高重建圖像的空間分辨率，上述方法在CPU上處理大規(guī)模的LBM計算速度較慢，不足以滿足快速成像的需求。由于LBM與GPU并行化的契合性十分高，二者的結(jié)合能解決空間分辨率和計算速度之間的矛盾問題。

因此我們提出了一種基于格子玻爾茲曼前向模型的GPU并行加速熒光擴(kuò)散斷層成像的方法。此方法通過構(gòu)建基于LBM的RET的前向模型，把LBM的碰撞、遷移、邊界處理過程在GPU上并行計算，可以大大提高求解速度，之后利用代數(shù)重建技術(shù)（Algebraic Reconstruction Technique，ART）計算重建結(jié)果，最終把計算結(jié)果可視化。為了驗證本方法的可行性和加速性能，我們還進(jìn)行了數(shù)值仿真實驗和物理仿真實驗。通過實驗，從耗時和準(zhǔn)確性方面證明了所提方法的可行性。

1 基于格子玻爾茲曼方法的熒光擴(kuò)散斷層成像模型

首先，從BGK-Boltzmann方程的離散形式推導(dǎo)，對速度、時間、空間都進(jìn)行離散化。由于微觀粒子無時無刻不做著非規(guī)則的運動，因此微觀粒子的運動速度方向可以是無窮多個；同樣的，粒子的離散速度也可以是無窮多個。離散速度定義為一個有限的集合，即{e0,e1,e2,...,en}，為離散速度的個數(shù)。同樣的，離散分布函數(shù)集合則為{f0,f1,f2,...,fn}，α為當(dāng)前離散速度方向，α=0,1,...,n。離散BGK形式為：

其中r定義為在空間Ω中的位置，t定義為時間，△t定義為時間間隔，往往一個△t對應(yīng)LBM的一次迭代，為松弛時間，fa(r,t)為當(dāng)前離散速度方向的分布函數(shù)，為當(dāng)前離散速度方向的平衡的分布函數(shù)。同時離散的輻射傳輸方程為：

其中c為光子在生物組織中的傳播速度，表示方向。表示在t時刻、在空間r處的方向為s的輻射強(qiáng)度。μa(r)和μs(r)分別表示在位置r處的吸收系數(shù)和散射系數(shù)。表示從方向s到s'的光子傳播的散射相函數(shù)。為光源場的光強(qiáng)分布函數(shù)。

為了把格子玻爾茲曼方程和輻射傳輸方程構(gòu)建至光傳播的離散模型，首先對式（1）兩邊都乘以v，其中v表示光子的頻率，再把輻射傳輸模型中的光源與LBM中的平衡態(tài)分布函數(shù)做類比，可以得到：

再把輻射強(qiáng)度密度函數(shù)和速度分布函數(shù)做關(guān)聯(lián)，h為普朗克常數(shù)，可以得到：

通過上述的關(guān)聯(lián)，把式(3)和(4)代入式(1)，最終可以得到基于LBM的RTE離散化模型：

其中φ(r+ci△t,si,t+△t)表示在下一個時刻，在離散速度方向上的平均功率能量密度。S(r,被離散化為:用羅賓邊界條件處理邊界上的權(quán)值變化，最后通過ART重建計算結(jié)果。

綜上所述，本研究建立了基于LBM的FDOT的前向模型并介紹了重建方式。

2 結(jié)合GPU的格子玻爾茲曼前向模型

由于LBM節(jié)點的獨立任務(wù)的內(nèi)存需求相對較低，考慮LBM固有的并行性，這種方法非常適合在GPU上實現(xiàn)。通過使用GPU和CUDA通用并行計算架構(gòu)對LBM進(jìn)行編程，以加速基于格子玻爾茲曼前向模型的計算任務(wù)。本方法的框架如圖1所示。框架分為兩部分：CPU和GPU。與CPU相關(guān)的程序包括網(wǎng)格建立、坐標(biāo)離散化、邊界參數(shù)計算、CPU-GPU數(shù)據(jù)傳輸、代數(shù)重建和數(shù)據(jù)可視化。GPU相關(guān)程序主要由碰撞處理、遷移處理和邊界處理組成。

其中，wi表示相位函數(shù)的離散權(quán)重因子。在LBM中，通常使用DmQn來描述粒子的空間維度和離散速度方向，其中m和n分別代表空間維度和離散方向。為了簡化模型復(fù)雜度，在仿真和仿體實驗中選擇D3Q6模型，即三維六個離散速度方向模型，定義wi=1/6，i=1,2,3,4,5,6，并采用各向同性的介質(zhì)。生物組織的光學(xué)參數(shù)（如吸收系數(shù)和散射系數(shù)）和組織邊界被離散在均勻分布的立方網(wǎng)格上，從LBM的角度來看，光子被視作離散在網(wǎng)格中的粒子，而描述生物組織中光子傳播的LBM主要可以分為碰撞和遷移過程，分別是：

圖1 基于LBM的FDOT在GPU上并行化實現(xiàn)Fig.1 Parallel implementation of LBM-based FDOT on GPU

因此粒子在離散網(wǎng)格上不斷地進(jìn)行碰撞和遷移以模擬光子在生物組織中的傳播過程。通過對網(wǎng)格上每個節(jié)點的6個相鄰節(jié)點對應(yīng)的分布函數(shù)進(jìn)行計算，每次迭代后更新當(dāng)前節(jié)點的分布函數(shù)值和平衡分布函數(shù)，即模擬光子在傳播過程中當(dāng)前位置的能量權(quán)值隨時間的變化過程。本研究采

利用兩級并行將所有LBM的節(jié)點值分配給每個線程。在第一級中，每個線程網(wǎng)格由幾個線程塊組成。當(dāng)調(diào)用CUDA核函數(shù)時，線程塊的線程在一個多處理器（MP）上同步執(zhí)行。在第二級，每個線程塊由多個線程組成，每個單獨的LBM的節(jié)點值被分配給線程，每個線程將在流處理器（SP）上獨立執(zhí)行。在其實際操作中，核函數(shù)以線程網(wǎng)格的形式組織。為了避免數(shù)據(jù)競爭并確保同步，分別創(chuàng)建了三個主要的核函數(shù)并按順序運行每個核函數(shù)。在編譯程序之后，CUDA將計算任務(wù)映射到硬件從而以動態(tài)調(diào)度的方式執(zhí)行線程，主要針對基于LBM的碰撞、遷移、邊界處理過程進(jìn)行并行加速計算。式（7）表明了碰撞過程只涉及當(dāng)前位置節(jié)點的相關(guān)參數(shù)的計算，通過分別啟動多個線程對單個節(jié)點的計算所需的參數(shù)讀取，從而完成并行計算碰撞過程；式（8）表明了遷移過程只涉及當(dāng)前位置節(jié)點的相應(yīng)離散方向上的數(shù)據(jù)轉(zhuǎn)移，通過啟動兩個不同的線程對單個節(jié)點的數(shù)據(jù)根據(jù)遷移方向覆蓋至下一個節(jié)點，即可完成遷移過程；而對于邊界處理過程，我們選擇了掩模處理方法。對已經(jīng)離散化的網(wǎng)格模型，完成邊界定義，生成一個與生物組織網(wǎng)格模型大小一樣的掩模矩陣，分別用0、1區(qū)分邊界和非邊界區(qū)域。當(dāng)啟動GPU線程時，僅對掩模標(biāo)記為1的區(qū)域進(jìn)行邊界化處理。

CUDA的硬件層面分為主機(jī)（Host）端和設(shè)備（Device）端。基于LBM的FDOT的完整CUDA算法流程如下：

（1）確定合適DmQn的形式和離散權(quán)重；

（2）對三維生物組織模型進(jìn)行網(wǎng)格化，包括吸收系數(shù)、散射系數(shù)等；

（3）對探測器采集的數(shù)據(jù)進(jìn)行相同離散尺寸的網(wǎng)格化；

（4）將計算涉及的數(shù)據(jù)由Host端傳輸?shù)紻evice端；

（5）對n個離散方向的粒子進(jìn)行LBM計算：包括粒子的碰撞、遷移、邊界處理、更新，直到模型達(dá)到穩(wěn)定狀態(tài)或者迭代次數(shù)上限的計算；

（6）當(dāng)完成核函數(shù)的計算任務(wù)，Device端的計算結(jié)果傳輸回Host端；

（7）對計算結(jié)果進(jìn)行ART重建；

（8）對重建結(jié)果進(jìn)行可視化操作。

綜上所述為結(jié)合GPU的并行化實現(xiàn)過程。

3 實驗與結(jié)果

為了驗證所提方法，我們選擇8核、16 GB內(nèi)存、64位、裝有NVIDIA GTX660顯卡的個人電腦進(jìn)行實驗。操作系統(tǒng)為Windows 7。本研究使用CUDA 6.5，OpenMP 4.0和Visual Studio 2012進(jìn)行編譯。

在單個GPU上進(jìn)行實驗測試，本研究不涉及GPU集群和消息傳遞接口（MPI）。實驗由仿真實驗和仿體實驗兩部分組成。為了簡化模型，均勻介質(zhì)圓柱體的吸收系數(shù)μa和散射系數(shù)μs分別選擇0.3 cm-1和10 cm-1。前者旨在評估在相同編譯環(huán)境下的加速度和準(zhǔn)確性，在仿真實驗中，分別實現(xiàn)了C++串行，GPU并行和OpenMP三種不同的編程方法，并具體分析了當(dāng)LBM的網(wǎng)格尺寸分別為15×15×26，31×31×51，41×41×61，51×51×81，51×51×101，61×61×81時的性能指標(biāo)；后者通過與傳統(tǒng)擴(kuò)散方程相比來驗證所提方法的實用性?；跀U(kuò)散方程方法通過有限元方法求解，有限元方法在Matlab中使用Gmsh[8]和Toast++[9]完成求解，其中Gmsh用于FEM的網(wǎng)格生成，Toast++用于FEM的求解。

3.1 仿真實驗

實驗采用一個里面嵌有兩個圓柱體的大圓柱體模擬均勻介質(zhì)下的光傳播情況。采用半徑為1.5 cm、高5 cm的大圓柱體模擬整個物理模型，采用兩個半徑為0.1 cm、高0.2 cm的小圓柱體分別模擬兩個不同位置的熒光標(biāo)記團(tuán)，記為FT1和FT2。假設(shè)大圓柱體的底部圓心作為三維坐標(biāo)系（x，y，z）的原點，布置不同位置的FT1和FT2以模擬不同情況，從而證明方法的穩(wěn)定性。本研究實現(xiàn)了3種模擬情況，分別命名為Case1、Case2和Case3，如圖2所示。FT1和FT2的具體坐標(biāo)如表1所示。

在保證同C++的計算結(jié)果一致、精度誤差極小的前提下，使用所提方法可以得到Case1至Case3的、離散尺度分別為15×15×26和61×61×81的重建結(jié)果，如圖3所示，圖3（a）代表離散尺度為15×15×26的重建結(jié)果；同理圖3（b）代表離散尺度為61×61×81的重建結(jié)果。其中黑色粗體圓分別表示真實熒光團(tuán)（圓柱體形狀）所在位置。從圖3可以發(fā)現(xiàn)，我們所提方法具有較好的準(zhǔn)確性。網(wǎng)格模型離散尺度為61×61×81的圖像重建質(zhì)量遠(yuǎn)好于離散尺度為15×15×26的重建結(jié)果。

為了對所提方法的加速度進(jìn)行評估，我們統(tǒng)計了C++、openMP、GPU三者方法的對于Case1至Case3的多次平均計算耗時情況，如表2所示。從表2中可以發(fā)現(xiàn)，當(dāng)處理離散尺度為15×15×26的網(wǎng)格模型時，三者的耗時情況在同一個數(shù)量級；但當(dāng)處理離散尺度不斷增大時，基于GPU的并行化的優(yōu)勢逐漸明顯；這是由于C++串行算法無法利用好資源，受到算法流程的限制，只能做等待，計算資源的利用率很低；OpenMP處理小尺寸的網(wǎng)格模型時，系統(tǒng)內(nèi)存資源豐富，可以有充足的空間進(jìn)行數(shù)據(jù)復(fù)制和線程切換，但是當(dāng)處理大尺寸的網(wǎng)格模型時，在線程之間切換需要額外的內(nèi)存空間以釋放資源，從而導(dǎo)致處理速度變慢；而GPU與LBM的契合性較高，在處理較大離散尺寸的網(wǎng)格模型時擁有并行性，提升了計算效率。

圖2 數(shù)值仿真實驗的Case1至Case3的實驗設(shè)計示意圖Fig.2 Diagram of experimental design of Case1 to Case3 in numerical simulation experiments

表1 三種數(shù)值仿真實驗的熒光團(tuán)具體坐標(biāo)Tab.1 Specific coordinates of fluorophores in three numerical phantom experiments

圖3 基于不同離散尺寸的網(wǎng)格模型的FDOT重建結(jié)果Fig.3 FDOT reconstruction results based on mesh models with different discrete sizes

表2 基于C++、OpenMP、GPU的三種不同編程方式的計算耗時統(tǒng)計Tab.2 Computing time-consuming statistics based on C++,OpenMP and GPU

3.2 仿體實驗

采用FDOT/XCT混合成像系統(tǒng)，將待成像物體放置在360o旋轉(zhuǎn)臺上緩慢旋轉(zhuǎn)，通過CCD相機(jī)測量從仿體表面溢出的光子的輻射強(qiáng)度，細(xì)節(jié)見GUO[10]的文章。本研究將直徑為0.4 cm的熒光團(tuán)管狀物嵌入高度為4.4 cm、半徑為1.5 cm的圓柱體中作為仿體模型。

通過對比基于DE方法的求解結(jié)果從而驗證我們所提方法的實用性。基于所提方法的LBM網(wǎng)格節(jié)點和基于DE方法的網(wǎng)格元素個數(shù)分別采用49 011離散節(jié)點和50 469個節(jié)點從而保證二者方法的節(jié)點個數(shù)相近以公平對比。本研究可視化了計算尺寸為31×31×51的計算結(jié)果，如圖4所示，（a）為真實仿體的示意圖；（b）為基于LBM的重建結(jié)果；（c）為基于DE方法的重建結(jié)果。LBM的重建結(jié)果有明顯的偽影，但中間熒光團(tuán)的形狀幾乎相同，與基于DE方法的計算結(jié)果保持一致。

圖4 仿體實驗的重建結(jié)果對比Fig.4 Comparison of reconstruction results of physical phantom experiments

表3為基于LBM方法和基于DE方法的耗時對比。當(dāng)計算尺寸為15×15×26時，本研究的方法計算速度比基于DE的方法快25倍；當(dāng)計算尺寸為31×31×51時，本研究的方法計算速度比基于DE的方法快118倍；當(dāng)離散尺度更大時，受到內(nèi)存資源的限制，基于DE的方法無法求解，而本研究的方法仍然可以求解。

表3 基于LBM方法和基于DE方法的計算耗時Tab.3 Time consumption obtained by the DE and LBM method

4 結(jié)論

綜上所述，本研究驗證了所提方法的準(zhǔn)確性和加速性。結(jié)果表明，本方法可以在保證準(zhǔn)確性的前提下，快速求解，有助于LBM的發(fā)展和FDOT的快速成像。但是本研究提出的方法仍然存在一些待改進(jìn)的地方：基于LBM的FDOT的光傳輸模型沒有解釋如何對離散時間和離散空間做耦合；D3Q6作為實現(xiàn)模型，過于簡化，帶來了偽影現(xiàn)象；生物模型做了較多的簡化，在實際應(yīng)用中仍然需要復(fù)雜化；基于GPU的結(jié)合方法目前仍然無法滿足實時的動態(tài)FDOT的技術(shù)需求。在今后的工作中，將重點研究非均勻介質(zhì)、不規(guī)則形狀和GPU集群方向。