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

伍紹佳陳 皓廖 麗桂建保

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

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

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

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

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

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

反投影濾波(Backprojection-Filter，BPF)算法憑借其可實(shí)現(xiàn)感興趣區(qū)域重建的優(yōu)點(diǎn)，近年來逐漸被應(yīng)用到錐束 CT 中。但是，由于算法的復(fù)雜性，實(shí)踐中存在耗時問題，同時其 GPU 加速的實(shí)現(xiàn)亦存在顯存不足等問題。因此，文章提出了一種基于 CUDA 的 BPF 并行加速算法。通過設(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 年第一臺 CT 安裝以來，CT 已經(jīng)經(jīng)歷了幾代的發(fā)展。近年來，錐束 CT 已經(jīng)成為醫(yī)療診斷和放療技術(shù)的重要工具，因此，如何提高錐束 CT 的性能是近年來研究的熱點(diǎn)之一。 針對 CT 的性能，劑量、圖像重建速度及圖像質(zhì)量是其重建的三個主要評判標(biāo)準(zhǔn)[1]。本文主要解決圖像重建的速度問題及其在 GPU 應(yīng)用過程中涉及的顯存問題。近十幾年來，解析法重建憑借其高效性得到大量的研究，與常規(guī)解析算法相比(如 FDK[2]等濾波反投影算法)，反投影濾波(Backprojection-Filter，BPF)算法[3,4]憑借利用較少量的投影數(shù)據(jù)重建出精確的圖像，得到學(xué)者的廣泛認(rèn)可。憑借感興趣區(qū)域掃描和重建，BPF 算法可實(shí)現(xiàn)降低劑量的目的，但其實(shí)現(xiàn)過程還存在運(yùn)算量大等問題。故本文提出一個基于 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ì)一個應(yīng)用于GPU 中通用 Pi 線選擇的方案，但是并沒有詳細(xì)地介紹如何利用 GPU 進(jìn)行優(yōu)化計(jì)算。此外，Zou等[8]曾提出基于 GPU 的 BPF 算法優(yōu)化方案，盡管時間得到大幅度的提升，但是其架構(gòu)并不是最優(yōu)，如在大數(shù)據(jù)問題上并不能很好地應(yīng)用或在實(shí)際應(yīng)用中存在顯存使用過多的問題；同時，由于其框架基于多 GPU，因而難以在商業(yè)項(xiàng)目中進(jìn)行應(yīng)用。除此之外，大部分相關(guān)研究也僅僅關(guān)注重建速度的優(yōu)化，并沒有考慮顯存的合理利用，因此，本文引入顯存池的概念去優(yōu)化顯存的管理。根據(jù) BPF 算法的特性，本文首先提出適合實(shí)際應(yīng)用的單 GPU 優(yōu)化框架；其次對 CUDA 實(shí)現(xiàn)過程中采用的優(yōu)化策略進(jìn)行總結(jié)；最后，本文將闡述如何利用 GPU 進(jìn)行 Shepp-Logon 模體投影數(shù)據(jù)的快速生成。

2 BPF 算法介紹

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

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

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

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

反投影部分可以描述為：

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

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

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

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

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

