伍紹佳陳 皓廖 麗桂建保

1(肇慶廣播電視大學(xué) 肇慶 526060)

2(中國(guó)科學(xué)院深圳先進(jìn)技術(shù)研究院 深圳 518055)

BPF 重建算法的 CUDA 并行實(shí)現(xiàn)

伍紹佳1陳 皓2廖 麗1桂建保2

1(肇慶廣播電視大學(xué) 肇慶 526060)

2(中國(guó)科學(xué)院深圳先進(jìn)技術(shù)研究院 深圳 518055)

反投影濾波(Backprojection-Filter，BPF)算法憑借其可實(shí)現(xiàn)感興趣區(qū)域重建的優(yōu)點(diǎn)，近年來(lái)逐漸被應(yīng)用到錐束 CT 中。但是，由于算法的復(fù)雜性，實(shí)踐中存在耗時(shí)問(wèn)題，同時(shí)其 GPU 加速的實(shí)現(xiàn)亦存在顯存不足等問(wèn)題。因此，文章提出了一種基于 CUDA 的 BPF 并行加速算法。通過(guò)設(shè)計(jì)高效的算法框架，在保留其重建精度的前提下，有效地減少所需顯存。此外，總結(jié)了正投影算法及 BPF 算法中采用的加速策略，如利用算法特征加速等，并引入顯存池的概念優(yōu)化算法架構(gòu)。仿真實(shí)驗(yàn)結(jié)果表明，在精確重建的前提下，采用新框架重建 512×512×512 數(shù)據(jù)只需 8.055 s，感興趣區(qū)域重建只需 4.566 s，只需 1.523 s 便可輸出第一部分?jǐn)?shù)據(jù)，且能把顯存占用從 2.5 GB 減少到 100 MB 以下，適用于大數(shù)據(jù)重建。

反投影濾波算法；錐束 CT；感興趣區(qū)域成像；圖形處理器；圖像重建；并行計(jì)算架構(gòu)

1 引 言

自 1971 年第一臺(tái) CT 安裝以來(lái)，CT 已經(jīng)經(jīng)歷了幾代的發(fā)展。近年來(lái)，錐束 CT 已經(jīng)成為醫(yī)療診斷和放療技術(shù)的重要工具，因此，如何提高錐束 CT 的性能是近年來(lái)研究的熱點(diǎn)之一。 針對(duì) CT 的性能，劑量、圖像重建速度及圖像質(zhì)量是其重建的三個(gè)主要評(píng)判標(biāo)準(zhǔn)[1]。本文主要解決圖像重建的速度問(wèn)題及其在 GPU 應(yīng)用過(guò)程中涉及的顯存問(wèn)題。近十幾年來(lái)，解析法重建憑借其高效性得到大量的研究，與常規(guī)解析算法相比(如 FDK[2]等濾波反投影算法)，反投影濾波(Backprojection-Filter，BPF)算法[3,4]憑借利用較少量的投影數(shù)據(jù)重建出精確的圖像，得到學(xué)者的廣泛認(rèn)可。憑借感興趣區(qū)域掃描和重建，BPF 算法可實(shí)現(xiàn)降低劑量的目的，但其實(shí)現(xiàn)過(guò)程還存在運(yùn)算量大等問(wèn)題。故本文提出一個(gè)基于 CUDA 的BPF 算法加速架構(gòu)，旨在更好地將 BPF 算法應(yīng)用到實(shí)際中。

