2017年高教社杯全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽題目 CT系統(tǒng)參數(shù)標(biāo)定及成像

在R2空間中的直線L。K/(D=f(x)\dx\Jl其基本思想為：radon變換可以理解為圖像在空間的投影，空間的每一點(diǎn)對(duì)應(yīng)一條直線，而radon變換是圖像像素點(diǎn)在每一條直線上的積分。因此，圖像中高灰度值的直線會(huì)在空間形成亮點(diǎn)，而低灰度值的線段在空間形成暗點(diǎn)。對(duì)直線的檢測(cè)轉(zhuǎn)化為在變換區(qū)域?qū)α咙c(diǎn)、暗點(diǎn)的檢測(cè)。Radon變換是一幅圖像在一個(gè)特定的角度下的徑向線方向的投影，一幅圖像的radon變換是每一個(gè)像素radon變換的集合。對(duì)于MATLAB中語(yǔ)句R=radon(l,theta)，如果theta是一個(gè)標(biāo)量，R則是一個(gè)包含在theta的列向量。如果theta是一個(gè)向量，R則是一個(gè)矩陣，矩陣的每一列是對(duì)應(yīng)其中一個(gè)theta的radon變換。而radon變換的逆運(yùn)算，就可以將CT系統(tǒng)對(duì)于每一直線上的X射線吸收率數(shù)據(jù)還原回未知介質(zhì)的物理性質(zhì)。radon反變換的公式是：—yg(xcos6+1/sjn.0,ff)dB該反變換操作比較簡(jiǎn)單，思路清晰，可以借助數(shù)學(xué)運(yùn)算軟件計(jì)算。4.2.2模型求解在MATLAB上將附件2和附件3的數(shù)據(jù)導(dǎo)入，利用radon逆變換將附件中數(shù)據(jù)重建未知介質(zhì)信息，得到圖像如圖4.2.2.1。圖4.2.2.1radon逆變換重建介質(zhì)圖像(1)為使逆拉東變換后得到的圖像與原圖像大小相等，在逆radon變換公式中取362條射線，每?jī)蓷l射線相距0.2770，所以得到的圖形長(zhǎng)寬也為100，由第一問(wèn)求解得知，X射線從約29度位置開(kāi)始旋轉(zhuǎn)，為消除重建模型與真實(shí)模型的旋轉(zhuǎn)角度差異，將從29度逆radon變換得到的數(shù)據(jù)導(dǎo)入Excel表格可得變換后圖像坐標(biāo)水平方向平移了31X0.2770=8.587，豎直方向平移了22X0.2770=6.094，所以將原圖中的坐標(biāo)按上述數(shù)值平移就得到了變換后的坐標(biāo)。從模板的逆radon變換產(chǎn)生的矩陣中可以發(fā)現(xiàn)，模板中所有點(diǎn)對(duì)應(yīng)的灰度都近似為0.5,又由于模板的吸收率為1,所以相對(duì)比例近似為2。據(jù)此利用MATLAB編程可以算出圖中對(duì)應(yīng)十個(gè)點(diǎn)的吸收率如下表所示：表4.2.2.1所求E十點(diǎn)的吸收率（1）序號(hào)12345吸收率0.01850.9701-0.00011.17031.0287序號(hào)678910吸收率1.44251.2815-0.0074-0.00250.0194然后根據(jù)逆radon變換的結(jié)果，將得的數(shù)據(jù)導(dǎo)入Excel表格，在不同范圍內(nèi)的數(shù)值填充成不同顏色對(duì)不同吸收率的部分進(jìn)行色塊填充，得到結(jié)果如圖4.2.2.2所示。再利用Excel表格中尋找各橢圓定點(diǎn)的坐標(biāo)，由于拉東變換中x射線之間的距離都是0.2770，所以表格中兩組數(shù)據(jù)在實(shí)際物體上的距離也是0.2770。由100一0.2770?362得知表格取362組數(shù)據(jù)時(shí)，總長(zhǎng)度與原圖基本相同，此時(shí)可以長(zhǎng)度為基礎(chǔ)，計(jì)算坐標(biāo)變換公式。根據(jù)前文坐標(biāo)變換的逆變換，用MATLAB編寫(xiě)程序運(yùn)算各色塊（由運(yùn)算可得各色塊均為橢圓形）的坐標(biāo)及長(zhǎng)短軸相應(yīng)數(shù)據(jù)。將問(wèn)題二中數(shù)據(jù)進(jìn)行radon逆變換后的圖像最低點(diǎn)在Excel行數(shù)和列數(shù)，將行數(shù)和列數(shù)乘以倍率0.2770即距離圖像邊界的距離，由于第一個(gè)圖距離兩邊界的距離已知，可以得到平移的方向和距離，具體結(jié)果如表4.2.2.2所示。D，以每個(gè)橢圓中心的吸收率代表整個(gè)橢圓的吸收率則A，B，C，D，E，F(xiàn)的吸收率分別為：0、1.1870、1.2914、0、0.9877、1.0632。表4.222不同吸收率介質(zhì)在托盤(pán)中的位置、幾何形狀與吸收率橢圓編號(hào)ABCDEFX063.017546.259054.707540.580551.245048.3365y033.101575.067071.189030.331553.184055.5385

