謝睿誠，劉楚欣，潘玉媚，李 健

（汕頭大學(xué)理學(xué)院數(shù)學(xué)系，廣東 汕頭 515063）

0 引言

CT，即計(jì)算機(jī)斷層成像技術(shù)，是一種依據(jù)外部投影數(shù)據(jù)重建物體內(nèi)部結(jié)構(gòu)圖像的無損檢測技術(shù).CT系統(tǒng)示意圖見圖1.

王召巴[1]在2001年從濾波反投影算法的基本原理著手，分析了CT系統(tǒng)旋轉(zhuǎn)中心偏移對圖像重構(gòu)所造成的影響；Olander等[2]利用Radon變換性質(zhì)把角度相差180°的投影經(jīng)過平移與翻轉(zhuǎn)，計(jì)算出旋轉(zhuǎn)中心的偏移量；李增云等[3]再次證明了物體質(zhì)心與其投影質(zhì)心關(guān)系定理，并提出了利用部分投影快速校正CT系統(tǒng)旋轉(zhuǎn)中心的方法.這些方法理論上可以準(zhǔn)確地確定CT系統(tǒng)旋轉(zhuǎn)中心的位置再進(jìn)一步進(jìn)行圖像的重構(gòu)，但現(xiàn)實(shí)中往往因?yàn)槠溆?jì)算復(fù)雜，實(shí)操性不強(qiáng)等原因而難以在醫(yī)學(xué)和工業(yè)方面被廣泛應(yīng)用.

由于CT系統(tǒng)對系統(tǒng)旋轉(zhuǎn)中心在托盤中的位置、探測器單元之間的距離以及X射線的180個方向等參數(shù)有極高的要求，這極大地增加了CT系統(tǒng)制造和安裝的難度.本文運(yùn)用濾波反投影和Radon變換等，針對CT系統(tǒng)參數(shù)標(biāo)定及成像展開研究.

1）建立了直接求解模型，求解CT系統(tǒng)旋轉(zhuǎn)中心在正方形托盤中的位置、探測器單元之間的距離以及該CT系統(tǒng)使用的X射線的180個方向，然后利用重構(gòu)法驗(yàn)證其結(jié)果.

2）建立濾波反投影模型計(jì)算未知介質(zhì)的衰減系數(shù)，然后通過Radon變換和Beer-Lambert定律導(dǎo)出吸收率的計(jì)算公式，并由此計(jì)算未知介質(zhì)吸收率.并根據(jù)此前確定的標(biāo)定參數(shù)對吸收率數(shù)據(jù)進(jìn)行校正，得出指定位置的吸收率，并重構(gòu)出未知介質(zhì)的圖像.

3）定義模板標(biāo)定參數(shù)與實(shí)際參數(shù)的誤差為精度度量指標(biāo)，定義重構(gòu)前后吸收率的變化為穩(wěn)定性度量指標(biāo).然后以矩形和正方形代替橢圓和圓建立新的標(biāo)定模板，減少濾波反投影中插值方法產(chǎn)生的誤差.最后，通過對比兩個模板重構(gòu)前后吸收率差距和重構(gòu)前后吸收率變化量，認(rèn)為方形模板的精度和穩(wěn)定性比橢圓形模板更好.

本文的原文獲得2017年全國大學(xué)生數(shù)學(xué)建模競賽一等獎，所采用的數(shù)據(jù)亦來自2017年全國大學(xué)生數(shù)學(xué)建模競賽（http://www.mcm.edu.cn/html_cn/node/460baf68ab0ed0e1e 557a0c79b1c4648.html）.

1 問題提出

1.1 問題提出

問題1：在正方形托盤上放置兩個均勻固體介質(zhì)組成的標(biāo)定模板，模板的幾何信息如圖2所示，其中每一點(diǎn)的數(shù)值反映了該點(diǎn)的吸收強(qiáng)度，這里稱為“吸收率”.根據(jù)這一模板及其接收信息，確定CT系統(tǒng)旋轉(zhuǎn)中心在正方形托盤中的位置，探測器單元之間的距離以及該CT系統(tǒng)使用的X射線的180個方向.

問題2：附件3是利用上述CT系統(tǒng)得到的某未知介質(zhì)的接收信息.利用問題1中得到的標(biāo)定參數(shù)，確定該未知介質(zhì)在正方形托盤中的位置，幾何形狀和吸收率等信息.另外，請具體給出圖3所給的10個位置處的吸收率.