GPU 憑借其強(qiáng)大的并行計(jì)算能力使不同的CT 重建算法得以發(fā)展。其中，傳統(tǒng)經(jīng)典的 CT 重建算法，如 FDK 算法和迭代算法的并行加速模型[5,6]已經(jīng)得到廣泛的研究。而有關(guān) BPF 加速的研究并不多，因此，本文著重于 BPF 的算法加速。根據(jù)文獻(xiàn)搜索，Zheng 等[7]設(shè)計(jì)一個(gè)應(yīng)用于GPU 中通用 Pi 線選擇的方案，但是并沒(méi)有詳細(xì)地介紹如何利用 GPU 進(jìn)行優(yōu)化計(jì)算。此外，Zou等[8]曾提出基于 GPU 的 BPF 算法優(yōu)化方案，盡管時(shí)間得到大幅度的提升，但是其架構(gòu)并不是最優(yōu)，如在大數(shù)據(jù)問(wèn)題上并不能很好地應(yīng)用或在實(shí)際應(yīng)用中存在顯存使用過(guò)多的問(wèn)題；同時(shí)，由于其框架基于多 GPU，因而難以在商業(yè)項(xiàng)目中進(jìn)行應(yīng)用。除此之外，大部分相關(guān)研究也僅僅關(guān)注重建速度的優(yōu)化，并沒(méi)有考慮顯存的合理利用，因此，本文引入顯存池的概念去優(yōu)化顯存的管理。根據(jù) BPF 算法的特性，本文首先提出適合實(shí)際應(yīng)用的單 GPU 優(yōu)化框架；其次對(duì) CUDA 實(shí)現(xiàn)過(guò)程中采用的優(yōu)化策略進(jìn)行總結(jié)；最后，本文將闡述如何利用 GPU 進(jìn)行 Shepp-Logon 模體投影數(shù)據(jù)的快速生成。

2 BPF 算法介紹

近年來(lái)，BPF 算法得到很大的推廣，其原型由潘曉川小組在 2004 年以三維重建的形式提出，后來(lái)又提出其扇束模式。Yu 等[9]在 2008 年提出基于圓周軌跡的 BPF 算法，進(jìn)一步地完善其在錐束 CT 的應(yīng)用。如圖1 所示，本文處理的投影數(shù)據(jù)亦為基于圓軌跡掃描獲得。

BPF 算法是一個(gè)基于弦的反投影濾波重建算法。如圖2 所示，弦被定義為射線源兩點(diǎn)間的連線，即和，因此，沿著弦的單位向量可以表述為：

該弦上所有點(diǎn)可以描述為：

圖1 BPF 算法幾何掃描圖Fig.1. Scan geometry of BPF

反投影部分可以描述為：

而邊界項(xiàng) C 則可以表達(dá)為：

P0是 X 射線經(jīng)過(guò) Pi 線中間點(diǎn)的投影數(shù)據(jù)平均值。

綜上所述，BPF 算法的實(shí)現(xiàn)可以分解為以下5 步：

(1)對(duì)加權(quán)的錐束數(shù)據(jù)進(jìn)行求導(dǎo)；

(2)沿著弦對(duì)加權(quán)后的投影數(shù)據(jù)進(jìn)行反投影；

(3)對(duì)反投影后的數(shù)據(jù)進(jìn)行濾波；

(4)計(jì)算邊界項(xiàng)進(jìn)行補(bǔ)償；

(5)累加邊界項(xiàng)和濾波后數(shù)據(jù)，并將結(jié)果乘以相關(guān)系數(shù)。

3 基于 CUDA 的 BPF 算法優(yōu)化方案

GPU 作為輔助處理器，專(zhuān)門(mén)為協(xié)助 CPU 進(jìn)行通用科學(xué)和工程計(jì)算而設(shè)計(jì)。接下來(lái)，本文首先針對(duì)模擬數(shù)據(jù)的生成進(jìn)行介紹并提及其加速方法，隨后詳細(xì)介紹基于 CUDA 的 BPF 算法新框架，最后介紹在 CUDA 實(shí)現(xiàn)過(guò)程中采用的加速策略。

3.1 正投影數(shù)據(jù)的生成

重建算法的魯棒性一般由其重建模體的精度來(lái)進(jìn)行驗(yàn)證，因此，獲得模體的投影數(shù)據(jù)是驗(yàn)證算法的必要條件。接下來(lái)，結(jié)合 GPU 簡(jiǎn)單地介紹本文采用的正投方案。