halfl11.19205.25577.08529.645422.68871.8478half26.11913.465312.88726.384539.53471.9439吸收率01.18701.291400.98771.0632其中XO、yO、halfl、half2分別表示橢圓中心橫縱坐標(biāo)和兩個(gè)半軸。數(shù)據(jù)均以左下角的點(diǎn)為坐標(biāo)原點(diǎn)建立坐標(biāo)系求得。4.3問(wèn)題三模型的建立與求解4.3.1建立模型與第二問(wèn)類似，仍利用radon變換的思想建立模型，將空間每一條射線所投影的點(diǎn)還原回一條直線，將數(shù)據(jù)合并后重建未知介質(zhì)的幾何性質(zhì)與吸收率。4.3.2模型求解將附件5的數(shù)據(jù)導(dǎo)入MATLAB,利用radon逆變換將附件中數(shù)據(jù)重建未知介質(zhì)信息，得到圖像如圖4.3.2.1。100:1001Q020030030040Q400500100200300100:1001Q020030030040Q4005001002003001002003D0400500W0200300400500100200300圖4.3.2.1radon逆變換重建介質(zhì)圖像（2）其中圖三圖四顯示完全，圖二圖四經(jīng)角度修正后得出。為具體算出所要求十個(gè)點(diǎn)的吸收率，將數(shù)據(jù)帶入MATLAB中運(yùn)算后結(jié)果如下表所示：表4.3.2.1所求十點(diǎn)的吸收率（2）序號(hào)12345吸收率0.07362.53057.13270.03320.8015序號(hào)678910吸收率3.14875.57140.05308.13620.0272根據(jù)radon逆變換的結(jié)果，將得的數(shù)據(jù)導(dǎo)入Excel表格，在不同范圍內(nèi)的數(shù)值填充成不同顏色對(duì)不同吸收率的部分進(jìn)行色塊填充，得到結(jié)果如圖4.3.2.2所示。圖4.3.2.2不同區(qū)域吸收率關(guān)系圖（2）其中無(wú)色處吸收率為0，黃色吸收率為0到2，紅色為2到4，淺藍(lán)色為4到6，綠色為6到8，紫色為8到9。從表格中大致取出圖案中心，調(diào)用程序得到圖案中心在（51.2450，47.9210）（若以橢圓中心為原點(diǎn)，則圖案中心在（1.2450，-2.0790））附近。4.4問(wèn)題四的分析與求解：為了便于求相鄰探測(cè)器之間的距離，考慮最好仍選擇圓形模板，但是要適當(dāng)增大圓形模板的半徑，從而減小偶然誤差；為了利用180次旋轉(zhuǎn)得到的投射圖像求出旋轉(zhuǎn)中心的坐標(biāo)，鑒于在原來(lái)的標(biāo)定模板中橢圓和圓的內(nèi)公切線的選擇有較大誤差，新的模板中應(yīng)盡量使得切線容易取得，且做出投射圖像后容易找到極值，并且圖像應(yīng)有一定的對(duì)稱性；圖像掃描后應(yīng)盡量少地得出重復(fù)數(shù)據(jù)，因此兩個(gè)圖像差別應(yīng)較大；考慮到方便地識(shí)別投射位置以確定旋轉(zhuǎn)中心，應(yīng)至少設(shè)置兩個(gè)模板且相互隔開(kāi)?；谝陨戏治觯O(shè)計(jì)出如下新模板：圖4.4.1新設(shè)計(jì)標(biāo)定模板標(biāo)定方法：首先借助于圓形模板很容易求得相鄰探測(cè)器間距，利用對(duì)稱性更容易求得射線水平、豎直和一條傾斜方向的位置，因此用和第1問(wèn)相同的思路，此模板相對(duì)來(lái)說(shuō)更能準(zhǔn)確地確定CT系統(tǒng)的參數(shù)。結(jié)果的分析與檢驗(yàn)5.1問(wèn)題一結(jié)果分析與檢驗(yàn)第一問(wèn)求得的單元間距，旋轉(zhuǎn)中心與旋轉(zhuǎn)角度與所搜集資料的實(shí)際值相差不

