ct系統(tǒng)參數(shù)標定及成像模型

1、<p>　　CT系統(tǒng)參數(shù)標定及成像模型</p><p><b>　　摘要</b></p><p>　　本文在CT工作原理的基礎(chǔ)上通過建立參數(shù)標定及成像模型，實現(xiàn)了對安裝好的CT系統(tǒng)標定參數(shù)，以及對未知結(jié)構(gòu)的樣品進行成像。</p><p>　　首先，根據(jù)尺寸、位置已知的模板的吸收率數(shù)據(jù)確定其材料的均勻性。將已知的吸收率和吸收信息代入La

2、mbert-Beers吸收定律，可建立吸收信息與被測物體厚度的函數(shù)關(guān)系。另外，吸收信息中非零值越多，則單元數(shù)越多，投影越長，物體寬度越大。由此可用物體尺寸求得探測器單元間距。篩選數(shù)據(jù)得到模板的最長和最短投影，利用模板在探測器的投影位置找到二者的幾何關(guān)系，可解出旋轉(zhuǎn)中心的坐標。最后利用一定次數(shù)旋轉(zhuǎn)的角度和旋轉(zhuǎn)次數(shù)作商，可求得每次旋轉(zhuǎn)的 角度，結(jié)合初始角，可知X射線的180個方向。</p><p>　　然后，因為衰減

3、系數(shù)在不均勻介質(zhì)中會隨著位置而變化，沿照射路徑作線積分即為接收信息。該過程與投影過程類似，因此運用直接傅里葉反變換法反解出衰減系數(shù)。因圖像較大，可調(diào)用imresize函數(shù)縮放。附件3物體形狀較為規(guī)則，可當作橢圓處理。由正方形托盤和該物體的幾何關(guān)系可求出位置。將圖形拆分開后，分別采用多項式擬合、冪指數(shù)擬合，以此描述形狀。吸收率和衰減系數(shù)有比例關(guān)系，可由模板求得，帶入被測件的衰減系數(shù)即可知其吸收率。通過坐標變換，將10個點的位置解出，則可知

4、吸收率。</p><p><b>　　一、問題重述</b></p><p>　?。茫约夹g(shù)通過Ｘ射線穿過被檢測物體產(chǎn)生衰減的原理，能夠在不破壞物體內(nèi)部結(jié)構(gòu)的前提下，對被掃描物體進行可視化成像[1]。具體的工作過程為從發(fā)射器射出平行射線，具有512個單元的探測器在被測物體之后接收信息。發(fā)射器和探測器繞某固定的旋轉(zhuǎn)中心逆時針旋轉(zhuǎn)180次，且二者相對位置不變，進而記錄下180

5、組數(shù)據(jù)，經(jīng)成像過程生成具體圖像。然而CT系統(tǒng)安裝時存在的誤差會影響成像質(zhì)量，因此需要對安裝好的CT系統(tǒng)進行參數(shù)標定。</p><p>　　問題1： 現(xiàn)已給定標定模板的位置、幾何信息、該介質(zhì)的吸收率和經(jīng)過掃描過程探測器的接收信息。需要根據(jù)已知條件求解出CT系統(tǒng)旋轉(zhuǎn)中心的位置、探測器每一個單元之間的距離和X射線入射的180個方向。</p><p>　　問題2：根據(jù)問題1的解，可以完全確定CT系

6、統(tǒng)的工作狀態(tài)，現(xiàn)再提供某未知介質(zhì)的接收信息。憑借CT系統(tǒng)各參數(shù)的值和該未知介質(zhì)最后得到的接收信息結(jié)果，確定出該未知介質(zhì)的位置、形狀、吸收率等信息。并在該吸收率情況下找出10個位置的吸收率。 </p><p>　　問題3：同問題2，利用CT系統(tǒng)各參數(shù)的值和接收信息，確定該未知介質(zhì)的相關(guān)信息。并求出10個位置的該未知介質(zhì)下的吸收率。</p><p>　　問題4：分析問題1求解的參數(shù)標定的精度和

7、穩(wěn)定性。并自行設(shè)計新模板、建立對應(yīng)的標定模型改進標定精度和穩(wěn)定性。</p><p><b>　　二、問題分析</b></p><p>　　本題需要解決的問題可分為兩部分：對CT系統(tǒng)進行完整的參數(shù)標定，繼而用CT系統(tǒng)實現(xiàn)成像功能。 </p><p>　　2.1系統(tǒng)參數(shù)標定模型</p><p>　　參數(shù)標定從本質(zhì)上講是成像過

8、程的逆推過程。即已知模板幾何尺寸，位置信息，接收信息，根據(jù)已知信息可以確定發(fā)射-接收系統(tǒng)的工作狀態(tài)。接收器的單元間距未知，但是從已知信息可以推算出模板投射在探測器的投影單元數(shù)，又知道模板的實際尺寸，就可知二者間換算關(guān)系，從而確定單元間距。旋轉(zhuǎn)中心難以直接計算，但因為吸收率一定，所以探測器顯示的吸收信息蘊含著物體的長度和寬度，再輔以旋轉(zhuǎn)180次得到的不同角度的信息，結(jié)合幾何關(guān)系，可以求解。角度問題可以直接對應(yīng)探測器吸收信息的變化情況求得。

