高教社杯全国大学生数学建模竞赛题目CT系统参数标定及成像

1、2018西交数模第一次模拟赛数学建模论文首页选题 A 队伍编号 696班级学号姓名队长钱学森 73 班孔令辉队员 1 钱学森 73 班杨峥队员 2 电气 64 李佳航2018 年 6 月 30 日摘要本题针对一种二维CT获取样品内部结构信息的工作方式及成像原理，意在通过借助已知结构样品进行参数标定，消除系统误差， 而后对未知结构的样品进行成像，得出该未知介质的相关信息。问题一中，要求通过标定模板的相关，确定CT系统的旋转中心、探测器单元之间的距离以及该CT 系统使用的 X 射线的 180 个方向。本文利用标定模板几何参数，通过条数、 以及探测器单元间距相等等信息，首先计算出探测器单元之间的距离

2、为 0.2778mm。根据 X 射线与标定模板的几何特性，以椭圆短轴所在直线为 x 轴，长轴所在直线为 y 轴建立直角坐标系， 得旋转中心在所建立的坐标系中的坐标为7060.5,0396.9。最后通过近似确定X射线 180 个方向的旋转角度具有高度线性相关性， 根据部分确定数据拟合出整体旋转角度。前十五组旋转角度为：29.649331.002731.557232.647633.679534.647835.648336.649037.647538.648639.648840.649141.648742.649743.6492 问题二中，要求通过已求得的标定参数，确定未知介质在正方形托盘中的位置、

3、几何形状和吸收率等信息， 并具体给出所要求十个位置的吸收率。本文依据CT 层析成像原理，利用逆拉东变换作出重构图像，并利用Excel中数据分布计算出原位置介质的相关性质。所求十个位置的吸收率为：序号1 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 问题三中，与问题二计算方法相似得十个位置的吸收率为：序号1 2 3 4 5 吸收率0.0736 2.5305 7.1327 0.0332 0.8015 序号6 7 8 9 10 吸收率3.1487

4、5.5714 0.0530 8.1362 0.0272 问题四中，要求对问题一中的标定模型进行改进以减小误差并增加稳定性，本文利用在问题一求解过程中遇到的问题进行思考，首先适当增大模板减小偶然误差，其次做出投射图像后应容易找到极值，并且图像应有一定的对称性；图像扫描后应尽量少地得出重复数据。关键词： CT层析成像 Radon变换与逆变换吸收率MATLAB算法1.问题重述1.1 问题背景CT(ComputedTomography)可以在不破坏样品的情况下，利用样品对射线能量的吸收特性对生物组织和工程材料的样品进行断层成像，由此获取样品内部的结构信息。1.2 相关信息本题介绍了一种二维CT 系统，

5、平行入射的 X 射线垂直于探测器平面，每个探测器单元看成一个接收点， 且等距排列。 X 射线的发射器和探测器相对位置固定不变，整个发射 -接收系统绕某固定的旋转中心逆时针旋转180 次。对每一个X 射线方向，在具有 512 个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量， 并经过增益等处理后得到180 组接收信息。并且，为消除安装误差，需要对安装好的CT 系统进行参数标定，即借助于已知结构的样品（称为模板）标定CT 系统的参数，并据此对未知结构的样品进行成像。题目附件中提供了标定模板的几何信息，接受信息，待测介质的接收信息以及图3所给位置的相应数据。1.3 需要解决

6、的问题在正方形托盘上放置两个均匀固体介质组成的标定模板，模板的几何信息如图 2 所示，请根据这一模板及其接收信息，确定CT 系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该CT 系统使用的 X 射线的 180 个方向。附件 3 是利用上述 CT 系统得到的某未知介质的接收信息。利用第一问中得到的标定参数， 确定该未知介质在正方形托盘中的位置、几何形状和吸收率等信息。另外，请具体给出图3 所给的 10 个位置处的吸收率。利用给出上述 CT 系统得到的另一个未知介质的接收信息。利用第一问中得到的标定参数，给出该未知介质的相关信息。另外，请具体给出图3 所给的 10个位置处的吸收率。分析

