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

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

0 引 言

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

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

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

1 ART算法基本原理

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

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

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

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

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

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

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

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

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

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

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

從圖2中可以看出，射線(xiàn)與網(wǎng)格的交點(diǎn)主要有兩類(lèi)，一類(lèi)是與Y軸的交點(diǎn)，一類(lèi)是與X軸的交點(diǎn).圖中射線(xià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），所有穿過(guò)圖像的射線(xiàn)條數(shù)總計(jì)為M，則總的射線(xiàn)投影方程為

圖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)值，代表交線(xiàn)長(zhǎng)度；PYcy表示與Y軸相交的第c點(diǎn)的縱坐標(biāo)值，代表網(wǎng)格編號(hào)；與X軸交點(diǎn)集合用數(shù)組PX（PXcx，PXcy）表示，PXcx表示與X軸相交的第c點(diǎn)的橫坐標(biāo)值，代表交線(xiàn)長(zhǎng)度；PXcy表示與X軸相交的第c點(diǎn)的縱坐標(biāo)值，代表網(wǎng)格編號(hào).

依據(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中交線(xiàn)長(zhǎng)度的坐標(biāo)點(diǎn)表達(dá)式

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

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

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

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

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

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

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

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

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

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

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

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

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

圖3 三維射線(xiàn)投影示意圖 Fig.3 3D projection schematic diagram

1）M′N(xiāo)′與XOY投影系數(shù)的計(jì)算

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

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

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

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

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

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

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

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

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

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

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

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

3 實(shí)驗(yàn)結(jié)果與分析

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

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

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

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

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

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

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

4 結(jié) 論

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

隨著迭代算法重建速度的提高，鑒于其在不完備投影數(shù)據(jù)情況下的優(yōu)異表現(xiàn)，必定會(huì)在實(shí)時(shí)的工業(yè)成像檢測(cè)系統(tǒng)中發(fā)揮重要的作用.下一步的工作應(yīng)重點(diǎn)將放在對(duì)算法的硬件加速改進(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］.計(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é)報(bào)，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ù)計(jì)算的快速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］肖波.基于離散化模型對(duì)稱(chēng)結(jié)構(gòu)改進(jìn)的EM算法研究［D］.北京：北京交通大學(xué)，2009.

［14］雷德川，許州，陳浩.基于CUDA的GPU加速代數(shù)迭代重建算法［J］.無(wú)損檢測(cè)，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）