2017國賽國家一等獎a題優(yōu)秀論文3_W

1、CT 系統參數標定及反投影重建成像摘要CT 技術在當代社會已廣泛應用在臨床醫(yī)學、工業(yè)工程等領域。在本題中首先通過建立離散模型并對其簡化從而對CT 儀器參數進行了標定；建立平行束濾波反投影重建模型，通過 Radon 變換及 R-L 濾波器解決了成像問題，并通過對平行束濾波反投影重建模型進行優(yōu)化，使模型具有降噪能力，從而得到了更精確符合實際的圖形與吸收率值。對問題一，由于附件 2 中的數據存在系統增益，通過對附件 2 數據作消除增益（增益系數為 1.776974）處理后，建立離散化模型，由于方程組模型難于求解，故建立一種簡化離散模型，且通過 MATLAB、Mathmatica 等軟件對方程進行求解

2、，解得 180 個方向的精確角度（結果見表 5.1）。通過 Excel 對附件 2 數據使用條件約束進行色階處理，可以一定程度上將小圓與橢圓的掃描情況清晰地顯示。我們利用附件 2 中的數據與模板存在的幾何關系，求得探測器單元間距范圍為 0.2759，0.2796 ，以橢圓幾何中心為原點建立的坐標系 xoy 中，解得旋轉中心坐標為 以.圓以以中 h . 圓 。 對問題二，由于傳統反投影重建算 引入星狀偽跡，我們決定使用基于吸收率的Radon 變化及 R-L 濾波器構成的平行束濾波反投影重建模型。由于系統旋轉中心和幾何中心不重合，使用 Radon 變換時會造成數據缺失，故我們對正方形托盤進行了“鑲

3、邊” 處理。在對數據除去增益（此處增益系數為 2.0033）處理并且在濾波結束后對“鑲邊” 進行去除后，我們對被檢測物體實現了重建并且得到題目所需十個點的吸收率分別為： 0.000、0.9722、0.0036、1.1761、1.0426、1.4652、1.2849、0.0007、0.0000、0.0175，其 對應的坐標見表 5.2。 對問題三，我們對平行束濾波反投影重建模型進行了優(yōu)化，加入了降噪處理函數。使用問題二中模型對附件 5 進行重建后發(fā)現所得重建圖像存在邊界模糊不清等問題。通過對模板進行數據反演后與標準值對比，發(fā)現吸收率數據會受非系統因素影響，存在一定波動。所以本文中通過多項對比選取

4、自適應濾波函數這一最優(yōu)降噪方法對數據進行降噪處理后，減小了波動范圍，并擬合出了邊界更為清晰的圖形（見圖 5.17），實現了重建并且得到題目所需十個點的吸收率分別為：0.0126、2.2902、5.9159、0.0163、0.0823、3.1336、6.0333、0.0000、7.7184、0.0861，其對應的坐標見表 5.3。 對問題四，首先我們基于標定 CT 的各項參數，建立歸一化均方差評價模型從而對精度和穩(wěn)定性進行分析。此外基于原模板存在的精度與穩(wěn)定性的問題，我們建立了新的模板（見圖 5.20、5.21），新模板具有更高的精度與穩(wěn)定性，能夠對該 CT 進行更精確地參數標定。 本題中我們所

5、建立的一系列模型均能夠滿足題目要求，且層層遞進、環(huán)環(huán)相扣，具有較高的精度與穩(wěn)定性。其中本文建立的平行束濾波反投影重建模型還可以由二維推廣至三維，可以用于太空垃圾的形狀確定與分類，航空器在太空中對未知天體的形狀確定等領域。 關鍵詞：離散模型、二維平行波反投影重建、Radon 變化、R-L 濾波器、降噪處理 25 一、問題重述CT 可以在不破壞樣品的情況下，利用樣品對射線能量的吸收特性對生物組織和工程材料的樣品進行斷層成像，由此獲取樣品內部的結構信息。本題 X 射線的發(fā)射器和探測器的相對位置固定不變，整個發(fā)射-接收系統繞位于正方形托盤下方某處旋轉中心逆時針旋轉 180 次。對每一個 X 射線方向，

6、發(fā)射接收裝置裝有 512 個等距單元探測器，用于測量位置固定不動的二維待測介質吸收衰減后的射線能量，并且通過增益等處理方式得到 180 組接收信息。然而由于存在系統誤差，所以需要對安裝好的 CT 系統進行參數標定，通過已知模板對CT 系統的參數進行標定，并根據標定的參數對未知結構的樣品進行成像。 具體問題重述為下： (1) 在正方形托盤上放置兩個均勻固體介質組成標定模板，模板的幾何信息如圖 2 給出， 相應的數據文件見附件 1，其中每一點的數值反映了該點的吸收強度“吸收率”。對應用于模板的接收信息見附件 2。問題一要求根據模板及其接收信息對 CT 系統進行參數標定，確定出此 CT 系統實際的旋

