基于射線驅(qū)動(dòng)的Siddon正投影算法并行計(jì)算研究

陳云斌，劉清華，李壽濤，王 遠(yuǎn)

(1.中國工程物理研究院 應(yīng)用電子學(xué)研究所，綿陽621900；2.國家X射線數(shù)字化成像儀器中心，綿陽 621000)

0 引言

計(jì)算機(jī)斷層掃描成像(Computed Tomography，CT)的實(shí)際執(zhí)行過程分為正投影和反投影兩大步驟。正投影指計(jì)算給定函數(shù)沿指定射線方向的線積分過程。反投影指由投影數(shù)據(jù)進(jìn)行圖像重建，獲取切片體數(shù)據(jù)。

迭代重建算法是CT重建算法的重要分支，多用于有限角、稀疏投影重建，包括SART[1](Simultaneous Algebraic Reconstruction Techniques)、EM[2](Expectation Maximization)和TV[3](Total Variation)等算法。迭代重建算法重復(fù)了CT成像的兩大步驟，即正投影和反投影。與實(shí)際物理成像的區(qū)別在于，迭代重建算法的正投影步驟是在數(shù)字圖像域，沿指定路徑對圖像矩陣的像素灰度值進(jìn)行線積分。

計(jì)算機(jī)模擬仿真是CT成像研究的重要技術(shù)手段，廣泛用于重建算法驗(yàn)證、掃描成像參數(shù)驗(yàn)證和物理現(xiàn)象的跟蹤分析等領(lǐng)域。計(jì)算機(jī)模擬仿真方法大致分為兩類：統(tǒng)計(jì)仿真方法和解析仿真方法。統(tǒng)計(jì)仿真方法的典型例子就是蒙特卡羅方法。解析仿真方法借助生成的數(shù)字模體，通過解析公式計(jì)算相關(guān)參量，其中典型的例子就是給定數(shù)字模體的圖像矩陣，通過正投影算法模擬生成投影數(shù)據(jù)。

綜上所述，正投影算法和計(jì)算策略是迭代重建算法和計(jì)算機(jī)模擬仿真的關(guān)鍵技術(shù)手段，在CT成像研究領(lǐng)域具有重要作用。

1 正投影算法概述

常見的正投影算法大體包括如下幾類，即Siddon算法[4]、Joseph算法[5]和等距采樣法[6]，下面逐一進(jìn)行簡要介紹。

Siddon算法的核心思想是計(jì)算射線L與單位像素f(i,j),0≤i

圖5-1 Siddon正投影算法示意圖

(1)

Joseph算法的核心思想是通過射線的端點(diǎn)坐標(biāo)(x1,y1)、(x2,y2)，并比較分量|x1-x2|、|y1-y2| 的大小確定該射線的主傳播方向。比如，|x1-x2|>|y1-y2|，則X方向?yàn)樯渚€的主傳播方向。Joseph算法只需計(jì)算射線在主傳播方向上的采樣點(diǎn)，如圖2所示。采用Joseph算法得到的正投影值為：

圖2 Joseph正投影算法示意圖

(2)

式中，f′(i,j)為采樣點(diǎn)經(jīng)過雙線性插值得到的像素灰度值。Joseph算法計(jì)算量小，但計(jì)算精度有所損失。

等距采樣法，顧名思義，沿著射線傳播方向等間距采樣，如圖3所示。逐漸累加得到射線穿過整個(gè)圖像矩陣的正投影值，如式(3)所示。

圖3 等距采樣法正投影算法示意圖

(3)

式中，Δl為采樣間距；f′k(i,j)為采樣點(diǎn)經(jīng)過雙線性插值得到的像素灰度值。

綜上所述，Joseph算法計(jì)算量小，但計(jì)算精度存在損失。等距采樣法計(jì)算精度由采樣間距決定，采樣間距大，能夠降低計(jì)算量，但計(jì)算精度也會下降。根據(jù)正投影算法的原理，Siddon算法是一種精確算法，計(jì)算量較大。實(shí)際應(yīng)用中，等距采樣法和Siddon算法應(yīng)用廣泛，是迭代重建領(lǐng)域主流的正投影技術(shù)手段。為了提高Siddon算法的計(jì)算效率，本文主要研究基于GPU(Graphics Processing Unit)的Siddon算法并行計(jì)算策略。

2 Siddon算法的串行計(jì)算策略