9、</p><p><b>　　2.2成像模型</b></p><p>　　標定各參數(shù)之后，整個發(fā)射-接收系統(tǒng)的工作狀態(tài)就可以完全確定。CT系統(tǒng)的工作原理為CT掃描被測物體后得到吸收信息，進而反解出被測物體的其他相關(guān)信息?？梢园l(fā)現(xiàn)該工作過程可以對應(yīng)直接傅里葉反變換法的計算過程，因此，可以用此方法方便快捷地求得正確的衰減系數(shù)。得到被測物體的衰減系數(shù)，就可以根據(jù)各參數(shù)的實際

10、意義與衰減系數(shù)的實際意義之間的聯(lián)系推算出其余信息，包括位置、形狀、吸收率等參數(shù)。已知整個未知介質(zhì)的吸收率后，即可求得同坐標系上任一點的吸收率。</p><p>　　2.3精度和穩(wěn)定性分析</p><p>　　由于建立模型和求解模型中將部分模型理想化，因此需要在最后做一個精度和穩(wěn)定性的分析，以避免誤差過大，影響模型的使用和推廣。</p><p><b>　　

11、三、模型假設(shè)</b></p><p>　　1.橢圓投影長度最長（短），對應(yīng)X射線與橢圓相切于長軸（短軸）。</p><p>　　2.投影長度相等的若干組數(shù)據(jù)，取中間一組為極值。</p><p>　　3.求未知介質(zhì)位置時，可看作規(guī)則圖形求其中心位置。</p><p><b>　　四、定義與符號說明</b><

12、;/p><p>　　五、模型的建立與求解</p><p>　　5.1 CT系統(tǒng)參數(shù)標定模型</p><p>　　5.1.1吸收率與接收信息</p><p>　　借助模板標定CT系統(tǒng)的參數(shù)，首先需要正確理解各參數(shù)的實際意義并明確相互之間的關(guān)系。</p><p>　　由題可知，附件1中的數(shù)據(jù)反映的是該點的吸收強度，即吸收率。通

13、過觀察可以發(fā)現(xiàn)，吸收率的數(shù)值分別為0.0000和1.0000，且1.0000對應(yīng)的點的分布具有一定的規(guī)律性。為了更直觀地觀察其分布情況，利用MTALAB讀取表格數(shù)據(jù)并繪制圖像（程序見附錄1），如圖5.1.1所示。</p><p>　　圖5.1.1 吸收率</p><p>　　圖中高亮的部分為吸收率為1.0000的點所在區(qū)域，其余部分是數(shù)值為0的區(qū)域。將圖5.1.1與模板的幾何信息對比可以發(fā)

14、現(xiàn)，吸收率非0的區(qū)域形狀、尺寸、位置均與模板匹配，即模板確為均勻介質(zhì)，且該介質(zhì)的吸收率為1.0000。</p><p>　　接收信息定義為X射線穿過模板的過程中，被模板所吸收的信息，反映在附件2，為512×180的表格。由CT的工作原理可知，某一列的512組數(shù)據(jù)即為X射線從某一個角度發(fā)射，探測器512個單元接收到的信息。某一行的180組數(shù)據(jù)代表該行對應(yīng)的單元，隨著反射器轉(zhuǎn)動180次而分別得到的數(shù)據(jù)。&l

15、t;/p><p>　　由CT的圖像重建原理可知，接收信息來源于X射線照射在被測物體上，部分能量被物體吸收而導(dǎo)致的射線強度衰減。這種衰減呈指數(shù)變化，如果發(fā)射器發(fā)射出的射線是單能的，并且照射在了均勻材料的物體上，則遵循Lambert-Beers吸收定律[2]：</p><p><b>　?。?）</b></p><p>　　μ為衰減系數(shù)，d為被測物體的

16、厚度，為X射線穿過被測物體前的初始強度，I為X射線穿過被測物體后的射出強度。</p><p>　　由吸收率的含義可知，吸收率a和衰減系數(shù)μ保持一定的比例關(guān)系：</p><p><b>　?。?）</b></p><p>　　若定義P為接收信息量，且滿足式（3）：</p><p><b>　　（3）</b&

17、gt;</p><p>　　根據(jù)吸收率和接收信息的實際意義，可知接收信息P的數(shù)值大小與吸收率a以及被測物體厚度d有關(guān)。聯(lián)立（1）、（2）、（3）式，可得三者之間的函數(shù)關(guān)系：</p><p><b>　?。?）</b></p><p>　　分析可知，吸收率一定的情況下，被測物體越厚，則接受信息的數(shù)值越大。</p><p>

18、　　5.1.2探測器單元之間的距離</p><p>　　求解探測器單元之間的距離，需要將模板的真實尺寸與其投影在探測器的長度進行比較。投影在探測器上的長度與附件2的數(shù)據(jù)相關(guān)。每一列的數(shù)據(jù)為X射線從某一個角度發(fā)射，探測器的512個單元接收到的信息。由接收信息的含義可知，數(shù)值為0的單元的X射線是從發(fā)射器直接發(fā)射到探測器的，不經(jīng)過模板的任意部位。而數(shù)值不為0的單元對應(yīng)的射線，不論數(shù)值大小，均穿過模板。因此，只要確定出每

19、一列的非零數(shù)據(jù)的數(shù)量，就知道有多少個單元的射線穿過模板。</p><p>　　由于使用Excel對每一列接收信息數(shù)據(jù)進行觀察的過程比較繁瑣，而且容易出現(xiàn)錯誤，因此，可以運用MATLAB導(dǎo)入數(shù)據(jù)，編寫程序得到非零數(shù)據(jù)的數(shù)量。由于某些列的非零數(shù)據(jù)分成了兩部分，所以所編寫的程序?qū)刹糠植痖_計數(shù)，最終得到的結(jié)果為兩組數(shù)據(jù)（程序見附錄2）。得到數(shù)據(jù)后進行觀察，發(fā)現(xiàn)其中一組數(shù)據(jù)絕大多數(shù)為29，因此可以確定這組數(shù)據(jù)是射線穿過圓

