2017年数学建模省一等奖高月_陈磊_饶运良-基于滤波反投影重建法的-CT-系统参数标定及成像的研究

1基于滤波反投影重建法的CT系统参数标定及成像的研究摘要本文是根据滤波投影和反投影重建原理，基于平行PI线段的平行射束的反投影滤波算法(parallelbeamBPF,PBBPF)，利用shepp-Logan头模型进行了滤波投影和卷积反投影，得出CT系统成像模型，使其能够通过接收数据进行参数标定。针对问题一：本文首先运用统计学方法对附件一的吸收率文件数据进行统计分析并做出三维立体图和MATLAB进行数据模拟生成图像，分析出该数据对应的是密度均匀的一个模板。再对附件二中该模板的接收信息进行数据分析，之后用MATLAB进行图像模拟，然后根据几何关系，计算得出CT系统的旋转中心、探测器之间的距离，以及CT系统使用X射线的180个方向，并用shepp-Logan头模型进行图像模拟验证投影的建立和重建。针对问题二：首先我们运用问题一的模型根据题中所给的未知介质的接收信息，在问题一得到的标定参数的基础上，通过几何算法和投影反重建原理，在MATLAB中用radon（）函数还原出该未知介质的投影形状，进而得出该未知介质在正方形托盘中的位置、几何信息和吸收率等数据。针对问题三：问题三与问题二问法相同，只是把接收信息的数据该变了，同样利用问题一得出的标定参数，利用MATLAB中radon()函数投影反重建出未知介质的投影图形，然后对照问题二的解题思路求解即可。针对问题四：我们对问题一中所给的模板介质是否均匀，旋转中心和起始角度的标定是否影响标定精度和稳定性，我们设计了一个新的模板，该模板由更好的精度和稳定性。关键词：CT系统成像；参数标定；滤波投影；PBBPF2一、问题重述与提出1.1阐述背景CT(ComputedTomography)可以在不破坏样品的情况下，利用样品对射线能量的吸收特性对生物组织和工程材料的样品进行断层成像，由此获取样品内部的结构信息。一种典型的二维CT系统如图1所示，平行入射的X射线垂直于探测器平面，每个探测器单元看成一个接收点，且等距排列。X射线的发射器和探测器相对位置固定不变，整个发射-接收系统绕某固定的旋转中心逆时针旋转180次。对每一个X射线方向，在具有512个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信息。CT系统安装时往往存在误差，从而影响成像质量，因此需要对安装好的CT系统进行参数标定，即借助于已知结构的样品（称为模板）标定CT系统的参数，并据此对未知结构的样品进行成像。1.2提出问题问题一：在正方形托盘上放置两个均匀固体介质组成的标定模板，模板的几何信息如图2所示，相应的数据文件见附件1，其中每一点的数值反映了该点的吸收强度，这里称为“吸收率”。对应于该模板的接收信息见附件2。根据这一模板及其接收信息，确定CT系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该CT系统使用的X射线的180个方向。问题二：附件3是利用上述CT系统得到的某未知介质的接收信息。利用(1)中得到的标定参数，确定该未知介质在正方形托盘中的位置、几何形状和吸收率等信息。另外，具体给出图3所给的10个位置处的吸收率，相应的数据文件见附件4。问题三：附件5是利用上述CT系统得到的另一个未知介质的接收信息。利用(1)中得到的标定参数，给出该未知介质的相关信息。另外，请具体给出图3所给的10个位置处的吸收率。问题四：分析(1)中参数标定的精度和稳定性。在此基础上自行设计新模板、建立对应的标定模型，以改进标定精度和稳定性，并说明理由。以上所有题目中的所有数值结果均保留4位小数。同时提供(2)和(3)重建得到的介质吸收率的数据文件（大小为256×256，格式同附件1，文件名分别为problem2.xls和problem3.xls）。3二、问题分析2.1问题一的分析由题意可知，本问的目的就是根据标定模板校准该CT系统的参数，并且求得其他未知物体的相关CT扫描信息。根据所测的标定模板的相关信息，如吸收强度，利用MATLAB进行数据分析，做iradon()函数变换，得出相关CT系统扫描图像，建立合适的坐标系，正确描述CT系统的各种参数，根据测量数据计算该CT系统的参数值，即给出CT系统的标定：CT系统的旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该CT系统使用的X射线的180个方向。2.2问题二的分析根据CT系统得到的某未知介质的接收信息,在问题一中得到的标定参数基础上，先用MATLAB进行数据分析，图像拟合，确定该未知介质在正方形托盘中的位置、几何形状。再结合吸收强度函数，可得出该物体的吸收率和已知点的吸收率。2.3问题三的分析需要测算另一未知物体的相关信息。分析思路同问题二。2.4问题四的分析对问题一所求的CT系统的标定参数进行精度和稳定性分析，采用误差分析，以设计新的标定模版，改进以提高精度和稳定性。三、模型假设1．假设每个探测器单元一个模型质点，不考虑其大小形状;2．假设模板的接收信息不受其他客观因素影响，只与其形状和内部结构影响;3．假设平行光源发出的所有光强度一样，X射线密度均匀，并且不受外界因素干扰;4．假设模板为一个平面图形，但能正常吸收X射线的能量，研究中不考虑所有模板的立体性;5．假设实验所收集的数据能客观反映实际情况。四、符号说明表4-1符号说明xo旋转中心的横坐标4yo旋转中心的纵坐标d探测器单元之间的距离d1小圆的直径θ1光源的初始位置∆θ每次旋转的角度x′变换后的横坐标y′变换后的纵坐标实际吸收率数据吸收率增益系数I光通量P0投影值五、模型的建立与求解5.1通用模型的建立根据物理学中磁通量的定义，我们在这里定义光通量，其定义是单位长度通过的光线数量，这里用I表示。定义：miKI，其中K为常数，i为平行光线数量(我们假定光线是按数量来统计的，例如：“n条光线”就对应i=n，m为与平行光线垂直的线段长度)。注：本文用能接受模板信息的探测器的数量大小来表征I的大小。我们再定义p0为一条光线重合地通过一条有限长度的线段时对应的投影值，并规定：dlpx00，其中是线段的线密度，x为线段的长度。5.2问题一5.2.1首先我们要了解投影的原理实际的投影数据及真实的投影数据是由检测器得到的，在计算机模拟时，则取自仿真模型。仿真模型常由一系列大小不同、位置各异，具有给定密度的椭圆组成，以模拟介质的图形[1]。5.2.2根据投影原理建立模型一建立并求解为了方便分析，先来讨论椭圆中心在直角坐标原点，长轴与y轴重合时，穿过椭圆5的某一射线投影p0的求法。图1求穿过椭圆状断面的投影图由图1可以看出：投影p点的值为QPdlpQP式中是椭圆模型的密度值，椭圆的方程为：12222byax，(1)射线QP的法线式为：dysinxcos，(2)各变量的定义见图4，由(1)和(2)解得射线和椭圆焦点P、Q的坐标：22222,1sincosxrdrabda22222,1cossyrdrabindb式中22222sincosrba因此：212212)()(xQPyyx222rab2rd则：2220rd-r2abQPp一般情况下，若与椭圆中心在（G，H），长轴与y轴的交角为（图2）则：6图2任意位置、任意倾角椭圆222rd-r2abQP式中sincosdHG，）-（sinb）-（cosar22222此时，射线投影值为2220rd-r2abp）-（sinb）-（cosa）sincos（-）-（sinb）-（cosa2ab222222222HG结论一：由上述分析可知，投影0p的数值大小不仅和椭圆的密度有关，而且还和光线与椭圆的交线线段长度也有关，而且0p的大小与密度成正比，与交线的线段长成正相关。5.2.3对问题一的求解（1）对附件一的分析首先对附件一的数据进行MATLAB分析，得出以下图形（图4），与题中给的模板示意图（图3）比较，发现其形状一致，我们先假设附件一的数据按比例反映了模板的几何信息。7图3模板示意图（单位：mm）图4附件一数据拟合图然后根据附件二的数据分析，并结合模板的几何信息，得出该模板再正方形托盘中的位置，探测器单元之间的距离，以及CT系统使用的X射线的180个方向（即解决了问题一）。（2）对附件二进行MATLAB数据拟合在MATLAB环境下进行数据拟合得到模板的接收数据示意图（图5），其与模板对照如下：图5附件二的图像拟合生成的分布图（左）与模板图（右）的对照对比分析不难看出长轴对应的即为左图x=-90所对应的最长列，小球的直径为左图x=0所对应的最短列。由此我们进行一下计算。步骤一：建立坐标系先在正方形圆盘中以其中心为原点，以椭圆长轴方向为y轴，以短轴方向为x，建立平面直角坐标系，并规定探测器在低部时（光源在正顶部），从右边起第一个探测器编号为1，往左边的探测器编号数逐一增加，直至编号512。步骤二：拟合数据图并分析找出与几何图形的关系用MATLAB对附件二的数据拟合生成数据分布图，并在图中自动生成合适的坐标系。分析附件二的数据，找出其椭圆光通量I最大的那一列，也就是图6中是横轴为-90时所对应的那一列，计录此时的旋转次数为61，本列中投影值p0最大对应的探测器编8号为235。找出一列的光通量I最小的那一列，也就是图6中是横轴为0时所对应的那一列，计录此时的旋转次数为151，本列中投影值p0最大对应的探测器编号为223。步骤三：计算探测器单元之间的距离同理对于小圆，各探测器在不同角度就收到小圆的数据如图6中的小带所示，由于圆的位置直径不唯一，所以无论CT系统旋转到什么角度，通过小圆的光通量I是一定的，总是和其直径的长度对应，接收到小圆信息的探测器个数是一定的。在对大量列数统计能接受到小圆信息的探测器个数后，发现每一列能接收到的小圆信息的探测器个数是29个（只有极少的几列是28个，这里我们取29个探测器为准），这就说明这29个探测器对应的是小圆直径1d，因此每个探测单元之间的距离为291d。代入数据求解的探测器单元之间的距离mmdd2759.0298291（3）步骤四：计算旋转中心编号为235的探测器与编号为256的距离即为探测器旋转中心在y方向偏离椭圆中心的距离0y;编号为223的探测器也编号为256的距离即为探测器旋转中心在x方向偏离椭圆中心的距离0x，从而计算出其旋转中心（00,yx）。代入数据求解，dx)223256(0（4）dy)235256(0（5）式中d为探测器单元之间的距离。将（3）式带入（4）（5）两式，解得1047.90x，7939.50y.所以CT系统的旋转中心为（-9.1047,5.7939）.步骤五：计算光源的初始位置及180个方向现已知两个方向上观测数据以及两个光源的位置（角度差）,那么根据这两个就可以计算出每次旋转的角度每次旋转角度：∆θ=。。161151)90(0(将其他一般位置代入检验满足这个结果)。光源的初始位置：1-150°或130°。95.2.4模型总结分析以上分析全都是建立在结论一与前面建立的通用模型的基础上分析的，正是因为各个探测器接受的投影值0p与光线和椭圆（前面已论证其密度均匀且相等）的交线断长度成正相关，才能准确找出那两列（即横轴为-90和横轴为0）的特殊位置以及确定长轴短轴对应的数值，最后才能确定CT系统的旋转中心。5.3问题二5.3.1.用MATLAB对附件三进行数据拟合得出图6（左）图6附件三的图像拟合生成的分布图（左）与附件二分布图（右）的对照根据图6我们很容易看出附件三的问题二中未知介质的大致形状即为模板图除去小圆后的图形，也就是一个椭圆。再进一步分析其颜色分布发现这个椭圆与模板的椭圆内部结构又有所不同。对此我们进行一下探索：我们认为CT系统对未知介质的接收信息数据值的大小与模型内部结构有关。对此我们对附件三的数据用MATLAB进行拉登反变换，处理数据得到介质各点的吸收率数