7、轉中心在正方形托盤中的位置、探測器單元間距以及該 CT 系統使用的 X 射線的 180 個方向。 (2) 利用問題一所標定的 CT 系統相關參數以及所建立模型，使用附件 3 所給未知介質的接收信息，確定出該未知介質在正方形托盤中的位置、幾何形狀以及吸收率等信息。另外利用附件 4 中數據給出圖 3 中所給 10 個位置處的吸收率。 (3) 附件 5 為利用該 CT 系統得到了另一未知介質的接收信息。同樣利用問題 1 中的標定參數與模型，給出未知介質的系列相關信息并給出圖 3 中 10 個位置處的吸收率。 (4) 對問題 1 中參數標定的精度以及穩(wěn)定性進行分析，并在此基礎上建立新模板、建立對應的標

8、定模型，以改進標定精度和穩(wěn)定性，并說明理由。 圖 1.1CT 系統示意圖圖 1.2 模板示意圖（單位：mm）圖 1.3 10 個位置示意圖 二、 問題分析2.1 題目整體分析 本題為典型的 X 射線斷層成像1.問題。通篇分析下來，解決問題一，對題設 CT 系統進行準確的參數標定是求解整個問題的關鍵。題目中已經給出用于標定的模型，以及使用該 CT 對模型進行掃描后得到的參數?；舅悸窞榻㈦x散化模型，并利用數學軟件對其進行求解。在對該CT 完成參數標定后，即可使用所得數據及模型對后面的問題進行分析和求解。 2.2 數據的分析與算法處理 由于本題數據量龐大，首先我們要對與數據進行預處理。特別是附件

9、 2 中找出題目所給數據與被測模型之間的關系。利用幾何關系及數學原理求出該CT 機的各項參數： 旋轉中心位置、探測器單元間距及該儀器使用 X 射線的 180 個方向，并對數據進行優(yōu)化， 對該 CT 機實現參數標定，標定完成后，使用多層面重建2方法，對附件中給出的參數進行二維重建。 2.2.1 附件 2 數據處理 由于附件 2 中數據量龐大，首先需要對數據進行預處理。本文將采用 Excel 條件格式對數據進行處理，使用色階將數據層次化，通過圖形特征及顏色變化進行分析。 由題意可以推知附件 2 所給數據為 X 射線的吸收強度，即 X 射線發(fā)出與接收時的強度差。根據相關資料可知，X 射線的衰減規(guī)律滿

10、足衰減系數公式3： I = I0e L(2.1)可以看出衰減系數與穿過物體長度存在一定關系，即附件 2 中的數據為對物體長度進行一確定法則的增益得到的。所以在數據處理中首先應先分析出增益函數，找出附件 2 中數據與實際待測物體長度之間的關系，將附件 2 中的數據增益還原，得到長度數據， 再進行后面的計算。由于探測器間距未知，本文將在探測器間距求得之后對增益函數進行確定。 2.2.2 平行束濾波反投影重建模型建立 根據題目分析，本文擬建立平行束濾波反投影模型。模型將由 Radon 變換與 R-L 濾波器兩部分構成。根據問題一我們可以推測探測器旋轉中心與正方形托盤的位置并不重合， 勢必會造成一定的

11、誤差，所以我們將尋求有效的方法，對正方形托盤進行處理，使得處理后的虛擬大托盤能夠滿足旋轉中心與托盤幾何中心重合，從而得到較為準確的重建圖形。 2.2.3 減噪處理-對模型進行優(yōu)化 通過題目分析，問題二與問題三的題目問法基本相同，我們由此推測問題三可能存在由于非測試系統引起的數據波動導致重建圖形較為不準確。若根據 2.2.2 中平行束濾波反投影重建模型得到的重建圖形較為不準確，我們猜想數據的波動可能由于噪聲干擾產生，由此我們可以對模型進行優(yōu)化，在模型中加入減噪處理函數，使得重建圖像更加精確 2.2.4 模型的精度分析與評價 第四問首先要求對第一問中所建立的模型進行評價，即對模型一得到的CT 參數

12、標定結果進行精度與穩(wěn)定性的分析，即對于 180 個方向的確定、探測器間距、旋轉中心坐標的精度進行分析。此外，在對模型一進行分析之后，需要提供一種新型模板，從而使得精度與穩(wěn)定性大大提高。 三、模型假設1、題目 CT 機正常運行，質量完好，題中所給數據不存在因機器故障而造成的錯誤； 2、旋轉中心出現偏差的來源在于安裝誤差，而系統本身不存在偏差，又因系統為對稱系統，故旋轉中心位于中垂線上某點； 3、X 射線僅有長度，不存在寬度； 4、X 射線強度只在穿透被測物體時發(fā)生衰減，空氣衰減系數為零； 5、X 射線在傳播到接收過程中發(fā)生的干涉與衍射忽略不計； 6、發(fā)射與接收裝置一一對應。 四、符號系統序號 符