Siddon 算法是 CT 算法中正投影的經(jīng)典算法。盡管該算法極大地保證其準(zhǔn)確率，但是該算法耗時(shí)長(zhǎng)，在實(shí)際應(yīng)用中的效率并不高。因此，本文結(jié)合 Nvidia 公司的 sdk 例子[10,11]，對(duì)原有的投影算法進(jìn)行優(yōu)化，得出本文的正投影算法。

Sidon 算法把將投影物體看作一個(gè)三維網(wǎng)格，通過(guò)連接射線源與探測(cè)器相應(yīng)像素點(diǎn)，求出該連線與其經(jīng)過(guò)每個(gè)小方格上的入射點(diǎn)與出射點(diǎn)，計(jì)算其穿過(guò)的距離，整理為一個(gè)系統(tǒng)矩陣，然后再通過(guò)插值求出網(wǎng)格的灰度值，最后累加并獲得探測(cè)器上每一點(diǎn)的投影值。但是，由于此方法需要計(jì)算每一個(gè)網(wǎng)格的穿過(guò)距離，并不能很好地進(jìn)行加速。因此，如圖3 所示，本文選擇只計(jì)算穿過(guò)整個(gè)大網(wǎng)格的入射點(diǎn)和出射點(diǎn)，然后根據(jù)需要，以特定的步長(zhǎng)把兩點(diǎn)間的連線進(jìn)行劃分，并利用三維線性插值計(jì)算每段的灰度值，最后累加獲得平板該點(diǎn)的投影值。

圖3 正投影解釋圖Fig.3. Introduction of projection algorithm

在 GPU 端，數(shù)據(jù)讀入后，以每條射線為一個(gè)線程，如探測(cè)器像素個(gè)數(shù)為 512×512，可以設(shè)512×512 的線程。不同射線劃分的段數(shù)并不一致，每個(gè)線程計(jì)算該射線線段組的累加值，而線段的累加值則通過(guò)把物體映射到三維紋理內(nèi)存上獲得。選擇三維紋理顯存有兩個(gè)原因：(1)由于此類(lèi)顯存獲取為隨機(jī)存取，紋理顯存能優(yōu)化讀取速度；(2)紋理顯存提供三維線性插值。加速后的正投影方案在保證投影計(jì)算質(zhì)量的前提下，其計(jì)算速度大幅地提升。此外，根據(jù)原始數(shù)據(jù)的不同大小，可以把上述原始數(shù)據(jù)劃分為若干塊，每次只讀取一塊進(jìn)行計(jì)算，并將計(jì)算值存于平板數(shù)組，進(jìn)而進(jìn)行下一數(shù)據(jù)塊的計(jì)算，這樣可以解決由于模體數(shù)據(jù)量大導(dǎo)致顯存不足的問(wèn)題。

3.2 BPF 算法實(shí)現(xiàn)

Zou 等[8]曾提出一個(gè)基于 GPU 的 BPF 算法框架，但其框架具有一定的局限性。由于醫(yī)生對(duì)圖像質(zhì)量要求的提高，512×512×512 的圖像重建往往不能滿足其需求。同時(shí) BPF 算法在原有的框架下需要很大的顯存，以致在實(shí)際應(yīng)用中要選用兩個(gè)GPU，這樣在實(shí)際應(yīng)用中勢(shì)必增加項(xiàng)目的成本。

根據(jù)上述算法討論，BPF 算法在實(shí)現(xiàn)中遇到兩個(gè)問(wèn)題：(1)計(jì)算反投影過(guò)程慢；(2)所需存儲(chǔ)空間大。而計(jì)算量大是 CT 重建一直存在的問(wèn)題，且由于 BPF 算法需要計(jì)算邊界項(xiàng)，其所需存儲(chǔ)空間為原本算法的兩倍，從而增加 BPF 算法實(shí)現(xiàn)的困難。

3.2.1 基于 CUDA 的 BPF 框架

