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

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

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

3、的吸收率。本文依据 CT层析成像原理，利用逆拉东变换作出重构图像，并利用Exce I中数据分布计算出原位置介质的相关性质。所求十个位置的吸收率为：序号12345吸收率0.0185P 0.9701-0.00011.17031.0287序号678910吸收率1.4425P 1.2815-0.0074-0.00250.0194问题三中，与问题二计算方法相似得十个位置的吸攵率为：序号12345吸收率0.07362.53057.13270.0332P 0.8015序号678910吸收率3.14875.57140.05308.1362:0.0272问题四中，要求对问题一中的标定模型进行改进以减小误差并增加

4、稳定性，本文利用在问题一求解过程中遇到的问题进行思考，首先适当增大模板减小偶然误差，其次做出投射图像后应容易找到极值， 并且图像应有一定的对称性；图像 扫描后应尽量少地得出重复数据。关键词：CT层析成像Radon变换与逆变换吸收率MATLAB算法1问题重述1.1问题背景CT(ComputedTomography)可以在不破坏样品的情况下，禾U用样品对射线能 量的吸收特性对生物组织和工程材料的样品进行断层成像，由此获取样品内部的结构信息。1.2相关信息本题介绍了一种二维CT系统，平行入射的X射线垂直于探测器平面，每个 探测器单元看成一个接收点，且等距排列。X射线的发射器和探测器相对位置固 定不变

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

6、信息，确定 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 个方向是相邻较密、不错位的。 （图中有色区域表示该点有吸收率，蓝色部分表示吸 收率大于 100）图 2.错误！未指定顺序。 .1附件 1的数据

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

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

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

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

12、方向照射后不会出现射线 分成两部分的情况；在II或III或IV区域，当吸收率出现最大值的时候射线也被 分成两部分，而不是像图1.2那样成为一部分，也排除；综合各因素可判断旋转 中心应该在V区域。4.132,发现此时得到两部分图像，说明两模板之间有一部分射线直接被探测器接收，直到第14列数据两部分图像结合在了一起，如图4.133。图4.1.3.2第一列数据分布图4.1.3.3第14列数据分布结合托盘的几何特征，在垂直于对称轴方向上应该会出现最大的吸收率，利 用MATLA求得出现最大吸收率的方向为第151个方向，在此方向结合增益率得 出的模板长度为79.9884 ： 80mm进一步验证了结果；平行

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

14、于x射线的发生装置是连续旋转的，所以在 512条射线中穿过模板的射线 条数应该是连续变化的，所以用MATLAB编写程序，计算每次旋穿过模板的x射 线条数，并画出图像如图4.134。图4.134穿过模板射线条数随旋转次数变化曲线图中横坐标为旋转次数，纵坐标为穿过模板的射线条数。从图中可以看出附 表2中的数据是按照射线发射装置旋转的顺序依次给出的，而且可以看出，180次旋转后装置共旋转了 180度，每次旋转的角度近乎相等。考虑到在前50次旋转中，穿过整个装置的x射线条数与穿过椭圆的条数相 等。设穿过椭圆的最边缘的射线到椭圆中心的距离为R,取椭圆中心为原点，模板的对称中心为x轴建立平面直角坐标系，由

15、于射线可以近似看成与椭圆和圆 忆十y2二R。都相切，可以得到以下方程组：图4.135前50次旋转角度变化曲线图4.135前50次旋转角度变化曲线化简以上方程组可得:402 - R2p£s - 15图4.135前50次旋转角度变化曲线图4.135前50次旋转角度变化曲线利用MATLAB®程得到图像如图4.135 所示图4.135前50次旋转角度变化曲线用同样方式得出最后面25组数据的图像为：图4.136后25次旋转角度变化曲线从图像中可以看出每次旋转角度近似为1度。为了更加精确地计算旋转角度以便据此得到估算其他角度地依据，本文设计了另外一种算法。从附表二中可以看出，前 15组