20、形模板得到的，29個單元對應(yīng)的就是圓形模板的直徑8mm。則探測器單元之間的距離L為： (5)</p><p>　　5.1.3 CT系統(tǒng)旋轉(zhuǎn)中心在正方形托盤中的位置</p><p>　　由于CT系統(tǒng)的旋轉(zhuǎn)中心一定在探測器的中垂線上，因此找

21、任意兩個探測器位置，分別做其中垂線，相交的點即為旋轉(zhuǎn)中心。</p><p>　　由5.1.2可知，探測器單元之間的距離已知，因此，取附件2的任意列，從上向下查到非零數(shù)據(jù)點的時候，有m組數(shù)據(jù)就對應(yīng)著m條射線沒有穿過模板。又由5.1.1可知模板投影到探測器對應(yīng)的單元個數(shù)n，且總單元數(shù)為512，由此可知模板數(shù)據(jù)下方的零數(shù)據(jù)點共有512-m-n個。</p><p>　　如圖5.1.2所示，AF為探

22、測器,B為正方形托盤左下角，也是坐標原點，C和E分別與橢圓相切，D點為探測器AF的中點。根據(jù)附錄2的數(shù)據(jù)，可以找到一組最小，一組最大數(shù)量的數(shù)據(jù)，最小數(shù)量為109組數(shù)據(jù)，為CE段，對應(yīng)上方零數(shù)據(jù)點有168個，為EF段，下方零數(shù)據(jù)點有235個，為AC段，對應(yīng)的是橢圓的短軸，即AC:CE:EF三者比值為235:109:168；最大數(shù)量為289組數(shù)據(jù)，上方零數(shù)據(jù)點有69個，下方零數(shù)據(jù)點有134個，對應(yīng)的是橢圓的長軸。</p>&l

23、t;p>　　以正方形左下角B為原點建立直角坐標系，則旋轉(zhuǎn)中心O的坐標可表示為：</p><p><b>　?。?）</b></p><p>　　同理，縱坐標可表示為：</p><p>　　綜上，旋轉(zhuǎn)中心坐標為（40.7931,56.0690），如圖5.1.2所示。</p><p>　　圖5.1.2 旋轉(zhuǎn)中心<

24、/p><p>　　5.1.4 X射線的180個方向</p><p>　　由5.1.3可知，最大數(shù)量的是第61組數(shù)據(jù)，最小數(shù)量的是第151組數(shù)據(jù)，并且這兩組數(shù)據(jù)分別對應(yīng)的是水平射線和豎直射線，即夾角為90度。</p><p>　　假設(shè)180次旋轉(zhuǎn)，每次轉(zhuǎn)動的角度相同，則：</p><p><b>　　(7)</b></p

25、><p>　　由于發(fā)射-接收系統(tǒng)繞逆時針旋轉(zhuǎn)，根據(jù)附件2的數(shù)據(jù)可推斷出，投影長度先增加后減小再增加，因此初始位置是從正方形托盤的右下方發(fā)射，左上方接收。又因為橢圓形長軸投影對應(yīng)的是第61組數(shù)據(jù)，即從初始位置旋轉(zhuǎn)了60次，因此，初始位置的射線與水平方向夾60度角。</p><p>　　綜上，X射線從與水平方向夾60度角開始旋轉(zhuǎn)，每次旋轉(zhuǎn)1度。</p><p>　　5.2

26、CT系統(tǒng)成像模型</p><p>　　5.2.1直接傅里葉反變換法（DFR）求衰減系數(shù)μ</p><p>　　直接傅里葉反變換法是由不同角度的投影的數(shù)據(jù)，恢復(fù)原始圖像的方法[3]。</p><p>　　由于衰減系數(shù)μ在不均勻介質(zhì)中會隨著位置的不同而變化，因此接收信息P與μ的關(guān)系滿足：</p><p><b>　　(8)</b&

27、gt;</p><p>　　其中d為X射線穿過物體的路徑，μ（x）表示d上某點x處的衰減系數(shù)。可以看出，投影過程和X射線照射過程均為線積分，因此由附件3中的接收信息P可以使用直接傅里葉反變換法方便地求解出衰減系數(shù)μ的數(shù)據(jù)。</p><p>　　具體過程為：首先原始圖像線積分得到投影，再由投影值求一維傅里葉變換，把變換所得當作二維傅里葉變換的一個切片。所有角度的切片數(shù)據(jù)放在一個極坐標系下，然

28、后將極坐標系柵格化，即使用二維插值求出笛卡爾坐標系中相對應(yīng)的值，再進行二維傅里葉逆變換，得到原始圖像。具體流程如圖5.2.1所示。</p><p>　　二維傅里葉逆變換 一維傅里葉變換</p><p>　　圖5.2.1 流程圖</p><p>　　分別導(dǎo)入附件3和附件

29、5，求出衰減系數(shù)μ的數(shù)據(jù)表格并生成圖像，發(fā)現(xiàn)數(shù)據(jù)量較大，因此調(diào)用MATLAB中的imresize函數(shù)運用插值法可進行一定程度的縮減，最終如圖5.2.2。（程序見附錄3）設(shè)定顏色越暗則衰減系數(shù)越小，反之，越亮則衰減系數(shù)越大。</p><p>　　圖5.2.2 附件3、附件5衰減系數(shù)圖</p><p>　　5.2.2位置和幾何形狀</p><p>　　根據(jù)圖5.2.2觀