如上所述，許多現(xiàn)行的框架都選擇把每一個(gè)重建點(diǎn)的計(jì)算分配給每一個(gè)線程，如為512×512×512 的數(shù)據(jù)開(kāi)辟 512×512×512 個(gè)線程去進(jìn)行反投影，但是此方法在 BPF 算法中存在嚴(yán)重的顯存問(wèn)題。因?yàn)槠洳粌H要保留512×512×512 的反投影數(shù)據(jù)，還需要保留相關(guān)的邊界項(xiàng)，因此，需要設(shè)計(jì)一種在一般 GPU 上能運(yùn)行并計(jì)算速度快的算法框架。

如圖4 所示，首先從 CPU 端讀入數(shù)據(jù)并存入 GPU 端，然后根據(jù)公式對(duì)投影數(shù)據(jù)進(jìn)行求導(dǎo)；接著對(duì)重建層進(jìn)行分塊，根據(jù)重建的不同塊所需的投影數(shù)據(jù)，導(dǎo)入相應(yīng)的求導(dǎo)數(shù)據(jù)，并對(duì)該塊進(jìn)行反投影(包括計(jì)算該塊的反投影值和其邊界項(xiàng)值)；反投影后再進(jìn)行濾波，濾波后將該層數(shù)據(jù)和邊界項(xiàng)數(shù)據(jù)進(jìn)行累加，這樣一段數(shù)據(jù)塊重建完成并輸出到終端，遍歷完所有重建塊后，則完成整個(gè)重建。此流程的優(yōu)點(diǎn)如下：

(1)可擴(kuò)展性。由于實(shí)際應(yīng)用中不僅進(jìn)行重建處理，還可能在其中加入降噪等校正操作，若重建處理占據(jù)所有顯存，則大大地降低算法的靈活性。

(2)顯存空間使用少。新框架中只需要反復(fù)利用同一塊顯存，不需要同時(shí)保留整個(gè)重建體數(shù)據(jù)及其邊界項(xiàng)數(shù)據(jù)，其分析見(jiàn)下節(jié)。

圖4 基于 CUDA 的 BPF 算法新框架Fig.4. Frame of BPF algorithm based on CUDA

(3)輸出快。重建一塊數(shù)據(jù)塊后便可輸出，并不需要等待所有數(shù)據(jù)重建完再進(jìn)行輸出，因此在實(shí)時(shí)性上，此框架比舊方法更為優(yōu)越。

此外，需要說(shuō)明以下兩點(diǎn)：首先，此框架并不需要擔(dān)心過(guò)度地進(jìn)行紋理綁定，因?yàn)榻壎ú⒉皇前褦?shù)據(jù)重新傳輸，只是綁定指針，所以內(nèi)部循環(huán)并不會(huì)大量提升計(jì)算時(shí)間。因此，算法實(shí)現(xiàn)過(guò)程中，雖然涉及紋理內(nèi)存的不斷映射，但由于其只把相應(yīng)指針進(jìn)行映射，并不會(huì)占據(jù)過(guò)多的時(shí)間；其次，雖然進(jìn)行一次 512×512×512 的數(shù)據(jù)重建處理的時(shí)間比目前分段處理的時(shí)間短(因?yàn)樾路椒ㄐ枰粩嗟亻_(kāi)啟 kenel 函數(shù))，但是，由于每一塊數(shù)據(jù)塊處理的數(shù)據(jù)量大，并已超過(guò)一次能同時(shí)處理的數(shù)量，即所有 SM 一次處理的線程數(shù)。所以，與一次處理相比，新方法所用的時(shí)間并不會(huì)大大的增加，而其降低顯存空間占用的優(yōu)勢(shì)則較舊方法優(yōu)越。

3.2.2 新框架內(nèi)存分析

根據(jù)圖4 所示，Zou 等[8]需要存儲(chǔ)整個(gè)重建空間的數(shù)據(jù)，并將其累加，其中包括圖像項(xiàng)和邊界項(xiàng)數(shù)據(jù)，過(guò)多地開(kāi)辟顯存容易導(dǎo)致其他操作無(wú)法進(jìn)行。新框架下只需要一塊重建塊的顯存空間，可以省略大部分圖像項(xiàng)和邊界項(xiàng)數(shù)據(jù)的空間。以下為簡(jiǎn)化的內(nèi)存分析：

