基于ART算法的投影系數(shù)快速計算方法

侯慧玲，王明泉，楊 娟，李世虎

（中北大學(xué) 儀器科學(xué)與動態(tài)測試教育部重點(diǎn)實驗室，山西 太原030051）

0 引 言

CT技術(shù)是無損檢測領(lǐng)域廣泛使用的檢測手段，通過設(shè)計和搭建CT成像系統(tǒng)完成對被檢工件的掃描［1］.對采集到的投影進(jìn)行重建的算法主要有兩大類：解析算法和迭代算法.濾波反投影法（Filtered Back Projection，F(xiàn)BP）是最基本的解析算法，運(yùn)算量小，重建速度快，但對投影數(shù)據(jù)的完備性要求高.代數(shù)重建技術(shù)（ART）是經(jīng)典的代數(shù)迭代算法，在投影數(shù)據(jù)部分缺失的情況下仍然適用，且重建圖像質(zhì)量高，但其計算量大，收斂速度慢，限制了其在實時工業(yè)檢測系統(tǒng)中的應(yīng)用［2-3］.

影響ART算法重建速度的因素有很多，如：投影數(shù)據(jù)的訪問方式、圖像初值的選擇、松弛因子的選擇、投影系數(shù)的計算等等［4-6］.其中，投影系數(shù)的計算在時間上占比最大［7］，因而有不少學(xué)者關(guān)注投影數(shù)據(jù)計算方法的研究，不斷有新的方法提出［8-12］，其中最為熟知的是Siddon于1985年提出的基于長度加權(quán)的基本計算方法，計算射線穿過像素網(wǎng)格的長度，計算耗時比較長［8］；2006年張順利提出投影系數(shù)快速計算方法，通過增量計算來確定射線穿過的網(wǎng)格編號及交線長度，重建速度比Siddon算法快了將近6倍，但引入頻繁的分支判斷［9］；2012年楊文良提出的快速計算方法，使用簡單的加減法及比較運(yùn)算來完成投影系數(shù)計算，重建效率進(jìn)一步提高［11］.

本文針對ART重建中投影系數(shù)計算時間冗長的問題，從研究射線穿過網(wǎng)格的相交規(guī)律出發(fā)，提出新的計算方法，力求在保證投影重建圖像質(zhì)量的前提下，減少運(yùn)算量及分支判斷，進(jìn)一步提高計算速度，加快重建效率.

1 ART算法基本原理

式中：i=1，2，3，…，M；j=1，2，3，…，N；pi為射線i的投影值；pij為像素i對像素j的投影貢獻(xiàn)；wij為投影系數(shù)；xj為像素j的像素值，即該像素對X射線的衰減系數(shù).直接求解該線性方程組極其困難.ART算法是通過迭代方式來求解該方程組的最優(yōu)解的過程，其迭代公式［2-3］為

式中：k是迭代次數(shù)；wij為投影系數(shù)；λk為松弛因子，λ∈（0，2）.

計算其中的投影系數(shù)wij有三種數(shù)學(xué)模型：點(diǎn)加權(quán)型，面積加權(quán)型以及長度加權(quán)型.

長度加權(quán)模型如圖1所示.在模型中，射線被看作是一條直線穿過網(wǎng)格，投影系數(shù)wij定義為射線i穿過像素j的交線長度.

就重建精度而言，長度加權(quán)模型低于面積加權(quán)型，優(yōu)于點(diǎn)加權(quán)型.在保證一定重建精度的情況下，計算速度快，計算復(fù)雜度低，所以在實際中運(yùn)用最多.本文便是基于長度加權(quán)模型研究投影系數(shù)的快速計算方法.

圖1 長度加權(quán)模型 Fig.1 Length-weighted model

2 投影系數(shù)的快速計算

2.1 二維情況下的投影系數(shù)快速計算

投影系數(shù)矩陣?yán)锓橇阒档膫€數(shù)較少，可用稀疏矩陣進(jìn)行描述，為節(jié)省存儲空間，運(yùn)算中僅保存非零值和對應(yīng)的空間坐標(biāo)（算法執(zhí)行中以記錄網(wǎng)格編號代替）.

圖2為射線投影對應(yīng)的某空間二維坐標(biāo)系，其重建區(qū)域面積為N×N，網(wǎng)格編號規(guī)則為從左下角開始，分別命名為0，1，2…N2-1，共有N2個網(wǎng)格，其中每個網(wǎng)格長度為L，文中為簡化分析，設(shè)L=1（即對射線源坐標(biāo)和探測器坐標(biāo)進(jìn)行歸一化處理）.