30、察可知，附件5的形狀不規(guī)則，外邊界近似為一個正方形，內(nèi)部有許多暗色區(qū)域，分布無規(guī)律，因此很難談?wù)撈湫螤钜约拔恢眯畔ⅰ?lt;/p><p>　　附件3的幾何形狀近似為橢圓，因此，可用類似求解橢圓中心的方法求解該未知介質(zhì)的位置問題。與附錄2程序功能類似，用MATLAB編寫程序?qū)⒏郊?的每一列數(shù)據(jù)中非零數(shù)據(jù)的數(shù)量算出，并以該數(shù)值為縱坐標，對應(yīng)旋轉(zhuǎn)次數(shù)為橫坐標建立直角坐標系，如圖5.2.3（程序見附錄4）。從圖5.2.3中可

31、以清晰地看出第58組的非零數(shù)最多，第146組的非零數(shù)最少，由于該組數(shù)非常接近第61和151組，因此可近似認為這二者對應(yīng)的投影長度分別為橢圓長軸和短軸投影長度。</p><p>　　圖5.2.3投影長度—旋轉(zhuǎn)次數(shù)</p><p>　　由于設(shè)備和旋轉(zhuǎn)中心均不變，因此探測器在第61和151組的位置不變。豎直放置的探測器非零數(shù)值的數(shù)據(jù)量為296，其上方的零數(shù)值數(shù)據(jù)量有96，下方的零數(shù)值數(shù)據(jù)量為12

32、0，由5.1.3的幾何關(guān)系可知，最下端到x軸的距離為10-69×=-9.0345mm；水平放置的探測器非零數(shù)值得數(shù)據(jù)量為156，其上方的零數(shù)值數(shù)據(jù)量有215，下方的零數(shù)值數(shù)據(jù)量為141，同理可求得最左端到y(tǒng)軸的距離為35-235×=29.8276mm。所以，該介質(zhì)的位置，即橢圓中心的坐標為：</p><p>　　x=215×-29.8276=50.9977mm

33、 (9)</p><p>　　y=+96×-9.0345=58.2783mm</p><p>　　綜上，未知介質(zhì)的位置，即中心坐標為（50.9977,58.2783）。</p><p>　　由5.2.1可得到衰減系數(shù)的數(shù)據(jù)，并可清晰地觀察分析出衰減系數(shù)在不同位置的大小。分析圖5.2.2可知，明亮部分衰減系數(shù)較大，說明此區(qū)域的形狀為該介質(zhì)的幾何形狀。但

34、具體的幾何形狀難以通過觀察分析得到，因此，可以取衰減系數(shù)值恰當，在圖像中處于明和暗的交界處的點(假設(shè)為此邊界衰減系數(shù)為0.15)，進行曲線擬合，得到具體的函數(shù)關(guān)系，以此來表達該介質(zhì)的幾何形狀。</p><p>　　但由于該圖形為封閉式圖形，無法直接用函數(shù)的擬合方式進行，故可運用MATLAB編寫程序?qū)⒁粋€閉環(huán)拆分成上下兩部分開環(huán)，分別進行擬合，用兩個式子來完整表達形狀（程序見附錄5）。另外通過對圖5.2.2的觀察可

35、發(fā)現(xiàn)，介質(zhì)內(nèi)部有兩塊區(qū)域衰減系數(shù)較小，甚至接近于外部區(qū)域，因此采用數(shù)值比較的方法將此兩塊區(qū)域和其他內(nèi)部區(qū)域分開考慮。而這兩塊區(qū)域的邊界也是封閉式的，同樣需要拆分開進行函數(shù)關(guān)系的擬合。綜上，可擬合出以下六式（圖5.2.3）：</p><p><b>　　（10）</b></p><p><b>　?。?1）</b></p><p

36、>　　圖5.2.4 外邊界擬合曲線</p><p><b>　?。?2）</b></p><p><b>　?。?3）</b></p><p><b>　?。?4）</b></p><p><b>　　（15）</b></p><

37、p>　　圖5.2.5 內(nèi)邊界擬合曲線</p><p>　　綜上所述，該未知介質(zhì)的幾何形狀即為（10）~（15）式聯(lián)立后圍住的區(qū)域，可由圖5.2.5觀察分析其形狀特征。</p><p><b>　　5.2.3吸收率</b></p><p>　　將5.1當中模板的吸收信息按照直接傅里葉反變換法進行計算，將會得到其對應(yīng)的衰減系數(shù)數(shù)據(jù)以及圖像（如

38、圖5.2.6）。</p><p>　　圖5.2.6 模板衰減系數(shù)圖</p><p>　　由5.1可知，衰減系數(shù)和吸收率呈正比關(guān)系，則通過現(xiàn)在模板的衰減系數(shù)和吸收率以及未知介質(zhì)的衰減系數(shù)，可以推測出未知介質(zhì)的吸收率。具體的做法是：從模板的衰減系數(shù)數(shù)據(jù)中判斷出橢圓位置，因為橢圓模板內(nèi)部的吸收率均為1.0000，因此可從橢圓內(nèi)部任取足夠多的點求其平均值，現(xiàn)取DF-EM列，76-109行數(shù)據(jù)，共1