設(shè)投影數(shù)據(jù)大小為 M×M，角度為 T，一般重建的大小與投影數(shù)據(jù)的大小相仿，設(shè)為 P×P×P，假設(shè) P 為 2 的冪次方(一般 P 為512)，由于含有圖像項(xiàng)，邊界項(xiàng)及合并項(xiàng)三部分，因此所占空間大小為 3×P×P×P。在濾波過(guò)程中會(huì)產(chǎn)生頻域數(shù)據(jù)，其長(zhǎng)度一般為其最接近的 2 的冪次方，以保證其精度，假設(shè)為2×P×P×P，因此，原方案為：

Memory＝M×N+5×P×P×P由于本架構(gòu)把重建層數(shù)進(jìn)行分塊，每一塊的層數(shù)為 H(其大小可以根據(jù)實(shí)際需要進(jìn)行控制)，因此總的顯存需求為：

Memory＝M×N+5×P×P×H由公式可以看出，在顯存的利用上，新框架有絕對(duì)的優(yōu)勢(shì)，新架構(gòu)所需空間僅為 Zou等[8]的 H/P，從而方便其他處理的進(jìn)行。若為512×512×512 的數(shù)據(jù)，以 16 層數(shù)據(jù)為一數(shù)據(jù)塊，所需顯存空間僅為原方案的 0.03125 倍，即需要不到 100 Mb 的顯存(考慮到其他輔助變量的存儲(chǔ))。因此，在大數(shù)據(jù)重建下，新框架的優(yōu)勢(shì)尤其突出，使其在重建數(shù)據(jù)量大的情況下能保證程序正確性的同時(shí)不影響其他操作。

3.2.3 顯存的統(tǒng)一管理

由于對(duì)重建層進(jìn)行分塊，程序會(huì)增加對(duì)顯存的申請(qǐng)和銷(xiāo)毀，這樣，由于申請(qǐng)內(nèi)存塊的大小不定，頻繁使用時(shí)會(huì)造成大量的顯存碎片進(jìn)而降低性能。本文對(duì)內(nèi)存池的思想進(jìn)行簡(jiǎn)化，并應(yīng)用至重建算法的顯存管理上。由于重建算法的顯存需求是固定的，因此，針對(duì)不同病人的圖像處理也是固定的，所以只需在程序開(kāi)始時(shí)分配固定大小的顯存空間即可。

其顯存管理機(jī)制如下：首先，在不影響其他操作的情況下，根據(jù)不同機(jī)器的顯存空間進(jìn)行一次性分配。申請(qǐng)整塊顯存空間后，根據(jù)不同的處理，選擇合適的計(jì)算單位進(jìn)行顯存的申請(qǐng)和釋放。如果每次申請(qǐng)空間過(guò)小，將對(duì)顯存進(jìn)行多次申請(qǐng)而影響算法的整體性能；若每次申請(qǐng)過(guò)大，則申請(qǐng)的顯存量會(huì)超出允許范圍。因此，本文根據(jù)上述內(nèi)存分析進(jìn)行分塊的自動(dòng)計(jì)算并自動(dòng)分配顯存空間。這樣，向顯存池申請(qǐng)顯存比直接向系統(tǒng)申請(qǐng)快，且不會(huì)產(chǎn)生碎片，在實(shí)際應(yīng)用中效果顯著。

3.2.4 CUDA 相關(guān)加速策略

由圖4 可看出，反投影濾波可以大致劃分為反投影、濾波、輸入輸出三部分，接下來(lái)根據(jù)不同部分介紹適用于 BPF 算法的加速策略。

(1)反投影部分