7、第一问中参数标定的精度和稳定性。在此基础上自行设计新模板、建立对应的标定模型，以改进标定精度和稳定性，并说明理由。图 1.CT 系统示意图图2.模板示意图（单位：mm）图 3.10 个位置示意图2.问题分析2.1 问题一分析结合图 2 和附件 1 表中数据，可以首先计算出CT系统探测器个数和模板长度度量的比值， 运用程序 1.1 可以得出模板的几何形状如图2.1.1，可大致认为它是对称的。对于附件2，由于对 180 个方向尚无清晰地认识，首先用同样的方式做出数据分布图 2.1.2，观察到图像比较平滑，因此认为按表格的顺序180 个方向是相邻较密、不错位的。 （图中有色区域表示该点有吸收率，蓝色

8、部分表示吸收率大于 100）图 2.错误！未指定顺序。 .1附件 1 的数据分布图图 2.1.2 附件 2 的数据分布图求出探测器单元之间的距离与增益比率之后，可以根据几何关系，自己设立坐标系并通过数学运算计算出旋转中心。最后在建立的坐标系内，将待求解的180 个旋转方向转换成 X射线的斜率进行数学运算。2.2 问题二分析题目给出了未知介质的接收信息，要求出介质的相关信息，可以搜索相应的数学模型，将附件中给出的按照射线条数与旋转角度列成的表格，一一对应为相应坐标点的吸收情况， 从而根据各坐标点的不同性质， 还原回该未知介质的几何信息与吸收率等信息。2.3 问题三分析问题三与问题二类似，可以大致

9、看出数据分布更具有一般性，不容易描述出未知介质的相关信息，可以通过图形大致描绘出介质相关信息。2.4 问题四分析问题四要求分析题目所给的二维CT系统，设计新的标定模板以提高原系统的精确度与稳定性。 可以搜索相关资料， 根据第一问的求解思路与求解过程，以规避求解过程中因标定模板自身性质而出现的误差为原则进行思路拓展，以设计高精度与稳定性的标定模板。3.模型假设假设在射线经过介质时能量只损失在介质中，及不考虑衍射等现象；假设附件中所给出数据是正确的、可以直接利用的；假设旋转中心在相邻两条射线的中间直线上；4.模型的建立与求解4.1 问题一模型的建立与求解从题目中可以得出，由于x 射线之间得间距相等

10、，不管x 射线怎么旋转，穿过托盘上圆的 x 射线条数应该是大致相等的。 可将穿过圆形标定模板的X射线与模板建立模型示意图如下：和。求解旋转中心时，以椭圆形标定模板短轴所在直线为X轴，长轴所在直线为Y轴建立直角坐标系， 将托盘进行划分。 由附件二的数据分布图的不对称性可知，旋转中心应在对称轴某一侧， 大致确定旋转中心方位后， 根据旋转中心两侧探测器个数不变且同一探测器接收的与距旋转中心的距离不变具体确定旋转中心的位置。求解 X 射线旋转角度时，设穿过椭圆的最边缘的射线到椭圆中心的距离为R，取椭圆中心为原点，模板的对称中心为x 轴建立平面直角坐标系， X射线所在直线的斜率为 k。4.1.2 符号说

11、明4.1.3 模型求解求探测器单元间距离：通过 MATLAB 编程求出计算出的条形带平均涉及探测器个数为8023.28N，沿直径方向上的平均吸收率1766.14；由图二圆形模板直径mm8，计算出探测器之间的平均距离：探测系统平均增益率：求旋转中心：分析图 2.1.2 可知，由于蓝色区域仅分布在后部分角度范围内，因此估计旋转中心的位置在对称轴的某一侧， 红色区域分为两部分时表示该角度下由两部分射线分别照射经过椭圆形和圆形； 为了讨论方便，现对正方形托盘做出如下划分：图 4.1.3.1 圆盘划分简单分析可知，如果旋转点在I 区域，则沿 180 个方向照射后不会出现射线分成两部分的情况； 在 II

