高教社杯全國大學生數(shù)學建模競賽題目 CT系統(tǒng)參數(shù)標定及成像

1、2018西交數(shù)模第一次模擬賽數(shù)學建模論文首頁選題 A 隊伍編號 696 班級學號姓名隊長錢學森73班2173714071孔令輝隊員1錢學森73班2171311098楊崢隊員2電氣642160400092李佳航 2018年 6月 30 日摘要本題針對一種二維CT獲取樣品內(nèi)部結構信息的工作方式及成像原理，意在通過借助已知結構樣品進行參數(shù)標定，消除系統(tǒng)誤差，而后對未知結構的樣品進行成像，得出該未知介質的相關信息。問題一中，要求通過標定模板的相關，確定CT系統(tǒng)的旋轉中心、探測器單元之間的距離以及該CT系統(tǒng)使用的X射線的180個方向。本文利用標定模板幾何參數(shù)，通過條數(shù)、以及探測器單元間距相等等信息，首先

2、計算出探測器單元之間的距離為0.2778mm。根據(jù)X射線與標定模板的幾何特性，以橢圓短軸所在直線為x軸，長軸所在直線為y軸建立直角坐標系，得旋轉中心在所建立的坐標系中的坐標為。最后通過近似確定X射線180個方向的旋轉角度具有高度線性相關性，根據(jù)部分確定數(shù)據(jù)擬合出整體旋轉角度。前十五組旋轉角度為：29.6493 31.0027 31.5572 32.6476 33.6795 34.6478 35.6483 36.6490 37.6475 38.6486 39.6488 40.6491 41.6487 42.6497 43.6492問題二中，要求通過已求得的標定參數(shù)，確定未知介質在正方形托盤中的位

3、置、幾何形狀和吸收率等信息，并具體給出所要求十個位置的吸收率。本文依據(jù)CT層析成像原理，利用逆拉東變換作出重構圖像，并利用Excel中數(shù)據(jù)分布計算出原位置介質的相關性質。所求十個位置的吸收率為：序號12345吸收率0.01850.9701-0.00011.17031.0287序號678910吸收率1.44251.2815-0.0074-0.00250.0194問題三中，與問題二計算方法相似得十個位置的吸收率為：序號12345吸收率0.07362.53057.13270.03320.8015序號678910吸收率3.14875.57140.05308.13620.0272問題四中，要求對問題一中

4、的標定模型進行改進以減小誤差并增加穩(wěn)定性，本文利用在問題一求解過程中遇到的問題進行思考，首先適當增大模板減小偶然誤差，其次做出投射圖像后應容易找到極值，并且圖像應有一定的對稱性；圖像掃描后應盡量少地得出重復數(shù)據(jù)。關鍵詞：CT層析成像 Radon變換與逆變換 吸收率 MATLAB算法1.問題重述1.1問題背景應用ng yong3構造ient in parabolic problemsJ. 33333333333333333333333333333333333333333333333333333333333333333CT(Computed Tomography)可以在不破壞樣品的情況下，利用樣品

5、對射線能量的吸收特性對生物組織和工程材料的樣品進行斷層成像，由此獲取樣品內(nèi)部的結構信息。1.2相關信息本題介紹了一種二維CT系統(tǒng)，平行入射的X射線垂直于探測器平面，每個探測器單元看成一個接收點，且等距排列。X射線的發(fā)射器和探測器相對位置固定不變，整個發(fā)射-接收系統(tǒng)繞某固定的旋轉中心逆時針旋轉180次。對每一個X射線方向，在具有512個等距單元的探測器上測量經(jīng)位置固定不動的二維待檢測介質吸收衰減后的射線能量，并經(jīng)過增益等處理后得到180組接收信息。并且，為消除安裝誤差，需要對安裝好的CT系統(tǒng)進行參數(shù)標定，即借助于已知結構的樣品（稱為模板）標定CT系統(tǒng)的參數(shù)，并據(jù)此對未知結構的樣品進行成像。題目附

