CT系统标定及成像

CT系统参数标定及成像问题的研究摘要本文通过建立数学模型，给出了由大量数据标定CT系统参数和数据转化为图像的方法。针对问题一，主要是对附件一、附件二中数据的分析和处理，标定CT系统的各个参数。根据附件一中模板的吸收率，对样品做连续化处理得到标定模板的方程，建立以椭圆中心为原点的坐标系，引入旋转中心等参数。用Matlab绘制出模板接收信息的图像，将其与标定模板示意图进行对比，又因为介质为均匀介质，分析得到：表格数据大小反映射线扫描到的介质长度，每一列不为零的数据个数反映了射线照射到有介质区域的单元探测器个数。通过模板在探测器上的投影占据的探测器单元个数和单个数据对应的探测器单元的相对位置，确定了CT系统的探测器单元间距为0.2767mm,旋转中心坐标为（-9.2233,6,0182）。每个方向每个单元探测器接收到信息的数值与表格数据相对应，列出512个方程，得到180个方程组，用最小二乘法算得误差最小时所对应的X射线旋转的180个方向。针对问题二，根据问题一得到的标定参数，利用上述CT系统和某未知介质的接收信息，模板接收信息矩阵作radon逆变换、平移、像素转换等操作后，得到模板原始图像的“吸收率”为模板接收信息的2.0751倍，由未知介质的接收信息、标定参数及真实吸收率与接收信息间的关系，确定该未知介质外部轮廓为一个椭圆，分布在正方形区域中心位置，其中第三和第四象限分别被挖空了一个小椭圆，一、二象限间填充了两个小椭圆。针对问题三，处理方法与问题二类似，由于图像受数字噪声影响，所以先降噪，再运用滤波反投影法（FBP），先修正、后投影重建图像的做法，得到原始图像的吸收率信息。针对问题四，首先对问题一中射线旋转方向、探测器单元间距和旋转中心坐标分别进行精度分析和稳定性分析，又根据附件中所提供的数据，分别算出三个参数的精度，得到精度为0.01左右。在此基础上建立新的模板，给出相关参量，最后对新模板的参数进行分析并作出评价。关键词：CT系统、最小二乘原理、Radon逆变换、滤波反投影重建法1一、问题重述CT(ComputedTomography)可以在不破坏样品的情况下，利用样品对射线能量的吸收特性对生物组织和工程材料的样品进行断层成像，由此获取样品内部的结构信息。一种典型的二维CT系统如图1所示，平行入射的X射线垂直于探测器平面，每个探测器单元看成一个接收点，且等距排列。X射线的发射器和探测器相对位置固定不变，整个发射-接收系统绕某固定的旋转中心逆时针旋转180次。对每一个X射线方向，在具有512个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信息。而CT系统安装时往往存在误差，从而影响成像质量，因此需要对安装好的CT系统进行参数标定，即借助于已知结构的样品（称为模板）标定CT系统的参数，并据此对未知结构的样品进行成像。请建立相应的数学模型和算法，解决以下问题：(1)在正方形托盘上放置两个均匀固体介质组成的标定模板，模板的几何信息如图2所示，相应的数据文件见附件1，其中每一点的数值反映了该点的吸收强度，这里称为“吸收率”。对应于该模板的接收信息见附件2。请根据这一模板及其接收信息，确定CT系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该CT系统使用的X射线的180个方向。(2)附件3是利用上述CT系统得到的某未知介质的接收信息。利用(1)中得到的标定参数，确定该未知介质在正方形托盘中的位置、几何形状和吸收率等信息。另外，请具体给出图3所给的10个位置处的吸收率，相应的数据文件见附件4。(3)附件5是利用上述CT系统得到的另一个未知介质的接收信息。利用(1)中得到的标定参数，给出该未知介质的相关信息。另外，请具体给出图3所给的10个位置处的吸收率。(4)分析(1)中参数标定的精度和稳定性。在此基础上自行设计新模板、建立对应的标定模型，以改进标定精度和稳定性，并说明理由。图1.CT系统示意图图2.模板示意图（单位：mm）图3.10个位置示意图2二、问题分析1、问题一该问题主要为数据分析处理类问题。首先对附件一进行分析，通过对模板每一点吸收率连续化处理分别得到椭圆和圆的方程。分析附件二中数据，并用MATLAB画出对应矩阵的图像，对比模板，𝑀𝑖反映探测器与模板的相对位置，𝑀𝑖最大时探测器平行于y轴，𝑀𝑖最时探测器平行于x轴；模板由均匀固体介质组成，则经过增益后的数据与实际数据成比例，附件中每个单元格中数据即为增益后的截模板弦长的量值。为减小误差，选取数据优先考虑数据量大处的数据，利用探测器平行于x轴和y轴两个特殊位置，选取表格中数据与对应的实际数据用近似最小二乘定理算得误差最小时的比例系数μ；同理，选取有数据点对应的探测器单元个数与实际数值对应求得单元探测器的间距d；选取探测器平行于x轴和y轴的位置的数据，取每一列表格最大值所对应的列数即椭圆中心，求均值得到与旋转中心的相对位置，代入d值分别求得旋转中心的横纵坐标；计算探测器的180个方向，首先通过椭圆和圆的方程以及d值和初末位置的表格数据算出探测器初末位置的方向，列出均含有512个方程的180个方程组，计算得到的弦长经比例放大后等于表格中数据，用最小二乘法求得误差最小时的180个方向。2、问题二该问题是问题一的逆向求解。根据问题一得到的CT系统的参数，以及题目所给的接收矩阵，需要反向得出未知介质的几何图形以及在正方形托盘的位置和吸收率。我们通过查阅资料了解到Radon变换以及Radon逆变换能够解决这类问题。Radon变换能够通过接收矩阵得到原始图形以及图形的吸收率，但是由于题目中的吸收率进行了增益处理，故还需要找到Radon变换与题目中增益的关系。3、问题三该问题与问题二类似，不过通过Radon变换后得到的原图像非常模糊，应该首先进行去噪处理。运用滤波反投影法（FBP），先修正、后投影重建图像的做法，可得到原始图像的吸收率信息。4、问题四新模板必须尽可能排除或减弱上述三个参数所带来的误差。旋转中心：新的模板尽可能避免出现重复过的像素的点的行或列；探测器单元间距：新的模板需尽可能缩短探测器间距，增大旋转次数，尽可3能的获得携带原始图像关键信息的数据；X射线的旋转方向：每次旋转的角度尽可能小，这样才可能是的误差的量级变小。而且，在旋转一定的次数后，暂定，微调，对比数据进行修正。在上述三个参数的分析基础上，我们可以进行新模板的建立。为了便于新模板对上述三个参数的求解及模板的实用性及合理性，我们决定决定采用“正三角切圆型”，即在正方形托盘中放置一个均匀固体物质构成的正三角形薄片，正三角形中心与正方形托盘中心重合，正三角形某一边的中垂线与正方形的与水平方向成135𝑜的对角线重合，且使得正三角形面积尽可能的大。在此基础上，再在上述中垂线上（即在正三角形右下方）放置一个与正三角形具有相同介质的圆薄片，且同样使得圆薄片面积尽可能大，此时圆薄片，刚好与正三角形薄片相切三、模型假设1、假设CT系统的旋转中心位于探测器中垂线上2、假设忽略光的衍射四、符号系统𝑀𝑖：第i列不为0的数据个数(i=0,1,2,„,180)x轴：与椭圆短半轴重合的数轴y轴：过椭圆中心且与长半轴重合的数轴𝑥0：旋转中心横坐标𝑦0：旋转中心纵坐标d：探测器单元之间的距离𝐵𝑖𝑗：第i次旋转第j个探测器单元射线的截距𝑃𝑖𝑗：表格中第i列第j行的数据(i=1,2„，180，j=1,2,„,512)𝐿𝑖𝑗：探测器第i次扫描时第j个单元探测器扫描到的介质弦长和(i=1,2„，180，j=1,2,„,512)μ：问题一中弦长与衰减后射线能量的比例系数𝑘𝑖：射线斜率(i=1,2,„,180)θi：旋转第i次与X正半轴的夹角(i=0,1,2,„,180)𝑛𝑗：第n个探测器单元（j=1,2,3，„，512）m：从模板接受信息经过逆变换到图形矩阵的转换系数xi：附件四中十个点的横坐标（i=1,„,10）4yi：附件四中十个点的纵坐标（i=1,„,10）𝑃𝑖𝑗1：问题二中附件一第i列第j行的数据(i=1,2„，180，j=1,2,„,512)𝑃𝑖𝑗2：问题二中附件二第i列第j行的数据(i=1,2„，180，j=1,2,„,512)五、模型建立1、问题一该问题主要为数据分析处理类问题。首先对附件一进行分析，通过对模板每一点吸收率连续化处理分别得到椭圆和圆的方程与理论一致分别为：𝑥2152+𝑦2402=1，(𝑥−45)2+𝑦2=16分析附件二中数据，并用MATLAB画出对应的图像（如图一），根据图像，探测器逆时针旋转，则探测器处于平行于y轴时，沿x负半轴方向分别为第1个到第512个探测单元。𝑀𝑖反映第i次扫描有介质区域的长度，反映探测器与模板的相对位置，对比模板实际位置，𝑀𝑖最大时探测器平行于y轴，𝑀𝑖最时探测器平行于x轴；模板由均匀固体介质组成，则经过增益后的数据与实际数据成比例，附件中每个单元格中数据即为增益后的截模板弦长的量值，即𝑃𝑖𝑗=μ𝐿𝑖𝑗。图一附件二投影图1）计算μ首先选取𝑀𝑖最大时的六列数据，可以有效减小系统误差，取每一列数据数值最大的即探测器位于平行于y轴的位置，最大值表示椭圆短轴和圆直径吸收衰减后的射线能量经增益处理的量值，取六个方向平均值L1=67.3419，对应的实际值为38；同理选取𝑀𝑖最小时的三组数据，此时探测器位于平行于x轴的位置，两段不为0数据中的最大值分别表示椭圆长半轴和圆直径吸收衰减后的射线能5量增益后的量值，取三个方向平均值分别得到L2=141.6893667,L3=14.17846667，对应的实际数值为80和8;为得到更加精确的结果，对这三个数据用近似最小二乘法算误差最小时取得的比例系数，最终算得μ=1.7713，LINGO程序见附录一。2）计算d选取数据时首先考虑数据量大的位置，可以有效避免系统误差，提高精度。计算探测器单元间距d时，选取探测器𝑀𝑖最大时的六组𝑀𝑖，即探测器位于平行于y轴的位置，对取平均值𝑀𝑖得到椭圆长轴所对应的探测器单元个数，并对应长轴实际长度算出探测器单元间距：𝑑1=80289=0.276816609选取Mi最小时的三组数据，即探测器平行于x轴的位置，由于这几组数据中间部分有探测器单元没有扫描到介质区域，所以不能直接选取选取Mi计算，分别选取模板中椭圆右侧端点与圆左侧端点、椭圆左侧端点与圆右侧端点及圆和椭圆两短轴端点间对应的探测器单元个数，分别对应实际的数据得到三组探测器单元间距：𝑑2=2694=0.276595744𝑑3=64∗3694=0.27665706𝑑4=38∗3412=0.276699029对于得到的四组探测器单元间距，为了减小系统误差，求平均值保留四位小数得到d=0.2767，详细数据见附件表一。3）计算旋转中心坐标计算探测器旋转中心坐标时，利用已经计算出的d值和两条不相交的直线确定一个点得到旋转中心坐标。取探测器平行于y轴和平行于x轴两个位置的数据，这些数据中量值最大的对应的即为椭圆的中点，在这两个位置将椭圆中心即坐标系原点与旋转中心之间的探测器单元数目差值分别确定，找到模板和探测器系统的相对位置，代入d值，分别求得纵坐标和横坐标。探测器平行于y轴时椭圆中心对应的探测器单元：𝑛𝑦=(239+236+233+229)4⁄=234.25𝑦0=(256−𝑛𝑦)∗d=6.0182探测器平行于x轴时椭圆中心对应的探测器单元：𝑛𝑥=(222+223+223)3⁄=22.66676𝑥0=(𝑛𝑥−256)∗d=−9.2233详细数据见附件表二，旋转中心相对坐标原点的位置示意图如图二。图二旋转中心坐标示意图4）CT系统射线的180个方向为确定探测器旋转的180个方向，设过旋转中心的探测器单元的射线方程为𝑦−𝑦0=𝑘(𝑥−𝑥0)，则任意一个探测器单元在任意旋转次数时的射线方程一般式为：𝑦−𝑦0=𝑘(𝑥−𝑥0)+(256−𝑛𝑗)∗𝑑𝑐𝑜𝑠θi(n=1