張 斌, 王李栓， 趙書俊

(1.鄭州大學(xué) 物理工程學(xué)院 河南 鄭州 450001； 2.鄭州市節(jié)能監(jiān)察中心 河南 鄭州 450006)

0 引言

正電子發(fā)射斷層掃描儀(positron emission tomography, PET)利用放射性核素標(biāo)記化合物作為分子探針，非侵入式地獲取生物活體的功能和解剖影像，其功能成像是當(dāng)前所有成像模式中靈敏度最高的，是獲得生物活體的生理生化信息的強(qiáng)有力工具，近年來在諸多科學(xué)研究領(lǐng)域得到廣泛應(yīng)用.在臨床醫(yī)學(xué)中，已經(jīng)證明PET在癌癥診斷和分期、神經(jīng)系統(tǒng)疾病評估、心臟病診療、放療計(jì)劃制定以及化療監(jiān)測等方面起到重要作用.在生物醫(yī)學(xué)工程中，已經(jīng)廣泛地利用小動物PET進(jìn)行基因傳送與表達(dá)、細(xì)胞販運(yùn)、信號傳輸通路中的蛋白質(zhì)相互作用等方面的研究[1].

在PET掃描過程中，不改變活體對象的生理狀態(tài)，向其注射放射性示蹤劑參與代謝.示蹤劑標(biāo)記物是構(gòu)成有機(jī)體的基本元素的正電子同位素，如11C、13N、15O、18F或其他核素.示蹤劑標(biāo)記物衰變發(fā)射正電子，發(fā)生正電子湮沒反應(yīng)，產(chǎn)生逆向發(fā)射的511 keV的γ光子對.使用符合探測技術(shù)探測成對出現(xiàn)的γ光子對，確定符合響應(yīng)線(line of response, LOR)，通過數(shù)據(jù)采集得到大量LOR，加以校正和濾波后，進(jìn)行圖像斷層重建，產(chǎn)生示蹤劑的時(shí)間序列濃度分布影像，最后得到活體代謝功能影像.圖像斷層重建算法是PET設(shè)備的關(guān)鍵組成部分，是決定PET設(shè)備性能與成像質(zhì)量的重要環(huán)節(jié).目前主要包括迭代重建算法與解析重建算法，兩者相比，迭代重建算法能夠提供更好的圖像重建質(zhì)量，已成為PET的標(biāo)準(zhǔn)重建方法[2].

迭代重建過程需要借助SRM，頻繁在圖像域和投影域之間進(jìn)行變換運(yùn)算，其執(zhí)行流程如圖1所示.當(dāng)變換運(yùn)算沿著一條射線路徑進(jìn)行時(shí)，稱為“線驅(qū)動的前向投影和背投影運(yùn)算”.前向投影需要計(jì)算每一個(gè)像素對投影值的貢獻(xiàn).背投影利用之前計(jì)算的校正因子，對估算的像素值進(jìn)行修正.這兩個(gè)計(jì)算過程都要反復(fù)用到像素的權(quán)重因子，需要占用大量的系統(tǒng)資源，是重建算法中的關(guān)鍵步驟和最耗時(shí)的部分.

圖1 圖像域和投影域之間的變換示意圖Fig.1 Diagram of the transformations between the image space and the projection space

SRM由諸多物理因素決定，包括掃描儀幾何效應(yīng)、隨機(jī)符合、散射、被掃描對象內(nèi)衰減、正電子行程、非共線性、晶體穿透、晶體內(nèi)散射等.任何情況下，幾何效應(yīng)都是主要因素.在PET中，巨量LOR線和像素空間中同樣數(shù)量龐大的像素組合所產(chǎn)生的SRM非常巨大，無論計(jì)算、存儲和調(diào)用都是迭代重建算法中的瓶頸.SRM可以離線計(jì)算并存儲在磁盤上，在重建過程中調(diào)入內(nèi)存使用.但由于過于龐大，通常不能一次全部裝入內(nèi)存，需要反復(fù)從磁盤上讀取，成為制約重建速度的一個(gè)主要因素.SRM也可以在重建過程中動態(tài)計(jì)算，即所謂on the fly模式，但計(jì)算強(qiáng)度非常大.PET中采集的數(shù)據(jù)主要有正弦圖數(shù)據(jù)和list-mode數(shù)據(jù)兩種組織模式.list-mode數(shù)據(jù)是以列表的模式逐一記錄每個(gè)湮沒事件的LOR線空間位置、發(fā)生時(shí)間、能量以及其他相關(guān)信息.正弦圖數(shù)據(jù)是通過對采集到的list-mode數(shù)據(jù)按照LOR線分類重組而獲取.