(3)對反投影后的數(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 作為輔助處理器，專門為協(xié)助 CPU 進(jìn)行通用科學(xué)和工程計(jì)算而設(shè)計(jì)。接下來，本文首先針對模擬數(shù)據(jù)的生成進(jìn)行介紹并提及其加速方法，隨后詳細(xì)介紹基于 CUDA 的 BPF 算法新框架，最后介紹在 CUDA 實(shí)現(xiàn)過程中采用的加速策略。

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

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

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

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

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

在 GPU 端，數(shù)據(jù)讀入后，以每條射線為一個線程，如探測器像素個數(shù)為 512×512，可以設(shè)512×512 的線程。不同射線劃分的段數(shù)并不一致，每個線程計(jì)算該射線線段組的累加值，而線段的累加值則通過把物體映射到三維紋理內(nèi)存上獲得。選擇三維紋理顯存有兩個原因：(1)由于此類顯存獲取為隨機(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)致顯存不足的問題。

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

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

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

3.2.1 基于 CUDA 的 BPF 框架

如上所述，許多現(xiàn)行的框架都選擇把每一個重建點(diǎn)的計(jì)算分配給每一個線程，如為512×512×512 的數(shù)據(jù)開辟 512×512×512 個線程去進(jìn)行反投影，但是此方法在 BPF 算法中存在嚴(yá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ù)公式對投影數(shù)據(jù)進(jìn)行求導(dǎo)；接著對重建層進(jìn)行分塊，根據(jù)重建的不同塊所需的投影數(shù)據(jù)，導(dǎo)入相應(yīng)的求導(dǎo)數(shù)據(jù)，并對該塊進(jìn)行反投影(包括計(jì)算該塊的反投影值和其邊界項(xiàng)值)；反投影后再進(jìn)行濾波，濾波后將該層數(shù)據(jù)和邊界項(xiàng)數(shù)據(jù)進(jìn)行累加，這樣一段數(shù)據(jù)塊重建完成并輸出到終端，遍歷完所有重建塊后，則完成整個重建。此流程的優(yōu)點(diǎn)如下：

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

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

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

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

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

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

根據(jù)圖4 所示，Zou 等[8]需要存儲整個重建空間的數(shù)據(jù)，并將其累加，其中包括圖像項(xiàng)和邊界項(xiàng)數(shù)據(jù)，過多地開辟顯存容易導(dǎo)致其他操作無法進(jìn)行。新框架下只需要一塊重建塊的顯存空間，可以省略大部分圖像項(xiàng)和邊界項(xiàng)數(shù)據(jù)的空間。以下為簡化的內(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。在濾波過程中會產(chǎn)生頻域數(shù)據(jù)，其長度一般為其最接近的 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由公式可以看出，在顯存的利用上，新框架有絕對的優(yōu)勢，新架構(gòu)所需空間僅為 Zou等[8]的 H/P，從而方便其他處理的進(jìn)行。若為512×512×512 的數(shù)據(jù)，以 16 層數(shù)據(jù)為一數(shù)據(jù)塊，所需顯存空間僅為原方案的 0.03125 倍，即需要不到 100 Mb 的顯存(考慮到其他輔助變量的存儲)。因此，在大數(shù)據(jù)重建下，新框架的優(yōu)勢尤其突出，使其在重建數(shù)據(jù)量大的情況下能保證程序正確性的同時不影響其他操作。

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

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

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

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

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

(1)反投影部分

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

(2)濾波部分

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

(3)輸入輸出部分

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

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

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

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

本文架構(gòu)可以將重建層劃分不同塊數(shù)，因此，第一塊重建完成后可以直接輸出到用戶界面上。在本次試驗(yàn)中，16 層為一組，僅需 1.523 s便可從界面上獲得重建數(shù)據(jù)。BPF 算法的另一個優(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ù)不同大小而改變重建策略。

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

5 結(jié) 語

隨著 GPU 在科學(xué)計(jì)算的廣泛使用，CT 重建算法計(jì)算耗時長的問題得以解決。本文提出一種更優(yōu)于以往研究結(jié)果的基于 CUDA 加速的優(yōu)化架構(gòu)。本文設(shè)計(jì)的新架構(gòu)采取以下策略：(1)合理的顯存利用，在充分利用 CUDA 核心的前提下，對顯存進(jìn)行復(fù)用；(2)利用 BPF 算法的特點(diǎn)進(jìn)行速度優(yōu)化；(3)對常用數(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)用提供一個高效的計(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] 張舒, 褚艷利, 趙開勇, 等. GPU 高性能運(yùn)算之CUDA [M]. 北京: 中國水利水電出版社, 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

國家科技支撐計(jì)劃(2012BAI13B04)；國家自然科學(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)用；桂建保(通訊作者)，博士，高級工程師，研究方向?yàn)?X 射線與 CT 成像技術(shù)，E-mail：jb.gui@siat.ac.cn。