【数学建模大赛】CT系统参数的标定与成像问题

PAGEPAGE21CT系统参数标定及成像作者：芋圆摘要本文主要研究CT系统参数的标定与成像问题，通过建立模型利用已知结构样品标定CT系统的参数，并对未知结构样品的图像进行重建。针对问题一，先对圆分析，无论射线的方向如何，所得到的投影长度总是等于圆的直径，由此关系可求得相邻探测器单元的距离为0.2775mm。再利用几何知识求出圆和椭圆在某个投影方向上的投影长度与该方向角度的关系式，并将其有较大跳跃或异常的点去除，利用傅里叶变换拟合出一个函数，再将附件2中180个角度对应的数据与拟合函数在同一角度得出的数据做差，用最小二乘法拟合出一个更加精确的函数，再将附件2中的数据代入进去求出180个角度，其中前三个以y轴正方向的夹角分别是，，。以托盘左下角为原点，水平方向为横坐标，竖直方向为纵坐标建立直角坐标系，同样利用最小二乘法得出旋转中心为（40.7543,56.2520）。针对问题二，首先找到投影坐标与原图像坐标之间的转化关系，再利用滤波反投影法得出图像的衰减系数分布函数，即所有经过滤波后的射线投影在度范围的衰减系数的叠加。考虑到CT系统是逆时针旋转180次，每一次对应一个角度，射线对物体进行平行束扫描，得到是180个离散投影数据，于是将衰减系数分布函数离散化，将180个角度得到的投影数据叠加，得到图像的像素点，从而重建出图像，并算出了附件四给出的10个位置的吸收率分别为：0.0442、0.9839、1.2682、0.9755、1.0099、0.9891、0.0033、0.9731、0.0476、0.0307。针对问题三，利用问题二中的模型重建出图像并解出10个位置的吸收率分别为：1.4813、0.0603、2.1298、0.0036、3.3900、0.0098、1.3358、2.6386、3.4746、0.0774。针对问题四，关于参数标定的精度与稳定性问题，利用CT系统的结构以及成像原理可知,射线对物体进行平行束扫描获得了180个离散数据，由于CT系统的探测器之间具有一定的间隔，所以物体的边界可能有一部分无法被收集到信息，从而造成投影长度缩小的情况，导致所求探测器的单元间隔的增大，并且造成旋转中心与射线旋转角的误差。通过分析椭圆投影到探测器上长度所对应的误差分析，得到投影长度越长，模板下图形面积越小，所对应的误差就越小，即可得出一个精度较高的新模板。关键词：最小二乘法；滤波反射法；衰变系数；傅里叶变换；一、问题的重述1.1问题的背景：计算机层析成像(ComputedTomography，CT)是通过对物体进行不同角度的射线投影测量而获取物体横截面信息的成像技术，涉及多个学科领域。CT技术不但给诊断医学带来革命性的影响．还成功地应用于无损检测、产品反求和材料组织分析等工业领域，其实质是由扫描所得到的投影数据反求出成像平面上每个点的衰减系数值。有一种典型的二维CT系统，平行入射的X射线垂直于探测器平面，每个探测器单元看成一个接收点，且等距排列。X射线的发射器和探测器相对位置固定不变，整个发射-接收系统绕某固定的旋转中心逆时针旋转180次。对每一个X射线方向，在具有512个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信息。1.2问题的目的：CT系统安装时往往存在误差，从而影响成像质量，因此需要对安装好的CT系统进行参数标定，即借助于已知结构的样品（称为模板）标定CT系统的参数，并据此对未知结构的样品进行成像。1.3所要解决的问题：问题一：在正方形托盘上放置两个均匀固体介质组成的标定模板，模板的几何信息已经给定，相应的数据文件见附件1，其中每一点的数值反映了该点的吸收强度，这里称为“吸收率”。对应于该模板的接收信息见附件2。请根据这一模板及其接收信息，确定CT系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该CT系统使用的X射线的180个方向。问题二：附件3是利用上述CT系统得到的某未知介质的接收信息。利用(1)中得到的标定参数，确定该未知介质在正方形托盘中的位置、几何形状和吸收率等信息。另外，请具体给出所给的10个位置处的吸收率，相应的数据文件见附件4。问题三：附件5是利用上述CT系统得到的另一个未知介质的接收信息。利用(1)中得到的标定参数，给出该未知介质的相关信息。另外，请具体给出所给出的10个位置处的吸收率。问题四：分析(1)中参数标定的精度和稳定性。在此基础上自行设计新模板、建立对应的标定模型，以改进标定精度和稳定性，并说明理由。二、问题的分析本文主要研究CT系统两方面的问题，一是利用已知结构的样品和它的投影信息得到CT系统的参数，二是用已知的标定CT系统的参数来对未知结构的样品进行成像。2.1问题一的分析：该问题要求利用附件1和附件2的信息来确定出CT系统的旋转中心、相邻探测单元之间的距离以及射线每次旋转的方向。针对这个问题，考虑到圆的图形比较规则且投影也是规则的，首先对圆单独进行考虑，对于任何一个方向的射线，圆对应的投影长度都等于圆的直径长，分析可知圆的投影长度与接收到该圆信息的探测器单元的个数有关，由于每个探测器单元是等间距，由此计算出相邻探测器单元间的距离。附件2中180列数据即是对应的每个射线方向的探测器上所接受到的信息。因为附件2中的每一列数据表示的都表示每次转角所对应的投影长度，用几何方法求解出的函数的那些突变值会对函数造成较大的影响。理论上用附件2的数据拟合和几何计算出的射线转角与投影长度的函数关系应该是吻合的，但由于去除特殊点所拟合的函数也具有一定的误差，用数据拟合的函数会与其有一定的差值，为了使数据所拟合的函数与几何方法拟合的函数更加接近，可以最小二乘法进行拟合。2.2问题二的分析：该问题实际上是CT图像重建问题，即利用CT系统扫描得到的投影数据解出投影像平面上各像素点的衰减系数，然后重建CT图像。由中心切片定理知道密度函数在某一方向上的投影函数的一维傅立叶变换函数是原密度函数的二维傅立叶变换函数平面上沿同一方向且过原点的直线上的值。首先找到投影坐标与原图形坐标之间的转化关系，然后求出某一角度下，在投影坐标系下射线方向上衰减系数的线积分，易知投影方向不同或射线离旋转中心的距离不同，射线穿过模板的厚度则不同，则衰减系数的线积分也会不同，所以衰减系数的线积分是与射线的角度以及投影射线到旋转中心的距离有关的一个函数。然后对其进行一维傅里叶变换，并利用坐标变换，转化为在原图像坐标系下的衰减系数的线积分函数，再对其进行二维的傅里叶变换，并进行二维傅里叶逆变换，然后经过滤波函数得到滤波后的投影像素点的衰减系数的一个分布函数，该分布函数对应的是经过滤波后的射线投影在度范围的衰减系数的叠加。但实际在CT扫描时，得到是180个离散投影数据，而不是连续的角度变化，于是将衰减系数分布函数离散化，将180个角度得到的投影数据叠加，得到图像的像素点，从而重建出图像。2.3问题三的分析：问题三的求解过程与问题二的相同，都是利用滤波反射模型进行求解，从重建出图像和10个位置的吸收率。2.4问题四的分析：要分析参数的精度与稳定性，首先找到造成参数误差的原因，即分别对问题一中所求的CT系统的旋转中心的位置、探测器单元间隔以及射线旋转方向进行误差分析。要进行误差分析，可以从CT系统的结构以及成像原理入手，由于探测器之间具有一定间隔以及投影得到的是离散数据，分析可得，当射线进行平行束扫描时，被投影物体的两侧的数据可能有一部分无法被接收器接受，从而导致投影的长度并不会等于对应原图像的长度，而是会比原像小一点，从而造成参数误差，再从这个误差来源入手，分析问题一的参数的精度问题。三、模型的假设1.假设射线不会因外界其它因数的干扰而衰落；2.假设每个探测器单元的长度都是相等的；3.假设每束射线都是一样的；4.射线的初始位置的方向是朝左上的；四、符号的说明符号符号说明圆的直径的大小第次旋转后，射线射向椭圆在探测器上对应非“0”值的探测器单元的个数第次旋转后，射线射向圆在探测器上对应非“0”值的探测器单元的个数每个探测器的长度椭圆和圆的投影总长度五、模型的建立与求解5.1对于问题一的求解5.1.1问题一的解题步骤：(1)根据附件2模板接受的球的信息来确定相邻两探测器单元之间的距离；(2)随射线旋转角度的变化，利用几何知识列出探测器上由椭圆和圆吸收射线所导致的非“0”值探测器单元的个数与射线方向的函数关系式。再用附件2中每一列非“0”值的探测器单元的个数，用最小二乘法来拟合出射线每个角度与对应的非“0”值探测器单元的个数的函数再与几何求解出的函数进行比较，从而确定射线方向与探测器上由椭圆和圆吸收射线所导致的非“0”值探测器单元的个数的函数关系式；(3)利用最小二乘法求解旋转中心的位置。5.1.2模型的建立与求解5.1.2.1射线入射方向模型（一）相邻两探测器单元之间距离的确定附件1中的数据中只有“0”和“1”的两种值，将其导入到MATLAB中可以发现，“1”值对应的为椭圆和圆的投影数据，“0”值则代表正方形托盘的其他的地方的投影数据，由此可以判断出问题给出的“吸收率”是指固体介质的对射线的吸收率，并且“吸收率”由于射线是平行照射到圆上，因此无论射线从任何方向照射圆，对应角度的探测器上的非“0”值的探测器的个数都是一个定值并且这些探测器单元的总长度都应该为球的直径的长度。附件2中的每列数据是对应每次旋转后的非“0”探测器单元的个数，理论上对于每列圆所对应于非“0”探测器单元的个数应该相等。但是附件2中每一列圆对应于非“0”探测器的个数不相等，导致这种现象的原因可能是有些时候探测器感应出现问题、射线在传播过程中发生了折射等因素导致的，如果用非“0”探测器的个数用最大或者最小值来代替会出现取的最大或最小值是由意外情况导致的，因此用平均值来进行计算得到每个探测器单元的长度为:其中为圆在探测器上投影的列数。（二）射线的180个方向的确定1.几何方法确定180个方向与投影总长度的关系根据图一所示的附件2中射线在探测器上的投影图（竖直方向的白色长度集合表示射线以不同入射角照射椭圆或圆投射在探测器上接受到非“0”值的所有探测器单元的总长度集合），为了更好的表述，定义椭圆投影长度为椭圆在某个投影方向上所对应的图一中的一段竖直的白色长度，圆投影长度为圆在某个投影方向上所对应的图一中所对应的一段竖直的白色长度。图一附件2中射线在探测器上的投影图以椭圆中心为原点，水平向圆心的方向为x轴的正方向，竖直向上的方向为y轴的正方向，建立直角坐标系，初始位置的射线入射方向向左上，射线的方向与y轴正方向的夹角为角。由初始位置（图像最左边）到A点椭圆的投影长度根据图一可以发现由图像最左边到A点的这段过程中在同一个竖直方向上有两段白色长度（即椭圆与圆的投影长度没有重叠），白色的长度即可认为是椭圆和圆投影长度的和。图二射线只与椭圆相切图当入射的射线只与椭圆相切（如图二所示）时假设入射的射线的方程为：y=x+；当入射的射线与椭圆相切时有：联立上述的两个方程组有：；因为入射的射线与椭圆相切，所以有：；；得到：；从而求解出入射的射线与椭圆相切的两个切点分别为：；；则这两条切线所对应的方程分别为：；；则这两条切线（椭圆投影长度）的距离为：（其中）；代入得：；根据两条切线的距离再结合附件2中给出的第一列数据中投影的长度，来求解出射线的初始入射角度大约为29.25。由A点到C点位置椭圆的投影长度根据图一可以发现竖直方向的两端白色带合为了一条，说明了椭圆和圆的投影长度有重叠，即图二所示，圆只有一部分吸收了射线，另一部分的射线被椭圆吸收了所以导致椭圆和圆的投影长度有重合。图三射线与圆相交于椭圆相切在图三椭圆的上半部分，设有直线和椭圆相切，直线和圆相切，且两直线互相平行，两直线的斜率分别，两直线与y轴正方向的夹角为。直线和直线的方程分别为：y=x+；y=x+；根据直线和椭圆相切得出下列方程：联立以上两个式子得：根据直线和圆相切得出下列方程：联立以上两个式子得：根据（3）和（4）式子可得：；；解得：或所以直线和直线的距离（即圆的投影长度）为：；椭圆的投影长度和①中的求法一致为：（）所以投影的总长度为：。在椭圆的下半部分求解一直线和椭圆相切，另一直线和圆相切，且两直线互相平行，这两直线的斜率和两直线与y轴正方向的夹角的求解和椭圆的上半部分的求解方法一样，解出的直线斜率为或，总的投影长度也是圆和椭圆的投影长度的和。③C点到边界的最右端由于C点到边界的最右端根据图一可以发现由图像最左边到A点的这段过程中在同一个竖直方向上有两段白色长度（即椭圆与圆的投影长度没有重叠），白色的长度即可认为是椭圆和圆投影长度的和为在（）上成立。综上所述：（其中），=8mm,表示180个射线方向角与圆和椭圆的投影总长度的关系)对上述函数先将有异常的点去除，再用傅里叶五级变换进行对其函数进行拟合，得到下列图像：图四傅里叶五级变换后的函数2.用附件2数据确定180个方向与投影总长度的关系对于附件2中所给出的512*180的数据，可以知道每一列即为对应的一个射线转动角度，每一列数据中非“0”值的个数即为圆和椭圆的投影总长度（射线被圆和椭圆吸收后在探测器上射线信号减落的探测器单元的总个数）。为了描述射线的入射角度的变化对圆和椭圆的投影总长度（一列上非“0”数值的个数）的变化关系，可以用每一列上的非“0”值的个数做为纵坐标，以列数（射线的每次转动的方向）做为横坐标，用最小二乘法来拟合出一个最佳的函数来。设表示第列数，表示的是第列数中非“0”值的个数，每一列（射线的每次转动的方向）所对应的（）则构成了含有180个点的点集，并且二者满足关系=……。现要根据附件2的数据来拟合出一个次数不超过次的多项式=，根据上述条件即求解:所解出的函数即为每次转动射线所对应的角度与投影总长度的函数关系式。3.射线方向与投影总长度关系的模型上文分别应用了几何方法求解函数关系和最小二乘法拟合函数的方法来确定的射线的方向与投影总长度之间的关系。这两种方法求解出的函数应该是同一个函数因为都是表示射线的方向与投影总长度之间的函数关系，只是这两种的方法在求解过程中会有一定的误差，可能导致两函数可能在相同的射线的方向上的投影总长度之间的函数值会有所偏差。因为用几何计算出的180个方向与投影长度的关系有一些异常的点，因此可以将那些异常的点去除对那些正常的点进行拟合，能拟合出函数为=；数据是用最小二乘法拟合，拟合出函数为。由于几何计算出的点有些是异常点，需要进行拟合，在拟合过程中会有很多不同的拟合方法，用数据拟合也是如此。因为上述两种方法求解的函数是同一个函数，理论上值应该相等，为了找出一个能更加准确地表示射线方向与投影总长度关系，在所有相同偏转角的情况下两种方法拟合的函数所对应的值的差值的平方应该尽肯能的小。因此可以建立射线方向与投影总长度关系的模型，模型如下所示：用MATALAB对进行模型求解后，得到了一个较优的射线方向与投影总长度的函数关系：图五优化后射线方向与投影总长度的函数再将附件2中每一列非“0”值的个数分别代入上述方程从而解得这180个角度方向如表一（详细的数据见附录）表一射线的180次的方向编号123．．．．．．178179180与y轴正方向的夹角29.066230.570231.5728．．．．．．207.0412208.0438209.04655.1.2.图六旋转图如图六中所示，以Y为坐标原点，射线射向的方向为水平方向的正方向建立直角坐标系。由旋转几何关系可以推断出直线的长度可以用随角度变化的函数来表示其长度：；设Y是一个常量：；所以这两点在一个角度为上的探测器上的的长度：；椭圆中心点在探测器上的投影长度可以通过附件2给出的数据来拟合出一个函数来表示它：；根据长度关系可得：；即：；因为那么可化为：，综上所述，可以建立旋转中心模型，模型如下所示：用MATALAB求解得：；；，即旋转中心坐标为（相对于正方形左下角建立的直角坐标系为参考系的）。5.2对于问题二的求解5.2.1问题二的求解步骤（1）求出投影坐标与原图像坐标之间的一个转化关系；（2）假设射线旋转角连续时，利用滤波反投影法求出经过滤波后的射线投影在度范围的衰减系数的叠加，求出衰减系数分布函数；（3）将衰减系数分布函数进行离散化，利用问题一求解出的180个角度进行累加，得出待重建图像的像素点，重建出CT图像；5.2.2问题二模型的建立与求解5.2.2.1投影坐标与原图像坐标之间的转化图七投影坐标与原图像坐标之间的转化关系图由投影与几何知识可知：；（1）；（2）其中为射线的方向（射线与原图像纵坐标的正方向的夹角）；5.2.2.2基于滤波反投影图像的重建模型（一）假设射线方向角度的变化是连续的以CT系统旋转中心为原点，射线方向为纵坐标（轴），垂直于射线方向的为横坐标（轴）建立直角坐标系，画出射线平行投影曲线，如图图八射线投影假设是CT系统扫描过程中得到的投影数据，，可知在数值上等于射线方向上衰减系数的线积分，由于该积分与射线穿过模板的厚度有关，对于不规则的图形来说，不同的投影方向或射线离旋转中心的距离不同，射线穿过模板的厚度则不同，所以的大小与射线的角度以及投影射线到旋转中心的距离有关。设是基于坐标（对应旋转后的坐标系）时，点处的衰减系数，基于坐标（对应原图形的坐标系）时，点处的衰减系数，则有：；（3）对进行一维傅里叶变换得：；（4）然后（3）代入到（4）得：；(5)结合（1）（2）可得；（6）再将进行二维傅里叶变换，可得：；（7）其中:（8）将（8）代入（7）中可得：；（9）对（7）进行傅里叶逆变换可得：；(10)结合（8）得：；；由对称性又可知：；；有此可得：；该函数为图像的衰减系数分布函数，表示在处，所有经过滤波后的射线投影在度范围的衰减系数的叠加。（二）将衰减系数分布函数进行离散化当假设射线旋转的角度是连续的，由（一）可知衰减系数分布函数，即可求出任意一点处，经过滤波后的射线投影在度范围的衰减系数的叠加，考虑到实际情况，X射线方向对应的角度并不是一个连续的，而是180个离散的值，故将该进行离散化。根据第一问的模型得出CT系统的旋转中心为（40.7543,56.2520），以及射线的180个方向。假设这180个方向对应的角度为，，，…,(见附录一)，则正方形托盘中任意一点的吸收率（也就是射线的衰减系数）为：；即可以求出在原图像坐标系下任一位置的吸收率。将加上未知介质的整个正方形托盘所对应图像的求出，重塑图像的每个像素点即可得重建出CT图像。利用MATLAB求解得出，该未知介质的形状如下图：图九问题二中未知介质的形状及在托盘中位置并求出附件四给出的10个位置处的吸收率见下表二：表二10个位置处的吸收率横坐标1034.543.54548.5505665.579.598.5纵坐标18253375.555.575.576.5371843.5吸收率0.04420.98391.26820.97551.00990.98910.00330.97310.04760.03075.3对于问题三的求解该题与问题二研究的问题一样，仍然是通过投影信息反求出图形的原像。首先利用问题二中的滤波反投影法得到任意一位置经过滤波后的射线投影在度范围的衰减系数的叠加，即问题三中模板的图像的衰减系数分布函数，再将其离散化，利用180个离散投影数据反投影求出问题三的原图形为：图十

