一種基于cuda的s-bpf重建算法加速方法
【技術(shù)領(lǐng)域】
[0001]該發(fā)明涉及一種CT成像的圖像重建算法，特別是涉及一種基于⑶DA的S-BPF重建 算法加速方法。
【背景技術(shù)】
[0002]計(jì)算機(jī)斷層成像(Computed Tomography,CT)可以在無損的條件下同時(shí)獲得物體 的外部和內(nèi)部幾何結(jié)構(gòu)特征。在CT成像過程中，圖像重建由于其時(shí)間復(fù)雜度高的缺點(diǎn)，已成 為限制CT成像速度的關(guān)鍵因素。在醫(yī)療上，為準(zhǔn)確地分析病情，醫(yī)院通常會(huì)借助CT對(duì)病人內(nèi) 部組織進(jìn)行成像。然而由于重建時(shí)間長(zhǎng)，病人往往會(huì)忍受長(zhǎng)時(shí)間的病痛后才能得到準(zhǔn)確的 治療;在手術(shù)導(dǎo)航中，為了準(zhǔn)確定位病灶位置，常使用C-arm CT對(duì)病人進(jìn)行成像，該應(yīng)用中 對(duì)重建速度具有嚴(yán)格的要求。因此圖像重建耗時(shí)長(zhǎng)已嚴(yán)重降低了 CT成像的使用性。
[0003] 目前，使用高性能的硬件設(shè)備對(duì)重建算法進(jìn)行并行加速是解決其計(jì)算時(shí)間過長(zhǎng)的 常用方法。在眾多硬件加速器(如計(jì)算機(jī)集群、DSP、FPGA、GPU等）中，GPU具有并行計(jì)算能力 強(qiáng)、性價(jià)比高、易于開發(fā)等優(yōu)勢(shì)，在圖像重建中得到了廣泛的研究和應(yīng)用。
[0004] 選擇性反投影濾波算法(Selective-BPF，S-BPF)是一種用于短掃描局部重建的快 速重建算法。S-BPF算法主要分為投影數(shù)據(jù)的加權(quán)、求導(dǎo)、反投影和DBP圖像的濾波四個(gè)步 驟，其并行性分析如下：
[0005] 步驟一:投影數(shù)據(jù)的加權(quán)。該操作是對(duì)平板探測(cè)器獲取的投影進(jìn)行幾何修正，其計(jì) 算表達(dá)式如下：
[0006]
[0007] 其中，R為光源到旋轉(zhuǎn)軸的距離，ρ(β，πι，η)表示在β掃描角度下平板探測(cè)器上坐標(biāo) 為(m，n)點(diǎn)的探測(cè)單元采集到的投影值對(duì)應(yīng)的衰減系數(shù)，^為余弦加權(quán)系數(shù)， Pw(P，m，n)為加權(quán)后的投影。由加權(quán)公式可知：1)加權(quán)與掃描角度β無關(guān)，即不同角度下的投 影加權(quán)操作可以并行計(jì)算;2)加權(quán)僅與ρ(β，πι，η)及幾何參數(shù)有關(guān)，即同一角度下不同位置 的投影可以并行計(jì)算。因此投影數(shù)據(jù)加權(quán)操作具有良好的并行性。
[0008] 步驟二:加權(quán)后投影數(shù)據(jù)的求導(dǎo)。該操作是將投影沿著掃描軌跡方向求導(dǎo)，在圓軌 跡中為沿探測(cè)器水平方向求導(dǎo)。該操作是對(duì)投影的一種濾波操作，在時(shí)域中求導(dǎo)對(duì)應(yīng)在頻 域中乘以j ω。在實(shí)際離散投影數(shù)據(jù)中，我們使用差分代替微分，其計(jì)算公式如下：
[0009]
[0010] 其中dpixel為探測(cè)器的各像素尺寸，ρ\(β，πι，η)為求導(dǎo)后的投影數(shù)據(jù)。與投影數(shù)據(jù) 加權(quán)相似，不同角度下的求導(dǎo)操作可以并行計(jì)算，同一角度下不同位置的投影也可并行計(jì) 算。
[0011] 步驟三:求導(dǎo)后投影數(shù)據(jù)的反投影。反投影操作是使用求導(dǎo)后投影數(shù)據(jù)計(jì)算差分 反投影(DBP)圖像的過程，其計(jì)算公式如下：
[0012]
[0013] MERGEFORMAT (3)
[0014] 其中 b (X，y，z)為DBP 圖像，Xmin，Xmax，Ymin，Ymax，Z min，Zmax 分別為X，y，z 的最小值和最 大值，Φ和:π-φ為掃描起始掃描角度和終止掃描角度，且it = arcsin(ymin/R) ^(β,χ,γ,ζ)和η (P，X，y，z)分別為空間一點(diǎn)(x，y，z)在掃描角度為β時(shí)的投影位置，其計(jì)算公式為：
[0015]
[0016]對(duì)于單個(gè)體素的反投影過程又可以分解為以下四個(gè)步驟：
[0017] 1)計(jì)算被重建物體中該體素(x，y，z)沿旋轉(zhuǎn)軸旋轉(zhuǎn)β度后的新坐標(biāo)(Υ 顯 然，計(jì)算各體素旋轉(zhuǎn)后的位置可以在不同或同一角度下并行執(zhí)行。
[0018] 么：-Ζ?
[0019] 2)計(jì)算該點(diǎn)在探測(cè)器上錐束投影位置的坐標(biāo)索引(π/，!/ )。
[0020]
[0021]其中，D為光源到探測(cè)器的距離，（πιο,ηο)是中心射束投射位置的坐標(biāo)索引。該過程 包括幾何放大以及坐標(biāo)-索引的轉(zhuǎn)換兩項(xiàng)操作，計(jì)算各體素點(diǎn)的投影索引僅與體素旋轉(zhuǎn)后 的坐標(biāo)和若干固定的幾何參數(shù)有關(guān)，因此該過程具有良好的并行性。
[0022] 3)計(jì)算投影中索引為(π/ ,γ/ )處的投影值。由于計(jì)算得到的(π/ ,γ/ )并不一定為正 整數(shù)，所以需要對(duì)其進(jìn)行插值得到(π/ ,η')處的投影值。通常采用三線性插值，其計(jì)算公式 如下：
[0023] p7 ,η7 ) = (η-[η]) ((m-[m] )p7 ?(β, [m]+l, [n]+l) + ([m]+l-m)p/ ?(β, [m], [n] +1)) +
[0024] ([n]+l-n)((m-[m])p/w(P,[m]+l,[n]) + ([mj+l-m)?7ν(β, [m],[n]))
[0025] \*MERGEF0RMAT(7)
[0026] 其中，[a]表示不大于a的最大整數(shù)。顯然，各點(diǎn)的三線性插值操作是可并行的。
[0027] 4)計(jì)算體素(x，y，z)在掃描角度β下的反投加權(quán)系數(shù)并對(duì)體素進(jìn)行累加。由公式可 知，反投加權(quán)系數(shù)為
[0028]
[0029] 該加權(quán)系數(shù)與步驟3)得到的投影值進(jìn)行相乘后累加到體素(x，y，z)上。該過程中 權(quán)值的計(jì)算在不同角度或同一角度都是可并行的。該過程中的累加操作在同角度下是可并 行的，但是對(duì)某一體素同時(shí)反投不同角度的投影時(shí)，存在讀寫沖突問題，因此該操作不能 在不同角度下實(shí)現(xiàn)并行。
[0030] 步驟四：DBP圖像的濾波。該步驟分別在DBP圖像上沿各PI線上進(jìn)行有限Hilbert逆 變換。為了計(jì)算方便，我們將PI線選擇成平行于y軸的線段。該操作可以分為以下5個(gè)步驟：
[0031] 1)對(duì)PI線上的DBP圖像進(jìn)行加權(quán)
[0032]
[0033] 2)對(duì)1^&，7,2)沿7方向進(jìn)行卩卩丁
[0034] fbw(x, ω ,z) = fft[bw(x,y,z)]\*MERGEFORMAT (10)
[0035] 3)在頻域中對(duì)fbw進(jìn)行希Hilbert濾波
[0036] fVw(x，ω ,z)=-jXfbw(x, ω ,z)\*MERGEF0RMAT (11)
[0037] 4)對(duì)作\沿7方向進(jìn)行正卩丁
[0038] b/w(x,y,z) = ifft[bw(x, ω ,z)]\*MERGEF0RMAT (12)
[0039] 5)對(duì)1^\進(jìn)行加權(quán)
[0040]
[0041] 通過并行性分析，很容易知道以上5個(gè)步驟均可實(shí)現(xiàn)步驟內(nèi)部并行執(zhí)行，但步驟之間 必須嚴(yán)格按照以上流程執(zhí)行。而各PI線之間相互獨(dú)立，因此各PI線之間的執(zhí)行是可并行的。