為了求解射線在重建區(qū)域不同位置的投影系數(shù)，這里設(shè)區(qū)域中點(diǎn)O為原點(diǎn)坐標(biāo)，點(diǎn)M和N分別為入射點(diǎn)和出射點(diǎn)位置，Ni中i為出射點(diǎn)個數(shù)，ki為不同射線對應(yīng)的斜率，PXc，PYc分別表示射線與X軸、Y軸的交點(diǎn)，下標(biāo)c代表交點(diǎn)的個數(shù).下面以圖2給出的投影示意圖為例，分析射線穿過重建區(qū)域的網(wǎng)格編號以及對應(yīng)的交線長度.

從圖2中可以看出，射線與網(wǎng)格的交點(diǎn)主要有兩類，一類是與Y軸的交點(diǎn)，一類是與X軸的交點(diǎn).圖中射線MN1與網(wǎng)格的交點(diǎn)不僅包含與Y軸的交點(diǎn)，也包含與X軸的交點(diǎn)，但MN2為其特例，僅包含與Y軸的交點(diǎn).為了具有普遍意義，本節(jié)主要以MN1為例，求解其與X軸以及Y軸的點(diǎn)坐標(biāo).

ART算法是將連續(xù)圖像離散化為J=N×N的像素單元，待求解圖像f（x，y），所有穿過圖像的射線條數(shù)總計為M，則總的射線投影方程為

圖2 投影系數(shù)分析原理圖 Fig.2 Schematic diagram of projection coefficients

設(shè)MN1的斜率0＜k1≤1，與Y軸交點(diǎn)集合用數(shù)組PY（PYcx，PYcy）表示，PYcx表示與X軸相交的第c點(diǎn)的橫坐標(biāo)值，代表交線長度；PYcy表示與Y軸相交的第c點(diǎn)的縱坐標(biāo)值，代表網(wǎng)格編號；與X軸交點(diǎn)集合用數(shù)組PX（PXcx，PXcy）表示，PXcx表示與X軸相交的第c點(diǎn)的橫坐標(biāo)值，代表交線長度；PXcy表示與X軸相交的第c點(diǎn)的縱坐標(biāo)值，代表網(wǎng)格編號.

依據(jù)MN1在圖2中的位置關(guān)系，式（4）給出了數(shù)組PY中坐標(biāo)點(diǎn)的表達(dá)式

式中：c∈［1，2n］，表示向下取整，（PY1x，PY1y）表示入射點(diǎn)M的起始坐標(biāo).

式（5）給出了數(shù)組PX中交線長度的坐標(biāo)點(diǎn)表達(dá)式

式中：（PX1x，PX1y）表示入射點(diǎn)M的起始坐標(biāo)，這里（PX1x，PX1y）=（PY1x，PY1y）.

式（5）中沒有給出數(shù)組PX中點(diǎn)的網(wǎng)格編號，依據(jù)圖2分析發(fā)現(xiàn)：A與B兩點(diǎn)的網(wǎng)格編號相差N.由于有許多類似B點(diǎn)的點(diǎn)出現(xiàn)，會導(dǎo)致式（4）中數(shù)組的交點(diǎn)數(shù)增加，這里定義新集合數(shù)組為PY′（PY′cx，PY′cy），式（4）的坐標(biāo)表達(dá)式更新為

數(shù)組PY′中點(diǎn)的集合包含射線MN1穿過重建區(qū)域所有點(diǎn)的交線長度及對應(yīng)的網(wǎng)格編號，由此即可得到射線在投影區(qū)域的所有投影系數(shù).

式（4）～式（6）給出了圖2中射線斜率0＜k1≤1時的投影系數(shù)計算過程，對于k1的其它幾種情況，下面給出了分析說明.

1）當(dāng)-1≤k1＜0時，從圖2可以看出-1≤k1＜0時的射線與0＜k1≤1時的射線關(guān)于X軸對稱，因此可將式（4）～式（6）中的縱坐標(biāo)進(jìn)行對稱變換即可得到-1≤k＜0時的投影系數(shù)，其結(jié)果按小到大的順序可排列為

2）當(dāng)k1＞1時，從圖2可以看出射線與X軸的交點(diǎn)較多，因此可按X軸方向進(jìn)行計算，結(jié)果即是將式（6）中的橫坐標(biāo)改為縱坐標(biāo)；