6、件中提供了標定模板的幾何信息，接受信息，待測介質的接收信息以及圖3所給位置的相應數(shù)據(jù)。1.3需要解決的問題問題一在正方形托盤上放置兩個均勻固體介質組成的標定模板，模板的幾何信息如圖2所示，請根據(jù)這一模板及其接收信息，確定CT系統(tǒng)旋轉中心在正方形托盤中的位置、探測器單元之間的距離以及該CT系統(tǒng)使用的X射線的180個方向。問題二附件3是利用上述CT系統(tǒng)得到的某未知介質的接收信息。利用第一問中得到的標定參數(shù)，確定該未知介質在正方形托盤中的位置、幾何形狀和吸收率等信息。另外，請具體給出圖3所給的10個位置處的吸收率。問題三利用給出上述CT系統(tǒng)得到的另一個未知介質的接收信息。利用第一問中得到的標定參數(shù)，

7、給出該未知介質的相關信息。另外，請具體給出圖3所給的10個位置處的吸收率。問題四分析第一問中參數(shù)標定的精度和穩(wěn)定性。在此基礎上自行設計新模板、建立對應的標定模型，以改進標定精度和穩(wěn)定性，并說明理由。 圖1. CT系統(tǒng)示意圖 圖2. 模板示意圖（單位：mm） 圖3. 10個位置示意圖2.問題分析2.1 問題一分析結合圖2和附件1表中數(shù)據(jù)，可以首先計算出CT系統(tǒng)探測器個數(shù)和模板長度度量的比值，運用程序1.1可以得出模板的幾何形狀如圖，可大致認為它是對稱的。對于附件2 ，由于對180個方向尚無清晰地認識，首先用同樣的方式做出數(shù)據(jù)分布圖，觀察到圖像比較平滑，因此認為按表格的順序180個方向是相鄰較密、

8、不錯位的。（圖中有色區(qū)域表示該點有吸收率，藍色部分表示吸收率大于100） 圖 2.1.1 附件1的數(shù)據(jù)分布圖 圖2.1.2 附件2的數(shù)據(jù)分布圖在圖中易觀查到有一紅色條形圖案，這是在不同的方向掃描到圓形時留下的，可以此為突破口首先求出探測器單元之間的距離。求出探測器單元之間的距離與增益比率之后，可以根據(jù)幾何關系，自己設立坐標系并通過數(shù)學運算計算出旋轉中心。最后在建立的坐標系內(nèi)，將待求解的180個旋轉方向轉換成X射線的斜率進行數(shù)學運算。2.2問題二分析題目給出了未知介質的接收信息，要求出介質的相關信息，可以搜索相應的數(shù)學模型，將附件中給出的按照射線條數(shù)與旋轉角度列成的表格，一一對應為相應坐標點的吸

9、收情況，從而根據(jù)各坐標點的不同性質，還原回該未知介質的幾何信息與吸收率等信息。2.3問題三分析問題三與問題二類似，可以大致看出數(shù)據(jù)分布更具有一般性，不容易描述出未知介質的相關信息，可以通過圖形大致描繪出介質相關信息。2.4問題四分析問題四要求分析題目所給的二維CT系統(tǒng)，設計新的標定模板以提高原系統(tǒng)的精確度與穩(wěn)定性。可以搜索相關資料，根據(jù)第一問的求解思路與求解過程，以規(guī)避求解過程中因標定模板自身性質而出現(xiàn)的誤差為原則進行思路拓展，以設計高精度與穩(wěn)定性的標定模板。3.模型假設u 假設在射線經(jīng)過介質時能量只損失在介質中，及不考慮衍射等現(xiàn)象；u 假設附件中所給出數(shù)據(jù)是正確的、可以直接利用的；u 假設旋

10、轉中心在相鄰兩條射線的中間直線上；4.模型的建立與求解4.1 問題一模型的建立與求解建立模型從題目中可以得出，由于x射線之間得間距相等，不管x射線怎么旋轉，穿過托盤上圓的x射線條數(shù)應該是大致相等的?？蓪⒋┻^圓形標定模板的X射線與模板建立模型示意圖如下：圖4.1.1.1X射線穿過圓形標定模板示意圖如圖4.1.1.1所示可以設穿過圓的在圓心兩旁的射線被圓所截的距離分別為x1,x2,x3和y1,y2,y3。求解旋轉中心時，以橢圓形標定模板短軸所在直線為X軸，長軸所在直線為Y軸建立直角坐標系，將托盤進行劃分。由附件二的數(shù)據(jù)分布圖的不對稱性可知，旋轉中心應在對稱軸某一側，大致確定旋轉中心方位后，根據(jù)旋轉