目前，基于正弦圖的迭代斷層重建算法在PET中得到廣泛應(yīng)用.基于list-mode數(shù)據(jù)的重建算法具有許多獨(dú)特優(yōu)勢，近年來，受到人們的高度關(guān)注，處于快速發(fā)展和完善過程中[3-5].對于正弦圖重建，SRM可以預(yù)先計(jì)算并存儲在硬盤上，也可以采用動態(tài)計(jì)算模式.而基于list-mode的斷層重建，迭代計(jì)算是一個(gè)事件接一個(gè)事件進(jìn)行，實(shí)時(shí)計(jì)算系統(tǒng)矩陣，可以有效降低數(shù)據(jù)存取時(shí)間和存儲負(fù)擔(dān)，提高斷層重建的速度.實(shí)時(shí)計(jì)算時(shí)，為了提高計(jì)算效率，主要考慮幾何效應(yīng)，并選用適當(dāng)?shù)母咚购?，圖像質(zhì)量幾乎不受影響[6-7].

本文面向list-mode斷層重建，實(shí)現(xiàn)了一種敏捷的實(shí)時(shí)計(jì)算系統(tǒng)矩陣的方法，并融合在前向投影過程中.

1 正交距離射線追蹤算法原理與實(shí)現(xiàn)

該算法以優(yōu)化Siddon算法為基礎(chǔ)，考慮PET探測特性，生成近似高斯形的LOR線，有效建模了探測器響應(yīng)函數(shù).

1.1 Siddon算法的基本原理

Siddon算法在圖像斷層重建中得到廣泛應(yīng)用，該算法是找到被LOR線穿越的全部像素并計(jì)算LOR線在每個(gè)像素中的截?cái)嚅L度，用截?cái)嚅L度代表像素相對于該LOR線的權(quán)重因子.它的時(shí)間復(fù)雜度為O(n)，而直接的射線追蹤方法時(shí)間復(fù)雜度為O(n3).它首先對視場空間像素進(jìn)行有序化標(biāo)記，再對LOR線進(jìn)行參數(shù)化表達(dá)，然后找到LOR線進(jìn)入視場空間的入射點(diǎn)和出射點(diǎn)，從而確定射線追蹤順序.這些操作構(gòu)成Siddon算法的實(shí)現(xiàn)基礎(chǔ)[8].

2D情況下，用2組分別垂直于x、y軸的平行等間距直線(3D時(shí)，用3組分別垂直于x、y、z軸的平行等間距平面)，把視場空間劃分成Nx×Ny的像素陣列，實(shí)現(xiàn)對視場空間的有序標(biāo)記，如圖2所示.對任一像素v(x,y)，可用(i,j)標(biāo)記，其中1≤i≤Nx，1≤j≤Ny.

假如P1(x1,y1)、P2(x2,y2)代表2個(gè)探測器的中心點(diǎn)，連接2中個(gè)心點(diǎn)的射線代表一條LOR線，可參數(shù)化表示為：x(α)=x1+α(x2-x1);y(α)=y1+α(y2-y1),其中α∈[0, 1]，表示LOR線上任意一點(diǎn)到起始點(diǎn)的距離與LOR線總長度的比值.參數(shù)方程實(shí)際上分別代表LOR線在x、y軸上的投影線.