大且符合現(xiàn)實(shí)認(rèn)知，在運(yùn)算旋轉(zhuǎn)角度時(shí)，由于線性關(guān)系良好，可以印證旋轉(zhuǎn)中心與單元間距計(jì)算誤差不大。5.2問(wèn)題二結(jié)果分析第二問(wèn)十個(gè)位置的吸收率計(jì)算，由于是直接由題目所給數(shù)據(jù)計(jì)算得出，結(jié)果較精確，所得結(jié)果與逆radon變換所得圖形相對(duì)應(yīng)。在求解具體坐標(biāo)與幾何關(guān)系時(shí)，通過(guò)Excel表格數(shù)據(jù)直接計(jì)算得出，可能存在誤差，但由具體結(jié)果運(yùn)用radon變換檢驗(yàn)后可以看出，基本符合題目所給數(shù)據(jù)。模型的優(yōu)缺點(diǎn)分析第一問(wèn)中求單元間距與解旋轉(zhuǎn)中心時(shí)利用了切線的特殊性質(zhì)，但切線的選擇不一定準(zhǔn)確，因?yàn)樯渚€寬度遠(yuǎn)小于探測(cè)器寬度且所用射線不一定恰好為標(biāo)定模板切線位置。從題目中可以得出，由于x射線之間得間距相等，不管x射線怎么旋轉(zhuǎn)，穿過(guò)托盤(pán)上圓的x射線條數(shù)應(yīng)該是大致相等的。根據(jù)附件二可以發(fā)現(xiàn)，前13組穿過(guò)圓的x射線條數(shù)都是29條，由此可以粗略的計(jì)算出x射線之間的距離為空=0.2758。29但是這種做法并不精確，因?yàn)樵?9條射線中最兩邊的x射線并不是與圓相切的，為盡量減小誤差，本文利用圖4.1.1.1進(jìn)行誤差檢驗(yàn)，則如下方程成立：JR2—X2—JR2_X2=dnn—1設(shè)模板的吸收率為p，第n個(gè)數(shù)據(jù)為每，上述方程可以轉(zhuǎn)化為:JR2—P2S2—JR2—P2S21=d=02^R2—P2S2—JR2—P2S2、—JR2—P2S￡=0利用MATLAB編程帶入多組數(shù)據(jù)求其平均值可得k=1.7725，進(jìn)而求得d=0.2770.從這里可以看出在圓兩邊的射線不管是否與圓完全相切，對(duì)結(jié)果的影響不是很大，所以該建模方式可以使用。模型的改進(jìn)簡(jiǎn)單地說(shuō)第2,3問(wèn)根據(jù)逆radon簡(jiǎn)單地說(shuō)第2,3問(wèn)根據(jù)逆radon矩陣求出坐標(biāo)值存在一定的偶然誤差，計(jì)算量大。而且本文的模型沒(méi)有考慮到噪聲等其他因素的影響，因此輸出圖像模糊有光暈。為了得到清晰的圖像，可以進(jìn)行頻域?yàn)V波。首先二維傅里葉變換對(duì)為：引入傅里葉切片定理，其中3是頻率分量:OCi6(?cos&-Fysiii—p)e仙兩曲dxdyIff/(^7cos+ysinQOCi6(?cos&-Fysiii—p)e仙兩曲dxdyoc-l了佃佰)oc-foopg=于匕切夕伽吩徳叭血旳血的J-QQ丿-OQ=F(ujcoe^jCJsiii^)這說(shuō)明一個(gè)投影的一維傅里葉變換，是二維投影矩陣的二維傅里葉變換的一個(gè)切片，執(zhí)行換元操作后，引入窗函數(shù)計(jì)算積分計(jì)算式并濾波，從而得到一個(gè)相對(duì)較好的結(jié)果。參考文獻(xiàn)：工業(yè)CT技術(shù)劉豐林工業(yè)CT系統(tǒng)旋轉(zhuǎn)中心定位方法研究劉明進(jìn)附錄%程序1.1作附件1的數(shù)據(jù)分布圖axisequalfori=1:256forj=1:256if(A(i,j)>0)plot(j,257-'r*')holonendendend%程序1.2作附件2的數(shù)據(jù)分布圖fori=1:512forj=1:180if(AS(i,j)>0)if(AS(i,j)>100)plot(j,513-'b*-')holonelseplot(j,5'r*-')hoonendendendendholdoff%程序1.3計(jì)算探測(cè)系統(tǒng)的增益率及探測(cè)器的平均距離ticd=ones(86,1);%統(tǒng)計(jì)圓模板對(duì)應(yīng)探測(cè)器的平均個(gè)數(shù)m=zeros(86,1);%統(tǒng)計(jì)圓模板對(duì)應(yīng)探測(cè)器的最大吸收率平均值fori=1:14k=0;x=0;forj=374:430if(AS(j,i)>0)k=k+1;endif(AS(j,i)>x)x=AS(j,i);endendd(i)=k;m(i)=x;endfori=109:180k=0;x=0;forj=45:110if(AS(j,i)>0)k=k+1;endif(AS(j,i)>x)x=AS(j,i);endendd(i-94)=k;m(i-94)=x;enda=mean(d)b=mean(m)toc%程序1.4求出豎直射線近似方向d=zeros(180,1);k=0forj=1:180if(k<max(AS(:,j)))d(j)=max(AS(:,j));endenddfori=1:180if(k<d(i))k=d(i)iendend%程序1.5推出水平射線近似方向d=zeros(180,1);u=1.7721;forj=14:109d(j)=max(AS(:,j));d(j)jend%程序1.6根據(jù)你和公式用matlab編寫(xiě)程序th=zeros(1,180);fori=1:1:180th(i)=0.9938*i+28.718;