問題3：附件5是利用上述CT系統(tǒng)得到的另一個未知介質(zhì)的接收信息.利用問題1中得到的標(biāo)定參數(shù)，給出該未知介質(zhì)的相關(guān)信息.另外，請具體給出圖3所給的10個位置處的吸收率.

問題4：分析問題1中參數(shù)標(biāo)定的精度和穩(wěn)定性.在此基礎(chǔ)上自行設(shè)計(jì)新模板，建立對應(yīng)的標(biāo)定模型，以改進(jìn)標(biāo)定精度和穩(wěn)定性.

圖1 CT系統(tǒng)示意圖

圖2 模板示意圖（mm）

圖3 10個位置示意圖

2 CT圖像重構(gòu)模型推導(dǎo)

2.1 Radon變換和吸收率

根據(jù)文獻(xiàn)[4]，附件3中的接受信息數(shù)據(jù)即Radon變換的投影函數(shù)值，問題2、3是要求根據(jù)投影函數(shù)值來求吸收率，我們以衰減系數(shù)為橋梁構(gòu)建投影函數(shù)值與吸收率之間的關(guān)系.Radon變換示意圖見圖4.

圖4 Radon變換示意圖

2.1.1 Radon變換（投影函數(shù)值與衰減系數(shù)的關(guān)系）

設(shè)p軸與X射線垂直，p＝xcosφ＋ysinφ.設(shè)未知介質(zhì)的函數(shù)為f（x，y），未知介質(zhì)經(jīng)X射線投影在p軸上得到，即：

上述是Radon變換的原理，但在這個問題中，我們已知的數(shù)據(jù)是附件3中的接收信息數(shù)據(jù)即Radon變換的投影函數(shù)值，要通過來求 f（x，y）的函數(shù)值，得出物質(zhì)的吸收率來畫出物質(zhì)的圖像，因此本問題可以轉(zhuǎn)換為求Radon變換的逆過程，根據(jù)文獻(xiàn)[4]，這個過程就是圖像重構(gòu)的過程，我們將采取濾波反投影的方法來得到f（x，y）.

2.1.2 Beer-Lambert定律（衰減系數(shù)轉(zhuǎn)換為吸收率）

由Beer-Lambert定律[8]可知：

其中，A為物體的吸光率，K為吸收率（對于某種均勻的物質(zhì)，K為常數(shù)），l為介質(zhì)厚度，c為吸光物質(zhì)濃度（對均勻物體，c為常數(shù)）.

設(shè)μ為衰減系數(shù)，令

即衰減系數(shù)和吸收率之間的關(guān)系是線性關(guān)系.

據(jù)公式（1）和公式（4）有：

因此可以根據(jù)物體在不同角度下X射線的接收信息利用濾波反投影原理來重構(gòu)物體.

2.2 濾波反投影原理

2.2.1 投影切片定律：

F（ωcosφ，ωsinφ）見 2.1和 2.2.2中說明.

2.2.2 二維傅里葉變換及極坐標(biāo)變換：

定義

為f（x，y）的二維傅里葉變換.

令 u＝ωcosφ，v＝ωsinφ，則有函數(shù)

式中G（φω）為通過φ角的F切片，對固定的φ，G（φω）只有一個自變量.

那么上式可以表示為：

2.2.3 濾波器（文獻(xiàn)[5]）：

根據(jù)傅里葉變換的卷積定理：

2.2.4 傅里葉逆變換

但我們運(yùn)用上式的離散化形式：

2.3 濾波反投影原理在本題中的應(yīng)用

濾波反重構(gòu)原理要應(yīng)用于本題重構(gòu)過程中，還需要解決下面兩個問題：插值方法的選擇和濾波器的選擇.

2.3.1 插值方法的選擇

圖5 濾波反投影步驟圖

a）線性插值

線性插值中插值點(diǎn)xi，i＝1，2，…，511對應(yīng)的投影函數(shù)的計(jì)算公式如下所示：

b）三次樣條插值

三次樣條插值法中插值點(diǎn) xi，i＝1，2，…，511對應(yīng)的投影函數(shù)的計(jì)算公式如下所示：

其中ai，bi，ci，di為參數(shù)，由確定.

2.3.2 濾波器的選擇

我們選擇的濾波器是Cosine濾波器，Cosine濾波器是在Ram-Lak濾波器的基礎(chǔ)上乘上一個Cosine函數(shù).