反投影部分包括圖像項(xiàng)和邊界項(xiàng)的計(jì)算，是整個(gè)重建耗時(shí)最長(zhǎng)的部分。由于其存在隨機(jī)存取的現(xiàn)象，因此把全局內(nèi)存的數(shù)據(jù)映射到紋理內(nèi)存可以大幅減少讀取數(shù)據(jù)的耗時(shí)，而二維紋理的利用也可比一維紋理優(yōu)化 1/3 左右的時(shí)間；此外，如 Zou 等[8]提及的，在反投影中需采用預(yù)判策略，即由于 BPF 算法是采取感興趣區(qū)域重建的策略，如果只需重建部分區(qū)域，可以提前進(jìn)行檢測(cè)從而減少不必要的反投影，此方法能大幅度地減少重建時(shí)間；BPF 算法涉及較多輔助變量，這些變量數(shù)量較多且不同重建點(diǎn)具有不同值，因此不便存儲(chǔ)于常量存儲(chǔ)器，但其在不同層上的值一致，因此本文采取計(jì)算一次并存儲(chǔ)于全局內(nèi)存的策略，減少不必要的計(jì)算。

(2)濾波部分

BPF 算法的濾波為希爾伯特濾波，與 FDK 算法的 Ramp 濾波器優(yōu)化方案基本一致。本文利用CUFFT 庫(kù)對(duì)傅里葉變換進(jìn)行加速優(yōu)化，在新框架下，它的批處理特點(diǎn)不再是投影數(shù)據(jù)的批處理，而是圖像數(shù)據(jù)塊的批處理。在新框架下可以擴(kuò)展出許多模型，可以選擇一次濾波兩層或更多，這樣，批處理的特點(diǎn)也可進(jìn)行最大化的利用。

(3)輸入輸出部分

由于新框架是一塊塊輸出重建數(shù)據(jù)，本文把合并過(guò)程歸入輸入輸出部分中。從 GPU 拷貝到CPU 的時(shí)間基本可以完全隱藏。通過(guò)創(chuàng)建不同的流：一個(gè)進(jìn)行運(yùn)算，另一個(gè)則用來(lái)進(jìn)行數(shù)據(jù)輸出。這樣一層計(jì)算完成后，可以通過(guò)輸出流把傳輸時(shí)間隱藏，再通過(guò)在 CPU 中開(kāi)辟新的線程把獲得的數(shù)據(jù)輸出到硬盤(pán)，這樣可以最大限度的縮短輸出時(shí)間。

4 數(shù)據(jù)實(shí)驗(yàn)結(jié)果和分析

為了驗(yàn)證算法的有效性，本文設(shè)計(jì)了 2 個(gè)測(cè)試方向：(1)重建的準(zhǔn)確度；(2)重建的耗時(shí)。測(cè)試模體為三維的 Shepp-Logan 模體，重建規(guī)模為 512×512×512。測(cè)試平臺(tái)為 3.2 GHz Intel Core(TM) i5-3470 雙核 CPU，16 GB 內(nèi)存，Nvidia GeForce GTX660 Ti 顯卡；開(kāi)發(fā)環(huán)境為Visual Studio 2010，cuda5.5 runtime API；模擬放線的幾何環(huán)境分別為：射線源到探測(cè)器距離為700 mm，射線源到旋轉(zhuǎn)中心的距離為 350 mm，探測(cè)器間距為 0.1 mm、長(zhǎng)寬均為 51.2 mm，Shepp-Logan 模體的大小為 512×512×512，采樣間隔為 0.05 mm。

4.1 準(zhǔn)確性比較

圖5 采用 Shepp-Logan 模體重建的不同截面斷層圖像Fig.5. The reconstructed slice images of Shepp-Logan phantom from different sections

圖5(a)和(b)分別為原圖及 GPU 重建結(jié)果的切片方向圖像比較；而圖5(e)則為第 321 層中第 256 行的灰度值曲線比較，結(jié)果顯示重建數(shù)據(jù)與原數(shù)據(jù)變化規(guī)律基本上一致，誤差在容忍范圍內(nèi)。由于本實(shí)驗(yàn)為三維重建，圓周軌跡下的掃描結(jié)果都是近似重建，圖5(c)和(d)分別為重建物體冠狀面和矢狀面的結(jié)果，結(jié)果基本符合小錐角情況。綜上分析，新框架在精度上滿足需求。

4.2 時(shí)間比較