13、號 符號說明 1X 射線吸收率 2d探測器間距 3yij = kixij + mij以橢圓中心為直角坐標系時，編號 j 的探測器在第 i 個方向接收的 X 射線的直線方程 4y0旋轉中心縱坐標 5x0旋轉中心橫坐標 6Rij以橢圓中心為直角坐標系時，編號 j 的探測 器在第 i 個方向接收的 X 射線所截物體的實際長度 五、模型建立5.1 問題一的模型建立與求解 5.1.1 離散模型的建立 對題目分析可知，附件 1、2 均給出了大量的離散數據，根據數據 2，數據由于存在180 個不同方向而被分為 180 組，每組有 512 個數據，這些數據與射線穿過物體的長度有關，由此可以建立 180*512

14、 個方程。且在此題的求解中確認 180 個方向為確認所有參數的重要基礎，由此我們可以建立離散模型，找到 180 組方向間的關系從而進行方程的聯立與求解，求出方向參數。 首先選取橢圓中心點作為坐標原點，將橢圓短軸方向為 x 軸，長軸方向為 y 軸建立 xoy 平面坐標系。 yij = kixij + mijx2y2(i = 1h2hh180j=1,2, ,512)（5.1）152 + 402 = 1根據方程聯立可以將方程式的解 xij1，yij1 xij2，yij2 的分別使用相應的mij表示，并且我們可以得到由該光線截橢圓所得到弦長的表達式為： (xij1 xij2)2 (yij1 yij2)

15、2dij1 =xx22ij3 ij4yij3ij4 y同理可以求得該光線對截圓所得弦長的表達式： dij2 =(5.2)(5.3)其中 xij3，yij3 ， xij4，yij4 為直線與圓的兩個交點。 由附件二中信息可以求出射線截物體的長度，由此聯立得到方程組： Rij = dij1 + dij2(i = 1h2hh180j=1,2, ,512)(5.4)理論上，由這 180 組方程組即可求出所有ki i = 1h2hh180 ，此時還有三個未知量帶求，即旋轉中心坐標 x0，y0 和探測器單元間距 d。 5.1.2 模型求解 5.1.2.1 數據處理在 5.1.1 的分析中我們得到了求解問題

16、的理論模型，理論上通過對式（5.6）進行方程求解即可得到 180 個 X 射線方向，并且將這 180 個角度作為已知條件與附件 1、2 一 起得出設備旋轉中心及探測器單元間距。但由于本題的特殊性，存在龐大的數據量及未知量，實際上無法直接通過方程的求解得到題目所求信息，從而直接對該CT 系統進行參數標定。所以本文將對離散模型進行簡化，通過幾何分析及數據處理，尋求出能夠實際求解的有效簡化方法，并給出相應精確的結果。 首先，我們在 Excel 中使用條件格式對數據進行了整體分析，附件 2 中非零的數據進行標記，將數據中為零的數據與非零數據進行區(qū)分，縮放后得到了如下圖形： 5.1.2.2 計算探測器間

17、距 d圖 5. 1 附件 2 非零數據標記后縮略圖結合數據讀圖 5.1，可以明顯的看出在 1-22 列以及 108-180 列時，能清晰的看到兩條色帶，即在 1-22 組和 108-180 組所對應的探測器照射角度值下，不存在任何一條 X 射線同時穿過圓與橢圓。同時，根據附件 1 及圖 1.2 所給的模板尺寸可知寬帶處為 X 射線穿過橢圓，窄帶為 X 射線穿過小圓。此外，寬帶和窄帶之間為零的區(qū)域表示 X 射線通過橢圓與小圓間的空白區(qū)域。由于模型尺寸已知，根據穿過條數即可確定出探測器間距范圍?？梢韵旅鎸男A、橢圓、寬帶和窄帶間空白區(qū)域三方面共同對探測器間距 d 進行確定，從而縮小范圍，使得得到