Ram-Lak濾波器：

3 問題1模型的建立與求解

3.1 確定探測器單元之間的距離d

探測器單元之間的距離d，可理解為探測器單元數(shù)密度的倒數(shù)，即該問題可通過求解探測器單元數(shù)密度來解決.

a）根據(jù)附件提供的第151次旋轉(zhuǎn)數(shù)據(jù)，光源平行于y軸時，接收到信息的探測器單元個數(shù)n為108個.結(jié)合模板的幾何信息可知，108個探測器單元個數(shù)對應(yīng)的實(shí)際距離l為30 mm，即探測器單元數(shù)密度ρ為：

b）根據(jù)附件中第61次旋轉(zhuǎn)的數(shù)據(jù)，光源平行于x軸時，接收到信息的探測器單位個數(shù)為288個.圖2模板示意圖提供的長軸長度為80 mm，類似地，根據(jù)a）可求出探測器單位數(shù)密度同樣為3.6個/mm，探測器單元之間的距離d為0.277 8 mm.

經(jīng)以上分析可得探測器單元之間的距離為0.277 8 mm.

3.2 確定CT系統(tǒng)使用的X射線的180個方向

利用重構(gòu)原理（詳見2.1 Radon變換和吸收率），將附件2數(shù)據(jù)進(jìn)行重構(gòu)圖像，得到的重構(gòu)圖像如圖6所示：

觀察圖6發(fā)現(xiàn)，附件2數(shù)據(jù)重構(gòu)圖像發(fā)生了明顯的角度偏轉(zhuǎn).利用橢圓與圓圓心連線與水平方向的夾角即重構(gòu)所得圖像的偏轉(zhuǎn)角度為29.649 3°.因此，掃描光源初始方向?yàn)榕cx軸負(fù)半軸夾角為60.350 7°，從第四象限指向第二象限的方向，掃描光源最終方向?yàn)槌跏挤较虻哪娣较颍看蜟T系統(tǒng)逆時針旋轉(zhuǎn)1°.

圖6 附件2數(shù)據(jù)重構(gòu)得到的圖像

3.3 求解CT系統(tǒng)旋轉(zhuǎn)中心

利用3.2中得到的初始旋轉(zhuǎn)角度，對重構(gòu)圖像進(jìn)行修正結(jié)果如圖7所示.

由于橢圓模板在標(biāo)定時放置在托盤的正中央，其重構(gòu)圖像也應(yīng)該在正中央，可是圖像中橢圓形的中心與圖像中心并不重合，因而可以利用橢圓中心與圖像中心的偏差來確定CT系統(tǒng)的旋轉(zhuǎn)中心.

以重構(gòu)圖像的橢圓圓心為原點(diǎn)，建立圖8所示坐標(biāo)系，圖像中心位于直角坐標(biāo)坐標(biāo)系的第二象限，求得圖像中心的位置為（-8.998 5，5.999 1）.

圖7 調(diào)整角度后重構(gòu)得到的圖像

圖8 置于坐標(biāo)系的重構(gòu)圖像

4 問題2，3模型的建立與求解

4.1 吸收率校正模型

利用濾波反投影原理對標(biāo)定模板的180個角度的接收信息進(jìn)行重構(gòu)，并利用3.2和3.3中所得角度和旋轉(zhuǎn)中心位置，對重構(gòu)圖像的方向和位置進(jìn)行校正.可得圖像標(biāo)定模板中吸收率為1.000 0處的衰減系數(shù)的觀測值中位數(shù)為0.490 4，標(biāo)定模板中吸收率為0.000 0處的吸收率中位數(shù)為0.000 1.由Beer-Lambert定律，可以得吸收率與衰減系數(shù)關(guān)系為：

4.2 問題2圖像重構(gòu)

依據(jù)2中建立的模型，利用濾波反投影原理對附件3所給的接收信息進(jìn)行重構(gòu)，并利用公式（21）進(jìn)行校正，得到重構(gòu)圖像如圖9所示，問題中指定的10個位置處的吸收率結(jié)果如表1所示.

圖9 原始重構(gòu)圖像（標(biāo)記指定位置）

表1 圖9所給10個位置處的吸收率結(jié)果

4.3 問題3圖像重構(gòu)

