2017年大学生数学建模A题CT系统标定成像论文

摘要：本文运用MATLAB等工具对已给出的数据进行分析和处理，通过反射投影算法，等比例转换法，radon变换和iradon变换，还原180次扫描信息和图形信息。

对于问题1，通过radon变换法，在MATLAB中得出该介质以正方形托盘左上角为原点的坐标系下的位置分布图，然后根据题目中已经给出的介质物体实际图，以椭圆圆心为原点建立直角坐标系，得出两个坐标系之间的比例关系，通过位置与长度的等比例变换得出旋转中心在正方形托盘中的坐标为（-8.7755, 6.1697），通过观察附件2发现存在探测器接收到的非零信号个数稳定在28个，对比小圆的直径得出探测器的间距为0.2857mm，探测器接收到的非零信号个数与角度曲线没有发生突变，且最高点与最低点横坐标相差90次，可以认为每次旋转1度，初始位置与坐标系X轴正方向夹角为29度。

对于问题2，通过使用iradon变换，得出了投影重建结构的解，对附件3中某未知介质的投影数据进行滤波反投影重建运算，实现从其它空间向图像空间进行转换的过程，最终通过MATLAB运行结果获得该未知介质模型的重建图像，得出该未知介质在正方形托盘中的几何形状和位置信息，然后采用比例变换的方式，根据10个点的位置和相对于实物图位置，得出这10个位置介质点的吸收率结果。

对于问题3，采用与问题2相似的方式，利用MATLAB中的iradon算法，根据附件5中提供的另一未知介质的吸收信息，通过反投影重建可以得到该未知介质的位置，形状和吸收率等信息，同样采用等比例变换的方式，根据点的位置和相对于实物图的位置，得出这10个位置点的吸收率结果。

对于问题4，通过对已经给定的数据进行分析，用iradon验证扫描次数对成像质量的影响，在不同滤波环境下比较成像质量，分别对18,36,90,180个角度投影进行观察和分析，能够得出随着投影角度个数的增加，图像的重影越来越少，也即是稳定性和精确度越来越高。运用shepp-lagon模型重新优化模型。

关键词：反射投影重建；MATLAB软件；radon变换；iradon变换；比例变换；成像质量；

1.问题的重述：

CT是一种利用样品对射线能量的吸收特性对生物组织和工程材料的样品进行进行断层成像的技术，由此可以获取样品的内部结构信息。本次使用的二维CT系统中平行入射的X射线垂直于探测器的平面，并且探测器单元可视为等距离排列的接收点，X射线的发射源和接收源的相对位置不变，整个系统绕着某个固定的旋转中心逆时针旋转180次，每个X射线方向，具有的512个等距离单位元探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信息。但是CT系统在安装时往往存在一定的误差而给测量质量带来一定的影响，本次借助于已知结构的样品标定CT系统的参数，题目要求建立相应的数学模型和算法解决以下四个问题：

（1）根据正方形托盘上的两个均匀介质模板的几何信息以及其所反应的吸收率和接收信息，确定CT系统旋转中心在正方形托盘中的具体位置，探测器单元之间的距离以及该CT系统所使用的X射线的180个方向；

（2）利用CT系统得到的某未知介质的接收信息以及1中的标定参数，确定未知介质在正方形托盘中的位置，几何形状和接收率，以及图3给定的10个位置处的吸收率；

（3）利用CT系统得到的另一未知介质的接收信息以及1中得到的标定参数，确定 该未知介质的信息，以及图3给定的10个未知处的吸收率；

（4）分析1中参数标定的精度和稳定性，自行设计新模板、建立对应的标定模型，改进标定精度和稳定性，给出理由；

2.问题的分析：

CT技术是一种依据外部投影数据重建被探测物体的内部结构的无损检测技术，具有高精度，高效率，无损检测等优点，近年来被广泛应用于生物医学，航空航天，材料加工，组织探伤等领域，图像重建算法是CT技术的核心之一，重建算法的优劣直接影响了CT检测性能的好坏，而对于未进行图像校准的CT检