问题三中未知介质的原图形同理得到附件四给出10个位置对应的吸收系数，如表三所示：表三附件四给出10个位置对应的吸收系数：横坐标1034.543.54548.5505665.579.598.5纵坐标18253375.555.575.576.5371843.5吸收率1.48130.06032.12980.00363.39000.00981.33582.63863.47460.07745.4对于问题四的求解5.4.1探测器单元之间的距离的误差分析图十一射线与圆两侧相切图十二射线与圆一侧相切图十三射线与圆两侧相交其中：为圆的直径，由题目可知；表示探测到信息的探测器的个数；表示探测器之间的距离；问题一求解探测器之间的距离的时候，假设的是图十这种情况即射线与圆的两侧相切的情况，此时：；但探测器之间是有间距的，并且得到的是512个离散的数据，即在实际旋转的过程中会出现图十一和图十二这种情况，此时对于圆的一侧或者两侧，探测器并没有接受到信息，导致计算的时候投影长度变大了，因此存在误差。对于图二和图三这两种情况，分析易知投影长度最大误差的时候也只是与圆的实际直径相差一个探测器的长度，可知：；化简的：；易知当：；即：；此时误差为：；所以探测器单元之间距离的误差范围为（0，）。5.4.2射线180个投影方向角的误差分析1.对于射线180个投影方向角来说，误差来源与探测器单元之间的距离误差来源一样，不同的地方是，当椭圆与圆投影没有重叠的部分时，其投影长度与实际长度相差值在范围内，而当椭圆与圆投影有重叠的部分时，其投影长度与实际长度相差值在范围内。结合第一问求射线180个投影方向角的模型易知：当和以及不存在时即射线与轴垂直时，椭圆与圆投影无重叠的地方，此时：；当时，，椭圆与圆投影有重叠的地方，此时：；2.相对误差的分析模型：椭圆在探测器上的投影长度角度的变化而变化，具体的关系式为所以椭圆长度d角度变化的曲线如下图图十四椭圆长度d角度变化的曲线图而通过探测器个数求解出对应椭圆的投影长度是有误差，如下图图十五投影长度与探测器单元个数的关系图所以相对误差=；即可得出相对误差随角度的变化关系曲线图十六相对误差波动随角度的变化关系图可以得出椭圆投影到探测器上的长度越长，相对误差就越小，长度越短，误差越大；所以为了提高精度，那么模板在每个角度下投影的长度尽量大，且实际模板下的图像的面积足够小，根据以上两个关系，就可以得出以下的新模板：图十七新模板成像图六、模型的检验滤波反射模型的检验：附件2是附件1对应图像的CT图像，因此通过CT图像进行还原得到一个吸收系数的矩阵，通过判别用本模型求解出来的结果与实际结果对应整个像素点值之间的差值情况，来辨别本模型的优略。利用滤波反射模型对附件2进行还原，还原出来的结果与实际数字图片显示后结果如下：通过观察可以看出实际图像与本模型求解出的图像非常接近，为了进一步刻画二者图像的相似情况，引入一个区别系数P：其中：；是还原图像矩阵第行第列对应像素点的吸收系数大小；计算得出：116.5686即两张图片对应每个像素点之间吸收率的误差仅为0.177%，由此可证明本模型可靠性较高。七、模型的评价7.1模型的优点：1.滤波反射模型能够将未知的CT参数进行成像，并且成像的效果显著，具有很强的推广性。2.射线入射模型能够通过利用一些CT参数数据就能确定出各个旋转方向的转角，结果准确性较高。7.2模型的缺点：射线入射模型是通过最小二乘拟合出来的转向的角度与探测器上非“0”值的探测单元的个数的函数关系，是通过拟合确定的，还是具有一定的误差性。八、模型的推广与应用滤波反射模型对于处理给定对于已知结构样品的所标定CT系统的参数来对未知结构样品来进行成像的效果显著，可以对不同的已知结构的样品的参数，并根据其参数进行成像，且还原精度较高，可以应用到生物科学或者考古等领域对样品进行断层成像。参考文献[1]赵静，但琦.数学建模与数学实验[M].北京:高等教育出版社,2014.[2]盖云英，包革军.复变函数[M].哈尔滨:哈尔滨工业大学出版社,1998.附录附录一：求解180个角度MATLAB程序clc,clearc=0;cc=[];%存储shuju=xlsread('附件二.xls');juli_xsd=0.2775;%接受信息点之间的距离fori=1:180shuju_1=shuju(:,i)';L=length(shuju_1);t_1=1:L;t_2=1:0.5:L;shuju_1_chazhi=interp1(t_1,shuju_1,t_2,'linear');forj=1:length(t_2)ifshuju_1_chazhi(j)~=0c=c+1;endendc=juli_xsd*c/2;cc=[cc,c];c=0;endliang=200;pppp=[];forp=0:0.0001:pi/4;k=linspace(0,pi,180);a_1=(2.*sqrt(225.*(tan(k+p+pi/2)).^2+1600))./(sqrt((tan(k+pi/2+p)).^2+1));sum(abs(a_1-cc));plot(p,(sum(abs(a_1-cc))),'r');holdonpppp=[pppp,sum(abs(a_1-cc))];end[a,b]=min(pppp);t=0:0.0001:pi/4;jd=t(b);%%%%%第一个角jiaodu_cc=[];jiaodu_cc(1)=jd;fori=2:180cc_x=cc(i);q=jd+0.0175*(i-1)-0.0175/2:0.0001:jd+0.0175*(i-1)+0.0175/2;jieguo=(abs((2.*sqrt(225.*(tan(jd+pi/2+q)).^2+1600))./(sqrt((tan(pi/2+jd+q).^2+1))-cc_x)));[x_1,x_2]=min(jieguo);jiaodu_cc=[jiaodu_cc,q(x_2)];endxlabel('初始旋转角度');ylabel('180个角度对应高度想相减之和');title('旋转角度和图像拟合程度曲线变化图');jiaodu_cc=jiaodu_cc*180/pi;disp('对应的180个角度');disp(jiaodu_cc);附件二：求解旋转中心程序%%%%%%%本程序采用暴力求解精度较高，但方法求解时间较长%%%%%%clc,clearcsjd=29.06615;%初始角度shuju=xlsread('附件2.xls');%juli_xsd=0.2775;%接受信息点之间的距离cc=[];%用来存储fori=1:180%求椭圆中心投影高度变化曲线aa1=shuju(:,i);cd=0;%长度forjj=1:512ifaa1(jj)~=0cd=cd+1;endendgd=0;forjjj=1:512ifaa1(513-jjj)==0gd=gd+1;endifaa1(513-jjj)==0&&aa1(512-jjj)~=0break;endendcc=[cc,(gd+cd/2)*0.2775];endbb=linspace(0,pi,180);c=20000;x=0;y=0;forii=39:0.0001:41%X的值forjj=55:0.0001:57%Y的值L=cos(atan((50-jj)/(50-ii))-(29.06615+180)*pi/180-bb)*sqrt((50-jj)^2+(50-ii)^2);ifvar(cc-L)<cc=var(cc-L);x=ii;y=jj;ddddd=abs(cc-L);%plot3(ii,jj,var(cc-L),'ro');holdonendendenddisp('旋转中心横坐标');disp(x);disp('旋转中心纵坐标');disp(y);disp('旋转中心点在投影处据向量原点的固定距离');disp(mean(ddddd));plot(ddddd);xlabel('旋转角度');axis([01804080]);ylabel('旋转中心投影的高度');title('旋转角度与旋转中心投影的关系图像');附录三：求解问题2程序clc,clear;shuju=xlsread('qiujie.xls');%%%将投影数据导入jiaodu=(xlsread('180个角度值.xls'));%%%将对应角度数据导入x1=40.7543;%旋转中心横坐标y1=56.2520;%旋转中心纵坐标%imshow(shuju);显示图片h_xzd=71.1127;%旋转中心点在投影处据向量原点的固定距离L_xsd=0.2775;jcd=floor((h_xzd*2-512*L_xsd)/L_xsd);%将原有矩阵加长度shuju1=[zeros(jcd,180);shuju];jiaodu=jiaodu+180;%%%%%对应角度%%与X轴jiaodu=round(jiaodu);[x1_1,y1_1]=find(jiaodu==360);%做反向fxjz=zeros(512,180);foriii=1:180;forjjj=1:512fxjz(jjj,iii)=shuju(513-jjj,iii);endendzh=[shuju,fxjz];shuju=zh(:,y1_1:y1_1+179);%imshow(shuju1);%显示图片L=360;P_1=zeros(L,L);%用于存放重构后的图像H_2_H=size(shuju,1);H_3_H=zeros((H_2_H*2-1),1);fori=0:H_2_H-1%%创造滤波ifi==0H_3_H(H_2_H-i)=1/4;elseifrem(i,2)==0H_3_H(H_2_H-i)=0;H_3_H(H_2_H+i)=0;elseH_3_H(H_2_H-i)=-1/(i*pi)^2;H_3_H(H_2_H+i)=-1/(i*pi)^2;endendN=180;%%%换成极坐标表示x=zeros(H_2_H,N);fori=1:N%%%%%%%%%%滤波%%%%%%%%%%%s=shuju(:,i);xx=conv(s',H_3_H');x(:,i)=xx(H_2_H:2*H_2_H-1);end%%%%%%%%%%%简单反投影算法%%%%%%%%MOBAN=zeros(L,L);fori=1:L%%x轴forj=1:L%%%y轴fork=1:180%%%旋转角度THE=(k)/180*pi;t=(j-L/2-0.5)*cos(THE)+(L/2+0.5-i)*sin(THE)+(H_2_H+1)/2;t1=floor(t);t2=floor(t+1);MOBAN(i,j)=MOBAN(i,j)+(t2-t)*x(t1,k)+(t-t1)*x(t2,k);endendendMOBAN=pi/N*MOBAN;xianshi=[];forii=1:360forjj=1:360xianshi(ii,jj)=MOBAN(361-ii,361-jj);endendfigure;imshow(xianshi,[]);title('滤波反投影法(未移动)');%%%图像转换%%坐标原点x1_1_1=round(abs((x1-50)/L_xsd));y1_1_1=round(abs((y1-50)/L_xsd));bzmb=MOBAN(:,1:360-x1_1_1);bzmb=bzmb(1:360-y1_1_1,:);[k1,k2]=size(bzmb);jia_1=zeros(360-k1,k2);jia_2=zeros(360,360-k2);bzmb=[jia_1;bzmb];bzmb=[jia_2,bzmb];bzmb=2*bzmb;%figure;imshow(bzmb,[]);title('');%%%%%求对应10个坐标的值%%%以及对应像素点值大小zb_10ge=xlsread('附件4.xls');%%%将投影数据导入zb_10ge=round(zb_10ge./L_xsd);%转换成像素点个数cc_xs=[];%用来存储吸收系数fori=1:10meige=zb_10ge(i,:);cc_xs=[cc_xs;bzmb(360-meige(2),meige(1))];enddisp('对应10个坐标的吸收强度');disp(abs(cc_xs));V=bzmb;[X,Y]=meshgrid(1:360);[Xq,Yq]=meshgrid(1:0.2:360);Vq=interp2(X,Y,V,Xq,Yq);jieguo=[];fori=1:256forj=1:256jieguo(i,j)=Vq(7*i,7*j);endendjieguo=abs(jieguo);xianshi1=[];forii=1:256forjj=1:256xianshi1(ii,jj)=jieguo(257-ii,257-jj);endendfigure;imshow(xianshi1,[]);title('滤波反投影最终图片')%%%%%答案结果disp(‘xianshi1矩阵是对应像素点的吸收强度’);附录四：求解第三问滤波反射法模型的MTALB程序clc,clear;shuju=xlsread('附件5.xls');%%%将投影数据导入jiaodu=(xlsread('180个角度值.xls'));%%%将对应角度数据导入x1=40.7543;%旋转中心横坐标y1=56.2530;%旋转中心纵坐标h_xzd=71.1127;%旋转中心点在投影处据向量原点的固定距离L_xsd=0.2775;jcd=floor((h_xzd*2-512*L_xsd)/L_xsd);%将原有矩阵加长度shuju1=[zeros(jcd,180);shuju];jiaodu=jiaodu+180;%%%%%对应角度%%与X轴jiaodu=round(jiaodu);[x1_1,y1_1]=find(jiaodu==360);%做反向fxjz=zeros(512,180);foriii=1:180;forjjj=1:512fxjz(jjj,iii)=shuju(513-jjj,iii);endendzh=[shuju,fxjz];shuju=zh(:,y1_1:y1_1+179);HHHHHHH=360;P_1=zeros(HHHHHHH,HHHHHHH);%用于存放重构后的图像H_1_1=size(shuju,1);h_1_1=zeros((H_1_1*2-1),1);foriii=0:H_1_1-1%%%创建滤波ifiii==0h_1_1(H_1_1-iii)=1/4;elseifrem(iii,2)==0h_1_1(H_1_1-iii)=0;h_1_1(H_1_1+iii)=0;elseh_1_1(H_1_1-iii)=-1/(iii*pi)^2;h_1_1(H_1_1+iii)=-1/(iii*pi)^2;endendN=180;%%%换成极坐标表示x=zeros(H_1_1,N);foriii=1:N%%%%%%卷积求%%%s=shuju(:,iii);xx=conv(s',h_1_1');x(:,iii)=xx(H_1_1:2*H_1_1-1);end%%%反投影到矩阵中去X_3_3=zeros(HHHHHHH,HHHHHHH);foriii=1:HHHHHHH%%x轴forj=1:HHHHHHH%%%y轴fork=1:180%%%旋转角度THE=(k)/180*pi;t=(j-HHHHHHH/2-0.5)*cos(THE)+(HHHHHHH/2+0.5-iii)*sin(THE)+(H_1_1+1)/2;t1=floor(t);t2=floor(t+1);X_3_3(iii,j)=X_3_3(iii,j)+(t2-t)*x(t1,k)+(t-t1)*x(t2,k);endendendX_3_3=pi/N*X_3_3;x1_1_1=round(abs((x1-50)/L_xsd));y1_1_1=round(abs((y1-50)/L_xsd));bzmb=X_3_3(:,1:360-x1_1_1);bzmb=bzmb(1:360-y1_1_1,:);[k1,k2]=size(bzmb);jia_1=zeros(360-k1,k2);jia_2=zeros(360,360-k2);bzmb=[jia_1;bzmb];bzmb=[jia_2,bzmb];bzmb=2*bzmb;%%%修正系数xianshi=[];forii=1:360forjj=1:360xianshi(ii,jj)=bzmb(361-ii,361-jj);endendfigure;imshow(xianshi,[]);title('滤波反射法形成的图');%%%%%求对应10个坐标的值zb_10ge=xlsread('附件4.xls');%%%将投影数据导入zb_10ge=round(zb_10ge./L_xsd);%转换成像素点个数cc_xs=[];%用来存储吸收系数foriii=1:10meige=zb_10ge(iii,:);cc_xs=[cc_xs;bzmb(360-meige(2),meige(1))];enddisp('对应10个坐标点的吸收率');disp(abs(cc_xs));V=bzmb;[X,Y]=meshgrid(1:360);[Xq,Yq]=meshgrid(1:0.2:360);Vq=interp2(X,Y,V,Xq,Yq);jieguo=[];fori=1:256forj=1:256jieguo(i,j)=Vq(7*i,7*j);endendjieguo=abs(jieguo);xianshi1=[];forii=1:256forjj=1:256xianshi1(ii,jj)=jieguo(257-ii,257-jj);endendimshow(xianshi1,[]);title('滤波反投影第三问图片')%%%%%答案结果disp(‘xianshi1矩阵是对应像素点的吸收强度’);附录五：滤波反射模型检验MATLAB程序clc,clear;shuju=xlsread('附件2.xls');%%%将投影数据导入jiaodu=(xlsread('180个角度值.xls'));%%%将对应角度数据导入x1=40.7543;%旋转中心横坐标y1=56.2530;%旋转中心纵坐标h_xzd=71.1127;%旋转中心点在投影处据向量原点的固定距离L_xsd=0.2775;jcd=floor((h_xzd*2-512*L_xsd)/L_xsd);%将原有矩阵加长度shuju1=[zeros(jcd,180);shuju];jiaodu=jiaodu+180;%%%%%对应角度%%与X轴jiaodu=round(jiaodu);[x1_1,y1_1]=find(jiaodu==360);%做反向fxjz=zeros(512,180);foriii=1:180;forjjj=1:512fxjz(jjj,iii)=shuju(513-jjj,iii);endendzh=[shuju,fxjz];shuju=zh(:,y1_1:y1_1+179);HHHHHHH=360;P_1=zeros(HHHHHHH,HHHHHHH);%用于存放重构后的图像H_1_1=size(shuju,1);h_1_1=zeros((H_1_1*2-1),1);foriii=0:H_1_1-1%%%创建滤波ifiii==0h_1_1(H_1_1-iii)=1/4;elseifrem(iii,2)==0h_1_1(H_1_1-iii)=0;h_1_1(H_1_1+iii)=0;elseh_1_1(H_1_1-iii)=-1/(iii*pi)^2;h_1_1(H_1_1+iii)=-1/(iii*pi)^2;endend%h_1_1=h_1_1+0.001;N=180;%%%换成极坐标表示x=zeros(H_1_1,N);foriii=1:N%%%%%%卷积求%%%s=shuju(:,iii);xx=conv(s',h_1_1');x(:,iii)=xx(H_1_1:2*H_1_1-1);end%%%反投影到矩阵中去X_3_3=zeros(HHHHHHH,HHHHHHH);foriii=1:HHHHHHH%%x轴forj=1:HHHHHHH%%%y轴fork=1:180%%%旋转角度THE=(k)/180*pi;t=(j-HHHHHHH/2-0.5)*cos(THE)+(HHHHHHH/2+0.5-iii)*sin(THE)+(H_1_1+1)/2;t1=floor(t);t2=floor(t+1);X_3_3(iii,j)=X_3_3(iii,j)+(t2-t)*x(t1,k)+(t-t1)*x(t2,k);endendendX_3_3=pi/N*X_3_3;%figure;imshow(X_3_3,[]);title('滤波反投影法');%%%图像转换%%坐标原点x1_1_1=round(abs((x1-50)/L_xsd));y1_1_1=round(abs((y1-50)/L_xsd));bzmb=X_3_3(:,1:360-x1_1_1);bzmb=bzmb(1:360-y1_1_1,:);%figure;imshow(bzmb,[]);[k1,k2]=size(bzmb);jia_1=zeros(360-k1,k2);jia_2=zeros(360,360-k2);bzmb=[jia_1;bzmb];bzmb=[jia_2,bzmb];bzmb=2*bzmb;fori=1:360forj=1:360if(i-342)^2+(j-180)^2<=225bzmb(j,i)=1;endendend%figure;imshow(bzmb,[]);title('滤波反射法形成的图');%%%%%求对应10个坐标的值zb_10ge=xlsread('附件4.xls');%%%将投影数据导入zb_10ge=round(zb_10ge./L_xsd);%转换成像素点个数cc_xs=[];%用来存储吸收系数foriii=1:10meige=zb_10ge(iii,:);cc_xs=[cc_xs;bzmb(360-meige(2),meige(1))];end%disp('对应10个坐标点的吸收率');disp(cc_xs);V=bzmb;[X,Y]=meshgrid(1:360);[Xq,Yq]=meshgrid(1:0.2:360);Vq=interp2(X,Y,V,Xq,Yq);jieguo=[];fori=1:256forj=1:256jieguo(i,j)=Vq(7*i,7*j);endendimshow(jieguo,[]);title('附件2还原图像')%%%%%答案结果jieguo3=xlsread('附件1