39、156個，得到。由于=1.0000是定值，因此，未知介質(zhì)的吸收率：</p><p><b>　?。?6）</b></p><p>　　用Excel軟件打開附件3、附件5對應(yīng)的衰減函數(shù)圖，每一個衰減系數(shù)乘一個固定的數(shù)值1.6040即可得到最終需要獲得的該介質(zhì)的吸收率數(shù)據(jù)表，具體數(shù)據(jù)見附件problem2.xls及problem3.xls。</p><

40、p>　　5.2.4 “空點”的篩除</p><p>　　由于探測器的單元間隔較大，投影次數(shù)太少等因數(shù)的影響，會導(dǎo)致探測器探測到的吸收信息有限，從而使得利用這些信息不足以完全還原二維情形。所以才會出現(xiàn)吸收信息矩陣中有0值，但是變換后的系數(shù)矩陣沒有0值的情況。對應(yīng)的實際情況即沒有物體的區(qū)域衰減系數(shù)也有非零值。所以為了精確得出二維真實情形可設(shè)計一個閾值來篩除‘’空點‘’。</p><p>

41、　　利用Excel表格來確定衰減系數(shù)矩陣的最大值M、最小值m。將左右端點值分別為m、M的區(qū)間等分為N個小區(qū)間，第i個小區(qū)間記為 ，則第i個小區(qū)間的右端點值為。再以落在某個小區(qū)間內(nèi)元素的個數(shù)作為該小區(qū)間的函數(shù)值。運用MATLAB繪出將橫坐標區(qū)間等分成10分的圖像。（如圖5.2.7）</p><p>　　圖像大致分為三段，第一段函數(shù)隨著橫坐標值增加，函數(shù)值迅速減小。該段表示此區(qū)間內(nèi)的點絕大部分的吸收系數(shù)接近0，并且整

42、體變化一致。故第一段應(yīng)該對應(yīng)“空點”即不是實體上的點。第二段前半部分曲線緩慢減小，且函數(shù)值很低，后半部分幾乎不變。這段圖像前半部分說明這類點數(shù)目少但吸收系數(shù)分布區(qū)間長，這類點出現(xiàn)的原因可從圖5.2.2看出，在靠近橢圓的黑色區(qū)域有“白線”存在，說明這些區(qū)域的點較黑色區(qū)域有較高的值。但這些點是由于radon逆變換時產(chǎn)生的干擾，它不反應(yīng)真實物體情況，所以設(shè)計的閾值應(yīng)該排除掉這些點。第二段曲線后半部分可以認為真實點開始變多，從而曲線變平緩。第三

43、段為上升段，該區(qū)域的點應(yīng)為實體上的點，應(yīng)該保留。綜上所述，我們選擇當?shù)诙吻€的變化剛接近平緩的點0.413為閾值。在Excel中用快速分析功能的色階分析表格可以大致得出類似于圖5.2.2的圖像如下，而后在其基礎(chǔ)上令值大于閾值的格子填充為粉色。另外再作一個圖，先使值大于閾值的格子填充為粉色，再用快速分析功能的色階描述表格，使后兩圖相互印證。</p><p>　　未填充的圖

44、 先色階分析再用粉色填充的圖 先用粉色填充再用色階分析的圖</p><p>　　如圖可看出粉色區(qū)域幾乎覆蓋了綠色區(qū)域，但邊緣還有一些綠色，也就說明閾值有點偏高。在篩選時可能篩掉實體點，但可以近似認為閾值是合適的。用閾值對衰減系數(shù)矩陣進行篩選。</p><p>　　對于衰減矩陣中的值小于閾值令其為0，大于等于0則保留原值。</p><p>　　

45、第三題閾值尋找方法如第二題所述，并作出N=10的曲線。</p><p>　　圖線大致分為兩端，第一段如第二題中所述對應(yīng)“空點”第二段連續(xù)而緩慢上升到最高，說明此時圖像干擾點不多可以忽略，此時可以只用把第一段對應(yīng)點排除。綜上所述以0.067為閾值。先在Excel中用快速分析功能的色階分析表格可以大致得出類似于圖5.2.2 的圖像如下第一個圖所示，而后在其基礎(chǔ)上令值大于閾值的格子填充為黃色如下第二個圖所示 。另外再作

46、一個圖，先使值大于閾值的格子填充為黃色，再用快速分析功能的色階描述表格，使后兩圖相互印證。</p><p>　　未填充的圖 先用色階分析再用黃色填充的圖 先用黃色填充再用色階分析的圖</p><p>　　從第二張看出黃色基本覆蓋了綠色，第三張除了一些邊緣綠色基本覆蓋了黃色，這說明該閾值有點低，可能會在

47、篩除時漏過一些“空點”但基本還適合。用與第二題同樣的法則對衰減系數(shù)矩陣進行篩選</p><p>　　5.2.5十個位置處的吸收率</p><p>　　由于X射線的初始位置與豎直方向夾30度角，因此經(jīng)過直接傅里葉反變換產(chǎn)生的圖像并非實際方向，而是與豎直方向夾30度角的。因此，想要求出10個位置的吸收率，首先需要把這些點的坐標轉(zhuǎn)換到Excel表格對應(yīng)的坐標系——衰減率坐標系中。</p&g

48、t;<p>　　已知模板的形狀和位置，即橢圓長軸平行于y軸，因此，坐標系可共用一個模板，即圍繞這個模板作兩個坐標系。但由于兩坐標系的原點不在同一位置，故需要采取一個中間坐標系進行轉(zhuǎn)換。因為在衰減率坐標系中難以找到橢圓中心，因此可選取圓形的圓心作為中間變換坐標系的原點。</p><p>　　利用Excel快速分析中的色階分析命令分析，可以找到圓形模板在表中所對應(yīng)的位置，如圖5.2.7所示。觀察圓的邊緣