由于不同研究小組所采用的機(jī)器不一致，導(dǎo)致獲得的重建時(shí)間的可比性不高，如 Zou 等[8]的配置為專(zhuān)門(mén)進(jìn)行科學(xué)運(yùn)算的 Tesla 雙顯卡。因此，本文結(jié)合 Zou 等[8]的架構(gòu)進(jìn)行改寫(xiě)并作為比較。通過(guò)比較各部分的計(jì)算時(shí)間，分析不同架構(gòu)下的性能，每個(gè)測(cè)試數(shù)據(jù)通過(guò) NVIDIA 公司提供的 Visual Pro fi ler 獲得，每個(gè)數(shù)據(jù)在測(cè)試 10 次后取平均值。如表1 所示，Zou 等[8]的架構(gòu)與本文的架構(gòu)時(shí)間上相差不大，但本文所用架構(gòu)顯存的消耗不到 100 MB，大大地減少顯存的消耗。因此，在實(shí)際應(yīng)用上提供極大的方便，如實(shí)際應(yīng)用中往往需要額外的顯存空間進(jìn)行校正處理等。

本文架構(gòu)可以將重建層劃分不同塊數(shù)，因此，第一塊重建完成后可以直接輸出到用戶界面上。在本次試驗(yàn)中，16 層為一組，僅需 1.523 s便可從界面上獲得重建數(shù)據(jù)。BPF 算法的另一個(gè)優(yōu)點(diǎn)在于其感興趣重建，若重建 512×512×256數(shù)據(jù)則只需要 5.609 s。

表1 本文結(jié)果與最新研究結(jié)果比較Table1. Result comparison between this paper and the latest outcomes

表2 不同大小的重建結(jié)果比較Table2. Comparison of reconstruction results based on different size

此外，表2 為新架構(gòu)在不同重建大小下的性能比較，其中 1024×1024×128 與512×512×512 是重建空間大小一致的比較。根據(jù)結(jié)果所示，本方案能在不改變架構(gòu)下適用于不同大小的圖像重建，不需要根據(jù)不同大小而改變重建策略。

最后，針對(duì) 3.2.2 所述，實(shí)驗(yàn)測(cè)得新方案能把重建 512×512×512 的顯存消耗由 2.5 GB 減少到 100 MB 以下，證明其適合大數(shù)據(jù)的重建。

5 結(jié) 語(yǔ)

隨著 GPU 在科學(xué)計(jì)算的廣泛使用，CT 重建算法計(jì)算耗時(shí)長(zhǎng)的問(wèn)題得以解決。本文提出一種更優(yōu)于以往研究結(jié)果的基于 CUDA 加速的優(yōu)化架構(gòu)。本文設(shè)計(jì)的新架構(gòu)采取以下策略：(1)合理的顯存利用，在充分利用 CUDA 核心的前提下，對(duì)顯存進(jìn)行復(fù)用；(2)利用 BPF 算法的特點(diǎn)進(jìn)行速度優(yōu)化；(3)對(duì)常用數(shù)據(jù)進(jìn)行預(yù)先計(jì)算。實(shí)驗(yàn)結(jié)果顯示，只需要不到 100 MB 的顯存，實(shí)現(xiàn)重建 512×512×512 數(shù)據(jù)只需 8.055 s，第一塊數(shù)據(jù)塊輸出僅需 1.523 s，而感興趣區(qū)域重建為 4.566 s。實(shí)驗(yàn)證明，新的架構(gòu)在不損失 BPF 算法原有精度和保證速度大幅度降低的前提下，只需要較少的顯存空間進(jìn)行重建，為實(shí)際應(yīng)用提供一個(gè)高效的計(jì)算框架。

[1] Jung KJ, Lee KS, Kim SY, et al. Low-dose, volumetric helical CT: image quality, radiation dose, and usefulness for evaluation of bronchiectasis [J]. Investigative Radiology, 2000, 35(9): 557-563.

[2] Feldkamp LA, Davis LC, Kress JW. Practival conebeam algorithm [J]. Journal of the Optical Society of America A: Optics, Image Science, and Visions, 1984, 1(6): 612-619.