利用2中建立的模型，根據(jù)濾波反投影模型將CT系統(tǒng)得到的未知介質(zhì)接收信息，即附件5數(shù)據(jù)進(jìn)行圖像重構(gòu)，并利用公式（21）進(jìn)行優(yōu)化，得到的重構(gòu)圖像如圖10所示，問題指定的10個位置處的吸收率結(jié)果如表2所示.

圖10 原始重構(gòu)圖像（標(biāo)記指定位置）

表2 圖3所給10個位置處的吸收率結(jié)果

5 問題4模型的建立與求解

5.1 模型的建立

5.1.1 精度評價(jià)

1）模板位置對校正精度的影響

通過平移，旋轉(zhuǎn)變換，改變模板的位置并對模板做Radon變換，得到模板的接收信息，并利用3中確定旋轉(zhuǎn)中心和初始投影角度的方法，求解旋轉(zhuǎn)中心和初始投影角度.與實(shí)際位置進(jìn)行比較，重復(fù)多次并求平均值，得到平均誤差作為精度估計(jì)依據(jù).

2）Radon變換再重構(gòu)的吸收率差距

對模板作Radon變換得出模板的投影數(shù)據(jù)，再用投影數(shù)據(jù)運(yùn)用2中建立的重構(gòu)方法還原圖像，在模板和重構(gòu)得到的圖像中建立相同的坐標(biāo)系，比較相同坐標(biāo)下兩者的吸收率數(shù)值差距，設(shè)前者為fo（x，y），后者為fl（x，y）.

用L2范數(shù)來衡量fo（x，y）和fl（x，y）之間的差距R，設(shè)已知吸收率的點(diǎn)構(gòu)成m×m的矩陣：

用均方誤差來衡量吸收率的差距大小，均方誤差計(jì)算公式為：

5.1.2 穩(wěn)定性評價(jià)

我們用各參數(shù)對精度的影響來衡量模板的穩(wěn)定性，在兩種模板的精度足夠好的情況下，通過改變反投影方向數(shù)、旋轉(zhuǎn)中心平移距離、初始旋轉(zhuǎn)角度大小以及計(jì)算精度指標(biāo)中的重構(gòu)前后吸收率差距的變化量來說明其穩(wěn)定性，即吸收率差距變化不大就認(rèn)為模板的參數(shù)標(biāo)定精度穩(wěn)定，反之則不穩(wěn)定.

5.2 新標(biāo)定模板設(shè)計(jì)

5.2.1 原始標(biāo)定模板誤差成因分析

a）角度變化引起的探測器數(shù)目變化量

附件2中被判斷為X射線掃描光源方向平行于x軸的數(shù)據(jù)有7組，即第58次掃描至第64次掃描得到的接收信息的探測器單元個數(shù)相同.考慮造成這樣現(xiàn)象的原因是當(dāng)X射線掃描光源方向旋轉(zhuǎn)至與x軸夾角較小的區(qū)域內(nèi)時，每逆時針旋轉(zhuǎn)一次引起標(biāo)定模板投影的接收數(shù)據(jù)變化量極小，且CT系統(tǒng)測量的精度不高，不足以反映出如此微小的變化量.

b）插值方法的適應(yīng)性

由于重構(gòu)方法中應(yīng)用的是線性插值方法，橢圓的邊緣為圓弧狀，線性插值并不能完全貼合橢圓邊緣，在線性插值中會引起誤差.即使用其它插值方法如多項(xiàng)式插值方法也不可以完全貼合，因此原始模板會在重構(gòu)圖形時產(chǎn)生誤差.

5.2.2 新模板設(shè)計(jì)

針對原始模板的不足之處，我們構(gòu)建新的模板（方形模板）如圖11所示：

圖11 新模板示意圖

圖12 新模板吸收率

如圖11所示，新模板（方形模板）由矩形和正方形組成，矩形中心在正方形托盤的中心，矩形和正方形均關(guān)于正方形托盤的水平對稱軸對稱，其幾何信息如圖11所示.

5.3 模板測試

我們對橢圓模板（原始模板）以及方形模板（新模板）做Radon變換以及重構(gòu)，并進(jìn)行精度分析和穩(wěn)定性分析.

5.3.1 精確度評價(jià)對比

對橢圓模板和方形模板做20次不同的平移變換和30次不同角度的旋轉(zhuǎn)變化，并計(jì)算坐標(biāo)和旋轉(zhuǎn)角度與實(shí)際值的平均誤差（見表3）.

表3 橢圓模板和方形模板的校正精度對比