沿著P1到P2，或者P2到P1的方向，分別計(jì)算2個(gè)投影線與2組平行線的交點(diǎn)，包括入射點(diǎn)和出射點(diǎn)，剔除重復(fù)項(xiàng)，并按照升序排列，得到參數(shù)化表示的與LOR線相交的一系列像素{α(0),α(1),…,α(n)}.

被LOR線穿越的第m個(gè)像素在像素空間中的索引，可通過i(m)=?x1+αmid(x2-x1)」和j(m)=?y1+αmid(y2-y1)」計(jì)算求出.

需要注意的是，在實(shí)踐中可以根據(jù)需要，靈活建立參考坐標(biāo)系，不會影響最終計(jì)算結(jié)果，例如，可以將坐標(biāo)原點(diǎn)設(shè)定在視場中心位置.像素空間中像素的排序以坐標(biāo)系為參照，程序自定義.

1.2 優(yōu)化的Siddon算法

文[9]通過對Siddon算法的優(yōu)化，將射線追蹤的計(jì)算時(shí)間減少2/3.在優(yōu)化算法中，首先求第一個(gè)被射線擊中的像素的索引坐標(biāo)，然后通過加、減運(yùn)算，依次求出射線路徑上其余的像素索引坐標(biāo)，把大量乘法運(yùn)算簡化為加減運(yùn)算.引入一個(gè)新的參量λ替代α，參數(shù)化公式轉(zhuǎn)化為：x(λ)=x1+(λ/L)·(x2-x1),y(λ)=y1+(λ/L)·(y2-y1),其中，λ表示當(dāng)前點(diǎn)到LOR線起點(diǎn)的距離.LOR線在像素中的截?cái)嚅L度l(m)=λ(m)-λ(m-1).

優(yōu)化算法的主要計(jì)算步驟如下：

1)初始化.初始化計(jì)算LOR線穿越像素空間時(shí)λ的取值范圍，即最小值λmin和最大值λmax；初始化計(jì)算LOR射線進(jìn)入圖像空間后被x、y方向軸對齊等間距平行線的截?cái)喑踔郸藊、λy；初始化計(jì)算被射線擊中的第一個(gè)像素的索引坐標(biāo)(ix,jy)，其中1≤ix≤Nx，1≤jy≤Ny.

2)找到λx和λy中的最小者，設(shè)λξ=min(λx，λy).如果λξ>λmax，則射線追蹤結(jié)束.否則計(jì)算LOR線在當(dāng)前像素中的橫斷長度l=λξ-λ.

3)按照單位像素大小，步進(jìn)增加λξ值，λξ=λξ+L/(ξ2-ξ1),ξ或者是x或者是y.

4)根據(jù)上一像素索引計(jì)算下一像素索引.然后返回第2)步，循環(huán)執(zhí)行，直到滿足結(jié)束射線追蹤條件.

1.3 正交距離射線追蹤算法的實(shí)現(xiàn)

用Siddon優(yōu)化算法生成的LOR射線較細(xì)，不能夠很好匹配有一定寬度的探測器.因此，具有一定寬度的LOR線(3D時(shí)，用響應(yīng)管(tube of response, TOR))建模探測器之間的響應(yīng)會更加精確.正交距離射線追蹤法以Siddon優(yōu)化算法為基礎(chǔ)，生成近似高斯形的LOR線，可有效建模探測器響應(yīng)函數(shù)，計(jì)算公式為aij=1-dij/(cδ),其中dij表示像素幾何中心到LOR線的距離，cδ表示歸一化因子.δ與PET系統(tǒng)點(diǎn)擴(kuò)展函數(shù)的半高寬相關(guān)，c是經(jīng)驗(yàn)調(diào)整系數(shù),aij表示權(quán)重因子.

在實(shí)時(shí)計(jì)算系統(tǒng)矩陣時(shí)，首先用Siddon算法，找到被LOR線穿越的像素.然后，設(shè)置aij的最小閾值，并結(jié)合cδ值，計(jì)算LOR穿過影像空間時(shí)，被LOR線直接穿越的像素及其鄰接像素的dij，求出全部關(guān)聯(lián)像素的權(quán)重因子.正確設(shè)置aij的最小閾值非常關(guān)鍵，需要考慮在計(jì)算時(shí)間和重建圖像品質(zhì)因數(shù)之間進(jìn)行折中.

可用矢量公式計(jì)算各像素的dij值[10].假設(shè)LOR線兩端晶體條中心點(diǎn)坐標(biāo)分別為P1(x1,y1)和P2(x2,y2)，求矢量r在v上的投影，也就是像素幾何中心P0到LOR線的距離，如圖3所示.

生成的LOR線效果如圖2所示.圖像空間左下角坐標(biāo)為(bx,by)；dx,dy分別表示x軸平行線和y軸平行線間的距離；Nx,Ny分別表示x軸平行線和y軸平行線的個(gè)數(shù)；dij表示像素幾何中心到LOR線的距離；P1P2確定一條LOR線，陰影部分表示用正交距離法求出的具有一定寬度的LOR線，每個(gè)像素的灰度表示其權(quán)重.

圖2 正交距離射線追蹤法原理示意圖Fig.2 Schematic of OD-RT method

圖3 矢量法計(jì)算像素中心到直線距離示意圖Fig.3 Schematic of calculating the distance from geometric centre of pixel to LOR using vector method

2 實(shí)驗(yàn)裝置與材料

Eplus-166是中國科學(xué)院高能物理研究所研制的我國第一臺擁有自主知識產(chǎn)權(quán)的小動物PET.該掃描儀由16個(gè)模塊組成，呈正十六邊型，每個(gè)模塊由2個(gè)block沿著軸向排列組成，每個(gè)block包含16×16個(gè)晶體條.晶體條的大小為1.9 mm×1.9 mm×10 mm，為了提高閃爍光收集效率，晶體條之間加入了反光材料，使得探測器單元沿軸向和切向斷面長度都為2.0 mm.系統(tǒng)總共形成32個(gè)探測器環(huán)，軸向長度為64 mm，其結(jié)構(gòu)如圖4所示[11].

實(shí)驗(yàn)采用自制的Derenzo仿體模型，模型內(nèi)各空心柱直徑分別為 1.4,1.6,1.9,2.2,2.5,3.0 mm，空心柱圓心間的距離為其直徑的2倍.空心柱內(nèi)均勻注入18F-FDG溶液，在Eplus-166中掃描，對采集的數(shù)據(jù)進(jìn)行歸一化、衰減以及散射校后，進(jìn)行l(wèi)ist-mode迭代重建.

3 List-mode重建算法

圖像重建采用S-LMEM (ordinary subsetized list-mode EM algorithm)算法，迭代公式為

4 結(jié)果分析與結(jié)論

本文以2D重建進(jìn)行算法驗(yàn)證，采用S-MLEM算法，使用正交距離射線追蹤方法以on-the-fly模式實(shí)時(shí)計(jì)算系統(tǒng)矩陣.把list-mode數(shù)據(jù)劃分為32個(gè)子集，每個(gè)子集約50 000個(gè)事例.Eplus-166中心視場的點(diǎn)擴(kuò)展函數(shù)半高寬約1.67 mm[11]，分別設(shè)cδ=3, 2, 1.5，當(dāng)cδ=2時(shí)，達(dá)到最佳效果，能夠清晰地區(qū)分出1.9 mm熱區(qū)，結(jié)果如圖5所示.

圖4 Eplus-166結(jié)構(gòu)示意圖Fig.4 Schematic for Eplus-166 structure

圖5 重建結(jié)果圖Fig.5 Comparison of image reconstruction result using S-MLEM

從重建結(jié)果可以看出，當(dāng)設(shè)置cδ為合適值時(shí)，內(nèi)含的線性核能夠更好地近似高斯模糊函數(shù)，有效建模探測器響應(yīng)函數(shù)，提高系統(tǒng)矩陣的計(jì)算精度，對于分辨率恢復(fù)起到?jīng)Q定性作用.實(shí)踐中發(fā)現(xiàn)，aij的下限越小，概率計(jì)算精度越高，但計(jì)算強(qiáng)度迅速增加.因此，aij并非越小越好，要根據(jù)重建目標(biāo)，合理選擇aij的下限值，達(dá)到圖像重建質(zhì)量和時(shí)間之間的平衡.

無論采用何種算法，PET斷層重建的計(jì)算強(qiáng)度都非常大.從list-mode算法的原理與實(shí)現(xiàn)來看，在一次迭代中某個(gè)像素值的更新不受其他像素更新的影響，整個(gè)計(jì)算過程特別適合采用多線程并行計(jì)算，進(jìn)行實(shí)時(shí)重建.CUDA技術(shù)非常適宜于實(shí)現(xiàn)多線程list-mode斷層重建[12-14]，在本研究中進(jìn)行了初步嘗試，重建速度明顯提高.

本文的研究表明，利用正交距離射線追蹤方法以on-the-fly模式實(shí)時(shí)計(jì)算系統(tǒng)響應(yīng)矩陣，進(jìn)行PET動態(tài)重建是可行的.這對于進(jìn)一步研究面向list-mode數(shù)據(jù)的低統(tǒng)計(jì)計(jì)數(shù)率(低放射性活度)情況下的斷層重建以及實(shí)現(xiàn)包含TOF信息的斷層重建都有非常重要的借鑒意義.

參考文獻(xiàn)：

[1] Hao Peng , Craig S L. Recent developments in PET instrumentation[J]. Current Pharmaceutical Biotechnology, 2010, 11(6):555-571.

[2] Reader A J, Zaidi H. The promise of new PET image reconstruction[J]. Physica Medica, 2008, 24:49-56.

[3] Barrett H H, White T, Parra L C. List-mode likelihood[J]. J Opt Soc Am A, 1997, 14(11):2914-2923.

[4] Parra L , Barrett H H. List-mode likelihood: EM algorithm and image quality estimation demonstrated on 2-D PET[J]. IEEE Trans Med Imaging, 1998, 17(2):228-235.

[5] Rahmin A, Lenox M, Reader A J, et al. Statistical list-mode image reconstruction for the high resolution research tomograph[J]. Phys Med Biol, 2004, 49:4239-4258.

[6] Aguiar P, Rafecas M, Ortuo J E, et al. Geometrical and Monte Carlo projectors in 3D PET reconstruction[J]. Med Phys, 2010, 37(11):5691-5702.

[7] Colas S. A fast tube of response ray-tracer[J]. Med Phys, 2006, 33(12):4744-4748.

[8] Jacobs F, Sundermann E, Sutter B D, et al. A fast algorithm to calculate the exact radiological path through a pixel or voxel space[J]. Journal of Computing and Information Technology, 1998, 6(1):89-94.

[9] Han Gouping, Liang Zhengrong, You Jiangsheng. A fast ray tracing technique for TCT and ECT studies[C]// IEEE Nuclear Science Symposium Conference Record. Washington, 1999:1515-1518.

[10] Weisstein E W. Point-line distance 2-dimensional[EB/OL].[2011-12-11]. http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html.

[11] 贠明凱．正電子發(fā)射斷層成像中的三維迭代重建研究[D]．北京：中國科學(xué)院高能物理研究所，2009．

[12] 劉琳，何劍鋒，王紅玲．GPU加速數(shù)據(jù)挖掘算法的研究[J]．鄭州大學(xué)學(xué)報(bào)：理學(xué)版，2010，42(2)：31-34．

[13] Pratx G, Chinn G, Habte F, et al. Fully 3-D list-mode OSEM accelerated by graphics processing units[C]// IEEE Nuclear Science Symposium Conference Record. San Diego CA, 2006:2196-2202.

[14] 王君，任永功.基于FP-tree 最大頻繁模式超集挖掘算法[J].鄭州大學(xué)學(xué)報(bào)：理學(xué)版，2011，43(1)：33-36.