3）當(dāng)k＜-1時，即是將k1＞1中的網(wǎng)格編號做關(guān)于X軸的對稱運(yùn)算即可；

4）當(dāng)k=0或k→∞時，射線與Y軸或X軸平行，只要確定了M點(diǎn)的坐標(biāo)，之后所有點(diǎn)的坐標(biāo)經(jīng)過簡單疊加即可.

由于ART算法主要是被射線驅(qū)動，一次迭代僅記錄一條射線的信息，實際處理中可分次處理不同射線的信息，這有利于節(jié)省內(nèi)存空間，提高算法重建速度.

2.2 三維情況下的投影系數(shù)快速計算

上面分析了二維情況下的投影系數(shù)計算方法，實際中被投影物體往往是三維狀態(tài)，為此需將圖2中的二維坐標(biāo)系擴(kuò)展為三維坐標(biāo)系，其中O為物體中心，也是坐標(biāo)原點(diǎn)，三維空間可將物體劃分為N層，每層包含N×N個立方體（1個立方體即可視為物體的一個體素單元），N層共N×N×N個立方體，每個正方體的具體編號可根據(jù)坐標(biāo)軸正方向進(jìn)行排列.

三維投影系數(shù)求解的基本思想為：首先將射線在三維坐標(biāo)下的投影映射到XOY和Z兩個二維平面，其次分別求出射線在XOY平面、Z平面的投影系數(shù)，最后將兩者獲得的結(jié)果進(jìn)行綜合分析，即可得到射線的三維投影系數(shù).其中XOY平面的投影系數(shù)為層內(nèi)索引，具體求解可采用2.1節(jié)方法；Z平面的投影系數(shù)為層間索引，層間索引（Z方向索引）主要是計算射線穿過被測物體的層信息.實際處理中，層信息通過計算射線在XOZ平面或YOZ平面的投影系數(shù)獲得，當(dāng)射線在X軸投影較長，選擇XOZ平面；反之，當(dāng)射線在Y軸投影較長，選擇YOZ平面.

如圖3所示，將射線MN的三維投影系數(shù)求解轉(zhuǎn)化為：①求解M′N′與XOY平面相交的投影系數(shù)；②求解M″N″與YOZ平面相交的投影系數(shù)；③將①和②的結(jié)果綜合處理.下面對這三個過程分別進(jìn)行說明.

圖3 三維射線投影示意圖 Fig.3 3D projection schematic diagram

1）M′N′與XOY投影系數(shù)的計算

分析圖3可以得出，射線M′N′投影系數(shù)的計算與圖2中射線MN1的計算完全相同，因此可采用式（4）～式（7）完成.

2）M″N″與YOZ投影系數(shù)的計算

分析圖3可以得出，射線M″N″與水平網(wǎng)格線l（圖3中的點(diǎn)劃線）相交的信息即可反映層索引的變化，所以這里僅分析M″N″與l相交的情況.

設(shè)M″N″的斜率0＜k′≤1，與Z軸交點(diǎn)集合用數(shù)組PZ（PZcy，PZcz）表示，PZcy表示與Y軸平行線相交的第c點(diǎn)的橫坐標(biāo)值，代表交線長度；PZcz表示與Y軸平行線相交的第c點(diǎn)的縱坐標(biāo)值，代表網(wǎng)格編號.按照式（4）的分析方法可得數(shù)組PZ中坐標(biāo)點(diǎn)的表達(dá)式為

式中：c∈［1，2n］，（PZ1y，PX1z）表示入射點(diǎn)M的起始坐標(biāo).最后，數(shù)組PZ內(nèi)的元素即為層間索引結(jié)果，包含層信息以及對應(yīng)的長度信息.

3）綜合1）和2）結(jié)果的分析

射線MN的完整投影系數(shù)是將1）結(jié)果和2）結(jié)果進(jìn)行綜合，其中1）得出了射線M′N′在XOY平面內(nèi)的網(wǎng)格編號和交線長度；2）得出了射線M″N″在YOZ內(nèi)的網(wǎng)格編號和交線長度；由于兩者都含有相同的橫坐標(biāo)軸Y，因此可利用其將XOY平面內(nèi)網(wǎng)格編號和YOZ平面內(nèi)編號相融合，具體融合方法為：