從表3可以看出，橢圓模板和方形模板的標(biāo)定參數(shù)的誤差都比較小，精度較高，因而兩個模板都適合參數(shù)標(biāo)定.

5.3.2 穩(wěn)定性評價(jià)對比

由于CT系統(tǒng)是以等角度間隔掃描一次方式離散進(jìn)行的，我們定義CT系統(tǒng)掃描1周所做投影的數(shù)量，也即重構(gòu)圖像時可以利用的輪廓數(shù)量為反投影方向數(shù).為探究在反投影方向數(shù)減少，投影方向的角度間隔增大時，兩種模板的穩(wěn)定性.以反投影方向數(shù)為橫坐標(biāo)，以模板原吸收率與重構(gòu)后模板吸收率的歐幾里得距離和均方誤差為縱坐標(biāo)，建立平面直角坐標(biāo)以呈現(xiàn)橢圓模板和方形模板在吸收率差距這一精確度指標(biāo)上的差別.

由圖13和圖14可以看出，在反投影方向數(shù)范圍為（0，180]時，方形模板與橢圓模板相比，重構(gòu)前后吸收率差距都比較小，并且在減少反投影方向數(shù)時，方形模板吸收率差距的變化量比橢圓模板的要小，這說明方形模板的穩(wěn)定性比橢圓模板的要高.

進(jìn)一步探究模板擺放時旋轉(zhuǎn)中心位置及旋轉(zhuǎn)角度對精度的影響，分別以旋轉(zhuǎn)中心平移距離和旋轉(zhuǎn)角度為橫坐標(biāo)，以模板原吸收率與重構(gòu)后模板吸收率的歐氏距離為縱坐標(biāo)，建立平面直角坐標(biāo)系得到圖15，圖16（平移距離）以及圖17圖18（旋轉(zhuǎn)坐標(biāo)），橢圓模板和方形模板在穩(wěn)定性指標(biāo)上的差別見圖15-圖18.

由圖15圖16可以看出，在調(diào)節(jié)旋轉(zhuǎn)中心的位置變化時，橢圓模板的誤差要稍小于方形模板，并且二者的誤差都很?。╩se＜0.001），很難區(qū)分方形模板與橢圓模板的精度和穩(wěn)定性差距.進(jìn)一步的，從旋轉(zhuǎn)角度變化的影響來探究（見圖17圖18）橢圓模板和方形模板的差距.

圖13 不同反投影方向數(shù)下兩種模板重構(gòu)前后吸收率差距對比（以歐幾里得距離衡量）

圖14 不同反投影方向數(shù)下兩種模板重構(gòu)前后吸收率差距對比（以均方誤差衡量）

圖15 不同平移距離下兩種模板重構(gòu)前后吸收率差距對比（以歐幾里得距離衡量）

圖16 不同平移距離下兩種模板重構(gòu)前后吸收率差距對比（以均方誤差衡量）

圖17 不同旋轉(zhuǎn)角度下兩種模板重構(gòu)前后吸收率差距對比（以歐幾里得距離衡量）

圖18 不同旋轉(zhuǎn)角度下兩種模板重構(gòu)前后吸收率差距對比（以均方誤差衡量）

由圖17和圖18可以看出，在改變旋轉(zhuǎn)角度（橫坐標(biāo)）時，在旋轉(zhuǎn)角度大于10°時，方形模板重構(gòu)前后的吸收率差距比橢圓模板的小，并且在旋轉(zhuǎn)角度小于10°時，二者誤差都很小，因而認(rèn)為模板的擺放角度對標(biāo)定精度影響不大.

綜上所述，改變反投影方向數(shù)時，方形模板的穩(wěn)定性優(yōu)于橢圓模板；平移，旋轉(zhuǎn)時方形模板和橢圓模板重構(gòu)的誤差都極小.故使用方形模板作為標(biāo)定模板要優(yōu)于橢圓模板.

6 模型評價(jià)

本文建立的濾波反投影模型實(shí)現(xiàn)簡單，精度高，運(yùn)算量小，能獲得較好的重建圖像，對實(shí)際生產(chǎn)中的CT系統(tǒng)參數(shù)的修正有一定的幫助，但重構(gòu)所得圖像邊緣出現(xiàn)偽影，雖然我們使用的中值濾波有效平滑去除圖像的噪聲點(diǎn)，但未能完全消除其影響，這也是后續(xù)深入研究的一個方向.