endth計(jì)算得出所有的方向?yàn)?至15列29.711830.705631.699432.693233.687034.680835.674636.668437.662238.656039.649840.643641.637442.631243.625016至30列44.618845.612646.606447.600248.594049.587850.581651.575452.569253.563054.556855.550656.544457.538258.532031至45列59.525860.519661.513462.507263.501064.494865.488666.482467.476268.470069.463870.457671.451472.445273.439046至60列74.432875.426676.420477.414278.408079.401880.395681.389482.383283.377084.370885.364686.358487.352288.346061至75列89.339890.333691.327492.321293.315094.308895.302696.296497.290298.284099.2778100.2716101.2654102.2592103.253076至90列104.2468105.2406106.2344107.2282108.2220109.2158110.2096111.2034112.1972113.1910114.1848115.1786116.1724117.1662118.160091至105列119.1538120.1476121.1414122.1352123.1290124.1228125.1166126.1104127.1042128.0980129.0918130.0856131.0794132.0732133.0670106至120列134.0608135.0546136.0484137.0422138.0360139.0298140.0236141.0174142.0112143.0050143.9988144.9926145.9864146.9802147.9740

121至135列148.9678149.9616150.9554151.9492152.9430153.9368154.9306155.9244156.9182157.9120158.9058159.8996160.8934161.8872162.8810136至150列163.8748164.8686165.8624166.8562167.8500168.8438169.8376170.8314171.8252172.8190173.8128174.8066175.8004176.7942177.7880151至165列178.7818179.7756180.7694181.7632182.7570183.7508184.7446185.7384186.7322187.7260188.7198189.7136190.7074191.7012192.6950166至180列199.6516207.6020193.6888194.6826195.6764196.6702197.6640198.6578199.6516207.6020200.6454201.6392202.6330203.6268204.6206205.6144206.6082%程序2.1畫(huà)出未知介質(zhì)1的幾何形狀I(lǐng)1=iradon(AX,0:179);I2=iradon(AX,0:179,'linear','Hann');I3=iradon(AX,0:179,'nearest','Ram-Lak');I4=iradon(AX,0:179,'spline','Cosine');I5=iradon(AX,0:179,'pchip','Hamming');subplot(2,3,1),imshow(AX),title('附寸件3');subplot(2,3,2),imshow(I1),title('odinary');subplot(2,3,3),imshow(I2),title('linear.Hann');subplot(2,3,4),imshow(I3),title('nearest.Ram-Lak');subplot(2,3,5),imshow(I4),title('spline.Cosine');subplot(2,3,6),imshow(I5),title('pchip.Hamming');%程序2.2求十個(gè)位置吸收率x0=[10,34.5,43.5,45,48.5,50,56,65.5,79.5,98.5];y0=[18,25,33,75.5,55.5,75.5,76.5