首先，對掃描物體進(jìn)行離散化，形成n行n列的圖像矩陣，像素灰度為f(i,j),0≤i

圖4 Siddon算法原理

(4)

為了便于描述，對圖像矩陣中各單位像素的坐標(biāo)進(jìn)行離散化，如下式所示：

(5)

式中，xi+1-xi=1，yi+1-yi=1。

2.1 確定射線穿過圖像矩陣的入射點(diǎn)和出射點(diǎn)

(6)

(7)

αstart=max(0,min(αx0,αxn),min(αy0,αyn))

(8)

出射點(diǎn)距射線起點(diǎn)S的距離占據(jù)射線總長|ST|的比例因子為αend，有：

αend=min(1,max(αx0,αxn),max(αy0,αyn))

(9)

2.2 確定射線穿過單位像素的入射點(diǎn)和出射點(diǎn)

(10)

αstart≤αxj,αyi≤αend

(11)

如圖4所示，射線與像素f(0,0)相交的入射點(diǎn)為A，該點(diǎn)為射線與邊界直線y=y0的交點(diǎn)，該點(diǎn)距射線起點(diǎn)S的距離占據(jù)射線總長|ST|的比例因子為αy0；出射點(diǎn)為B，該點(diǎn)為射線與邊界直線x=x1的交點(diǎn)，比例因子為αx1。截距|AB|的實(shí)際長度為：

|AB|=(αx1-αy0)|ST|

(12)

射線經(jīng)過像素f(0,0)的區(qū)間線積分為：

pf(0,0)=f(0,0)·(αx1-αy0)·|ST|

(13)

3 Siddon算法的并行計(jì)算策略

Siddon算法的計(jì)算量較大。對于n×n的圖像矩陣，每條射線的正投影值需要計(jì)算(n+1)×(n+1)個(gè)交點(diǎn)坐標(biāo)，且需對交點(diǎn)距射線起點(diǎn)的距離占據(jù)射線總長的比例因子進(jìn)行排序計(jì)算；如果探元數(shù)量為m，投影角度為v，則射線數(shù)量為v×m。合計(jì)需要計(jì)算v×m×(n+1)×(n+1)個(gè)交點(diǎn)坐標(biāo)，進(jìn)行v×m次排序計(jì)算。如果擴(kuò)展到錐束掃描的情況，射線數(shù)量和交點(diǎn)數(shù)量還會急劇增加，串行計(jì)算策略難以滿足實(shí)際工程需要。

根據(jù)Siddon算法原理，射線與射線之間的正投影計(jì)算獨(dú)立進(jìn)行，具備并行計(jì)算的前提條件。本文基于CUDA架構(gòu)，將圖像矩陣綁定至紋理存儲器，每個(gè)線程負(fù)責(zé)一條射線的正投影計(jì)算，實(shí)現(xiàn)Siddon正投影算法的GPU并行加速運(yùn)算。

原Siddon算法需要計(jì)算射線與各單位像素邊界直線的交點(diǎn)，并存儲交點(diǎn)距射線起點(diǎn)的距離占據(jù)射線總長的比例因子，然后進(jìn)行排序計(jì)算。然而，上述計(jì)算策略存在兩大缺點(diǎn)，不利于GPU并行計(jì)算，具體表現(xiàn)為：

(1)由于每個(gè)線程對應(yīng)一條射線的正投影計(jì)算，因此正投影計(jì)算所需的比例因子{αxj,αyi}作為私有變量需要由每個(gè)線程開辟存儲空間進(jìn)行獨(dú)立維護(hù)。一般情況下，線程私有變量會被分配到延遲訪問極低的寄存器中。但是，比例因子{αxj,αyi}數(shù)據(jù)量大，寄存器空間不足，最終會被分配到延遲訪問較高的顯存中。例如，圖像矩陣為1024×1024，單條射線維護(hù)的交點(diǎn)數(shù)量為1025×1025，需要的存儲空間為1 M。若同時(shí)計(jì)算的射線數(shù)量為1024，則需要的存儲空間為1 G。因此，原Siddon算法在進(jìn)行并行計(jì)算時(shí)需耗費(fèi)大量的存儲空間。

(2)每個(gè)線程需要對比例因子{αxj,αyi}進(jìn)行排序計(jì)算，增加額外的計(jì)算量。

基于上述原因，本文對Siddon算法進(jìn)行計(jì)算策略優(yōu)化，使其利于GPU并行計(jì)算，提高算法執(zhí)行效率。

