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

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入射光强d ij点i→j的距离I0出射光强λ两直角坐标系间的比例系数P探测器探测到的能量h（x,y）待变换目标函数μ衰减系数q傅立叶变换原函数Δl单位方格宽度w吸收强度f（x,y）衰减系数函数T灰度值n扫过圆的探测器个数βw与t间的比值p i第i个点的坐标4.模型的建立与求解：4.1问题1的模型建立与求解：4.1.1：求初始位置角度和探测器单元间距离：由朗伯比尔定律可知：用一束强度为I 0的x 光照射物体，如果物体内部均匀，已知物体对于x 光的衰减系数为μ，那么出射光强I 与入射光强I 0存在如下关系：I =I 0ⅇ−μc探测器探测到的能量：p =In (I 0∕I ),可知，p =−μc在沿某个方向入射的X 射线，探测器探测的能量与X 射线经过的物体的长度成正比。

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

图4-1 朗伯比尔定律示意图 探测器接收到的能量p =∑f i ∆l=ln I0I ,当∆l 趋向于无穷小时，得到线积分p =∫f i dl针对假设，我们得知在某个方向上，p 可以看做在该方向上所有射线的投影集合，通过求解一定角度的p 来求解物体内部平面各点的衰减系数函数f(x,y)。

图4-2-1 模板示意图（单位：mm ） 图4-2-2 反射投影模板结构图· 利用MATLAB 绘制的附件中模板图中椭圆和圆的方程：椭圆方程为：x 2225+y 21600=1 圆方程为：（x −45）2+y 2=42选取椭圆短轴所在直线为x 轴，长轴所在直线方向为y 轴平行束反投影重建算法(Filter back projeCTion,abbr. FBP)傅立叶中心切片定理：假设f(x,y)是待重建物体的密度函数，p(t i ,α)为f(x,y)在角度α时的平行束投影，有傅立叶中心切片定理的数学表达式为：F I =F(ρ,α)F 1表示一维傅立叶变换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 )∞∞δ(xcosα+ysinα−x )dxdy, 其中0≤α≤π我们构建接收器接收到的信号和经过Radon变换得到的图形投影变换模型：P=-μ*p(x′,α)接收到衰减信号的接收器个数t(t i)正比于穿过图形的个数，如图所示：图4-3 接收到衰减信号的接收器个数从A题附件得知在DE列之后，图像分为两部分；分析可知，对于模型有一个椭圆和一个小圆，根据提出的模型，接收器接收到的信号p>0的区域也分为两个部分，因此，假设两者之间存在对应关系，选取从A列至FX列、从111行数据利用matlab统计存在p>0的个数，如图所示：图4-4 111行能接受到非零信号的探测器个数与转动次数示意图可知从DE列到最后72列探测器信号p>0的个数在一定范围内趋于一致，利用MATLAB求其平均值，n为个数，n̅=∑n i(p>0)求得n̅=28.8333，取整n̅=29，对应于小圆宽度为2*r=0.2857mmd=2∗rn̅−1对A题附件2进行n(p>0)统计，利用MATLAB得出如下图形：图4-5 能接受到非零信号的探测器个数与转动次数图通过图形能够得出最大值为289，对应的转动次数分别为58、59、60、61、62、63、64、64，取其均值作为存在最大值的转动次数，均值为61，最小值为137，对应的转动次数为150、152，取其均值作为存在最小值的转动次数，均值为151，根据图像是平滑的曲线，不存在任何突变，且不存在两组完全相同的探测器数据，因此可以假设CT为均匀旋转，且每次旋转1度，则初始位置为29度。