12、或 III 或 IV 区域，当吸收率出现最大值的时候射线也被分成两部分，而不是像图1.2 那样成为一部分，也排除；综合各因素可判断旋转中心应该在 V 区域。4.1.3.2，发现此时得到两部分图像，说明两模板之间有一部d 相邻探测器之间的距离N 对应于圆形模板的射线条数或探测器个数标定模板的吸收率处理数据时的增益率圆形模板的直径分射线直接被探测器接收，直到第14 列数据两部分图像结合在了一起，如图4.1.3.3。图 4.1.3.2第一列数据分布图4.1.3.3第 14 列数据分布结合托盘的几何特征，在垂直于对称轴方向上应该会出现最大的吸收率，利用 MATLAB 求得出现最大吸收率的方向为第151

13、 个方向，在此方向结合增益率得出的模板长度为80mm79.9884进一步验证了结果；平行于对称轴方向上最大的吸收率出现在沿对称轴的直线上，计算得出为第61 个方向，此方向的模板长度8mm3037.9988，也验证了结果。根据极近似水平方向为第61 方向、最大吸收率出现在 235 号探测器，极近似竖直方向为151方向、最大吸收率出现在223号探测器，设旋转中心到竖直轴、水平轴的距离分别为x,y 。为求解还需要另一个方向的等量关系，选取椭圆和圆的一条外公切线的方向，经计算得出经椭圆与圆公切线所在直线的X射线为第 372号射线。则在所建立坐标系内，切线方程为：进而列出二元方程组：解之得 : 因此旋转

14、中心在所建立的坐标系中的坐标为7060. 5,0396. 9。求 180 个旋转方向：由于 x射线的发生装置是连续旋转的，所以在512 条射线中穿过模板的射线条数应该是连续变化的， 所以用 MATLAB编写程序，计算每次旋穿过模板的x 射线条数，并画出图像如图4.1.3.4。图 4.1.3.4穿过模板射线条数随旋转次数变化曲线图中横坐标为旋转次数，纵坐标为穿过模板的射线条数。从图中可以看出附表 2 中的数据是按照射线发射装置旋转的顺序依次给出的，而且可以看出，180次旋转后装置共旋转了180 度，每次旋转的角度近乎相等。考虑到在前 50次旋转中，穿过整个装置的x 射线条数与穿过椭圆的条数相等。

15、设穿过椭圆的最边缘的射线到椭圆中心的距离为R ，取椭圆中心为原点，模板的对称中心为 x 轴建立平面直角坐标系，由于射线可以近似看成与椭圆和圆都相切，可以得到以下方程组：化简以上方程组可得：利用 MATLAB 编程得到图像如图 4.1.3.5所示图 4.1.3.5前 50 次旋转角度变化曲线用同样方式得出最后面 25组数据的图像为：图 4.1.3.6后 25 次旋转角度变化曲线从图像中可以看出每次旋转角度近似为1度。为了更加精确地计算旋转角度以便据此得到估算其他角度地依据，本文设计了另外一种算法。从附表二中可以看出，前15组数据中穿过椭圆的射线与穿过圆的射线没有交叉， 所以数据中换算出来的最大吸

16、收距离就是近似穿过椭圆中心的射线被椭圆所截的距离，弦长公式为：且这个距离由直线的斜率k唯一确定，已知 k、m时可解得：所以用 MATLAB编程计算得出了较为精确的前15 次旋转角度的结果如下：29.649331.002731.557232.647633.679534.647835.648336.649037.647538.648639.648840.649141.648742.649743.6492 可以看出，这些结果的线性相关性非常好，据此利用Excel 进行拟合，得到如下图表：图4.1.3.7旋转角度与旋转次数拟合公式从图表中可以得出拟合公式为：, 相关系数为： 0.9996 根据拟合公式

17、利用 MATLAB 编写程序计算得出所有方向角度。4.2 问题二模型的建立与求解通过分析问题与附件数据可以发现，题目所给数据与介质相关性质的二维分布具有对应关系， X射线将介质一条线上的性质投影为一点。据此，本文利用radon 变换与 radon 逆变换进行运算，以通过投影后的讯号重建原始未知介质相关性质的二维分布。该变换的定义为：令密度函数 f(X)=f(x,y)是一个的定义域为2R的紧致台 (compactsupport)。令 R为 radon 变换的运算子， 则 Rf(x,y)是一个定义在2R 空间中的直线 L。其基本思想为： radon 变换可以理解为图像在空间的投影, 空间的每一点对

