高教社杯全国大学生数学建模竞赛题目 CT系统参数标定及成像Word下载.docx
《高教社杯全国大学生数学建模竞赛题目 CT系统参数标定及成像Word下载.docx》由会员分享，可在线阅读，更多相关《高教社杯全国大学生数学建模竞赛题目 CT系统参数标定及成像Word下载.docx（29页珍藏版）》请在冰点文库上搜索。
2
3
4
5
吸收率
0.0185
0.9701
-0.0001
1.1703
1.0287
6
7
8
9
10
1.4425
1.2815
-0.0074
-0.0025
0.0194
问题三中，与问题二计算方法相似得十个位置的吸收率为：
0.0736
2.5305
7.1327
0.0332
0.8015
3.1487
5.5714
0.0530
8.1362
0.0272
问题四中，要求对问题一中的标定模型进行改进以减小误差并增加稳定性，本文利用在问题一求解过程中遇到的问题进行思考，首先适当增大模板减小偶然误差，其次做出投射图像后应容易找到极值，并且图像应有一定的对称性；
图像扫描后应尽量少地得出重复数据。
关键词：
CT层析成像Radon变换与逆变换吸收率MATLAB算法
1.问题重述
1.1问题背景
应用ngyong﷽﷽﷽﷽﷽﷽﷽﷽﷽﷽﷽﷽﷽﷽﷽﷽﷽﷽﷽﷽构造ientinparabolicproblems[J].CT（ComputedTomography）可以在不破坏样品的情况下，利用样品对射线能量的吸收特性对生物组织和工程材料的样品进行断层成像，由此获取样品内部的结构信息。
1.2相关信息
本题介绍了一种二维CT系统，平行入射的X射线垂直于探测器平面，每个探测器单元看成一个接收点，且等距排列。
X射线的发射器和探测器相对位置固定不变，整个发射-接收系统绕某固定的旋转中心逆时针旋转180次。
对每一个X射线方向，在具有512个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信息。
并且，为消除安装误差，需要对安装好的CT系统进行参数标定，即借助于已知结构的样品（称为模板）标定CT系统的参数，并据此对未知结构的样品进行成像。
题目附件中提供了标定模板的几何信息，接受信息，待测介质的接收信息以及图3所给位置的相应数据。
1.3需要解决的问题
1.3.1问题一
在正方形托盘上放置两个均匀固体介质组成的标定模板，模板的几何信息如图2所示，请根据这一模板及其接收信息，确定CT系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该CT系统使用的X射线的180个方向。
1.3.2问题二
附件3是利用上述CT系统得到的某未知介质的接收信息。
利用第一问中得到的标定参数，确定该未知介质在正方形托盘中的位置、几何形状和吸收率等信息。
另外，请具体给出图3所给的10个位置处的吸收率。
1.3.3问题三
利用给出上述CT系统得到的另一个未知介质的接收信息。
利用第一问中得到的标定参数，给出该未知介质的相关信息。
1.3.4问题四
分析第一问中参数标定的精度和稳定性。
在此基础上自行设计新模板、建立对应的标定模型，以改进标定精度和稳定性，并说明理由。
图1.CT系统示意图图2.模板示意图（单位：
mm）图3.10个位置示意图
2.问题分析
2.1问题一分析
结合图2和附件1表中数据，可以首先计算出CT系统探测器个数和模板长度度量的比值，运用程序1.1可以得出模板的几何形状如图2.1.1，可大致认为它是对称的。
对于附件2，由于对180个方向尚无清晰地认识，首先用同样的方式做出数据分布图2.1.2，观察到图像比较平滑，因此认为按表格的顺序180个方向是相邻较密、不错位的。
（图中有色区域表示该点有吸收率，蓝色部分表示吸收率大于100）
图2.1.1附件1的数据分布图图2.1.2附件2的数据分布图
在图2.1.2中易观查到有一红色条形图案，这是在不同的方向扫描到圆形时留下的，可以此为突破口首先求出探测器单元之间的距离。
求出探测器单元之间的距离与增益比率之后，可以根据几何关系，自己设立坐标系并通过数学运算计算出旋转中心。
最后在建立的坐标系内，将待求解的180个旋转方向转换成X射线的斜率进行数学运算。
2.2问题二分析
题目给出了未知介质的接收信息，要求出介质的相关信息，可以搜索相应的数学模型，将附件中给出的按照射线条数与旋转角度列成的表格，一一对应为相应坐标点的吸收情况，从而根据各坐标点的不同性质，还原回该未知介质的几何信息与吸收率等信息。
2.3问题三分析
问题三与问题二类似，可以大致看出数据分布更具有一般性，不容易描述出未知介质的相关信息，可以通过图形大致描绘出介质相关信息。
2.4问题四分析
问题四要求分析题目所给的二维CT系统，设计新的标定模板以提高原系统的精确度与稳定性。
可以搜索相关资料，根据第一问的求解思路与求解过程，以规避求解过程中因标定模板自身性质而出现的误差为原则进行思路拓展，以设计高精度与稳定性的标定模板。
3.模型假设
◆假设在射线经过介质时能量只损失在介质中，及不考虑衍射等现象；
◆假设附件中所给出数据是正确的、可以直接利用的；
◆假设旋转中心在相邻两条射线的中间直线上；
4.模型的建立与求解
4.1问题一模型的建立与求解
4.1.1建立模型
从题目中可以得出，由于x射线之间得间距相等，不管x射线怎么旋转，穿过托盘上圆的x射线条数应该是大致相等的。
可将穿过圆形标定模板的X射线与模板建立模型示意图如下：
图4.1.1.1X射线穿过圆形标定模板示意图
如图4.1.1.1所示可以设穿过圆的在圆心两旁的射线被圆所截的距离分别为
和
求解旋转中心时，以椭圆形标定模板短轴所在直线为X轴，长轴所在直线为Y轴建立直角坐标系，将托盘进行划分。
由附件二的数据分布图的不对称性可知，旋转中心应在对称轴某一侧，大致确定旋转中心方位后，根据旋转中心两侧探测器个数不变且同一探测器接收的与距旋转中心的距离不变具体确定旋转中心的位置。
求解X射线旋转角度时，设穿过椭圆的最边缘的射线到椭圆中心的距离为R，取椭圆中心为原点，模板的对称中心为x轴建立平面直角坐标系，X射线所在直线的斜率为k。
d
相邻探测器之间的距离
N
对应于圆形模板的射线条数或探测器个数
μ
标定模板的吸收率
λ
处理数据时的增益率
φ
圆形模板的直径
4.1.2符号说明
4.1.3模型求解
求探测器单元间距离：
通过MATLAB编程求出计算出的条形带平均涉及探测器个数为
，沿直径方向上的平均吸收率
；
由图二圆形模板直径
，计算出探测器之间的平均距离：
探测系统平均增益率：
求旋转中心：
分析图2.1.2可知，由于蓝色区域仅分布在后部分角度范围内，因此估计旋转中心的位置在对称轴的某一侧，红色区域分为两部分时表示该角度下由两部分射线分别照射经过椭圆形和圆形；
为了讨论方便，现对正方形托盘做出如下划分：
图4.1.3.1圆盘划分
简单分析可知，如果旋转点在I区域，则沿180个方向照射后不会出现射线分成两部分的情况；
在II或III或IV区域，当吸收率出现最大值的时候射线也被分成两部分，而不是像图1.2那样成为一部分，也排除；
综合各因素可判断旋转中心应该在V区域。
因为发射-接收系统逆时针旋转，且图2.1.2中两部分红色区域，间距缩小，说明旋转之后在垂直于发射-接收方向上二者的距离是缩短的，从而确定旋转中心在短轴的上半侧，先运用Excel对附件二第1列数据作图（即画出第一个方向上的扫描图像）如图4.1.3.2，发现此时得到两部分图像，说明两模板之间有一部分射线直接被探测器接收，直到第14列数据两部分图像结合在了一起，如图4.1.3.3。
图4.1.3.2第一列数据分布图4.1.3.3第14列数据分布
结合托盘的几何特征，在垂直于对称轴方向上应该会出现最大的吸收率，利用MATLAB求得出现最大吸收率的方向为第151个方向，在此方向结合增益率得出的模板长度为
进一步验证了结果；
平行于对称轴方向上最大的吸收率出现在沿对称轴的直线上，计算得出为第61个方向，此方向的模板长度
，也验证了结果。
根据极近似水平方向为第61方向、最大吸收率出现在235号探测器，极近似竖直方向为151方向、最大吸收率出现在223号探测器，设旋转中心到竖直轴、水平轴的距离分别为x,y。
为求解还需要另一个方向的等量关系，选取椭圆和圆的一条外公切线的方向
，经计算得出经椭圆与圆公切线所在直线的X射线为第372号射线。
则在所建立坐标系内，切线方程为：
进而列出二元方程组：
解之得:
因此旋转中心在所建立的坐标系中的坐标为
求180个旋转方向：
由于x射线的发生装置是连续旋转的，所以在512条射线中穿过模板的射线条数应该是连续变化的，所以用MATLAB编写程序，计算每次旋穿过模板的x射线条数，并画出图像如图4.1.3.4。
图4.1.3.4穿过模板射线条数随旋转次数变化曲线
图中横坐标为旋转次数，纵坐标为穿过模板的射线条数。
从图中可以看出附表2中的数据是按照射线发射装置旋转的顺序依次给出的，而且可以看出，180次旋转后装置共旋转了180度，每次旋转的角度近乎相等。
考虑到在前50次旋转中，穿过整个装置的x射线条数与穿过椭圆的条数相等。
设穿过椭圆的最边缘的射线到椭圆中心的距离为R，取椭圆中心为原点，模板的对称中心为x轴建立平面直角坐标系，由于射线可以近似看成与椭圆和圆
都相切，可以得到以下方程组：
化简以上方程组可得：
利用MATLAB编程得到图像如图4.1.3.5所示
图4.1.3.5前50次旋转角度变化曲线
用同样方式得出最后面25组数据的图像为：
图4.1.3.6后25次旋转角度变化曲线
从图像中可以看出每次旋转角度近似为1度。
为了更加精确地计算旋转角度以便据此得到估算其他角度地依据，本文设计了另外一种算法。
从附表二中可以看出，前15组数据中穿过椭圆的射线与穿过圆的射线没有交叉，所以数据中换算出来的最大吸收距离就是近似穿过椭圆中心的射线被椭圆所截的距离，弦长公式为：
且这个距离由直线的斜率k唯一确定，已知k、m时可解得：
所以用MATLAB编程计算得出了较为精确的前15次旋转角度的结果如下：
可以看出，这些结果的线性相关性非常好，据此利用Excel进行拟合，得到如下图表：
图4.1.3.7旋转角度与旋转次数拟合公式
从图表中可以得出拟合公式为：
相关系数为：
0.9996
根据拟合公式利用MATLAB编写程序计算得出所有方向角度。
4.2问题二模型的建立与求解
4.2.1建立模型
通过分析问题与附件数据可以发现，题目所给数据与介质相关性质的二维分布具有对应关系，X射线将介质一条线上的性质投影为一点。
据此，本文利用radon变换与radon逆变换进行运算，以通过投影后的讯号重建原始未知介质相关性质的二维分布。
该变换的定义为：
令密度函数f（X）=f（x,y）是一个的定义域为
的紧致台（compactsupport）。
令R为radon变换的运算子，则Rf（x,y）是一个定义在
空间中的直线L。
其基本思想为：
radon变换可以理解为图像在空间的投影,空间的每一点对应一条直线，而radon变换是图像像素点在每一条直线上的积分。
因此，图像中高灰度值的直线会在空间形成亮点，而低灰度值的线段在空间形成暗点。
对直线的检测转化为在变换区域对亮点、暗点的检测。
Radon变换是一幅图像在一个特定的角度下的径向线方向的投影，一幅图像的radon变换是每一个像素radon变换的集合。
对于MATLAB中语句R=radon（I,theta），如果theta是一个标量，R则是一个包含在theta的列向量。
如果theta是一个向量，R则是一个矩阵，矩阵的每一列是对应其中一个theta的radon变换。
而radon变换的逆运算，就可以将CT系统对于每一直线上的X射线吸收率数据还原回未知介质的物理性质。
radon反变换的公式是：
该反变换操作比较简单，思路清晰，可以借助数学运算软件计算。
4.2.2模型求解
在MATLAB上将附件2和附件3的数据导入，利用radon逆变换将附件中数据重建未知介质信息，得到图像如图4.2.2.1。
图4.2.2.1radon逆变换重建介质图像
（1）
为使逆拉东变换后得到的图像与原图像大小相等，在逆radon变换公式中取362条射线，每两条射线相距0.2770，所以得到的图形长宽也为100，由第一问求解得知，X射线从约29度位置开始旋转，为消除重建模型与真实模型的旋转角度差异，将从29度逆radon变换得到的数据导入Excel表格可得变换后图像坐标水平方向平移了
，竖直方向平移了
，所以将原图中的坐标按上述数值平移就得到了变换后的坐标。
从模板的逆radon变换产生的矩阵中可以发现，模板中所有点对应的灰度都近似为0.5，又由于模板的吸收率为1，所以相对比例近似为2。
据此利用MATLAB编程可以算出图中对应十个点的吸收率如下表所示：
表4.2.2.1所求十点的吸收率
（1）
然后根据逆radon变换的结果，将得的数据导入Excel表格，在不同范围内的数值填充成不同颜色对不同吸收率的部分进行色块填充，得到结果如图4.2.2.2所示。
再利用Excel表格中寻找各椭圆定点的坐标，由于拉东变换中x射线之间的距离都是0.2770，所以表格中两组数据在实际物体上的距离也是0.2770。
由
得知表格取362组数据时，总长度与原图基本相同，此时可以长度为基础，计算坐标变换公式。
根据前文坐标变换的逆变换，用MATLAB编写程序运算各色块（由运算可得各色块均为椭圆形）的坐标及长短轴相应数据。
将问题二中数据进行radon逆变换后的图像最低点在Excel行数和列数，将行数和列数乘以倍率0.2770即距离图像边界的距离，由于第一个图距离两边界的距离已知，可以得到平移的方向和距离，具体结果如表4.2.2.2所示。
图4.2.2.2不同区域吸收率关系图
（1）
以每个椭圆中心的吸收率代表整个椭圆的吸收率则A，B，C，D，E，F的吸收率分别为：
0、1.1870、1.2914、0、0.9877、1.0632。
表4.2.2.2不同吸收率介质在托盘中的位置、几何形状与吸收率
椭圆编号
A
B
C
D
E
F
x0
63.0175
46.2590
54.7075
40.5805
51.2450
48.3365
y0
33.1015
75.0670
71.1890
30.3315
53.1840
55.5385
half1
11.1920
5.2557
7.0852
9.6454
22.6887
1.8478
half2
6.1191
3.4653
12.8872
6.3845
39.5347
1.9439
1.1870
1.2914
0.9877
1.0632
其中X0、y0、half1、half2分别表示椭圆中心横纵坐标和两个半轴。
数据均以左下角的点为坐标原点建立坐标系求得。
4.3问题三模型的建立与求解
4.3.1建立模型
与第二问类似，仍利用radon变换的思想建立模型，将空间每一条射线所投影的点还原回一条直线，将数据合并后重建未知介质的几何性质与吸收率。
4.3.2模型求解
将附件5的数据导入MATLAB，利用radon逆变换将附件中数据重建未知介质信息，得到图像如图4.3.2.1。
图4.3.2.1radon逆变换重建介质图像
（2）
其中图三图四显示完全，图二图四经角度修正后得出。
为具体算出所要求十个点的吸收率，将数据带入MATLAB中运算后结果如下表所示：
表4.3.2.1所求十点的吸收率
（2）
根据radon逆变换的结果，将得的数据导入Excel表格，在不同范围内的数值填充成不同颜色对不同吸收率的部分进行色块填充，得到结果如图4.3.2.2所示。
图4.3.2.2不同区域吸收率关系图
（2）
其中无色处吸收率为0，黄色吸收率为0到2，红色为2到4，浅蓝色为4到6，绿色为6到8，紫色为8到9。
从表格中大致取出图案中心，调用程序得到图案中心在（51.2450，47.9210）（若以椭圆中心为原点，则图案中心在（1.2450，-2.0790））附近。
4.4问题四的分析与求解：
为了便于求相邻探测器之间的距离，考虑最好仍选择圆形模板，但是要适当增大圆形模板的半径，从而减小偶然误差；
为了利用180次旋转得到的投射图像求出旋转中心的坐标，鉴于在原来的标定模板中椭圆和圆的内公切线的选择有较大误差，新的模板中应尽量使得切线容易取得，且做出投射图像后容易找到极值，并且图像应有一定的对称性；
图像扫描后应尽量少地得出重复数据，因此两个图像差别应较大；
考虑到方便地识别投射位置以确定旋转中心，应至少设置两个模板且相互隔开。
基于以上分析，设计出如下新模板：
图4.4.1新设计标定模板
标定方法：
首先借助于圆形模板很容易求得相邻探测器间距，利用对称性，更容易求得射线水平、竖直和一条倾斜方向的位置，因此用和第1问相同的思路，此模板相对来说更能准确地确定CT系统的参数。
5.结果的分析与检验
5.1问题一结果分析与检验
第一问求得的单元间距，旋转中心与旋转角度与所搜集资料的实际值相差不大且符合现实认知，在运算旋转角度时，由于线性关系良好，可以印证旋转中心与单元间距计算误差不大。
5.2问题二结果分析
第二问十个位置的吸收率计算，由于是直接由题目所给数据计算得出，结果较精确，所得结果与逆radon变换所得图形相对应。
在求解具体坐标与几何关系时，通过Excel表格数据直接计算得出，可能存在误差，但由具体结果运用radon变换检验后可以看出，基本符合题目所给数据。
6.模型的优缺点分析
第一问中求单元间距与解旋转中心时利用了切线的特殊性质，但切线的选择不一定准确，因为射线宽度远小于探测器宽度且所用射线不一定恰好为标定模板切线位置。
根据附件二可以发现，前13组穿过圆的x射线条数都是29条，由此可以粗略的计算出x射线之间的距离为
但是这种做法并不精确，因为在29条射线中最两边的x射线并不是与圆相切的，为尽量减小误差，本文利用图4.1.1.1进行误差检验，则如下方程成立：
设模板的吸收率为p，第n个数据为
，上述方程可以转化为：
即
利用MATLAB编程带入多组数据求其平均值可得
，进而求得
.从这里可以看出在圆两边的射线不管是否与圆完全相切，对结果的影响不是很大，所以该建模方式可以使用。
7.模型的改进
简单地说第2，3问根据逆radon矩阵求出坐标值存在一定的偶然误差，计算量大。
而且本文的模型没有考虑到噪声等其他因素的影响，因此输出图像模糊有光晕。
为了得到清晰的图像，可以进行频域滤波。
首先二维傅里叶变换对为：
引入傅里叶切片定理，其中ω是频率分量：
这说明一个投影的一维傅里叶变换，是二维投影矩阵的二维傅里叶变换的一个切片，执行换元操作后，引入窗函数计算积分计算式并滤波，从而得到一个相对较好的结果。
参考文献：
工业CT技术刘丰林
工业CT系统旋转中心定位方法研究
刘明进
附录
%程序1.1作附件1的数据分布图
axisequal
fori=1:
256
forj=1:
if（A（i,j）>
0）
plot（j,257-i,'
r*'
）
holdon
end
end
%程序1.2作附件2的数据分布图
512
180
if（AS（i,j）>
100）
plot（j,513-i,'
b*-'
else
plot（j,513,'
r*-'
holdoff
%程序1.3计算探测系统的增益率及探测器的平均距离
tic
d=ones（86,1）;
%统计圆模板对应探测器的平均个数
m=zeros（86,1）;
%统计圆模板对应探测器的最大吸收率平均值
14
k=0;
x=0;
forj=374:
430
if（AS（j,i）>
k=k+1;
x）
x=AS（j,i）;
d（i）=k;
m（i）=x;
fori=109
升级会员 