設(shè)a，b，c，e，f，h，i為射線M′N′在Y軸上的映射點(diǎn)，其網(wǎng)格編號分別對應(yīng)為p1，p2，p3，p4，p5，p6，p7，點(diǎn)d，g，j為射線M″N″在Y軸上的映射點(diǎn)，其分別對應(yīng)的層間網(wǎng)格編號為q1，q2，q3.根據(jù)圖3中各點(diǎn)在Y軸投影關(guān)系可知相鄰點(diǎn)間的線段網(wǎng)格編號規(guī)則為：

1）當(dāng)相鄰點(diǎn)在Y軸坐標(biāo)都小于d坐標(biāo)，則相鄰點(diǎn)間線段網(wǎng)格編號為其原有網(wǎng)格編號和d的網(wǎng)格編號相加之和.如：［a，b］間線段網(wǎng)格編號為p1+q1，［b，c］間線段網(wǎng)格編為p2+q1.

2）當(dāng)相鄰點(diǎn)間在Y軸坐標(biāo)一個小于d坐標(biāo)，另一個大于d坐標(biāo)，則表明這兩點(diǎn)不在同一層，需要分別進(jìn)行網(wǎng)格編號.如［c，e］間線段網(wǎng)格編號分為p3+q1和p3+q2.

依據(jù)上述規(guī)則，可得［e，f］，［f，g］，［g，h］，［h，i］間線段的網(wǎng)格編號分別為：p4+q2，p5+q2，p5+q3，p6+q3.另外，對于這些點(diǎn)的交線長度可由距離公式獲得.

考慮到系統(tǒng)設(shè)計時，即將射線源、被測物體中心以及探測器中心置于同一平面內(nèi)的同一直線上，所以探測器每一列對應(yīng)的射線在XOY平面的層內(nèi)索引結(jié)果是一樣的，只是層間索引結(jié)果不同，因此，在計算一個投影角度下射線的投影系數(shù)時，可以采取列優(yōu)先的方法，一個列上的投影射線只需計算一次層內(nèi)索引.而且，利用存在的對稱關(guān)系［13］，在列優(yōu)先的基礎(chǔ)上，可以只計算Z軸正半軸的射線的投影系數(shù)，負(fù)半軸方向的射線的投影系數(shù)則直接根據(jù)對稱關(guān)系由正半軸得到.由此，又可減少一半的運(yùn)算量.

3 實驗結(jié)果與分析

以Sheep＿Logan模型為例進(jìn)行算法仿真實驗，重建圖像大小為128×128×128，體素尺寸為0.256 mm，探測器矩陣大小為128×128，像元尺寸0.512 mm.射線源距旋轉(zhuǎn)中心的距離以及探測器中心距旋轉(zhuǎn)中心的距離均設(shè)定為780 mm.每繞旋轉(zhuǎn)中心旋轉(zhuǎn)1°采集一幅投影，一共得到360幅投影.采用ART算法進(jìn)行重建，設(shè)定圖像初值x（0）=0，松弛因子λ=0.025，順序訪問方式，采用本文提出的投影系數(shù)快速計算方法.迭代3次的重建結(jié)果切片與原始模型切片做對比，如圖4所示.圖5為對應(yīng)第64行重建切片與原始模型切片灰度曲線的對比圖.

圖4 原始模型與重建結(jié)果圖 Fig.4 Original model and reconstructed result

圖5 切片的灰度對比圖 Fig.5 Grayscale contrast of original model and reconstructed result

由圖4，圖5對比可得，經(jīng)過3次迭代后，重建圖像切片與原始模型切片的灰度曲線相近，說明在ART算法中采用本文提出的投影系數(shù)計算方法，重建圖像質(zhì)量能夠得到保障.

下面從投影重建時間的角度列表分析，將本文方法與傳統(tǒng)的Siddon算法以及張順利所提的投影系數(shù)計算方法進(jìn)行比較.重建時間對比如表1所示.

表1 不同投影系數(shù)據(jù)算方法重建時間對比 Tab.1 Reconstruction time contrast of different projection coefficient computation method

從表1可以看出，本文所提方法較傳統(tǒng)Siddon算法重建速度提高約13倍；相比張順利所提的方法，重建速度提高1.26倍.并且由于所提算法中分支判斷大為減少，所以后期在進(jìn)行CUDA算法加速時，有利于并行加速的實現(xiàn)［14］.利用存在的對稱關(guān)系，進(jìn)行對稱優(yōu)化后，速度可再提高約3倍，加速效果更加明顯.

4 結(jié) 論