49、，當某一行的數(shù)據(jù)出現(xiàn)突變時，可以認為該行所在直線是圓的切線，記錄行數(shù)。同理可以找到另一條切線并記錄該切線所在行數(shù)，由此可知圓的直徑大小為兩行行數(shù)之差。同樣，可以找到與圓相切的兩列，得到直徑大小。最后以二者的平均值作為圓的最終直徑，兩列序號的中間值為圓心的縱坐標，兩行序號的中間值為圓心的橫坐標。</p><p>　　圖5.2.7 圓形模板位置</p><p>　　以正方形托盤的左下角為原點o

50、，水平向右為x軸，豎直向上為y軸建立平面xoy坐標系。模板中的圓形模板在此坐標系下的圓心位置為（95,50）。再以Excel表格的左上角（1，A）為原點o2，水平向右為v軸，豎直向下為u軸建立平面uo2v坐標系。圓心在此坐標系下的坐標為（，）。</p><p>　　為了能實現(xiàn)10個位置坐標的變換，需要以圓形模板的圓心o1為原點，平行于x軸的軸為x1軸，平行于y軸的軸為y1軸建立x1o1y1坐標系。同理建立u1o1

51、v1坐標系。坐標系的建立如圖5.2.8所示。</p><p><b>　　圖5.2.8 坐標</b></p><p>　　假設(shè)A點為10個點中的任一點，在xoy坐標系下的坐標為（），在x1o1y1坐標系下的坐標為（），在uo2v坐標系下的坐標為（），在u1o1v1坐標系下的坐標為（）。</p><p>　　已知圓形圓心o1在xoy坐標系下的坐標

52、為（95,50），由向量相減準則可知：</p><p><b>　　整理得：</b></p><p><b>　?。?7）</b></p><p>　　因為x1o1y1和u1o1v1坐標系原點均為o1，夾角為-120度，故可用坐標變換公式求得（）：</p><p>　　整理得：

53、 </p><p><b>　?。?8）</b></p><p>　　再運用向量相加準則：</p><p><b>　?。?9）</b></p><p>　　經(jīng)過（17）、（18）、（19）式，可將10個點的坐標變換到u02v坐標系中，取整后查找Excel表格中對應(yīng)坐標的數(shù)值，即為該點

54、的吸收率。最后，得到坐標變換前后坐標及10個點的吸收率。如表1所示。</p><p>　　表一 10個點的坐標變換及吸收率</p><p><b>　　六、模型評價</b></p><p>　　首先，第一題有3個小問題，我們分別分析3個問題的精確度</p><p>　　第一個是求探測器單元之間的距離，我們是通過小球的長

55、度除以投影的數(shù)量(射線通過小球投影到探測器單元的數(shù)量)即8/29.但這個數(shù)就存在誤差。如圖所示，只有如圖一，單元格對應(yīng)的射線左側(cè)和小球相切時小球的直徑才對應(yīng)射線通過小球投影到探測器單元的數(shù)量，而正常情況下如圖2，顯然單元格射線線左側(cè)和和小球之間存在距離，這就是是誤差。而當圖三所示探測器單元的右側(cè)和小球相切時即為臨界值，即為誤差最大值，此時射線通過小球投影到探測器單元的數(shù)量其實為27個，探測器單元之間的距離為8/27。綜上所述探測器單元之

56、間的距離x為8/29<=x<8/27，所以它的誤差為8/27-8/29=0.0204.</p><p>　　圖一圖二</p><p><b>　　圖三</b></p><p>　　第二個是求X射線的180個方向，我們是通過求射線水平射入橢圓的組數(shù)和垂直射入的組數(shù)差即為90度，即投影最大時(80)是水平射入，投影為30

57、時是垂直切入,但是無法準確找到投影為80和30的點。只能找到79.7241和30.0690這兩個點，這是誤差一，而這2個點本身就是通過射線通過橢圓投影到探測器單元的數(shù)量除以探測器單元之間的距離得到的，而探測器單元之間的距離就有誤差，這是誤差二。而且有8組投影是79.7241,分別為58-65組，4組是30.0690，分別為149,151,153,154組。雖然我們是通過2組數(shù)據(jù)的中間值相減，但還是存在誤差，這是誤差三。通過誤差三可以求出

58、每次旋轉(zhuǎn)的方向角a范圍為0.9375<a<1.0714, 誤差為1.0714-0.9375=0.1339。入射角（和水平方向的夾角）的范圍為57.1875~65.3554，誤差為8.1679。</p><p>　　第三個是求CT系統(tǒng)旋轉(zhuǎn)中心在正方形托盤中的位置。我們是通過求射線水平射入探測器時和垂直射入時中垂線的交點即為旋轉(zhuǎn)中心。這里用到探測器單元之間的距離，因為這個值的范圍為8/29~8/27，存在

59、誤差。且有8組水平投影最接近80，4組垂直投影接近30，找到4組中2個臨界點，得出5.1的算法的AC:CE:EF=236:109:167和234:109:169，和8組中2個臨界點比例為132:289:91和136:289:87。最后可以得出最后旋轉(zhuǎn)中心的橫坐標的范圍為40.5172~41.5185和縱坐標的范圍為55.5172~60.0741，誤差分別為1.0013和4.5569。</p><p>　　求穩(wěn)定度

