2017国赛A题省二等奖论文

1CT系统参数标定及成像摘要CT（ComputedTomography）可在不破坏样品的情况下利用样品对射线能量的吸收特性对生物组织和工程材料样品进行断层成像，获取样品内部的结构信息。针对问题一，首先将附件1模板的几何信息数据及附件2中模板的接收信息利用MATLAB软件进行处理，通过Radon变换及基于傅里叶变换的中心切片定理进行X射线CT系统的二维重建；其次根据重建的二维图像，找出旋转中心与几何中心之间的距离关系，利用MATLAB编程求出旋转中心与几何中心的距离10.8003m，旋转中心的坐标00(,)xy为(9.2734,5.5363)；然后将求探测器单元之间的距离问题转化为求小圆的模板直径长度与小圆直径所占探测器的个数之比，得出探测器单元之间的距离为0.2667d,对该问题进行模型改进后，求得10.2768d，为最优解；最后通过投影最宽与最窄的位置之间的关系，解出每个旋转方向的角度为0.9783，180个方向起始点为-56.74°，终止点为119.35°。针对问题二及问题三，通过Radon反变换，利用MATLAB将附件3及附件5的不同未知介质的接收信息进行处理，分别得到两介质各自在正方形托盘的位置图像。由图以及第一问，知未知介质的初始位置均偏离几何中心，对其校正使其旋转中心与几何中心重合，利用MATLAB软件编程，得到某未知介质的吸收率三维图像。通过取边界点求标准方程法确认问题二中的未知介质形状为椭圆；问题三中的图像明显为不规则图形。最后将附件4的十个位置坐标带入各自的程序中，用MATLAB求出了问题二、问题三中给定的十个位置的吸收率k。针对问题四，依据对问题一中模板的精度及稳定性分析，重建长方形模板，使CT系统的旋转中心与托盘的几何中心重合，观察放射角度与投影长度的关系，最终整理出的图像数据与实际探测数据基本吻合，说明其精确度较高，从图像显示出的线性关系也证明了重建模板的良好稳定性。关键词：Radon变换傅里叶变换Radon反变换二维图像重建吸收率2一、问题重述CT可以在不破坏样品的情况下，利用样品对射线能量的吸收特性对生物组织和工程材料的样品进行断层成像，由此获取样品内部的结构信息。整个发射-接收系统绕某固定的旋转中心逆时针旋转180次，对每一个X射线方向，在具有512个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信息。CT系统安装时往往存在误差，因此需要对安装好的CT系统进行参数标定，即借助于已知结构的样品（称为模板）标定CT系统的参数，并据此对未知结构的样品进行成像。建立相应的数学模型和算法，解决以下问题：（1）在正方形托盘上放置了两个由均匀固体介质组成的标定模板，根据附件1中该模板的数据、及附件2中的接收信息，确定CT系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该CT系统使用的X射线的180个方向。（2）利用附件3中某未知介质的接收信息和（1）中得到的标定参数，确定该未知介质在正方形托盘中的位置、几何形状和吸收率等信息，并利用附件4相应的数据文件具体给出图3所给的10个位置处的吸收率。（3）利用附件5中另一个未知介质的接收信息及（1）中得到的标定参数，给出该未知介质的相关信息，并具体给出图3所给的10个位置处的吸收率。（4）分析（1）中参数标定的精度和稳定性。在此基础上自行设计新模板、建立对应的标定模型，以改进标定精度和稳定性，并说明理由。二、问题分析2.1问题一本问要求我们通过给定的均匀固体介质组成的标定模板的信息，求解关于CT系统的一系列标定参数。而题中附件二中给出的是某模板的接收信息，我们需要将其扫描信息转化为模板的几何信息，并建立平行投影重建的系统模型，即可确定CT系统旋转中心在正方形托盘中的位置等参数。因此我们首先引进3Radon变换及基于傅里叶变换的中心切片定理，解决X射线CT系统的二维重建问题，得到标定模板的二维重建图。其次，进行CT系统投影旋转中心位置的确定。但由于CT系统在安装时往往存在误差等因素，使得旋转中心与几何中心不重合，找出两者之间的距离关系并作图，利用MATLAB编程求出旋转中心与几何中心的距离从而求出旋转中心的坐标。然后选取图中所得到的小圆形结合附件1中相应的数据，利用比例求解法，求出探测器单元之间的距离，但这种人工手动查询的方法存在测量误差，得到的结果不准确，因此在模型的改进中进一步提出相应的改进模型，使得出的结果精度提高。2.2问题二及问题三：题中给出了附件3利用上述CT系统得到的某未知介质的接收信息，要利用这些数据，则需要对这些数据进行处理变换，本问通过Radon反变换，通过MATLAB编程求解，即可以得到某介质在正方形托盘的位置。由于其旋转中心与几何中心并不重合，则要对其进行校正，对校正后的图像进行三维转化，建立某未知介质的吸收率三维图像，对图像进行坐标分析，可得出该图形的大体形状。并将题中所给定的十个具体点的位置坐标带入相应代码中，直接得出这十个位置的吸收率。第二问与第三问所要解决的问题本质是一样的，该问针对附件5的数据信息进行分析求解。对比第二问与第三问，这两个不同的未知介质的形状及吸收率都存在较大差异。2.3问题四经过问题一问题二的分析，我们知道了产生误差的主要原因。题目中要分析其精度和稳定性，并自己设计一套方案。因此我们将CT系统的旋转中心尽可能的做到与托盘的几何中心重合，它是重建集合坐标的原点位置，其误差会影响到CT图像上的伪影。重新建立的模板是个长方形模板，有利于我们去观察在不同角度下与投影长度的关系，分析数据的结果与CT系统照射方向相吻合进而证明了其精度，从得到图像的线性关系可以看出该次检测没有波动并全部覆盖，因此稳定性也得到了保证。4三、模型假设1、假设可将标定模板的质心看作质点；2、X射线在照射过程中忽略光的波动及光的衍射对实验的影响；3、假设每一个CT系统使用的X射线均匀变动；4、模板及模型重建过程中，假设所选模板的旋转中心及几何中心是重合的；5、假设实验所收集的数据能客观反映实际情况。四、符号说明序号符号符号说明1平面内直线到原点的距离2原点到直线的垂线与x轴的夹角3D整个图像平面4n探测器单元的个数5d探测器单元之间的距离6S像素坐标7Y叠加平均后得到的像素值8m旋转中心与几何中心的距离五、模型的建立与求解5.1问题一5.1.1平行投影二维重建（Radon变换及傅里叶变换重建二维图像）Radon变换是指将(,)xy平面空间中的一条直线l:=xcossiny做线积分，映射成Radon空间的点(,)的变换。即只要知道了一个未知二维分布函数的所有方向上的线积分，即可求出连续图像的Radon变换。根据简单的狄拉克函数（Delta函数）：1,0()0,0ttt并结合直线方程可求出Delta函数的另一种表示：51,cossin0(cossin)0,cossin0xyxyxy即在线l上的点(,)xy满足()1x，其他非l上的点则为()0x。综上，可将Radon变换表示为：(,)(,)(cossin)DRfxyxydxdyD为整个图像平面，(,)fxy为图像上某一点(,)xy的像素点灰度值，特征函数为狄拉克函数，为(,)xy平面内直线到原点的距离，为原点到直线的垂线与x轴的夹角。Radon变换可以理解为图像在变换域(,)空间的投影，(,)空间的每个点对应图像空间的一条直线上的积分，也可以理解为图像顺时针旋转角度之后在水平轴上的投影，特征函数使得积分沿着直线=xcossiny进行。(,)在这里实际上就是直线的Hough变换[1]。若求Radon变换的逆变换,即可求出被测物截面函数每一点的(,)fxy值。Radon变换的正变换和逆变换恰好对应了CT技术的射线投影过程和图像二维分布函数的重建过程（程序见附录1）。图1对附件1的0-90度采5个方向进行radon变换而在实际重建的过程中，可依据中心切片定理，将所知的投影数据通过简单6的傅里叶反变换得到标定模板截面的评估，即平行投影的一维傅里叶变换等同于原始标定模板的二维傅里叶变换的一个切片[2]。首先定义代表截面函数的二维傅里叶变换：2()(,)(,)juxvyFuvfxyedxdy则某角度下一条投影()Pt的傅里叶变换为：2()()jtFtedt若假设沿着0的直线方向，从平行投影的定义来看，可得到：0()(,)xfxydy将(,)xy坐标系统旋转得到''(,)xy坐标系统，如上图1，并将以上两式合并，借助旋转坐标定义可得：2(cos()sin())()(,)jxyFfxyedxdy综上所示，对(,)fxy在与x轴成12,,k角度的直线上的一系列投影逐个进行傅里叶变换，可得傅里叶变换后''(,)Fxy在与x轴成12,,k角度的直线上的值，按理论来说，当投影线束无限多，即可求得''(,)Fxy在频域上对应线束的对应直线的值，从而可利用傅里叶逆变换得到原图(,)fxy，如图2（程序见附录1）。图2几何信息的图像实现75.1.2CT系统投影旋转中心位置的确定1、据中国专利技术提供的一种CT系统投影旋转中心定位方法[3]，给出其详细步骤：a、采集360度投影数据，得到投影正弦图；b、设定合适阈值，对投影正弦图进行图像分割，将空气背景部分的像素值设为0，其余不变；c、对于B步骤得到的分割后的正弦图，将各个投影角度的投影数据按照对应的像素叠加平均，即投影正弦图的各行叠加求平均得到一维信号Y=F(S)，其中，S为像素坐标，Y为叠加平均后得到的像素值；d、利用二次多项式Y=As2+Bs+C采用最小二乘法拟合Y=F(S)；e、计算抛物线Y=As2+Bs+C的对称中心横坐标值S=|，此即为投影旋转中心坐标-2a值。2、按原理来说，探测器的中心即为正方形的几何中心，CT系统的中心也应该与正方形的几何中心相同，但由于CT系统在安装时往往存在误差、实验仪器对测量的影响等也能产生误差影响，使得CT系统的中心与正方形中心存在位置偏差。用中心探测器即第256个探测器所旋转的路径交点作为旋转装置的旋转中心，给出旋转中心与几何中心的距离关系。假设取模板的几何中心为坐标原点，分别以水平、竖直方向为横纵坐标建立如图3的二维平面直角坐标系，图3旋转中心与几何中心位置图像8利用MATLAB编程（程序见附录1）求得：旋转中心与几何中心的距离10.8003m，旋转中心的坐标00(,)xy为(9.2734,5.5363)。5.1.3探测器单元距离及X射线的方向1、分别将附件1及附件2中的数据通过MATLAB进行处理(程序见附录3)，可分别得到题中所给的两个均匀固体介质组成的标定模板的几何信息示意图及投影图，将两个图进行对比，并用X射线从不同视角进行照射，可得到两介质旋转过程中任意方向上的边缘曲线图，图4中左右两图分别给出的是采用1度、135度方向照射所得图像。图4分别采用1度及135度方向X射线照射投影图利用经过部分照射的小圆的表现来求距离，采用不同视角进行多次验证，并结合图4，进行游标移动，可得小圆直径方向上占用的探测器个数为30。由题可知，小圆的半径为4mm,则其直径为8mm。由于探测器单元之间的距离可等价转化为求小圆的模板直径长度与小圆直径所占探测器的个数之比，得出探测器单元之间的距离为80.266730d。2、要确定该CT系统使用的X射线的180个方向，首先假设每一个CT系统使用的X射线均匀变动，投影最宽的位置与投影最窄的位置相差90度，且旋转角度相差92个方向，则每个旋转方向的角度为900.978392。选择以x正方向为0°（第58个方向是0°，每个方向转0.9783°）逆时针转动，180个方向起始点为-56.74°，终止点为119.35°95.2问题二5.2.1Radon反变换原理Radon反