【發(fā)明內(nèi)容】

[0042] 本發(fā)明克服了現(xiàn)有技術(shù)中，CT成像的圖像重建算法時(shí)間長(zhǎng)的問題，提供一種效果 明顯的基于CUDA的S-BPF重建算法加速方法。
[0043]含有如下步驟：
[0044]步驟1:從硬盤上讀取若干投影到內(nèi)存，再傳入顯存中并在CPU中計(jì)算用于有限 Hilbert逆變換的常數(shù)C;
[0045]步驟2:GPU使用已傳入的若干投影做求導(dǎo)-反投影操作，得到DBP圖像；
[0046] 步驟3 :GPU對(duì)步驟2計(jì)算得到的DBP圖像做有限Hilbert逆變換，同時(shí)將計(jì)算好的結(jié) 果傳回內(nèi)存，供應(yīng)用處理或使用。
[0047]所述步驟2中GPU執(zhí)行求導(dǎo)-反投影的特征包含以下步驟：
[0048] 2.1、對(duì)投影數(shù)據(jù)進(jìn)行加權(quán)和求導(dǎo)；
[0049] 2.2、將結(jié)果從全局存儲(chǔ)器拷貝到紋理存儲(chǔ)器中；
[0050] 2.3、執(zhí)行反投影操作。
[0051]所述步驟3中的有限Hilbert逆變換，需要將數(shù)據(jù)裝載到一個(gè)特定的數(shù)據(jù)類型 cufftComplex中，具體的執(zhí)行流程為：
[0052] 3.1、對(duì)PI線上的DBP數(shù)據(jù)加權(quán)并賦值到cufftComplex中；
[0053] 3.2、使用CUFFT實(shí)現(xiàn)快速傅里葉變換；
[0054] 3.3、在頻域進(jìn)行希爾伯特濾波；
[0055] 3.4、使用(：服?1'庫做逆傅里葉變換；
[0056] 3.5、對(duì)其進(jìn)行加權(quán)并賦值DBP圖像；
[0057] 3.6、將結(jié)果拷貝到CPU中。
[0058]所述步驟1的通過使用兩倍投影大小的內(nèi)存空間解除讀取投影操作和計(jì)算常數(shù)C 操作之間的耦合，再利用多線程技術(shù)實(shí)現(xiàn)兩者的并行執(zhí)行。
[0059]所述步驟2通過使用步驟1中兩倍投影大小的顯存空間解除內(nèi)存-顯存的投影數(shù)據(jù) 傳輸和GPU中求導(dǎo)-反投影操作之間的耦合，再利用CUDA中的異步技術(shù)實(shí)現(xiàn)兩者的并行執(zhí) 行。
[0060] 所述GPU執(zhí)行求導(dǎo)-反投影方法通過使用兩倍于步驟1中投影大小的顯存空間來解 除步驟2.1、步驟2.2、步驟2.3中的耦合，再利用CUDA中的多流技術(shù)和異步技術(shù)實(shí)現(xiàn)GPU中數(shù) 據(jù)傳輸和計(jì)算的并行執(zhí)行。
[0061] 所述GHJ執(zhí)行求導(dǎo)-反投影方法中，為了完全隱藏步驟2中數(shù)據(jù)傳輸時(shí)延，如CPU-GPU數(shù)據(jù)傳輸時(shí)延和GPU-GPU數(shù)據(jù)傳輸時(shí)延，根據(jù)GPU的性能和圖像重建的規(guī)模對(duì)反投影進(jìn) 行了分塊處理。
[0062]所述步驟3中的有限Hilbert逆變換方法中，為節(jié)約顯存，在一定圖像重建規(guī)模下， 根據(jù)GHJ計(jì)算和數(shù)據(jù)傳輸性能，對(duì)有限Hilbert逆變換進(jìn)行分塊處理，實(shí)現(xiàn)了完全隱藏?cái)?shù)據(jù) 傳輸時(shí)延。
[0063]所述步驟3中的有限Hilbert逆變換方法中，利用兩倍分塊大小的顯存空間保證兩 組有限Hilbert逆變換-回傳操作的獨(dú)立性，再利用CUDA中的多流技術(shù)和異步技術(shù)實(shí)現(xiàn)一組 有限Hilbert逆變換與另一組的結(jié)果回傳操作的并行執(zhí)行。
[0064]與現(xiàn)有技術(shù)相比，本發(fā)明基于CUDA的S-BPF重建算法加速方法具有以下優(yōu)點(diǎn)：目前 重建算法的GPU加速已經(jīng)取得了顯著的加速效果，通信時(shí)延成為了限制現(xiàn)有加速策略的瓶 頸。本發(fā)明基于選擇性反投影濾波算法的并行性分析及CPU-GP