3.1 確定射線穿過圖像矩陣的入射點(diǎn)和出射點(diǎn)

①考查邊界直線x=x0，該邊界線與射線的交點(diǎn)坐標(biāo)為：

(14)

如果x0≤x≤xn且y0≤y≤yn，則該交點(diǎn)作為射線穿過圖像矩陣的入射點(diǎn)，如圖5(1)。否則，繼續(xù)尋找入射點(diǎn)，如圖5(2)。

圖5 確定射線穿過圖像矩陣的入射點(diǎn)

②考查邊界直線x=xn，該邊界線與射線的交點(diǎn)坐標(biāo)為：

(15)

如果x0≤x≤xn且y0≤y≤yn，若入射點(diǎn)已存在，則該交點(diǎn)作為射線穿過圖像矩陣的出射點(diǎn)；如果x0≤x≤xn且y0≤y≤yn，若入射點(diǎn)不存在，則該交點(diǎn)作為射線穿過圖像矩陣的入射點(diǎn)。如果xxn或yyn，則繼續(xù)尋找入射點(diǎn)或出射點(diǎn)。

③同理，分別考查邊界直線y=y0、y=yn與射線的交點(diǎn)，確定入射點(diǎn)或出射點(diǎn)。

如果始終無法確定入射點(diǎn)，說明該射線與圖像矩陣不相交，沿該射線的正投影值為p=0。

3.2 確定射線穿過圖像矩陣的后繼點(diǎn)

如前文所述，Siddon算法的串行計(jì)算策略需要進(jìn)行排序計(jì)算，導(dǎo)致每條射線對應(yīng)的線程需要維護(hù)大量的私有變量，造成存儲資源的消耗和計(jì)算量的增加。為了克服上述缺點(diǎn)，首先，確定射線穿過某單位像素的起點(diǎn)后，應(yīng)直接定位后繼點(diǎn)，完成起點(diǎn)至后繼點(diǎn)區(qū)間內(nèi)的線積分運(yùn)算。然后，更新起點(diǎn)坐標(biāo)，重新定位后繼點(diǎn)，直到后繼點(diǎn)為射線穿過圖像矩陣的出射點(diǎn)。其中第一個(gè)起點(diǎn)為射線穿過圖像矩陣的入射點(diǎn)。

假設(shè)射線穿過某單位像素的起點(diǎn)坐標(biāo)為(xpre,ypre)，后繼點(diǎn)坐標(biāo)為射線穿過該單位像素的出射點(diǎn)(xnext,ynext)。根據(jù)Siddon算法原理，后繼點(diǎn)一定是射線與單位像素某邊界直線的交點(diǎn)。這里分別沿著射線方向的x分量和y分量搜索單位像素的邊界直線，確定后繼點(diǎn)。x分量方向的搜索起點(diǎn)為x=xj，其中xj滿足下述條件：

(16)

y分量方向的搜索起點(diǎn)為y=yi，其中yi滿足下述條件：

(17)

后繼點(diǎn)有兩種可能：①射線與單位像素邊界直線xnext=xj+sgn(xt-xs)·1的交點(diǎn)，其中sgn(·)為符號函數(shù)，該點(diǎn)距射線起點(diǎn)S的距離占據(jù)射線總長|ST|的比例因子αnext，其表達(dá)式為：

(18)

(19)

圖6 確定后繼點(diǎn)坐標(biāo)

從起點(diǎn)(xpre,ypre)到后繼點(diǎn)(xnext,ynext)，射線在該區(qū)間內(nèi)的線積分為：

f(xj,yi)·(αnext-αpre)·|ST|

(20)

更新起點(diǎn)坐標(biāo)，將當(dāng)前后繼點(diǎn)作為新的起點(diǎn)，繼續(xù)搜索下一個(gè)后繼點(diǎn)，直到搜索至射線穿過整個(gè)圖像矩陣的出射點(diǎn)，完成整個(gè)區(qū)間內(nèi)的線積分，得到該射線穿過圖像矩陣的正投影值。

Siddon算法的計(jì)算策略經(jīng)過改進(jìn)后，可以發(fā)現(xiàn)，每條射線不再維護(hù)龐大的私有變量數(shù)組，且算法迭代搜索后繼點(diǎn)，并完成局部區(qū)間內(nèi)的線積分計(jì)算，無需排序，便于采用GPU并行計(jì)算措施，提高算法執(zhí)行效率。

4 實(shí)驗(yàn)