60、時，我們分析附件2存在誤差，即在每一列中0和非0數(shù)據(jù)之間臨界值誤差。因此我們把每列中0和非0相鄰的0變成非0。即把每個角度射線通過圖像投影到探測器單元的數(shù)量加2（上面和下面分別加1），求出每個角度投影的長度，然后得出問題一的結(jié)果。分析這個結(jié)果和我們求得的結(jié)果是否相似，以此來確定它的穩(wěn)定性。</p><p>　　首先我們把每個角度射線通過圖像投影到探測器單元的數(shù)量加2后，通過小球投影到探測器單元的數(shù)量由29變成31

61、，所以探測器單元之間的距離為8/31。然后是求X射線的180個方向，垂直射入時的組數(shù)是不變的(因為是找投影最大值)，仍然為為58-65組，但水平射入時有變化，根據(jù)附錄程序可求出第159組投影最接近30，為29.9355。以61組垂直射入時的組數(shù)，因此可求得旋轉(zhuǎn)角度為90/ (159-61)=0.9184，入射角（和水平方向的夾角）為(61-1)*0.9184=55.1020。最后求CT系統(tǒng)旋轉(zhuǎn)中心在正方形托盤中的位置，以61組為水平射入

62、方向，159組為垂直射入方向，得出5.1的算法的AC:CE:EF分別為228:116:168和133:291:88。計算出旋轉(zhuǎn)中心為(42.2258，53.3548)。</p><p>　　第一題求得的4個數(shù)據(jù)分別為距離8/29、旋轉(zhuǎn)度為1度、入射角為60度、旋轉(zhuǎn)中心(40.7931,56.0690)，因此，探測器單元之間的距離的相對誤差為(8/29-8/31)/(8/29)=0.0645, 旋轉(zhuǎn)角度的相對誤差為

63、為（1-0.9184）/1=0.0816，入射角的相對誤差為(60-55.1020)/60=0.0816，旋轉(zhuǎn)中心橫坐標的相對誤差為(42.2258-40.7931)/ 40.7931=0.0351縱坐標的相對誤差為(53.3548-56.0690)/ 56.0690=-0.0484。因此，當初始數(shù)據(jù)有擾動時，參數(shù)依然可以標定相對穩(wěn)定。</p><p><b>　　七、總結(jié)</b></

64、p><p><b>　　八、參考文獻</b></p><p>　　[1]郭立倩. CT系統(tǒng)標定與有限角度CT重建方法的研究[D].大連理工大學(xué),2016.</p><p>　　[2] 盧彥斌. X射線CT成像技術(shù)與多模態(tài)層析成像技術(shù)研究[D].北京大學(xué),2012.</p><p>　　[3]莊天戈. CT原理與算法[M].

65、上海交通大學(xué)出版社, 1992.</p><p><b>　　九、附錄</b></p><p>　　附錄1：根據(jù)吸收信息求吸收率</p><p>　　data=xlsread('A題附件.xls','附件1');</p><p>　　imshow(data);</p><

66、;p>　　title('吸收率圖像');</p><p>　　附錄2：橢圓和圓投影所占的探測器單元數(shù)量</p><p>　　data=xlsread('A題附件.xls','附件2');</p><p>　　max=[];%定義放橢圓投影的數(shù)組</p><p>　　min=[];%定義放圓

67、投影的數(shù)組</p><p>　　for i = 1:size(data, 2) %180列循環(huán)</p><p>　　num1=0; %記錄2個圖之一的投影</p><p>　　num2=0; %

68、記錄2個圖之一的投影</p><p>　　flag=0; %標志位</p><p>　　for j = 1:size(data, 1) %512行循環(huán)</p><p>　　if data(j,i)~=0</p><p><b>　　i

69、f flag<2</b></p><p><b>　　flag=1;</b></p><p><b>　　end</b></p><p>　　if flag==1</p><p>　　num1 =num1+1; %如果flag為1的時候num1+1<

70、;/p><p>　　else %如果flag不為1的時候num2+1</p><p>　　num2 =num2+1;</p><p><b>　　end</b></p><p>　　elseif flag==1 </p><p><b>

71、　　flag=2;</b></p><p><b>　　else </b></p><p><b>　　end</b></p><p><b>　　end</b></p><p>　　if num1>num2

72、 %如果num1大于num2，num1屬于max，num2屬于min</p><p>　　max(i)=num1;</p><p>　　min(i)=num2;</p><p>　　else %反之，num2屬于max，num1屬于min</p><p>　　max(i)

73、=num2;</p><p>　　min(i)=num1;</p><p><b>　　end</b></p><p><b>　　end </b></p><p>　　disp(max) %分別輸出max和min數(shù)組</p&g

74、t;<p><b>　　disp(min)</b></p><p>　　附錄3：吸收信息轉(zhuǎn)換為衰變系數(shù)</p><p>　?。?）data=xlsread('A題附件.xls','附件3'); %讀入附件3的數(shù)據(jù)變成矩陣</p><p>　　show_img(data);

75、 %進入函數(shù)</p><p>　?。?）function [ output_args ] = show_img( x )</p><p>　　p_N=256; %圖像默認大小256</p><p>　　theta_N=180;