18、的結果更為精確。 1） 對小圓進行分析 圖 5. 2 窄帶對應小圓的掃描數據圖5. 3 窄帶部分放大圖不妨將小圓部分取出單獨分析（如圖 5.2 及 5.3），每一列黃色部分的數據不為零即表示該方向中某些射線穿過了小圓。通過編寫程序對數據分析可知，穿過小圓的 X 光條數范圍為 28-29 條，28 根時有可能在小圓兩邊存在兩條相切的光線，故可由以下兩種穿透情況確定探測器間距的上下界： 圖 5. 4 28 根射線穿過小圓且邊界相切圖5. 5 29 根射線穿過小圓由圖 5.4 和 5.5 可以得到探測器間距與小圓直徑 的關系如下： 27d829(5.5)化簡可得探測器間距 的范圍為： (0.2759

19、h0.2963(5.6)2） 對橢圓進行分析 通過對圖 5.1 分析可以直觀的看出圖中右側中部出現明顯的“頸縮”段。假設頸縮最短處 X 射線與橢圓長軸基本平行，對此情況做如小圓的類似分析，利用處理后的附件2 的數據得出 X 射線與長軸平行時穿過橢圓的射線條數與橢圓寬度的關系（由于圖形過大且上下對稱，此處僅做出橢圓的上半部分）： 圖 5. 6 X 射線穿過橢圓示意圖分析同小圓類似，可以得到探測器間距與橢圓短軸的關系如下： 107d30109(5.7)化簡可得探測器間距的范圍為： 0.2752h0.2803(5.8)3） 對大圓與小圓間空白距離分析 同樣分析“頸縮”段，寬帶與窄帶之間有一段全部為零

20、的數據，與模板中橢圓與小圓之間的空白距離相對應。由圖 1.2 可知橢圓與小圓中心線與橢圓短軸在同一直線上， 故當 X 射線與橢圓長軸平行時，也與橢圓與小圓中心線連線垂直。對此區(qū)域的數據與模板尺寸做如上相同處理，可以得到探測器間距與空白區(qū)域的關系如下： 93d2695(5.9)化簡可得探測器間距的范圍為： 0.2737，0.2796(5.10)綜上所述，求出式(5.6)、(5.8)、(5.10)的交集即為探測器間距的精確范圍： d 0.2759，0.2796(5.11)為了計算方便，本文取 d = 0.2776 進行后續(xù)計算。 5.1.2.3 增益函數的確定由所閱讀的文獻，CT 機通常對衰減系數

21、公式做如下處理3： P = ln 0 = L(5.12)對此我們猜想，附件 2 所給信息與實際長度成正比，即： Rij =t(5.13)其中Rij表示實際長度， t表示第 i 行 j 列的數據。 為了驗證我們的猜想，我們需要挑選出橢圓與圓分離的數據，為了選定僅通過橢圓而不通過小圓的射線，我們對數據進行了進一步的處理。在 Excel 中對附件 2 的數據進行了條件格式的限定，采用色階對數據進行整合，得到了如下圖像： 圖 5. 7 使用色階處理后的附件 2 縮略圖選取最窄的一列數據，即圖中畫黑線的一列，認為此時射線束剛好平行于橢圓長軸射入，找到這一列最大的數據，用該值除以橢圓長軸，得到比例系數 1

22、，再用上方屬于圓的數據中的最大值，再除以圓的直徑，得到比例系數 2，兩個比例系數非常接近，可近似認為猜想成立。 下面將對該比例系數進行精確求解，單取出圓的數據進行分析。由微積分的思想： S =t.(5.14)t其中t表示圓的數據所在列數 由于間距 d 非常小，令 = ，且Rij = t，上式即變?yōu)椋?S = t.(5.15)t對圓的多列數據進行求值，并取平均，即可得到較為精確的比例系數 的值。通過計算，我們求得增益系數的值 = 1.776974。 5.1.2.4 180 個角度的確定雖然理論上能夠對 180*512 的方程進行求解，但由于數據量過于龐大，實際電腦無法進行運算，所以我們對模型進行

23、了改進，提出了簡化離散模型：對于每一種角度，隨機選取 3 根 X 光線為一組數據作為研究對象，三根 X 射線滿足以下條件：1）跨度較大； 2） 等距分布；3）射線僅通過橢圓而不通過小圓。 為了選定僅通過橢圓而不通過小圓的射線，我們對數據進行了進一步的處理。在Excel 中對附件 2 的數據進行了條件格式的限定，采用色階對數據進行整合，得到了如下圖像： 圖 5. 8 使用色階處理后的附件 2 縮略圖如圖，利用色階可以明顯的看到色帶重合部分內存在一條顏色較深的色帶，此顏色所對應數據即為 X 射線既通過橢圓又通過小圓后吸收量數據，此部分數據不能夠作為數據組中的數據進行運算。 對于每一種角度所選數據組

24、進行如下運算： 假設： 其中： yi1 = kixi1 + mi1 yi2 = kixi2 + mi2 yi3 = kixi3 + mi3i = 1h2hh180(5.16)mi1 + mi3 = m(5.17)2i2將式(5.16)中的三式分別與橢圓方程聯立： yij = kixij + mijx2y2j = 1h2h3(5.18) 152 + 402 = 1根據方程聯立可以將方程式的解 xij1，yij1 (xij2，yij2)分別使用相應的mij表示，并且我們可以得到由該組光線截橢圓所得到弦長的表達式為： xx22ij1 ij2yij1ij2 ydij =(xij1 xij2)2 (yi

25、j1 yij2)2弦長可根據附件二中信息求出，記為Rij，由此聯立式(5.15)得到方程組： Rij = dij =同時聯系式(5.17)，可以得到 180 組方程式： (5.19)(5.20)Rij= dij =xij1 xij22yij1 yij22mi1 + mi3 = m 2i2(5.21)i = 1h2hh180j = 1h2h3對于每一個獨立方程組進行獨立求解，可以求得本組數據所對應的光線斜率k，進行反三角運算可得到該組光線的角度值，并且對于每一種斜率取多組數據組進行上述運算，將所得結果求平均值以減小隨機帶來的誤差。并將得到的 180 個數據進行擬合，擬合后的 180 個角度值見下