测系统来说，确定其旋转中心和每次旋转的角度又极为重要，下边的几个问题就是解决CT旋转系统的中心，每次旋转的角度，通过建立CT成像算法，来进行模型检验和优化求解问题。

（1） 对于问题1，通过对附件中所给出的相应数据进行分析，附件1中是将图2模板示意图利用建立二维坐标系，将原本连续的模板划分成独立的小单元，每个单元都能吸收X射线，我们可以把每个单元看成一个个像素，X射线从一个方向射向模板，，穿过模板后被探测器检测到，而每个探测器接收到的信号强弱各有不同，因此可以简单分析认为，探测器的信号和被接受的X射线经过模板的长度有关，CT围绕旋转中心旋转进行探测，每次旋转过一定的角度，而模板并不是一个完全对称的物体，在各个方向上的长度各有不同，X射线通过的长度不同， 2中已经显示每次旋转得到的探测其能量和旋转角度关系，通过分析某一固定位置的探测器接收到不同能量的变换来求得每次旋转的角度，探测器的间隔可以根据模板的长度和探测到X射线能量变化的接收器个数来确定。

（2）对于问题2，需要使用在问题1中所求出的CT系统旋转中心在正方形托盘中的位置，以及探测器单元之间的距离和CT系统中X射线的方向这些结果，附件3给定了未知介质的接收信息数据，通过对其中的数据进行分析，使用与问题1互逆的方式，通过使用算法求出正方形托盘上未知介质的相对位置，几何形状和吸收率，另外还给出了正方形托盘中的10个点的位置，来得出这10个点的吸收率，在解决问题1时，是使用吸收率得到相对位置，因而求10个点的吸收率也是使用进行问题1的逆向求解过程。

（3）对于问题3，可以采用与问题2相似的方法，也是对问题1的过程进行逆向求解，从而得出该未知介质的相关信息，以及另外10个位置点的吸收率的大小。

（4）对于问题4，通过结合前面3个问题的结论与实现方法，自行设计一种模板和相应的模型，对于处于不同位置的介质点，其吸收率也会不尽相同，可以仿照附件1中的模板介质的吸收率和附件2与附件3的介质相关信息，先设定一组吸收率或者介质的相关信息，使用问题1的逆向求解过程，得出模板的位置和几何形状信息，并求出问题1和问题4中相应的参数标定误差和稳定性的大小。

3.模型的假设和符号的说明

3.1模型的假设：

（1）假设X射线穿过介质时，由于X射线与介质的原子相互作用过程中的散射，衍射等对X射线强度的衰减的影响可以忽略不计，也即是X射线衰减只与介质吸收有关；

（2）假设整个CT发射-接收系统绕着旋转中心旋转的过程是匀速进行的； （3）假设入射的X射线的能量分布是均匀的，且每次发射源发射的X射线都是相同的；

（4）假设介质是等密度分布的，单位质量介质中的电子数目相同； 3.2 符号的说明

符号说明表

I I0

入射光强 出射光强

dij λ

点i→j的距离 两直角坐标系间的比例系数

P

探测器探测到的能量

μ

衰减系数

q h（x,y）

待变换目标函数

傅立叶变换原函数

Δl f（x,y） n

单位方格宽度 衰减系数函数 扫过圆的探测器个数

pi

第i个点的坐标

w T β

吸收强度 灰度值 w与t间的比值

4.模型的建立与求解：

4.1问题1的模型建立与求解：

4.1.1：求初始位置角度和探测器单元间距离：

由朗伯比尔定律可知：用一束强度为I0的x光照射物体，如果物体内部均匀，已知物体对于x光的衰减系数为μ，那么出射光强I与入射光强I0存在如下关系：

I=I0ⅇ−μc

探测器探测到的能量：

p=In(I0∕I),可知，p=−μc

在沿某个方向入射的X射线，探测器探测的能量与X射线经过的物体的长度成正比。

fi为x-y平面的衰减系数函数，然后将物体划分成大小为∆l×∆l均匀的方格，使其离散化。

I0 f1 ∆l 图4-1 朗伯比尔定律示意图

f2 f3 fn I 探测器接收到的能量