18、应一条直线，而 radon 变换是图像像素点在每一条直线上的积分。因此，图像中高灰度值的直线会在空间形成亮点， 而低灰度值的线段在空间形成暗点。对直线的检测转化为在变换区域对亮点、暗点的检测。Radon变换是一幅图像在一个特定的角度下的径向线方向的投影，一幅图像的 radon 变换是每一个像素radon 变换的集合。对于MATLAB中语句R=radon(I,theta)，如果 theta 是一个标量， R则是一个包含在theta 的列向量。如果 theta 是一个向量， R则是一个矩阵，矩阵的每一列是对应其中一个theta 的radon 变换。而 radon变换的逆运算，就可以将CT 系统对于

19、每一直线上的X 射线吸收率数据还原回未知介质的物理性质。radon反变换的公式是：该反变换操作比较简单，思路清晰，可以借助数学运算软件计算。radon 逆变换重建介质图像（1）为使逆拉东变换后得到的图像与原图像大小相等，在逆 radon 变换公式中取362 条射线，每两条射线相距0.2770，所以得到的图形长宽也为100，由第一问求解得知， X射线从约 29 度位置开始旋转，为消除重建模型与真实模型的旋转角度差异，将从 29 度逆 radon 变换得到的数据导入Excel 表格可得变换后图像坐标水平方向平移了，竖直方向平移了， 所以将原图中的坐标按上述数值平移就得到了变换后的坐标。从模板的逆

20、radon 变换产生的矩阵中可以发现，模板中所有点对应的灰度都近似为 0.5 ， 又由于模板的吸收率为1， 所以相对比例近似为2。 据此利用 MATLAB编程可以算出图中对应十个点的吸收率如下表所示：所求十点的吸收率（1）序号1 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 然后根据逆 radon 变换的结果，将得的数据导入Excel 表格，在不同范围内的数值填充成不同颜色由于拉东变换中x 射线之间的距离都是0.2770，所以表格中两组数据在

21、实际物体上的距离也是0.2770。由得知表格取362 组数据时，总长度与原图基本相同，此时可以长度为基础，计算坐标变换公式。根据前文坐标变换的逆变换，用MATLAB 编写程序运算各色块（由运算可得各色块均为椭圆形） 的坐标及长短轴相应数据。 将问题二中数据进行radon 逆变换后的图像最低点在Excel 行数和列数，将行数和列数乘以倍率0.2770 即距离图像边界的距离， 由于第一个图距离两边界的距离已知，可以得到平移的方向和距离，具体不同区域吸收率关系图（1）以每个椭圆中心的吸收率代表整个椭圆的吸收率则A，B，C，D，E，F的吸收率分别为： 0、1.1870、1.2914、0、0.9877、

22、1.0632。不同吸收率介质在托盘中的位置、几何形状与吸收率椭圆编号A B C D E F x063.017546.259054.707540.580551.245048.3365y033.101575.067071.189030.331553.184055.5385half111.19205.25577.08529.645422.68871.8478half26.11913.465312.88726.384539.53471.9439吸收率0 1.1870 1.2914 0 0.9877 1.0632 其中X0、y0、half1、half2分别表示椭圆中心横纵坐标和两个半轴。数据均以左下角的

23、点为坐标原点建立坐标系求得。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）序号1 2 3 4