11、中心兩側探測器個數(shù)不變且同一探測器接收的與距旋轉中心的距離不變具體確定旋轉中心的位置。求解X射線旋轉角度時，設穿過橢圓的最邊緣的射線到橢圓中心的距離為R，取橢圓中心為原點，模板的對稱中心為x軸建立平面直角坐標系，X射線所在直線的斜率為k。d相鄰探測器之間的距離N對應于圓形模板的射線條數(shù)或探測器個數(shù)標定模板的吸收率處理數(shù)據(jù)時的增益率圓形模板的直徑4.1.2 符號說明4.1.3模型求解求探測器單元間距離：通過MATLAB編程求出計算出的條形帶平均涉及探測器個數(shù)為，沿直徑方向上的平均吸收率；由圖二圓形模板直徑，計算出探測器之間的平均距離：探測系統(tǒng)平均增益率：求旋轉中心：分析圖2.1.2可知，由于藍色

12、區(qū)域僅分布在后部分角度范圍內(nèi)，因此估計旋轉中心的位置在對稱軸的某一側，紅色區(qū)域分為兩部分時表示該角度下由兩部分射線分別照射經(jīng)過橢圓形和圓形；為了討論方便，現(xiàn)對正方形托盤做出如下劃分：圖4.1.3.1 圓盤劃分簡單分析可知，如果旋轉點在I區(qū)域，則沿180個方向照射后不會出現(xiàn)射線分成兩部分的情況；在II或III或IV區(qū)域，當吸收率出現(xiàn)最大值的時候射線也被分成兩部分，而不是像圖1.2那樣成為一部分，也排除；綜合各因素可判斷旋轉中心應該在V區(qū)域。因為發(fā)射-接收系統(tǒng)逆時針旋轉，且圖中兩部分紅色區(qū)域，間距縮小，說明旋轉之后在垂直于發(fā)射-接收方向上二者的距離是縮短的，從而確定旋轉中心在短軸的上半側，先運用E

13、xcel對附件二第1列數(shù)據(jù)作圖（即畫出第一個方向上的掃描圖像）如圖4.1.3.2，發(fā)現(xiàn)此時得到兩部分圖像，說明兩模板之間有一部分射線直接被探測器接收，直到第14列數(shù)據(jù)兩部分圖像結合在了一起，如圖4.1.3.3。 圖4.1.3.2第一列數(shù)據(jù)分布 圖4.1.3.3第14列數(shù)據(jù)分布結合托盤的幾何特征，在垂直于對稱軸方向上應該會出現(xiàn)最大的吸收率，利用MATLAB求得出現(xiàn)最大吸收率的方向為第151個方向，在此方向結合增益率得出的模板長度為進一步驗證了結果；平行于對稱軸方向上最大的吸收率出現(xiàn)在沿對稱軸的直線上，計算得出為第61個方向，此方向的模板長度，也驗證了結果。根據(jù)極近似水平方向為第61方向、最大吸收

14、率出現(xiàn)在235號探測器，極近似豎直方向為151方向、最大吸收率出現(xiàn)在223號探測器，設旋轉中心到豎直軸、水平軸的距離分別為x,y。為求解還需要另一個方向的等量關系，選取橢圓和圓的一條外公切線的方向，經(jīng)計算得出經(jīng)橢圓與圓公切線所在直線的X射線為第372號射線。則在所建立坐標系內(nèi)，切線方程為：進而列出二元方程組：解之得:因此旋轉中心在所建立的坐標系中的坐標為。求180個旋轉方向：由于x射線的發(fā)生裝置是連續(xù)旋轉的，所以在512條射線中穿過模板的射線條數(shù)應該是連續(xù)變化的，所以用MATLAB編寫程序，計算每次旋穿過模板的x射線條數(shù)，并畫出圖像如圖4.1.3.4。圖4.1.3.4 穿過模板射線條數(shù)隨旋轉次