p=∑𝑓𝑖∆l=lnI0,

当∆l趋向于无穷小时，得到线积分

p=∫𝑓𝑖𝑑𝑙

针对假设，我们得知在某个方向上，p可以看做在该方向上所有射线的投影集合，通过求解一定角度的p来求解物体内部平面各点的衰减系数函数f(x,y)。

I

图4-2-1 模板示意图（单位：mm） 图4-2-2 反射投影模板结构图·

利用MATLAB绘制的附件中模板图中椭圆和圆的方程： 椭圆方程为：

x2𝑦2

+=1 2251600

圆方程为：

（x−45）+𝑦2=42

选取椭圆短轴所在直线为x轴，长轴所在直线方向为y轴 平行束反投影重建算法(Filter back projeCTion,abbr. FBP) 傅立叶中心切片定理：

假设f(x,y)是待重建物体的密度函数，p(ti,α)为f(x,y)在角度α时的平行束投影，有傅立叶中心切片定理的数学表达式为：

F𝐼=𝐹(𝜌,𝛼)

F1表示一维傅立叶变换F(ρ,α)是二维傅立叶变换的极坐标表示。

Radon变换就是将数字图像矩阵在某以指定角度射线方向上做投影变换，例如，二维函数的投影就是在其指定方向上的线积分，在垂直方向上的线积分就是在x轴上的投影；在水平方向上的二维线积分就是在y轴上的投影。 设直角坐标系（x,y）转动α角后得到旋转坐标系（x′,y′），由此得知：

x′=xcosα+ysinα

p(x′,𝛼)为原函数f(x′,y′)的投影f(x,y)沿着旋转坐标系中x′轴α方向的线积分， 根据定义公式知其表达式为：

p(x′,𝛼)=∬f(x,y)𝛿(𝑥𝑐𝑜𝑠𝛼+𝑦𝑠𝑖𝑛𝛼−𝑥)𝑑𝑥𝑑𝑦, 其中0≤𝛼≤π ∞

∞

2

我们构建接收器接收到的信号和经过Radon变换得到的图形投影变换模型：

P=-μ*p(x′,𝛼)

接收到衰减信号的接收器个数t(ti)正比于穿过图形的个数，如图所示：

图4-3 接收到衰减信号的接收器个数

从A题附件得知在DE列之后，图像分为两部分；分析可知，对于模型有一个椭圆和一个小圆，根据提出的模型，接收器接收到的信号p>0的区域也分为两个部分，因此，假设两者之间存在对应关系，选取从A列至FX列、从111行数据利用matlab统计存在p>0的个数，如图所示：

图4-4 111行能接受到非零信号的探测器个数与转动次数示意图

可知从DE列到最后72列探测器信号p>0的个数在一定范围内趋于一致，利用MATLAB求其平均值，n为个数，n̅=∑𝑛𝑖(𝑝>0) 求得n̅=28.8333， 取整n̅=29，对应于小圆宽度为2*r

d=n=0.2857mm ̅−1

对A题附件2进行n(p>0)统计，利用MATLAB得出如下图形：

2∗𝑟

图4-5 能接受到非零信号的探测器个数与转动次数图

通过图形能够得出最大值为289，对应的转动次数分别为58、59、60、61、62、63、64、64，取其均值作为存在最大值的转动次数，均值为61，最小值为137，对应的转动次数为150、152，取其均值作为存在最小值的转动次数，均值为151，根据图像是平滑的曲线，不存在任何突变，且不存在两组完全相同的探测器数据，因此可以假设CT为均匀旋转，且每次旋转1度，则初始位置为29度。

4.1.2 CT系统旋转中心位置：

判断本CT系统的旋转中心，我们通过radon变换法，在MATLAB中模拟出介质物体的中心点的位置如下图所示，并且能够得出相应点的坐标位置：

图4-6 介质物体位置分布示意图