76、 %180度平行投影</p><p>　　pad_N=1024; %投影后變?yōu)?67*180，每列數(shù)據(jù)就是當前角度下的投影值。367<512考慮到傅里葉變換是需要基2對其，故擴展為1024*180</p><p>　　theta=0:(

77、theta_N-1); %0~180度平行投影</p><p>　　output_args=0;</p><p>　　%% 根據(jù)每個投影角度下的投影數(shù)據(jù)，求其一維傅里葉變換</p><p>　　proj_sino=x;</p><p>　　if mod(length(pro

78、j_sino(:,1)),2)==1 %如果奇數(shù)個接受柵格，則需要擴展為偶數(shù)，填充0</p><p>　　proj_sino=[proj_sino:zeros(1,size(proj_sino,2))]; %填充一行，367*180變?yōu)?68*180</p><p><b>　　end</b></p><p&

79、gt;　　pad_row=(pad_N-size(proj_sino,1))/2; %（1024-368）/2=328 </p><p>　　proj_sino=padarray(proj_sino,[pad_row 0],0,'both');

80、 % 368*180,上下填充328行，變?yōu)榱?024*180</p><p>　　L_pad=pad_row+ceil(((p_N.*sqrt(2)+2)-p_N)/2)+1; %原始圖像大小為p_N,傅里葉逆變換后為pad_N,則需要截斷數(shù)據(jù)</p><p>　　proj_sino=ifftshift(proj_sino,1); %將實際數(shù)據(jù)移動到兩端，

81、中間的是填充的0</p><p>　　f_p=fft(proj_sino,[],1); %求取正弦圖的一維傅里葉變換，1024*180</p><p>　　f_p=fftshift(f_p,1); %反向處理，中間的數(shù)據(jù)移動到兩端，兩端的數(shù)據(jù)移動到中間</p><p>　　%%極坐標化為笛卡爾坐標系</

82、p><p>　　nfp=length(f_p(:,1));</p><p>　　omega_sino=(-(nfp-1)/2:(nfp-1)/2).*(2*pi/size(f_p,1)); </p><p>　　%極半徑范圍 -pi~pi，分成1024等分</p><p>　　theta=theta*pi/180;

83、 %角度轉(zhuǎn)化為弧度，范圍0~pi,180等分</p><p>　　[theta_grid,omega_grid]=meshgrid(theta,omega_sino); </p><p>　　%網(wǎng)格后，omega_grid和theta_grid都是1024*180，theta_grid_x每一列都是一樣的角度，omega_grid沒一行都是一樣的位置<

84、;/p><p>　　omega_image=omega_sino; %根據(jù)極半徑大小，建立笛卡爾坐標系</p><p>　　[omega_grid_x,omega_grid_y]=meshgrid(omega_image,omega_image); </p><p>　　%網(wǎng)格后，omega_grid和theta_g

85、rid都是1024*1024 omega_grid_x每一列的x坐標一樣，omega_grid_y每一行的y坐標一樣</p><p>　　[coo_th,coo_r]=cart2pol(omega_grid_x,omega_grid_y);</p><p>　　coo_r=coo_r.*sign(coo_th); %第一象限和第二象限內(nèi)，笛卡爾半徑為負

86、</p><p>　　coo_th(coo_th<0)=coo_th(coo_th<0)+pi; %第四象限和第二象限內(nèi)，笛卡爾角度為0~pi/2；第一象限和第三象限，笛卡爾角度為pi/2~pi</p><p>　　Fourier2=interp2(theta_grid,omega_grid,f_p,coo_th,coo_r,'cubic',(0+1i.*0)

87、);</p><p>　　%根據(jù)極坐標處的值，二維插值得出笛卡爾坐標處對應(yīng)的值，插值后也是1024*1024的復(fù)數(shù)矩陣</p><p>　　crop=L_pad;</p><p>　　target=fftshift(ifft2(ifftshift(Fourier2))); %二維傅里葉逆變換，得到的仍是1024*1024的矩陣</p>

88、<p>　　target=target(crop+1:end-crop,crop+1:end-crop); %原始數(shù)據(jù)大小256*256，因此需要對上面的數(shù)據(jù)截斷</p><p>　　I_a=abs(target); %復(fù)數(shù)的模值</p><p>　　I_a=(I_a-min(I_a(:)))./(

89、max(I_a(:))-min(I_a(:))); %歸一化0~1</p><p>　　Lg_I_a=log(1+I_a); %為了好的顯示效果，區(qū)對數(shù)</p><p>　　figure,subplot(1,2,1),</p><p>　　A=imresize(Lg_I_a,[256 256]

90、); %把400*400矩陣變成256*256</p><p>　　%xlswrite('心臟.xls',A);</p><p>　　imshow(A); %顯示圖像</p><p>　　% imshow(flipud(A

91、),[ ]);</p><p>　　title('圖像');</p><p>　　xlswrite('心臟灰度.xls',G); %生成表格</p><p><b>　　end</b></p><p>　　附錄4：光線旋轉(zhuǎn)0到180度圖形的

92、投影</p><p><b>　　clc;</b></p><p>　　data=xlsread('A題附件.xls','附件3'); %讀入附件3的數(shù)據(jù)變成180*512矩陣</p><p>　　x=[]; %定義一個

93、一維數(shù)組 記錄180組圖形的投影</p><p>　　for i = 1:size(data, 2) %180列循環(huán)</p><p>　　flag=0; %標志位</p><p>　　max=0;

94、 %記錄每一列有值的數(shù)量</p><p>　　for j = 1:size(data, 1) %512行循環(huán)</p><p>　　if data(j,i)~=0</p><p>　　max =max+1; %如果一個值不等于0 max+1</p><p><

95、b>　　end</b></p><p><b>　　end </b></p><p>　　x(i)=max*8/29; %根據(jù)單位數(shù)量和長度的比求出投影值，并賦給x數(shù)組</p><p><b>　　end</b></p><p><b>　　di