15、數(shù)變化曲線圖中橫坐標為旋轉次數(shù)，縱坐標為穿過模板的射線條數(shù)。從圖中可以看出附表2中的數(shù)據(jù)是按照射線發(fā)射裝置旋轉的順序依次給出的，而且可以看出，180次旋轉后裝置共旋轉了180度，每次旋轉的角度近乎相等?？紤]到在前50次旋轉中，穿過整個裝置的x射線條數(shù)與穿過橢圓的條數(shù)相等。設穿過橢圓的最邊緣的射線到橢圓中心的距離為R，取橢圓中心為原點，模板的對稱中心為x軸建立平面直角坐標系，由于射線可以近似看成與橢圓和圓x2+y2=R2都相切，可以得到以下方程組：x2+y2=R2y=kx+mx2152+y2402=1y=kx+m化簡以上方程組可得：k=402-R2R2-152=402-p2s2p2s2-152利

16、用MATLAB編程得到圖像如圖4.1.3.5所示圖4.1.3.5 前50次旋轉角度變化曲線用同樣方式得出最后面25組數(shù)據(jù)的圖像為：圖4.1.3.6后25次旋轉角度變化曲線從圖像中可以看出每次旋轉角度近似為1度。為了更加精確地計算旋轉角度以便據(jù)此得到估算其他角度地依據(jù)，本文設計了另外一種算法。從附表二中可以看出，前15組數(shù)據(jù)中穿過橢圓的射線與穿過圓的射線沒有交叉，所以數(shù)據(jù)中換算出來的最大吸收距離就是近似穿過橢圓中心的射線被橢圓所截的距離，弦長公式為：x2152+y2402=1y=kx+m且這個距離由直線的斜率k唯一確定，已知k、m時可解得：l=1200k2+1152k2+402152k2-m2+

17、402所以用MATLAB編程計算得出了較為精確的前15次旋轉角度的結果如下：29.6493 31.0027 31.5572 32.6476 33.6795 34.6478 35.6483 36.6490 37.6475 38.6486 39.6488 40.6491 41.6487 42.6497 43.6492可以看出，這些結果的線性相關性非常好，據(jù)此利用Excel進行擬合，得到如下圖表：圖4.1.3.7 旋轉角度與旋轉次數(shù)擬合公式從圖表中可以得出擬合公式為：y=0.9938+28.718,相關系數(shù)為：0.9996根據(jù)擬合公式利用MATLAB編寫程序計算得出所有方向角度。4.2 問題二模型的

18、建立與求解建立模型通過分析問題與附件數(shù)據(jù)可以發(fā)現(xiàn)，題目所給數(shù)據(jù)與介質相關性質的二維分布具有對應關系，X射線將介質一條線上的性質投影為一點。據(jù)此，本文利用radon變換與radon逆變換進行運算，以通過投影后的訊號重建原始未知介質相關性質的二維分布。該變換的定義為：令密度函數(shù)f(X)=f(x,y)是一個的定義域為的緊致臺(compact support)。令R為radon變換的運算子，則Rf(x,y)是一個定義在空間中的直線L。其基本思想為：radon變換可以理解為圖像在空間的投影,空間的每一點對應一條直線，而radon變換是圖像像素點在每一條直線上的積分。因此，圖像中高灰度值的直線會在空間形成

19、亮點，而低灰度值的線段在空間形成暗點。對直線的檢測轉化為在變換區(qū)域對亮點、暗點的檢測。Radon變換是一幅圖像在一個特定的角度下的徑向線方向的投影，一幅圖像的radon變換是每一個像素radon變換的集合。對于MATLAB中語句R = radon(I, theta)，如果theta是一個標量，R則是一個包含在theta的列向量。如果theta是一個向量，R則是一個矩陣，矩陣的每一列是對應其中一個theta的radon變換。而radon變換的逆運算，就可以將CT系統(tǒng)對于每一直線上的X射線吸收率數(shù)據(jù)還原回未知介質的物理性質。radon反變換的公式是：該反變換操作比較簡單， 思路清晰，可以借助數(shù)學運

20、算軟件計算。4.2.2模型求解在MATLAB上將附件2和附件3的數(shù)據(jù)導入，利用radon逆變換將附件中數(shù)據(jù)重建未知介質信息，得到圖像如圖4.2.2.1。圖4.2.2.1 radon逆變換重建介質圖像（1）為使逆拉東變換后得到的圖像與原圖像大小相等，在逆radon變換公式中取362條射線，每兩條射線相距0.2770，所以得到的圖形長寬也為100，由第一問求解得知，X射線從約29度位置開始旋轉，為消除重建模型與真實模型的旋轉角度差異，將從29度逆radon變換得到的數(shù)據(jù)導入Excel表格可得變換后圖像坐標水平方向平移了31×0.2770=8.587，豎直方向平移了22×0.27