26、表（與 x 正方向夾角）： 表 5. 1：180 個角度值-60.2152-58.8587-58.3008-57.2081-56.1751-55.2026-54.1989-53.1945-52.1913-51.1883-50.1841-49.1784-48.1748-47.1711-46.1670-45.0092-44.1559-43.1498-42.1464-41.1399-40.1346-39.1288-38.1194-37.1142-36.1026-35.0940-34.0861-33.0800-32.0708-31.0575-30.0481-29.1362-28.0209-27.010

27、4-25.9942-24.9819-23.9634-22.9444-21.9239-20.9011-19.8815-18.8563-17.8288-16.799-15.7555-14.7159-13.6685-12.6137-11.5409-10.4695-9.3730-8.2484-7.1066-5.8943-4.5786-3.0205-2.5839-2.1474-1.7108-1.2742-0.8377-0.40110.03540.47200.90851.34513.57474.88146.24087.45258.56439.7321110.797411.844012.917113.989

28、015.019616.035217.107818.114419.152520.191821.218122.224823.243924.261325.150226.295127.112528.324429.330030.340531.354532.364433.373934.382735.391236.399237.406738.413939.420840.427241.533942.439243.444944.450145.455346.460147.465048.469249.473450.477651.481352.485053.488654.491955.494956.498057.50

29、0758.503459.505960.508161.510262.512063.513864.515265.516666.517667.518368.518869.518870.518671.517972.516873.515174.512875.509776.505677.500478.493879.485380.474381.460182.441583.416284.381285.330186.250787.116287.858388.306191.779192.339193.128694.014494.945395.899896.868097.845098.827899.8147100.

30、8045101.7966102.7904103.7858104.7818105.7789106.7768107.7753108.7744109.7738110.7737111.7741112.7745113.7751114.7763115.7775116.7791117.7812118.77255.1.2.4 旋轉中心的確定 由題意我們可以知道，旋轉中心出現偏差的來源在于安裝誤差，而系統本身不存在偏差，所以本文將從此處入手求解旋轉中心位于方形托盤的位置。首先我們對 512 個探測器標號 1-512 號。由于儀器本身不存在偏差，為使得儀器能夠準確的掃描得出結果， 該 CT 肯定圍繞中垂線上某點旋

31、轉。所以將旋轉中心投影至上時對應點的應為編號為 256 或 257 的。橢圓的中心點即為正方形托盤的中心，所以接下來對 橢圓的掃描情況進行分析。由圖 5.7 可以看到橢圓帶有明顯的最寬處和最窄處。根據幾何分析可知，最寬處對應的掃描方式為：X 射線幾乎平行于橢圓短軸且從正方形托盤左方進行入射（如圖）；最窄處對應的掃描方式為：X 射線幾乎平行于橢圓長軸且從正方形托盤下方進行入射（如圖）。 512 1 圖 5. 9 研究對象選取位置5121 圖 5. 10 X 射線平行于短軸入射圖5. 11X 射線平行于長軸入射1） 旋轉中心 的確定 當 X 射線以幾乎平行于橢圓短軸掃描時，根據附件 2 的數據我們

32、可以根據最大值找到橢圓中心點投影至上對應的編號。如果在安裝過程中沒有出現誤差，理論上該投影應與旋轉中心投影重合。然由于存在安裝誤差，導致了橢圓中心點投影與接 收器中點不重合，而兩點之間的距離則為旋轉中心y0的值，即橢圓中心點投影所對編號與中心位置編號之差與探測器間距 d 的乘積。 2） 旋轉中心 的確定 當 X 射線以幾乎平行于橢圓長軸掃描時，根據附件 2 的數據我們可以根據最大值找到橢圓中心點投影至上對應的編號。如果在安裝過程中沒有出現誤差，理與中心位置編號只差與探測器間距 d 的乘積。 具體計算由 MATLAB 進行實現，程序見附錄。利用 MATLAB 求得旋轉中心的坐標為： 以.圓以以中

33、h . 圓綜上所述，我們對此 CT 完成了參數標定。探測器間距為 0.2776mm、旋轉中心坐標為（-9.2996,5.5520）。180 個角度值在表 5.1 中進行了全部的羅列，由于數據量過大再此不再重復敘述。 5.2 問題二的建模與求解 5.2.1 平行束濾波反投影重建模型的建立 本題為二維平行波反投影重建問題4，在問題一中已將該 CT 各項參數完成標定。為了得出準確的圖像，成功實現二維平行波反投影重建，應建立平行束反投影重建模型。傳統模型為反投影重建模型，但其存在嚴重的缺點是會引入星狀偽跡，對此需要引入濾波函數，建立平行束濾波反投影模型。建模過程分為兩部分：1)對數據進行基于吸收率的R

