数学建模国赛一等奖论文

A题：CT系统参数标定及成像ParametercalibrationandimagingofCTsystem报告人：陈淑珂队友：贾梦颖端启航CONTENT问题重述RESTATEMENT01第一问FIRSTQUESTION02二三问SECONDANDTHIRDQUESTION03第四问FOURTHQUESTION04问题重述RESTATEMENT01问题重述RESTATEMENT01CT系统介绍一种典型的二维CT系统如图1所示，平行入射的X射线垂直于探测器平面，每个探测器单元看成一个接收点，且等距排列。X射线的发射器和探测器相对位置固定不变，整个发射-接收系统绕某固定的旋转中心逆时针旋转180次。对每一个X射线方向，在具有512个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信息。CT系统安装时往往存在误差，从而影响成像质量，因此需要对安装好的CT系统进行参数标定，即借助于已知结构的样品（称为模板）标定CT系统的参数，并据此对未知结构的样品进行成像。问题重述RESTATEMENT01建立模型，解决以下问题(1)在正方形托盘上放置两个均匀固体介质组成的标定模板，模板的几何信息如图2所示，相应的数据文件见附件1，其中每一点的数值反映了该点的吸收强度，这里称为“吸收率”。对应于该模板的接收信息见附件2。请根据这一模板及其接收信息，确定CT系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该CT系统使用的X射线的180个方向。(2)附件3是利用上述CT系统得到的某未知介质的接收信息。利用(1)中得到的标定参数，确定该未知介质在正方形托盘中的位置、几何形状和吸收率等信息。另外，请具体给出图3所给的10个位置处的吸收率，相应的数据文件见附件4。(3)附件5是利用上述CT系统得到的另一个未知介质的接收信息。利用(1)中得到的标定参数，给出该未知介质的相关信息。另外，请具体给出图3所给的10个位置处的吸收率。(4)分析(1)中参数标定的精度和稳定性。在此基础上自行设计新模板、建立对应的标定模型，以改进标定精度和稳定性，并说明理由。图2：几何模版图3：10个位置点第一问FIRSTQUESTION02第一问FIRSTQUESTION02MathematicaMathematica是一款科学计算软件，很好地结合了数值和符号计算引擎、图形系统、编程语言、文本系统、和与其他应用程序的高级连接。很多功能在相应领域内处于世界领先地位，它也是使用最广泛的数学软件之一。Mathematica的发布标志着现代科技计算的开始。Mathematica是世界上通用计算系统中最强大的系统。自从1988发布以来，它已经对如何在科技和其它领域运用计算机产生了深刻的影响。第一问FIRSTQUESTION02探测器单元间的距离直径为8mm的圆遮挡了29个探测单元，于是可以计算出探测单元之间的距离为0.2759mm，我们还根据当探测器处于模版正上方和正左方时，椭圆的长轴和短轴的长度以及其投影在探测器单元上的个数，分别计算得到的距离为0.2759mm和0.2727mm。旋转中心在正方形托盘中的位置旋转中心在探测器上的投影位于第256个探测单元上。以探测器在正上方时最左边的点的下方512个单元的地方建立坐标系，记为坐标系1，那么旋转中心的坐标为(256，256)。第一问FIRSTQUESTION02旋转中心在正方形托盘中的位置射线穿过的长度越长，衰减的能量越多，我们用数据画出了CT探测器接收到的衰减能量的分布图如下。我们假设探测器在正上方时左边为第1个单元，右边为第512个单元。当探测器位于模版正左方时（如右上图所示），圆中心与椭圆中心的连线和CT探测器垂直，所以检测到的数值最大的传感器对应的横坐标是235，即表示椭圆中心对应检测器单元序号为235。当探测器处于正下方，即两个物体成像的中心相距最远时（如右下图所示），小圆中心与椭圆中心的连线与CT射线垂直，所以得到了四组椭圆中心对应的探测器单元序数222、223、224、225，取平均值为223.5。第一问FIRSTQUESTION02旋转中心在正方形托盘中的位置所以椭圆中心在坐标系1中坐标为（288.5,235），经过坐标变换将旋转中心坐标转换到坐标系2内，得出旋转中心坐标为（-32.5，21）。结合上一部分求出的探测器单元之间的距离，可以计算出旋转中心在距正方形托盘左边41.0345，距上边44.2069处。由图二中给出的各个长度数据我们可以看出旋转中心在椭圆的内部，右图一为任一角度时的旋转中心与椭圆的位置关系，右图二为旋转中心在坐标系2中的位置。第一问FIRSTQUESTION02CT系统使用的X射线的180个方向建立几何模型，在选定的直角三角形中利用距离计算，运用反三角函数计算出探测器与水平线之间的角度，然后再转换成射线与水平线的夹角，也就是系统中射线的180个方向。要找出两个圆心在接收器上投影的位置，我们采用的方法是寻找极大值，通过圆心的射线一定比未通过圆心的射线衰减能量更多，所以通过圆心的射线在探测器上的接收数值一定是一个极大值。第一问FIRSTQUESTION02CT系统使用的X射线的180个方向arccosds探测器能收到两个不同片区的信号，也即是椭圆和圆在某一角度范围内对X光的遮挡没有重叠作用。此时，我们已知椭圆中心和圆心之间的距离，并且可以计算出二者在探测器上对应的距离，由于探测器与X光垂直，所以在图1-5的三角形中，有其中是探测器与水平线间的夹角，那么X光与水平线夹角为：90第一问FIRSTQUESTION02CT系统使用的X射线的180个方向探测器只能接收到一个片区的信号，也即是二者对X光的照射有重叠作用。所以我们将这一段单独拿出来，对数据进行筛选，选择那些被椭圆遮挡而未被圆遮挡的地方，进行曲线拟合，补齐被圆遮挡的部分，然后按照上面的几何模型选取点计算角度。由于函数值域为，所以超过的度数得到的结果为负，所以我们人为的做了些矫正，用减去得到的负值得到在内的角度。arccos[0,]180180[0,2]第一问FIRSTQUESTION02CT系统使用的X射线的180个方向次数角度次数角度次数角度129.40566189.6517121149.1928231.48686290.5258122150.5944331.48686391.5583123151.3184432.80806492.6679124152.0596534.08366593.5032125152.8193634.70596694.196126154.4011735.31856795.6332127155.2271837.10366896.34007128156.962937.68246997.40224129156.9621038.81807098.11175130157.80315180.1087111139.5291171199.12865281.1775112140.0737172201.16965381.8882113141.182173202.12265482.9521114142.3176174203.0385583.8448115143.4831175203.92035684.7210116144.6815176204.77295785.9552117145.2941177206.40085886.8247118146.5488178206.40085987.7429119147.192179207.94046088.7738120148.5132180208.6816二三问SECONDANDTHIRDQUESTION03二三问SECONDANDTHIRD03拉东变换(,)(,):cossinLgRfxydLLRxy2(,)1(,)cos()2gRRfrdRpr逆拉东变换每个投影的变换对应物体变换中的一条线，我们根据投影得到物体的每条线的话，就能够重构整个物体。滤波反投影fourierfourier二三问SECONDANDTHIRD03几何图形的重建由附件三和附件五我们知道未知介质的接收信息，并且可以用拉动算法画出相对应的正弦图。标定模版、待测介质1、2的正弦图如下二三问SECONDANDTHIRD03我们利用上面得到的正弦图采用滤波反投影法反投影出模版原本的几何形状。采用Mathematica中逆拉东算法，将数据信息反投影，得到了初步的几何图形。几何图形的重建二三问SECONDANDTHIRD03随后我们采用同样的方法制作了标定模版的逆拉东CT图来检验图形位置和角度的正确性。但是得到的图像与题中已给的图像信息不符，位置和角度都发生了偏移。我们把180组的每组数据选取接收信息的极大值（衰减能量最多的时候是射线经过椭圆圆心），将它对应到256的中心位置，做到对旋转中心角度的修正。所以以下是我们对图的进行矫正的过程。首先是旋转中心的矫正，如下图（左边是几何图形，右边是正弦图）：几何图形的重建二三问SECONDANDTHIRD03旋转角度进行修正考虑到前面求出的旋转方向中初始位置的角度不是，而是，所以我们把超过的部分截取放到前面，进行旋转中心位置的调整。矫正之后我们绘制了下列图像：几何图形的重建029.4056180二三问SECONDANDTHIRD03修正处理后我们得到二三问中要求的几何图形（左二右三）：几何图形的重建二三问SECONDANDTHIRD03对附件2中均匀标定介质接收信息进行逆拉东变换后得到的几何图像不同颜色区域对应不同的测量值，结果为一个的矩阵，记为，记附件一中吸收率的数据矩阵为，标定模板中吸收率为1的部分测量值分布图如下：吸收率的计算256256256256A256256A二三问SECONDANDTHIRD03当时，，显示该被测点处无介质；当时，计算对应的平均值为0.133。对附件三和附件五做同样的操作，我们分别得到了和的表格，于是我们如下比式：吸收率的计算ij0aij0aij1aija256256B256256Cijijijij0.1331bcbc然后导出了两个的Excel表格，就是待测物体1、2的吸收率的表格，256256二三问SECONDANDTHIRD03我们进行相应的坐标转换，由于题中给出的是点与边界之间的距离，所以我们要把长度坐标轴中点的坐标转换成导出的表格中的具体位置。特定十个点的吸收率256256二三问SECONDANDTHIRD03我们进行相应的坐标转换，由于题中给出的是点与边界之间的距离，所以我们要把长度坐标轴中点的坐标转换成导出的表格中的具体位置。特定十个点的吸收率256256二三问SECONDANDTHIRD03然后分别在待测介质1、2吸收率表格中找到相对应的点，就能得到对应点的吸收率。结果如下表特定十个点的吸收率长度坐标系坐标系待测介质1吸收率待测介质2吸收率10.000018.000025.6000209.92000.00000.000034.500025.000088.3200192.00001.01941.687643.500033.0000111.3600171.52000.00007.086245.000075.5000115.200062.72001.21170.000048.500055.5000124.1600113.92001.07130.000050.000075.5000128.000062.72001.52463.094656.000076.5000143.360060.16001.29606.279865.500037.0000167.6800161.28000.00000.000079.500018.0000203.5200209.92000.00007.618998.500043.5000252.1600144.64000.00000.0000256256256