21、70=6.094，所以將原圖中的坐標按上述數(shù)值平移就得到了變換后的坐標。從模板的逆radon變換產(chǎn)生的矩陣中可以發(fā)現(xiàn)，模板中所有點對應的灰度都近似為0.5，又由于模板的吸收率為1，所以相對比例近似為2。據(jù)此利用MATLAB編程可以算出圖中對應十個點的吸收率如下表所示：表4.2.2.1 所求十點的吸收率（1）序號12345吸收率0.01850.9701-0.00011.17031.0287序號678910吸收率1.44251.2815-0.0074-0.00250.0194然后根據(jù)逆radon變換的結果，將得的數(shù)據(jù)導入Excel表格，在不同范圍內(nèi)的數(shù)值填充成不同顏色對不同吸收率的部分進行色塊填充

22、，得到結果如圖4.2.2.2所示。再利用Excel表格中尋找各橢圓定點的坐標，由于拉東變換中x射線之間的距離都是0.2770，所以表格中兩組數(shù)據(jù)在實際物體上的距離也是0.2770。由100÷0.2770362得知表格取362組數(shù)據(jù)時，總長度與原圖基本相同，此時可以長度為基礎，計算坐標變換公式。根據(jù)前文坐標變換的逆變換，用MATLAB編寫程序運算各色塊（由運算可得各色塊均為橢圓形）的坐標及長短軸相應數(shù)據(jù)。將問題二中數(shù)據(jù)進行radon逆變換后的圖像最低點在Excel行數(shù)和列數(shù)，將行數(shù)和列數(shù)乘以倍率0.2770即距離圖像邊界的距離，由于第一個圖距離兩邊界的距離已知，可以得到平移的方向和距離

23、，具體結果如表4.2.2.2所示。圖4.2.2.2 不同區(qū)域吸收率關系圖（1）以每個橢圓中心的吸收率代表整個橢圓的吸收率則A，B，C，D，E，F(xiàn)的吸收率分別為：0、1.1870、1.2914、0、0.9877、1.0632。表4.2.2.2 不同吸收率介質在托盤中的位置、幾何形狀與吸收率橢圓編號ABCDEFx063.017546.259054.707540.580551.245048.3365y033.101575.067071.189030.331553.184055.5385half111.19205.25577.08529.645422.68871.8478half26.11913.46

24、5312.88726.384539.53471.9439吸收率01.18701.291400.98771.0632其中X0、y0、half1、half2分別表示橢圓中心橫縱坐標和兩個半軸。數(shù)據(jù)均以左下角的點為坐標原點建立坐標系求得。4.3 問題三模型的建立與求解建立模型與第二問類似，仍利用radon變換的思想建立模型，將空間每一條射線所投影的點還原回一條直線，將數(shù)據(jù)合并后重建未知介質的幾何性質與吸收率。4.3.2模型求解將附件5的數(shù)據(jù)導入MATLAB，利用radon逆變換將附件中數(shù)據(jù)重建未知介質信息，得到圖像如圖4.3.2.1。圖4.3.2.1 radon逆變換重建介質圖像（2）其中圖三圖四顯

25、示完全，圖二圖四經(jīng)角度修正后得出。為具體算出所要求十個點的吸收率，將數(shù)據(jù)帶入MATLAB中運算后結果如下表所示：表4.3.2.1 所求十點的吸收率（2）序號12345吸收率0.07362.53057.13270.03320.8015序號678910吸收率3.14875.57140.05308.13620.0272根據(jù)radon逆變換的結果，將得的數(shù)據(jù)導入Excel表格，在不同范圍內(nèi)的數(shù)值填充成不同顏色對不同吸收率的部分進行色塊填充，得到結果如圖4.3.2.2所示。 圖4.3.2.2 不同區(qū)域吸收率關系圖（2）其中無色處吸收率為0，黃色吸收率為0到2，紅色為2到4，淺藍色為4到6，綠色為6到8，