34、adon 變換5；2)使用 R-L 濾波器4.進行濾波處理。 5.2.1.1 Radon 變換及其性質 假設 f(xhy)為待重建物體的密度函數，其 Radon 變換的定義5為沿一組平行 X 射線的第一類曲線積分： R hr =L hr f xhy ds(5.22)其中 R hr 是函數 f(xhy)的 Radon 變換。每一根射線 M 由和 r 兩個參數決定，其中是射線 M 的垂線和 x 軸的夾角，r 是射線 L 到原點的距離（如圖 5.12）。 圖 5. 12 Radon 變換基本參數示意圖對于本題，我們在第一問中已經將全部的給出。對于每一個給出二維成像物體 f(xhy)一維投影的全部集合

35、。利用變量 r = xcos + ysin和沖激函數的抽樣特性可得： R hr =(L2)f xhy (xcos + ysin y)dxdy(5.23)以上為對 Radon 變換定義的敘述。 在本題中題目要求對被掃描物體進行 1:1 重現，即不需要考慮對函數 f xhy 進行各項同性的縮放、旋轉或平移變換。 假設p(xr)為 f xhy 在角度 0時的平行束投影。理論上投影p(xr)和密度分布函數 f xhy 在時域上可以使用一維線積分的堆積來表述，但根據線積分尋找投影重建的方法 是非常困難的。而傅里葉切片定理5在頻域上提供了投影與圖像之間更簡單的數學關系。其數學表達式為： F1 p xr=

36、F h |= 0(5.24)投影圖像重建的問題，原則上可按如程求解：采集不同視角下的投影，求出各投影的一維傅里葉變換，由傅里葉定理即可得到圖像的二維傅里葉變換的各切片，然后匯集成圖像的 2D 傅里葉變換，最后求反變換得到重建圖像。 5.2.1.2 R-L 濾波器 基于上述原則，我們采用了 R-L 濾波器4。R-L 濾波器的離散形式由的學者 2BRamachandran 和 Lakshminarayanan 提出，其頻域函數為： 式中： HR L = W = rect12B1h = (5.25)rect=2d0h其他 (5.26)其中，是空間頻率，W 是窗函數。相應的時域卷積函數hR L xr

37、為： hR L xr = 2B2 sin c 2xhBB2sinc2 xhB(5.27)以xr = nd 帶入式（5.27）得到hR L xr 的離散形式： 14d2 0 1n22d2= 0hR L xr = 偶數 = 奇數 (5.28)R-L 濾波器形式簡單實用，重建效果較好，輪廓清楚。 5.2.1.3 平行束濾波反投影重建模型 設需要重建的圖像為 b(x,y),它的二維傅里葉變換為 B1h 2 。根據傅里葉切片定理， h 可通過 b(x,y)在不同視角下的投影 ( )的一維傅里葉變換求得，即： B1h 2 =h= F1 p xr = th(5.29)需要重建的圖像11t( + )h= b