16、数据中穿过椭圆的射线与穿过 圆的射线没有交叉，所以数据中换算出来的最大吸收距离就是近似穿过椭圆中心 的射线被椭圆所截的距离，弦长公式为：(X2/L = 1' 152402,y = kx + m且这个距离由直线的斜率k唯一确定，已知k、m时可解得：1200P 十1？?I =宀J152k2-mz + 40?152Jt2 十 402所以用MATLAB编程计算得出了较为精确的前15次旋转角度的结果如下：29.649331.002731.557232.647633.679534.647835.648336.64903 7.647538.648639.648840.649141.648742.64

17、9743.6492可以看出，这些结果的线性相关性非常好，据此利用Excel进行拟合，得到如下图表：图4.1.3.7旋转角度与旋转次数拟合公式从图表中可以得出拟合公式为1-:,相关系数为：0.9996根据拟合公式利用MATLAB编写程序计算得出所有方向角度。4.2问题二模型的建立与求解通过分析问题与附件数据可以发现， 题目所给数据与介质相关性质的二维分 布具有对应关系，X射线将介质一条线上的性质投影为一点。据此，本文利用 radon变换与radon逆变换进行运算，以通过投影后的讯号重建原始未知介质相 关性质的二维分布。该变换的定义为：令密度函数f(X)=f(x,y是一个的定义域为R2 的紧致台(

18、compactsupport)。令R为radon变换的运算子，贝U Rf(x,y是一个定义在 R2空间中的直线L。其基本思想为：radon变换可以理解为图像在空间的投影，空间的每一点对应一 条直线，而radon变换是图像像素点在每一条直线上的积分。因此，图像中高灰度值的直线会在空间形成亮点，而低灰度值的线段在空间形成暗点。对直线的检 测转化为在变换区域对亮点、暗点的检测。Radon变换是一幅图像在一个特定的角度下的径向线方向的投影，一幅图像 的radon变换是每一个像素radon变换的集合。对于MATLAB中语句R=radon（l,theta），如果theta是一个标量，R则是一个包含在the

19、ta的列向量。女口 果theta是一个向量，R则是一个矩阵，矩阵的每一列是对应其中一个 theta的 radon变换。而radon变换的逆运算，就可以将CT系统对于每一直线上的X射线吸收率数据还原回未知介质的物理性质。radon反变换的公式是：该反变换操作比较简单，思路清晰，可以借助数学运算软件计算。radon逆变换重建介质图像（1）为使逆拉东变换后得到的图像与原图像大小相等，在逆radon变换公式中取362条射线，每两条射线相距0.2770，所以得到的图形长宽也为100，由第一问 求解得知，X射线从约29度位置开始旋转，为消除重建模型与真实模型的旋转 角度差异，将从29度逆radon变换得到

20、的数据导入Excel表格可得变换后图像坐标水平方向平移了 7 一，-，竖直方向平移了' ''I，所以将原图中的坐标按上述数值平移就得到了变换后的坐标。从模板的逆radon变换产生的矩阵中可以发现，模板中所有点对应的灰度都 近似为0.5，又由于模板的吸收率为1,所以相对比例近似为2。据此利用MATLAB 编程可以算出图中对应十个点的吸收率如下表所示：所求十点的吸收率（1）序号12345吸收率0.01850.9701-0.00011.17031.0287序号678910吸收率1.44251.2815-0.0074-0.0025T 0.0194然后根据逆radon变换的结果，

21、将得的数据导入Excel表格，在不同范围内的数 值填充成不同颜色由于拉东变换中 x射线之间的距离都是0.2770，所以表格中 两组数据在实际物体上的距离也是 0.2770。由-'上得知表格取 362组数据时，总长度与原图基本相同，此时可以长度为基础，计算坐标变换公 式。根据前文坐标变换的逆变换，用 MATLABS写程序运算各色块（由运算可得 各色块均为椭圆形）的坐标及长短轴相应数据。将问题二中数据进行radon逆变 换后的图像最低点在Excel行数和列数，将行数和列数乘以倍率 0.2770即距离 图像边界的距离，由于第一个图距离两边界的距离已知，可以得到平移的方向和距离，具体不同区域吸

22、收率关系图（1）以每个椭圆中心的吸收率代表整个椭圆的吸收率则 A，B，C, D，E，F的吸收率 分别为：0、1.1870、1.2914、0、0.9877、1.0632。不同吸收率介质在托盘中的位置、几何形状与吸收率椭圆编号ABCDEFx063.017546.259054.707540.580551.245048.3365y033.101575.067071.189030.331553.184055.5385half111.19205.25577.08529.645422.68871.8478half26.11913.465312.88726.384539.53471.9439吸收率01.187

23、01.291400.98771.0632其中XO、yO、halfl、half2分别表示椭圆中心横纵坐标和两个半轴。数据均以 左下角的点为坐标原点建立坐标系求得。4.3问题三模型的建立与求解431建立模型与第二问类似，仍利用radon变换的思想建立模型，将空间每一条射线所投 影的点还原回一条直线，将数据合并后重建未知介质的几何性质与吸收率。4.3.2模型求解将附件5的数据导入MATLAB利用radon逆变换将附件中数据重建未知介质 信息，得到图像如图4.321。图4.3.2.1radon 逆变换重建介质图像（2）其中图三图四显示完全，图二图四经角度修正后得出。为具体算出所要求十个点的吸收率，将数

24、据带入MATLAB运算后结果如下表所示：表4.3.2.1所求十点的吸收率（2）序号12345吸收率0.07362.53057.13270.03320.8015序号678910吸收率3.14875.57140.05308.13620.0272根据radon逆变换的结果，将得的数据导入 Excel表格，在不同范围内的数 值填充成不同颜色对不同吸收率的部分进行色块填充，得到结果如图4.322所示。图不同区域吸收率关系图（2）其中无色处吸收率为0，黄色吸收率为0到2，红色为2到4，浅蓝色为4 到6，绿色为6到8，紫色为8到9。从表格中大致取出图案中心，调用程序得到图案中心在（51.2450，47.92

25、10） （若以椭圆中心为原点，贝U图案中心在（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射线之间的距离 2R为 一|丄匚；。但是这种做法并不精确，因为在29条射线中最两边的x射线并不是与圆相切的， 为尽量减小误差，本文利用图4.1.1.1

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

29、入傅里叶切片定理，其中co是频率分量：这说明一个投影的一维傅里叶变换， 是二维投影矩阵的二维傅里叶变换的一个切 片，执行换元操作后， 引入窗函数计算积分计算式并滤波， 从而得到一个相对较 好的结果。参考文献：工业CT技术刘丰林工业CT系统旋转中心定位方法研究？?刘明进附录%程序 1.1 作附件 1 的数据分布图axis equalfor i=1:256for j=1:256if (A(i,j)>0)plot(j,257-i, 'r*' )hold onendendend%程序 1.2 作附件 2 的数据分布图for i=1:512for j=1:180if (AS(i,j

30、)>0)if (AS(i,j)>100)plot(j,513-i, 'b*-' )hold onelseplot(j,513, 'r*-' )hold onendendendendhold off%程序 1.3 计算探测系统的增益率及探测器的平均距离ticd=ones(86,1);%统计圆模板对应探测器的平均个数m=zeros(86,1);%统计圆模板对应探测器的最大吸收率平均值for i=1:14k=0;x=0;for j=374:430if (AS(j,i)>0)k=k+1;endif (AS(j,i)>x)x=AS(j,i);end

31、end d(i)=k;m(i)=x; end for 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=0 for j=1:180 if (k<max(AS(:,j) d(j)=max(AS(:,j);end end d for i=1:180 if (k<d(i) k=d(i

32、) i end end%程序 1.5 推出水平射线近似方向 d=zeros(180,1);u=1.7721;编写程序：for j=14:109 d(j)=max(AS(:,j); d(j) j end%程序 1.6 根据你和公式用 matlabth=zeros(1,180);for i=1:1:180 th(i)=0.9938*i+28.718; endth计算得出所有的方向为：1 至 15 列29.711830.705631.699432.693233.687034.680835.674636.668437.662238.656039.649840.643641.637442.631243.

33、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

34、.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.15381

35、20.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.9740121 至 135 列148.9678149.9616150.9554151.9492152.9430153.936

36、8154.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.7

37、198189.7136190.7074191.7012192.6950166 至 180 列193.6888194.6826195.6764196.6702197.6640198.6578199.6516200.6454201.6392202.6330203.6268204.6206205.6144206.6082207.6020%程序 2.1 画出未知介质I1=iradon(AX,0:179);I2=iradon(AX,0:179,I3=iradon(AX,0:179,I4=iradon(AX,0:179,I5=iradon(AX,0:179,1 的几何形状'linear' , 'Hann' ); 'nearest' , 'Ram-Lak' ) 'spline' , 'Cosine' ); 'pchip'