26、紫色為8到9。從表格中大致取出圖案中心，調用程序得到圖案中心在（51.2450，47.9210）(若以橢圓中心為原點，則圖案中心在（1.2450，-2.0790）)附近。4.4 問題四的分析與求解：為了便于求相鄰探測器之間的距離，考慮最好仍選擇圓形模板，但是要適當增大圓形模板的半徑，從而減小偶然誤差；為了利用180次旋轉得到的投射圖像求出旋轉中心的坐標，鑒于在原來的標定模板中橢圓和圓的內(nèi)公切線的選擇有較大誤差，新的模板中應盡量使得切線容易取得，且做出投射圖像后容易找到極值，并且圖像應有一定的對稱性；圖像掃描后應盡量少地得出重復數(shù)據(jù)，因此兩個圖像差別應較大；考慮到方便地識別投射位置以確定旋轉中心

27、，應至少設置兩個模板且相互隔開。基于以上分析，設計出如下新模板：圖4.4.1 新設計標定模板 標定方法：首先借助于圓形模板很容易求得相鄰探測器間距，利用對稱性，更容易求得射線水平、豎直和一條傾斜方向的位置，因此用和第1問相同的思路，此模板相對來說更能準確地確定CT系統(tǒng)的參數(shù)。5.結果的分析與檢驗5.1問題一結果分析與檢驗第一問求得的單元間距，旋轉中心與旋轉角度與所搜集資料的實際值相差不大且符合現(xiàn)實認知，在運算旋轉角度時，由于線性關系良好，可以印證旋轉中心與單元間距計算誤差不大。5.2問題二結果分析第二問十個位置的吸收率計算，由于是直接由題目所給數(shù)據(jù)計算得出，結果較精確，所得結果與逆radon變

28、換所得圖形相對應。在求解具體坐標與幾何關系時，通過Excel表格數(shù)據(jù)直接計算得出，可能存在誤差，但由具體結果運用radon變換檢驗后可以看出，基本符合題目所給數(shù)據(jù)。6.模型的優(yōu)缺點分析第一問中求單元間距與解旋轉中心時利用了切線的特殊性質，但切線的選擇不一定準確，因為射線寬度遠小于探測器寬度且所用射線不一定恰好為標定模板切線位置。從題目中可以得出，由于x射線之間得間距相等，不管x射線怎么旋轉，穿過托盤上圓的x射線條數(shù)應該是大致相等的。根據(jù)附件二可以發(fā)現(xiàn)，前13組穿過圓的x射線條數(shù)都是29條，由此可以粗略的計算出x射線之間的距離為2R29=0.2758。但是這種做法并不精確，因為在29條射線中最兩

29、邊的x射線并不是與圓相切的，為盡量減小誤差，本文利用圖4進行誤差檢驗，則如下方程成立：R2-xn2-R2-xn-12=d設模板的吸收率為p，第n個數(shù)據(jù)為sn，上述方程可以轉化為：R2-p2sn2-R2-p2sn-12=d即2R2-p2sn2-R2-p2sn-12-R2-p2sn+12=0利用MATLAB編程帶入多組數(shù)據(jù)求其平均值可得k=1.7725，進而求得d=0.2770.從這里可以看出在圓兩邊的射線不管是否與圓完全相切，對結果的影響不是很大，所以該建模方式可以使用。7.模型的改進簡單地說第2，3問根據(jù)逆radon矩陣求出坐標值存在一定的偶然誤差，計算量大。而且本文的模型沒有考慮到噪聲等其他

30、因素的影響，因此輸出圖像模糊有光暈。為了得到清晰的圖像，可以進行頻域濾波。首先二維傅里葉變換對為：引入傅里葉切片定理，其中是頻率分量：這說明一個投影的一維傅里葉變換，是二維投影矩陣的二維傅里葉變換的一個切片，執(zhí)行換元操作后，引入窗函數(shù)計算積分計算式并濾波，從而得到一個相對較好的結果。參考文獻：工業(yè)CT技術 劉豐林工業(yè)CT系統(tǒng)旋轉中心定位方法研究   劉明進附錄%程序1.1 作附件1的數(shù)據(jù)分布圖axis equalfor i=1:256 for j=1:256 if(A(i,j)>0) plot(j,257-i,'r*') hold on en