38、xhy =1h 2=1h 2121 224 2=t( h) t2th ()(5.30)0上式中后半部分積分，可寫成空域變量為 的傅里葉反變換式：( = th ( )t( h) t2th ()=t( h) t2= hph= gh= gth ()h(5.31)式中， gh = hph(5.32)而 h=1，ph =1 t( h ) 。把式(5.30)代入式(5.29)后得到： 22h=thh(5.33)0上式的物理意義是：經過給定點(rh)的所有濾波后的射線投影在 0范圍內的累加， 反投影重建，得出(rh)點的像素值，這就是濾波反投影方程，可以集中表現出濾波反投 影算法的各個步驟。 綜上所述，我們

39、將應用平行束濾波反投影重建模型，對第二問進行求解。 5.2.2 平行束濾波反投影重建模型的實現與問題求解 由 5.2.1 已經分析出算法原理，而軟件 MATLAB 程序中 iradon 函數則正是基于 R-L 濾波器的濾波反投影函數，所以此模型可直接由 iradon 函數進行實現。 5.2.2.1 比例系數處理 由于題目所給數據經過了增益處理，所以在使用計算機進行反投影時需要建立計算機所求得數據與經過增益處理后的數據之間的關系。對此我們先對附件 2 進行反投影計算，求得各個像素點的吸收率，對比附件 1 所給出的吸收率，建立二者之間的關系函數， 分析可知計算機所求得的吸收率與附件 1 所給的吸收

40、率近似成正比，利用 MATLAB 求得比例系數 k=2.033。即： 2 =1(5.34)其中 1為電腦計算的吸收率， 2為附件 1 所給的吸收率。 5.2.2.2 “鑲邊”處理 根據第一問可知，由于安裝誤差，旋轉中心并不在正方形托盤的幾何中心，導致了直接使用 iradon 函數會出現較大的偏差。對此我們將正方形托盤進行了“鑲邊”處理（如圖 5.13 其中“”為鑲邊前正方形托盤的幾何中心，“”為旋轉中心），使得鑲邊后的正方形的幾何中心即為旋轉中心（即“”為大正方形幾何中心）。 圖 5. 13 經過鑲邊處理后大正方形及各組分幾何關系圖在對正方形托盤進行“鑲邊”處理之后，得到的“新大正方形托盤”滿

41、足旋轉中心 位于托盤中心，故消除了系統誤差，此時可以直接使用 iradon 函數對其進行濾波反投影計算。 5.2.2.3 濾波反投影計算 利用 MATLAB 軟件編寫對應程序，利用 iradon 函數對附件 3 數據進行濾波反投影計算，計算結果見 problem2。但得到的結果為大正方形下的運算結果（如圖 5.14）。 圖 5. 14“鑲邊”后附件 3 反演重建圖要求得真實正方形托盤上被掃描物體的幾何形狀及其確定位置，還需要進行“去邊” 處理。去邊處理并且乘上比例系數后的運行結果見圖 5.15，各個點吸收率的計算結果見文件 problem2.xls 。 圖 5. 15 去掉“鑲邊”后的反演圖十

42、個確定點的參數如下表： 表 5.2 問題二中十個確定點吸收率 XY吸收率10.000018.00000.000034.500025.00000.972243.500033.00000.003645.000075.50001.176148.500055.50001.042650.000075.50001.465256.000076.50001.284965.500037.00000.000779.500018.00000.000098.500043.50000.01755.3 問題三的模型與求解 5.3.1 平行束濾波反投影重建模型運算 由題設，要求基本與問題二相似，根據問題二求解的思想，我們同

43、樣進行鑲邊等處理，然后使用 iradon 函數編寫程序，對問題二的數據實現了濾波反投影重建，得到的圖像如下： 圖 5. 16 未降噪時的反演圖像十個確定點的參數如下： 表 5.3 問題三中十個確定點未降噪吸收率 XY未降噪吸收率10.000018.00000.011434.500025.00002.253043.500033.00005.930645.000075.50000.016948.500055.50000.041050.000075.50003.134356.000076.50005.950765.500037.00000.000079.500018.00007.764698.500

44、043.50000.05515.3.2 平行束濾波反投影重建模型優(yōu)化 由圖 5.16 我們可以看出，利用 5.2 中的平行束濾波反投影重建算法對附件 5 中的進行計算，得到的圖像邊界存在大量毛刺、邊緣不光滑、邊界模糊等問題。經過分析我們 推斷由于噪點的存在，影響了圖像的質量，對掃描數據及處理存在干擾。所以我們著手對數據進行降噪處理。 5.3.2.1 吸收率的波動范圍 對于問題三，已知數據只有受噪點影響的數據附件 5，以及根據附件 5 得出的相關 數據。僅有這些數據，在沒有標準值進行參考的情況下，我們對噪點造成的結果波動無法 定量的進行確定。但對于模板試件，我們已知其幾何參數及掃描后的得到數據，

45、故本文對附件 2 進行反演，利用 iradon 函數計算出反演矩陣，將反演矩陣和初始矩陣（附件1） 進行比對，從而確定噪點是否存在，及其對圖像的影響程度。 我們使用 5.2 中的算法對附件 2 進行了反演，將得到的反演矩陣與附件 1 標準矩陣相減，把得到的差值矩陣中的每個元素求絕對值以后求和，并將其平均分攤到 256*256 個像素點中，得到平均每個像素點的干擾值為 0.0174，即由 iradon 得到的吸收率數值平均波動為0.0174。 5.3.2.2 降噪以后吸收率波動范圍 根據噪聲的概率分情況來看，可分為高斯噪聲、瑞利噪聲、伽馬噪聲、指數噪聲和 均勻噪聲。MATLAB 軟件中有多種不同

46、的降噪工具函數，但針對不同的噪聲，不同的降噪工具函數的降噪效果也有很大區(qū)別。在本題中，我們無法確定造成噪點的噪聲種類， 所以將使用三種處理不同類型噪聲的降噪函數對已知參數的模板掃描結果進行降噪處理， 分別求出降噪后數據的波動范圍，選取最小的波動范圍所對應的工具函數作為本題的降噪（濾波）工具。我們初步選取的三種降噪工具函數6如下： 1） 均值濾波函數6均值濾波是典型的線性濾波算法。其采用的主要方法為領域平均法，即對待處理的某個像素點(xhy)，選擇一個模板，該模板由其近鄰的若干像素組成，求模板中所有像素的均值，再把該均值賦予當前像素點(xhy)。作為處理后圖像在該點上的灰度 g xhy ，即 g

47、 xhy= 1M(mhn)Sf(xhy)，s 為模板,M 為該模板中包含當前像素在內的像素總個數。 2） 自適應濾波函數6 它能根據圖象的局部方差來調整濾波器的輸出,局部方差越大,濾波器的平滑作用越強。它的最終目標是使恢復圖像 f(xhy)與原始圖像 f(xhy)的均方誤差e2 = E f xhy f xhy 2 最小。 3） 中值濾波函數6它是一種非線性平滑濾波函數,其基本原理是把數字圖像或數字序列中一點的值用該點的某個領域中各點值的中值代換，其主要功能是讓周圍象素灰度值的差比較大的像素改取與周圍的像素值接近的值,從而可以消除孤立的噪聲點,所以中值濾波對于濾除圖像的椒鹽噪聲非常有效。 我們分

48、別使用三種濾波工具，對反演矩陣進行降噪處理，經過同 5.3.2.1 對模板進行數據分析處理后，得到三種工具降噪后噪點造成的數據波動范圍為： 表 5.4 降噪工具與其處理后波動值 降噪工具 無 Filter2Wiener2Medfilt2波動值 0.01740.01590.01520.0158根據上表對三種降噪工具進行對比，我們最終將維納濾波法用于問題三中數據的去燥處理。 5.3.2.3 降噪后問題三的解 經過Wiener 去噪處理后，我們得出了問題三的優(yōu)化解。計算結果見文件 problem3.xls。我們對降噪前后得到的吸收率數據表進行了色階處理，通過對數據左上角區(qū)域進行分析 對比可以看出，經

49、過降噪處理后的噪點數量大大減小，圖像的邊界更為明確，降噪效果較好： 圖 5. 17 降噪前的數據圖 5. 18 降噪后的數據經過降噪處理后的圖像為： 十個確定點的參數為： 圖 5. 19 降噪處理后的重建圖像表 5.5 十個確定點降噪后吸收率 XY降噪后吸收率10.000018.00000.012634.500025.00002.290243.500033.00005.915945.000075.50000.016348.500055.50000.082350.000075.50003.133656.000076.50006.033365.500037.00000.000079.500018.

50、00007.718498.500043.50000.08615.4 模型的評價與優(yōu)化 5.4.1 模型精度和穩(wěn)定性評價 5.4.1.1 探測器單元之間的距離歸一化均方差評價模型： 歸一化是一種無量綱處理手段，使物理系統數值的絕對值變成某種相對值關系。簡化計算，縮小量值的有效辦法。我們通過計算重建圖像與原圖像的歸一化均方差，即： t t 1 htt 2 ht22tt 2 htdIF = 1其中 t 1表示重建圖像的像素矩陣， t 2表示原圖像的像素矩陣。 (5.35)計算所得 dIF 越接近于 1，表明重建圖像與原圖像差別越小，表明所用參數精度越高；反之，表明重建圖像與原圖像差別越大，表明所用參

51、數精度越低，以此建立歸一化均方差評價模型。 對精度分析，將探測器單元之間的距離進行歸一化均方差評價，用我們所采用的間距（d=0.2776）進行計算，求得 dIF = 0.9730,該數值較接近 1，認為該參數精度較高。 對穩(wěn)定性分析，在我們得到的間距范圍 0.2759h0.2796 ,以步長為 0.0001 進行搜索，求得一系列 值，其變化較為平緩，認為模型的穩(wěn)定性較好。 5.4.1.2 方向 由于原離散模型數據量過于龐大，我們采用一種簡化的離散模型算法，即隨機選取 一部分數據作為計算對象，但由于數據量小，算得的角度和遍歷所有數據組算出來的角度 存在著偏差，所以，得到的 180 個方向，與真實

52、值有差異。為了減小這種差異，我們認為， 可以適當提高選取數據的組數，對其求平均值，并且確定一個誤差評價的規(guī)則， 當新的一組數據計算時，得到的角度與前面的平均值的角度相差在一個范圍內的時候， 我們認為該平均值可以作為這個方向的角度值。 5.4.1.3 旋轉中心位置 在計算旋轉中心的時候，我們認為 X 射線分別平行于短軸、長軸入射，但計算得到的方向不完全平行。因此得到的旋轉中心 位于實際旋轉中心的下側，旋轉中心 位于實際旋轉中心的下方，左側。 5.4.2 新模板的設計 由 5.1.2.2 計算探測器間距，主要思路是對 180 個角度的生成的 512*180 的矩陣進行一個圖形構成、角度對應的分析，用以限制探測器間的間距。對此，我們提出這樣一個模板 圖 5. 20 新模板設計方案圖 5. 21 新模板偏心掃描大致圖像其對應的 512 個偏心掃描圖像大致形狀如圖 5.21。使用新模板，我們能夠找到更多可供分析的特征角度，由于能偶看到明顯的四條帶分層情況，我們可以更精確地對探測器間距數值的范圍進行縮放，可以在一定程度上減小誤差。 參數標定模型與 5.1.2.2 基本相同，而由于特定角度下可供分析的數據量大大增加， 在對探測器間距數值的范圍進行縮放后，探測器間距本身的精度得到提高，而通過探測 器間距而計算求得的旋轉中心點的坐標以及 1