為了驗(yàn)證Siddon算法GPU并行計(jì)算結(jié)果的有效性和準(zhǔn)確性，本文設(shè)計(jì)人體頜部數(shù)字模體作為測試對象，模體尺寸512×512，如圖7所示。

圖7 人體頜部數(shù)字模體

設(shè)計(jì)正投影掃描參數(shù)，如表1所示。

表1 掃描參數(shù)

采用前文提出的GPU并行計(jì)算策略，對頜部模體進(jìn)行正投影模擬計(jì)算，并得到圓周掃描范圍內(nèi)的扇束投影正弦圖，如圖8左圖。為了驗(yàn)證正投影模擬計(jì)算結(jié)果的有效性，采用標(biāo)準(zhǔn)扇束FBP重建算法，重建矩陣依然設(shè)計(jì)為512×512，重建體素尺寸為0.2 mm，得到模擬投影對應(yīng)的重建圖像，如圖8右圖?？梢园l(fā)現(xiàn)，由模擬投影能夠重建得到正確的斷層圖像，說明本文提出的GPU并行計(jì)算策略具備有效性。

圖8 頜部模體的扇束掃描正弦圖和重建圖像

為了驗(yàn)證GPU并行計(jì)算結(jié)果的準(zhǔn)確性，這里按照Siddon算法的串行計(jì)算策略采用CPU主機(jī)端進(jìn)行計(jì)算，按表1所示的掃描條件重新獲得頜部模體的扇束投影正弦圖。分別取同一投影角度下CPU串行計(jì)算結(jié)果和GPU并行計(jì)算結(jié)果的投影數(shù)據(jù)進(jìn)行比較，如圖9所示。

圖9 CPU串行計(jì)算結(jié)果和GPU并行計(jì)算結(jié)果比較

為了定量表征GPU并行計(jì)算結(jié)果的準(zhǔn)確性，以CPU串行計(jì)算結(jié)果為參考標(biāo)準(zhǔn)，引用歸一化平均絕對距離判據(jù)r對GPU并行計(jì)算結(jié)果進(jìn)行綜合評價(jià)。歸一化平均絕對距離判據(jù)r的表達(dá)式為：

(21)

經(jīng)計(jì)算，GPU并行計(jì)算結(jié)果r=0.31%。

綜上所述，CPU串行計(jì)算結(jié)果和GPU并行計(jì)算結(jié)果基本上保持一致，說明本文提出的GPU并行計(jì)算策略其計(jì)算結(jié)果具有準(zhǔn)確性。

CPU串行計(jì)算和GPU并行計(jì)算的計(jì)算效率如下：

計(jì)算機(jī)CPU型號為Intel i5-6300HQ，主頻2.3 GHz，內(nèi)存4 G，Windows 7 64位操作系統(tǒng)。GPU型號為NVIDIA GeForce GTX 950 M，獨(dú)立顯存2 G。采用CPU串行計(jì)算耗時(shí)230 s，GPU并行計(jì)算完成同樣的計(jì)算量耗時(shí)僅需9 s，運(yùn)算速度提高25倍，運(yùn)算效率得到大幅度提升。

5 結(jié)論

Siddon正投影算法原串行計(jì)算策略存在兩大缺點(diǎn)：①每條射線的正投影計(jì)算都需要維護(hù)龐大的私有變量數(shù)組，消耗存儲資源；②每條射線的正投影計(jì)算都需要進(jìn)行排序計(jì)算，增加額外的計(jì)算量。上述串行計(jì)算策略不利于基于多線程同步的并行計(jì)算。針對這種情況，本文提出Siddon正投影算法的并行計(jì)算策略。沿著射線方向順序搜索、更新射線穿過單位像素的起點(diǎn)和后繼點(diǎn)，逐漸累加射線穿過單位像素的線積分值，得到指定射線穿過圖像矩陣的正投影值。

Siddon算法經(jīng)過計(jì)算優(yōu)化后，采用GPU實(shí)現(xiàn)扇束CT投影的多線程并行計(jì)算。實(shí)驗(yàn)結(jié)果表明，Siddon算法的并行計(jì)算結(jié)果在保證計(jì)算準(zhǔn)確度的同時(shí)，獲得20倍以上的計(jì)算效率提升，顯著提升Siddon算法的工程實(shí)用價(jià)值。擴(kuò)展到三維錐束掃描及重建的情況，Siddon算法的GPU并行計(jì)算將顯得尤為必要。