迭代算法是基礎(chǔ)的CT重建算法，本文以ART算法為例，介紹了ART算法的基本原理，其中，投影系數(shù)的計算是重建速度和重建質(zhì)量的最大影響因素.本文分別針對二維情況和三維情況下，以Siddon算法和張順利所提算法為基礎(chǔ)，提出投影系數(shù)的快速計算方法，從射線穿過網(wǎng)格的相交規(guī)律出發(fā)，通過順序增量計算，快速推導(dǎo)出穿過的網(wǎng)格編號并計算其交線長度，極大地減少了運(yùn)算量，分支判斷也大為減少.仿真實驗證明本文所提方法在保證重建圖像重建質(zhì)量的基礎(chǔ)上，算法速度大幅提高，相比其他常用算法有明顯的速度優(yōu)勢.

隨著迭代算法重建速度的提高，鑒于其在不完備投影數(shù)據(jù)情況下的優(yōu)異表現(xiàn)，必定會在實時的工業(yè)成像檢測系統(tǒng)中發(fā)揮重要的作用.下一步的工作應(yīng)重點(diǎn)將放在對算法的硬件加速改進(jìn).

［1］Li J，Li LT，Cong P，et al.Rotating polar-coordinate ART applied in industrial CT image reconstruction［J］.NDTffE International，2007，40（4）：333-336.

［2］張朝宗，郭志平，張鵬，等.工業(yè)CT技術(shù)與原理［M］.北京：科學(xué)出版社，2009.

［3］莊天戈.CT原理與方法［M］.上海：上海交通大學(xué)出版社，1992.

［4］冷駿.ART算法中關(guān)于松弛因子的研究［J］.光學(xué)儀器，2014，36（1）：46-51.Leng Jun.Research on the relaxation factor in algebraic reconstruction technique［J］.Optical Instruments，2014，36（1）：46-51.（in Chinese）

［5］Yang Juan，Hou Huiling，Shi Lang.A method based on vector type to sparsely store and quickly access projection matrix［J］.Journal of Measurement Science and Instrumentation，2015，6（1）：53-56.

［6］Xue Z，Zhang L L.A new algorithm for calculating the radiological path in CT image reconstruction［C］.International conference on Electronic ff Mechanical Engineering and Information Technology，2011：4527-4530.

［7］Miao C.Comparative studies of different system models for iterative CT image reconstruction［D］.Master Dissertation.Winston-Salem：Wake Forest University，2013.

［8］Siddon R L.Fast calculation of the exact radiological path for a three-dimensional CT array［J］.Med.Phys.，1985，12：252-255.

［9］張順利，張定華，趙歆波.代數(shù)重建法中的一種快速投影系數(shù)計算方法［J］.計算機(jī)應(yīng)用研究，2007，24（5）：38-40.Zhang Shunli，Zhang Dinghua，Zhao Xinbo.Approach for fast projection coefficient computation in algeb reconstruction technique［J］.Application Research of Computers，2007，24（5）：38-40.（in Chinese）

［10］張順利，張定華，黃魁東，等.錐束ART算法快速圖像重建［J］.儀器儀表學(xué)報，2009，30（4）：887-892.Zhang Shunli，Zhang Dinghua，Huang Kuidong，et al.Fast image reconstruction with cone-beam ART algorithm［J］.Chinese Journal of Scientific Instrument，2009，30（4）：887-892.（in Chinese）

［11］楊文良，魏東波.一種改進(jìn)投影系數(shù)計算的快速ART算法［J］.CT理論與應(yīng)用研究，2012，21（2）：187-195.Yang Wenliang，Wei Dongbo.A fast algebraic reconstruction algorithm based on improved projection coefficient computation［J］.CT Theory and Applications，2012，21（2）：187-195.（in Chinese）

［12］Zhang S L，Zhang D H，Gong H，et al.Fast and accurate computation of system matrix for area integral model-based algebraic reconstruction technique［J］.Optical Engineering.2014，53（11）：1-9.

［13］肖波.基于離散化模型對稱結(jié)構(gòu)改進(jìn)的EM算法研究［D］.北京：北京交通大學(xué)，2009.

［14］雷德川，許州，陳浩.基于CUDA的GPU加速代數(shù)迭代重建算法［J］.無損檢測，2012，34（8）：5-9.Lei Dechuan，Xu Zhou，Chen Hao.Accelerating simultaneous algebraic reconstruction technique based on CUDA-enabled GPU［J］.NDT，2012，34（8）：5-9.（in Chinese）