24、 5 吸收率0.0736 2.5305 7.1327 0.0332 0.8015 序号6 7 8 9 10 吸收率3.1487 5.5714 0.0530 8.1362 0.0272 根据 radon 逆变换的结果，将得的数据导入Excel 表格，在不同范围内的数值填充成不同颜色对不同吸收率的部分进行色块填充，得到结果如图4.3.2.2所示。图不同区域吸收率关系图（2）其中无色处吸收率为0，黄色吸收率为 0 到 2，红色为 2 到 4，浅蓝色为 4到 6，绿色为 6 到 8，紫色为 8 到 9。从表格中大致取出图案中心， 调用程序得到图案中心在 （51.2450，47.9210）( 若以椭圆中

25、心为原点，则图案中心在（1.2450，-2.0790 ）) 附近。4.4 问题四的分析与求解：为了便于求相邻探测器之间的距离，考虑最好仍选择圆形模板，但是要适当增大圆形模板的半径， 从而减小偶然误差； 为了利用 180 次旋转得到的投射图像求出旋转中心的坐标， 鉴于在原来的标定模板中椭圆和圆的内公切线的选择有较大误差， 新的模板中应尽量使得切线容易取得， 且做出投射图像后容易找到极值，并且图像应有一定的对称性； 图像扫描后应尽量少地得出重复数据，因此两个图像差别应较大； 考虑到方便地识别投射位置以确定旋转中心，应至少设置两个模板且相互隔开。基于以上分析，设计出如下新模板：图 4.4.1新设计标

26、定模板标定方法：首先借助于圆形模板很容易求得相邻探测器间距，利用对称性，更容易求得射线水平、 竖直和一条倾斜方向的位置， 因此用和第 1 问相同的思路，此模板相对来说更能准确地确定CT系统的参数。5.结果的分析与检验5.1 问题一结果分析与检验第一问求得的单元间距， 旋转中心与旋转角度与所搜集资料的实际值相差不大且符合现实认知， 在运算旋转角度时， 由于线性关系良好， 可以印证旋转中心与单元间距计算误差不大。5.2 问题二结果分析第二问十个位置的吸收率计算，由于是直接由题目所给数据计算得出，结果较精确，所得结果与逆 radon 变换所得图形相对应。 在求解具体坐标与几何关系时，通过 Excel

27、表格数据直接计算得出， 可能存在误差， 但由具体结果运用 radon变换检验后可以看出，基本符合题目所给数据。6.模型的优缺点分析第一问中求单元间距与解旋转中心时利用了切线的特殊性质，但切线的选择不一定准确，因为射线宽度远小于探测器宽度且所用射线不一定恰好为标定模板切线位置。从题目中可以得出，由于x 射线之间得间距相等，不管x 射线怎么旋转，穿过托盘上圆的x 射线条数应该是大致相等的。根据附件二可以发现，前13 组穿过圆的 x 射线条数都是 29 条，由此可以粗略的计算出x 射线之间的距离为。但是这种做法并不精确， 因为在 29 条射线中最两边的x 射线并不是与圆相切的，为尽量减小误差，本文利

28、用图4.1.1.1 进行误差检验，则如下方程成立：设模板的吸收率为p，第 n 个数据为，上述方程可以转化为：即利用 MATLAB编程带入多组数据求其平均值可得，进而求得.从这里可以看出在圆两边的射线不管是否与圆完全相切，对结果的影响不是很大，所以该建模方式可以使用。7.模型的改进简单地说第 2，3 问根据逆 radon 矩阵求出坐标值存在一定的偶然误差，计算量大。而且本文的模型没有考虑到噪声等其他因素的影响，因此输出图像模糊有光晕。为了得到清晰的图像，可以进行频域滤波。首先二维傅里叶变换对为：引入傅里叶切片定理，其中 是频率分量：这说明一个投影的一维傅里叶变换， 是二维投影矩阵的二维傅里叶变换

29、的一个切片，执行换元操作后， 引入窗函数计算积分计算式并滤波，从而得到一个相对较好的结果。参考文献：工业CT技术刘丰林工业 CT系统旋转中心定位方法研究?刘明进附录%程序 1.1作附件 1 的数据分布图axisequalfori=1:256forj=1:256if(A(i,j)0)plot(j,257-i,r*)holdonendendend%程序 1.2作附件 2 的数据分布图fori=1:512forj=1:180if(AS(i,j)0)if(AS(i,j)100)plot(j,513-i,b*-)holdonelseplot(j,513,r*-)holdonendendendendholdoff%程序 1.3计算探测系统的增益率及探测器的平均距离ticd=ones(86,1);%统计圆模板对应探测器的平均个数m=zeros(86,1);% 统计圆模板对应探测器的最大吸收率平均值fori=1:14k=0;x=0;forj=374:430if(