通过结合题目中给出的椭圆形模板和圆形模板，进行坐标系之间的等比例变换，通过建立相应的平面直角坐标系，可以分别得到椭圆圆心，圆的圆心在该平面直角坐标系里的坐标分别为：P1 (273.6658,291.6983）和P（，2415.2547,372.0849）假设旋转中心点P3（256.0000,256.0000），

由平面内两点之间的距离公式：d²=(Xi-Xj)²+(Yi-Yj)²可得：

d12=162.8171, d13=39.8303, d23=197.0730

由P1和P2点的坐标可得通过点P1和P2的直线的方程为：

Y1=0.5677*X+136.3382，

设P3点与直线Y1相交于P4点，由平面直线之间的垂直关系可得P4（244.9793,275.4129）,可得通过P3和P4的直线方程为Y2=-1.7615*X+706.9440，从而可得出交点P4与椭圆圆心PI和旋转中心P3之间的距离分别为d14=31.7511，d34=22.3230，

在题目中所给的模板示意图中，建立以椭圆中心为原点，以椭圆圆心和圆圆心连线为x轴的另一个平面直角坐标系如下：

图4-7 以椭圆圆心为原点的正方形托盘平面直角坐标系

可得在此平面内椭圆圆心Q1(0,0)，圆的圆心Q2(45,0),二者之间的距离D12=45，由此可得两个平面直角坐标系之间的相对比例系数λ=D12、d12，，记在该坐标系下的旋转中心坐标为Q3(x0,y0),由坐标系之间的对应关系发现Q3位于第二象限之内，

其横纵坐标分别为：x0= d14λ=-8.7755，y0= d34*λ=6.1697，

可得该CT系统的旋转中心坐标为：Q3 (-8.7755, 6.1697)，以在正方形托盘平面内，以椭圆圆心为原点建立直角坐标系，可得旋转中心点Q3相对于椭圆和圆的位置如下图所示：

图4-8 旋转中心的相对位置示意图

4.2 问题2的模型建立与求解：

4.2.1 iradon变换求未知介质的形状与结构：

在图像处理中，我们一般会将iradon函数用到X射线断层扫描中，irodon变换是radon变换的逆过程，给出了投影重建结构的解，下面利用傅里叶反变换计算iradon函数的表达式：

h(x,y)为需要进行变换的目标函数，用2-D傅里叶反变换写成极坐标形式为：

h(x,y)=dθπqp)dp , qH(qt)exp(j20π其中方括号内的表达式为∣q∣H(q,t)的1-D傅里叶反变换的表达式，结合傅里叶变换的卷积定理，将q的1-D傅里叶反变换写成为：

H-1 {∣q∣}= H-1 {q *sgnq}=H-1 (j2πq)○× H-1 (sgn q、j2π);

利用微分性质可得F-1 (j2π)= δ(p)，经过整理可得radon反变换的表达式为：

πh(x,y)=1/（2π²）0dθRh(p,t) ○×(1/2π²) Γ（1/p）

根据上述的radon反变换的表达式，也即是iradon变换，通过使用MATLAB中的iradon函数实现反投影重建算法，对附件3中已经给出的某未知介质的投影数据进行滤波反投影重建运算，从而从其它空间向图像空间进行转换的过程，最终通过MATLAB运行结果获得该未知介质模型的重建图像如下：

图4-9 未知介质模型重建图像示意图

从上图中能够直看出该未知介质在正方形托盘中的位置。也能够能够较为直观看出该未知介质在某个空间角度下的结构图。

4.2.2求未知介质的相关信息：

通过使用iradon函数将附件3中的接收信息矩阵W1重建，得到图： 然后将该图片的值矩阵导入到到excel图表中， 基于此原理，对附件2采用相同的方式，得到各自的灰度值矩阵。对于附件二，其模型介质均匀，吸收强度均为1，其灰度值保持在0.5左右（干扰值除外，即数值接近0的数值），根据附件二的结果和关系，可以得到吸收强度W与灰度值T之间的相关比例为：β=W、T=1、0.5=2，即吸收强度：

W=灰度值*2=2T;

因此，在去除干扰点后，把上式应用到未知介质的模型的灰度值矩阵上，可得到对应点上的吸收率，选取部分数据如下：

表4-1 第87-95行，159-161列的吸收率数据表：

87 88 89 90 91 92 93 94 95

FC 0.5874 0.5879 0.5927 0.5920 0.5903 0.5865 0.5737 0.5834 0.5841

FD 0.5893 0.5931 0.5793 0.5863 0.5917 0.5941 0.5945 0.5779 0.5795

FE 0.5917 0.5885 0.6001 0.6027 0.5888 0.5906 0.5807 0.5820 0.5861

4.2.3求10个位置点的吸收率：

采用比例变换的方式，根据附件4给出的10个点的相对位置：

图4-10 附件4中10个点的位置分布图

相对于实物图，可在吸收强度矩阵的对应位置中找到这10个目标点的吸收率如下：

表4-2 问题2的10个位置的吸收率

位置 1 2 3 4 5 6 7 8 9 10

X 10.0000 34.5000 43.5000 45.0000 48.5000 50.0000 56.0000 65.5000 79.5000 98.5000

Y 18.0000 25.0000 33.0000 73.5000 55.5000 75.5000 76.5000 37.0000 18.0000 43.5000

吸收率 0 0 1 0 1 0 0 1 0 0

4.3问题3的模型建立与求解： 4.3.1求未知介质相关信息：

采用与问题2相同的方式，利用MATLAB中的iradon算法，根据附件5中提供的另一未知介质的吸收信息，通过反投影重建可以得到有滤波和无滤波条件下该未知介质的几何形状分别为：

图4-11 有滤波和无滤波条件下未知介质几何形状图

从上面的两图中能够看出，在有滤波和无滤波条件下的未知介质的几何形状无明显差别。

4.3.2求10个位置点的吸收率：

采用比例变换的方式，根据附件4给出的10个点的相对位置，得出这10个位置点的吸收率为：

表4-3 问题3的10个位置的吸收率

位置 1 2 3

X 10.0000 34.5000 43.5000

Y 18.0000 25.0000 33.0000

吸收率 0 0 1.5470

4 5 6 7 8 9 10

45.0000 48.5000 50.0000 56.0000 65.5000 79.5000 98.5000

73.5000 55.5000 75.5000 76.5000 37.0000 18.0000 43.5000

1.4459 1.7115 3.3970 0 2.4807 0 0

4.4问题4的模型建立与求解： 4.4.1建立新模板与对应的标定模型：

针对参数标定的精度和稳定性，首先对题目给出的数据进行分析，本题目是在180个扫描方向上得到的扫描信息，用iradon验证扫描次数对成像质量的影响，另外，根据滤波反射投影模型，在不同滤波环境下比较成像质量，进一步分析确定参数标定的精度和稳定性。新模板采用Shepp-Logan模型进行标定，分别得到4组不同个数角度投影对应的图形如下：

图4-12 不同个数角度下iradon变换图

通过对上面的图进行对比：18个角度投影图中存在较多的重影，说明其稳定性和精确度较低，36个角度投影图中的重影比18个角度投影图中的重影少，90个角度投影图中的重影与180个角度投影图中的重影更少，并且通过观察可得：随

着投影角度数目的逐渐增多，所得到的投影图的重影数目越来越少，对图形稳定性和精确性的影响也越来越小。

5.模型的评价

1. 将连续的二维物体碎片化，使其分散成像素点，提升了运算速度； 2. 模型基于在180个方向所得的探测器信号进行iradon变换，基于傅里

叶变换和反射投影算法，利用MATLAB处理数据，该大大提升了数据处理能力；

3. 对模型中增添不同的滤波函数，提升成像的品质，将iradon变换得到

的灰度图类比吸收率简化算法；

4. 对旋转中心求解方法利用通过定点的相交直线求得； 5. 计算效率较低，模型处理的数据类型比较单一。

参考文献

[1] 高飞.MATLAB图像处理375例. [M]. 北京人民邮电出版社 2015 186-190

[2] 马晨欣. CT图像重建关键技术研究[D]. 河南. 解放军信息工程大学. 2011

[3] 张顺利，李卫斌，唐高峰. 滤波反射图像重建算法研究[N]. 咸阳师范学院

学报， 2008-8（23）

[4] 刘明进. 工业CT系统旋转中心定位方法研究[D]. 重庆 . 重庆大学2014

[5] 毛小渊. 二维CT图像重建算法研究[D]. 江西. 南昌航空大. 2016

[6] 孟凡勇，李忠传，杨民，李静海.CT旋转中心的精确确定方法[J]. 中国体视

学与图像分析 2013（18.4）

附录：

问题2程序如下： （1）图4-2-2程序：

A=xlsread('d:、A题附件',2); [R,F]=iradon(A,29:208,512); imshow(R)

（2）图4-4程序： clc

d=xlsread('F:\\建模相关\\A2',1);

jishu2=zeros(180,1);%jishu2为小圆对应的个数 for i=109:180; for j=1:111; if d(j,i)<=0;

jishu2(i,1)=jishu2(i,1)+1; end end end

111-jishu2(109:180) plot(111- jishu2)

xlabel('转动次数'),ylabel('非零信号探测器个数'); mean(111-jishu2)%求出平均值

title('111行能接受到非零信号的探测器个数与转动次数图') [ma,i1]=max(jishu)%求出最大数值和所在位置 （3）图4-5程序：

data=xlsread('F:\\建模相关\\A2',1); jishu=zeros(180,1); for i=1:180;

for j=1:512; if data(j,i)<=0;

jishu(i,1)=jishu(i,1)+1; end end end

plot(512-jishu)

xlabel('转动次数'),ylabel('非零信号探测器个数'); title('能接受到非零信号的探测器个数与转动次数图') [a,b]=max(512-jishu)%最大值和位置 [a,b]=min(512-jishu)%最小值和位置 （4）图4-6程序： A=xlsread('d:、A题附件',2); [R,F]=iradon(A,1,512); bw = im2bw (R,0.1); L = bwlabel(bw);

stats = regionprops(L, 'Centroid'); imshow(L); hold on

for i = 1 : length(stats) temp = stats(i).Centroid; temp(1) temp(2)

plot(temp(1), temp(2), 'r*'); end

（5）图4-7程序： t=0:pi、50:2*pi; x1=15*sin(t); y1=40*cos(t);

x2=45+4*sin(t); y2=4*cos(t); x3=0; y3=0; x4=45; y4=0;

plot(x1,y1,x2,y2,x3,y3,'*',x4,y4,'*'); axis([-50 50 -50 50]); grid on;

title('以椭圆圆心为原点的正方形托盘平面直角坐标系')

（6）图4-8程序： t=0:pi、50:2*pi; x1=15*sin(t); y1=40*cos(t); x2=45+4*sin(t); y2=4*cos(t); x3=-8.7755; y3=6.1697;

plot(x1,y1,x2,y2,x3,y3,'*'); axis([-50 50 -50 50]); grid on;

title('旋转中心的相对位置')

（7）图4-9程序：

A=xlsread('d:、A题附件',3); [R,F]=iradon(A,1，512); imshow(R)

（8）图4-11程序：

A=xlsread('d:、A题附件',5); Subplot(1,2)

[R,F]=iradon(A,1，512); imshow(R)

（9）图4-12程序：

I = xlsread('F:\\建模相关\\A1',1); P=I; imshow(P)

theta1=0:10:170;[R1,xp]=radon(P,theta1);%存在18个角度投影 theta2=0:5:175;[R2,xp]=radon(P,theta2);%存在36个角度投影 theta3=0:2:179;[R3,xp]=radon(P,theta3);%存在90个角度投影 theta4=0:1:179;[R4,xp]=radon(P,theta4);%存在90个角度投影 figure,

xlabel('\heta');ylabel('x\\prime');%? 定义坐标轴 I1=iradon(R1,10); % 开始反变换，还原图像 I2=iradon(R2,5); I3=iradon(R3,2); I4=iradon(R4,1);

subplot(221);imshow(I1),title('I1,18个角度投影'); subplot(222);imshow(I2),title('I2,36个角度投影'); subplot(223);imshow(I3),title('I3,90个角度投影'); subplot(224);imshow(I4),title('I4,180个角度投影');