[3] Pan XC, Zou Y, Xia D. Image reconstruction in peripheral and central regions-of-interest and data redundancy [J]. Medical Physics, 2005, 32(3): 673-684.

[4] Cho S, Bian J, Pelizzari CA, et al. Region-ofinterest image reconstruction in circular cone-beam microCT [J]. Medical Physics, 2007, 34(12): 4923-4933.

[5] 韓玉, 閆鑌, 宇超群, 等. 錐束 CT FDK 重建算法的 GPU 并行實(shí)現(xiàn) [J]. 計(jì)算機(jī)應(yīng)用, 2012, 32(5): 1407-1410.

[6] 雷德川, 陳浩, 王遠(yuǎn), 等. 基于 CUDA 的多 GPU加速 SART 迭代重建算法 [J]. 強(qiáng)激光與粒子束, 2013, 25(9): 2418-2422.

[7] Zheng H, Yu YY, Kang Y, et al. Investigation on PI-line selecting method based on GPU accelerated backprojection fi ltered VOI reconstruction [C] // Medical Imaging 2010: Physics of Medical Imaging, 76222G, doi:10.1117/12.839831.

[8] Zou J, Chen H, Zhang QY, et al. Fast Cone-beam CT image reconstruction based on BPF algorithm: application to ortho-CT [J]. International Journal of Computational Methods, 2013, doi: 10.1142/ S0219876213500679.

[9] Yu LF, Zou Y, Sidky EY, et al. Region of interest reconstruction from truncated data in circular cone-beam CT [J]. IEEE Transactions on Medical Imaging, 2006, 25(7): 869-881.

[10] NVIDIA Corporation. NVIDIA CUDA Programming C Guide, Version 5.5. 2013 [OL]. http://www.nvidia.cn/object/cuda-cn.html.

[11] 張舒, 褚艷利, 趙開(kāi)勇, 等. GPU 高性能運(yùn)算之CUDA [M]. 北京: 中國(guó)水利水電出版社, 2009. Science, 2012, 336(6085): 1171-1174.

CUDA Based Parallel Implementation of BPF Algorithm

WU Shaojia1CHEN Hao2LIAO Li1GUI Jianbao2
1( Zhaoqing Radio & Television University, Zhaoqing 526060, China )
2(Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, Shenzhen 518055, China )

Based on its region-of-interest (ROI) reconstruction advantage, backprojection-filter algorithm has been used in cone-beam CT recently. However, because of its complexity and computation, there is memory insufficiency in the implementation of GPU acceleration. Hence, CUDA-based parallel implementation for BPF algorithm was proposed. Meanwhile, an accelerated projection scheme and other accelerated techniques were included as well, such as features of BPF for acceleration. Besides, video memory pool was introduced to optimize the implementation. With an ef fi cient structure, the simulation results show that it takes only 8.055 seconds by using the new structure to reconstruct 512×512×512 data and only 4.566 seconds for the ROI reconstruction. The output of fi rst block data takes only 1.523 seconds. With a great decrease of memory occupation from 2.5 GB to less than 100 MB, the new scheme is suitable for big data reconstruction.

backprojection-filter algorithm; cone-beam computed tomography; region of interest imaging; graphics processing unit; image reconstruction; parallel computing architecture

TP 391.41

A

2014-06-18

國(guó)家科技支撐計(jì)劃(2012BAI13B04)；國(guó)家自然科學(xué)基金(61102161)；深圳市項(xiàng)目(JCYJ20130401170306796, JC201105190923A)

伍紹佳，碩士，講師，研究方向?yàn)橛?jì)算機(jī)應(yīng)用與網(wǎng)絡(luò)技術(shù)；陳皓，碩士研究生，研究方向?yàn)?CT 重建與加速算法；廖麗，碩士，講師，研究方向?yàn)檐浖夹g(shù)應(yīng)用；桂建保(通訊作者)，博士，高級(jí)工程師，研究方向?yàn)?X 射線與 CT 成像技術(shù)，E-mail：jb.gui@siat.ac.cn。