31、d endend%程序1.2 作附件2的數(shù)據(jù)分布圖for i=1:512 for j=1:180 if(AS(i,j)>0) if(AS(i,j)>100) plot(j,513-i,'b*-') hold on else plot(j,513,'r*-') hold on end end endendhold off%程序1.3 計算探測系統(tǒng)的增益率及探測器的平均距離ticd=ones(86,1);%統(tǒng)計圓模板對應探測器的平均個數(shù)m=zeros(86,1);%統(tǒng)計圓模板對應探測器的最大吸收率平均值for i=1:14 k=0;x=0; for j=

32、374:430 if(AS(j,i)>0) k=k+1; end if(AS(j,i)>x) x=AS(j,i); end end d(i)=k;m(i)=x;endfor i=109:180 k=0;x=0; for j=45:110 if(AS(j,i)>0) k=k+1; end if(AS(j,i)>x) x=AS(j,i); end end d(i-94)=k;m(i-94)=x;end a=mean(d)b=mean(m)toc%程序1.4 求出豎直射線近似方向d=zeros(180,1);k=0for j=1:180 if(k<max(AS(:,j)

33、 d(j)=max(AS(:,j); endenddfor i=1:180 if(k<d(i) k=d(i) i endend%程序1.5 推出水平射線近似方向d=zeros(180,1);u=1.7721;for j=14:109 d(j)=max(AS(:,j); d(j) jend%程序1.6 根據(jù)你和公式用matlab編寫程序：th=zeros(1,180);for i=1:1:180 th(i)=0.9938*i+28.718;endth計算得出所有的方向為：1 至 15 列 29.7118 30.7056 31.6994 32.6932 33.6870 34.6808 35.

34、6746 36.6684 37.6622 38.6560 39.6498 40.6436 41.6374 42.6312 43.6250 16 至 30 列 44.6188 45.6126 46.6064 47.6002 48.5940 49.5878 50.5816 51.5754 52.5692 53.5630 54.5568 55.5506 56.5444 57.5382 58.5320 31 至 45 列 59.5258 60.5196 61.5134 62.5072 63.5010 64.4948 65.4886 66.4824 67.4762 68.4700 69.4638 70.

35、4576 71.4514 72.4452 73.4390 46 至 60 列 74.4328 75.4266 76.4204 77.4142 78.4080 79.4018 80.3956 81.3894 82.3832 83.3770 84.3708 85.3646 86.3584 87.3522 88.3460 61 至 75 列 89.3398 90.3336 91.3274 92.3212 93.3150 94.3088 95.3026 96.2964 97.2902 98.2840 99.2778 100.2716 101.2654 102.2592 103.2530 76 至 90

36、 列 104.2468 105.2406 106.2344 107.2282 108.2220 109.2158 110.2096 111.2034 112.1972 113.1910 114.1848 115.1786 116.1724 117.1662 118.1600 91 至 105 列 119.1538 120.1476 121.1414 122.1352 123.1290 124.1228 125.1166 126.1104 127.1042 128.0980 129.0918 130.0856 131.0794 132.0732 133.0670 106 至 120 列 134.

37、0608 135.0546 136.0484 137.0422 138.0360 139.0298 140.0236 141.0174 142.0112 143.0050 143.9988 144.9926 145.9864 146.9802 147.9740 121 至 135 列 148.9678 149.9616 150.9554 151.9492 152.9430 153.9368 154.9306 155.9244 156.9182 157.9120 158.9058 159.8996 160.8934 161.8872 162.8810 136 至 150 列 163.8748 1

38、64.8686 165.8624 166.8562 167.8500 168.8438 169.8376 170.8314 171.8252 172.8190 173.8128 174.8066 175.8004 176.7942 177.7880 151 至 165 列 178.7818 179.7756 180.7694 181.7632 182.7570 183.7508 184.7446 185.7384 186.7322 187.7260 188.7198 189.7136 190.7074 191.7012 192.6950 166 至 180 列 193.6888 194.6826 195.6764 196.6702 197.6640 198.6578 199.6516 200.6454 201.6392 202.6330 203.6268 204.6206 205.6144 206.6082 207.6020%程序2.1 畫出未知介質1的幾何形狀I1=iradon(AX,0:179);I2=iradon(AX,0:179,'linear','Hann');I3=iradon(AX,0:179,'nearest','Ram-Lak');I4=iradon(AX,0: