PET斷層重建中動態射線追蹤算法的實現

張 斌, 王李栓， 趙書俊

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

0 引言

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

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

迭代重建過程需要借助SRM，頻繁在圖像域和投影域之間進行變換運算，其執行流程如圖1所示.當變換運算沿著一條射線路徑進行時，稱為“線驅動的前向投影和背投影運算”.前向投影需要計算每一個像素對投影值的貢獻.背投影利用之前計算的校正因子，對估算的像素值進行修正.這兩個計算過程都要反復用到像素的權重因子，需要占用大量的系統資源，是重建算法中的關鍵步驟和最耗時的部分.

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

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

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

本文面向list-mode斷層重建，實現了一種敏捷的實時計算系統矩陣的方法，并融合在前向投影過程中.

1 正交距離射線追蹤算法原理與實現

該算法以優化Siddon算法為基礎，考慮PET探測特性，生成近似高斯形的LOR線，有效建模了探測器響應函數.

1.1 Siddon算法的基本原理

Siddon算法在圖像斷層重建中得到廣泛應用，該算法是找到被LOR線穿越的全部像素并計算LOR線在每個像素中的截斷長度，用截斷長度代表像素相對于該LOR線的權重因子.它的時間復雜度為O(n)，而直接的射線追蹤方法時間復雜度為O(n3).它首先對視場空間像素進行有序化標記，再對LOR線進行參數化表達，然后找到LOR線進入視場空間的入射點和出射點，從而確定射線追蹤順序.這些操作構成Siddon算法的實現基礎[8].

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

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

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

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

需要注意的是，在實踐中可以根據需要，靈活建立參考坐標系，不會影響最終計算結果，例如，可以將坐標原點設定在視場中心位置.像素空間中像素的排序以坐標系為參照，程序自定義.

1.2 優化的Siddon算法

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

優化算法的主要計算步驟如下：

1)初始化.初始化計算LOR線穿越像素空間時λ的取值范圍，即最小值λmin和最大值λmax；初始化計算LOR射線進入圖像空間后被x、y方向軸對齊等間距平行線的截斷初值λx、λy；初始化計算被射線擊中的第一個像素的索引坐標(ix,jy)，其中1≤ix≤Nx，1≤jy≤Ny.

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

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

4)根據上一像素索引計算下一像素索引.然后返回第2)步，循環執行，直到滿足結束射線追蹤條件.

1.3 正交距離射線追蹤算法的實現

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

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

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

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

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

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

2 實驗裝置與材料

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

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

3 List-mode重建算法

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

4 結果分析與結論

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

圖4 Eplus-166結構示意圖Fig.4 Schematic for Eplus-166 structure

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

從重建結果可以看出，當設置cδ為合適值時，內含的線性核能夠更好地近似高斯模糊函數，有效建模探測器響應函數，提高系統矩陣的計算精度，對于分辨率恢復起到決定性作用.實踐中發現，aij的下限越小，概率計算精度越高，但計算強度迅速增加.因此，aij并非越小越好，要根據重建目標，合理選擇aij的下限值，達到圖像重建質量和時間之間的平衡.

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

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

參考文獻：

[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] 贠明凱．正電子發射斷層成像中的三維迭代重建研究[D]．北京：中國科學院高能物理研究所，2009．

[12] 劉琳，何劍鋒，王紅玲．GPU加速數據挖掘算法的研究[J]．鄭州大學學報：理學版，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].鄭州大學學報：理學版，2011，43(1)：33-36.