CT系统参数标定及成像模型建模论文.docx

CT系统参数标定及成像模型
摘要
本文在CT工作原理的基础上通过建立参数标定及成像模型，实现了对安装好的CT系统标定参数，以及对未知结构的样品进行成像。
首先，根据尺寸、位置已知的模板的吸收率数据确定其材料的均匀性。将已知的吸收率和吸收信息代入Lambert-Beers吸收定律，可建立吸收信息与被测物体厚度的函数关系。另外，吸收信息中非零值越多，则单元数越多，投影越长，物体宽度越大。由此可用物体尺寸求得探测器单元间距。筛选数据得到模板的最长和最短投影，利用模板在探测器的投影位置找到二者的几何关系，可解出旋转中心的坐标。最后利用一定次数旋转的角度和旋转次数作商，可求得每次旋转的 角度，结合初始角，可知X射线的180个方向。
然后，因为衰减系数在不均匀介质中会随着位置而变化，沿照射路径作线积分即为接收信息。该过程与投影过程类似，因此运用直接傅里叶反变换法反解出衰减系数。因图像较大，可调用imresize函数缩放。附件3物体形状较为规则，可当作椭圆处理。由正方形托盘和该物体的几何关系可求出位置。将图形拆分开后，分别采用多项式拟合、幂指数拟合，以此描述形状。吸收率和衰减系数有比例关系，可由模板求得，带入被测件的衰减系数即可知其吸收率。通过坐标变换，将10个点的位置解出，则可知吸收率。

目录
CT系统参数标定及成像模型 1
摘要 1
一、问题重述 3
二、问题分析 3
2.1系统参数标定模型 3
2.2成像模型 3
2.3精度和稳定性分析 4
三、模型假设 4
四、定义与符号说明 4
五、模型的建立与求解 5
5.1 CT系统参数标定模型 5
5.1.1吸收率与接收信息 5
5.1.2探测器单元之间的距离 6
5.1.3 CT系统旋转中心在正方形托盘中的位置 6
5.1.4 X射线的180个方向 7
5.2 CT系统成像模型 8
5.2.1直接傅里叶反变换法（DFR）求衰减系数μ 8
5.2.2位置和几何形状 9
5.2.3吸收率 12
5.2.4十个位置处的吸收率 12
六、模型评价 14
七、总结 16
八、参考文献 16
九、附录 17

一、问题重述
ＣＴ技术通过Ｘ射线穿过被检测物体产生衰减的原理，能够在不破坏物体内部结构的前提下，对被扫描物体进行可视化成像[1]。具体的工作过程为从发射器射出平行射线，具有512个单元的探测器在被测物体之后接收信息。发射器和探测器绕某固定的旋转中心逆时针旋转180次，且二者相对位置不变，进而记录下180组数据，经成像过程生成具体图像。然而CT系统安装时存在的误差会影响成像质量，因此需要对安装好的CT系统进行参数标定。
问题1： 现已给定标定模板的位置、几何信息、该介质的吸收率和经过扫描过程探测器的接收信息。需要根据已知条件求解出CT系统旋转中心的位置、探测器每一个单元之间的距离和X射线入射的180个方向。
问题2：根据问题1的解，可以完全确定CT系统的工作状态，现再提供某未知介质的接收信息。凭借CT系统各参数的值和该未知介质最后得到的接收信息结果，确定出该未知介质的位置、形状、吸收率等信息。并在该吸收率情况下找出10个位置的吸收率。
问题3：同问题2，利用CT系统各参数的值和接收信息，确定该未知介质的相关信息。并求出10个位置的该未知介质下的吸收率。
问题4：分析问题1求解的参数标定的精度和稳定性。并自行设计新模板、建立对应的标定模型改进标定精度和稳定性。
二、问题分析
本题需要解决的问题可分为两部分：对CT系统进行完整的参数标定，继而用CT系统实现成像功能。
2.1系统参数标定模型
参数标定从本质上讲是成像过程的逆推过程。即已知模板几何尺寸，位置信息，接收信息，根据已知信息可以确定发射-接收系统的工作状态。接收器的单元间距未知，但是从已知信息可以推算出模板投射在探测器的投影单元数，又知道模板的实际尺寸，就可知二者间换算关系，从而确定单元间距。旋转中心难以直接计算，但因为吸收率一定，所以探测器显示的吸收信息蕴含着物体的长度和宽度，再辅以旋转180次得到的不同角度的信息，结合几何关系，可以求解。角度问题可以直接对应探测器吸收信息的变化情况求得。
2.2成像模型
标定各参数之后，整个发射-接收系统的工作状态就可以完全确定。CT系统的工作原理为CT扫描被测物体后得到吸收信息，进而反解出被测物体的其他相关信息。可以发现该工作过程可以对应直接傅里叶反变换法的计算过程，因此，可以用此方法方便快捷地求得正确的衰减系数。得到被测物体的衰减系数，就可以根据各参数的实际意义与衰减系数的实际意义之间的联系推算出其余信息，包括位置、形状、吸收率等参数。已知整个未知介质的吸收率后，即可求得同坐标系上任一点的吸收率。
2.3精度和稳定性分析
由于建立模型和求解模型中将部分模型理想化，因此需要在最后做一个精度和稳定性的分析，以避免误差过大，影响模型的使用和推广。
三、模型假设
1.椭圆投影长度最长（短），对应X射线与椭圆相切于长轴（短轴）。
2.投影长度相等的若干组数据，取中间一组为极值。
3.求未知介质位置时，可看作规则图形求其中心位置。

四、定义与符号说明

a
吸收率

μ
衰减系数

d
被测物体的厚度

初始强度

I
射出强度

P
接收信息

L
检测器单元之间的距离

θ
发射-接收系统每次转动的角度

五、模型的建立与求解
5.1 CT系统参数标定模型
5.1.1吸收率与接收信息
借助模板标定CT系统的参数，首先需要正确理解各参数的实际意义并明确相互之间的关系。
由题可知，附件1中的数据反映的是该点的吸收强度，即吸收率。通过观察可以发现，吸收率的数值分别为0.0000和1.0000，且1.0000对应的点的分布具有一定的规律性。为了更直观地观察其分布情况，利用MTALAB读取表格数据并绘制图像（程序见附录1），如图5.1.1所示。

图5.1.1 吸收率
图中高亮的部分为吸收率为1.0000的点所在区域，其余部分是数值为0的区域。将图5.1.1与模板的几何信息对比可以发现，吸收率非0的区域形状、尺寸、位置均与模板匹配，即模板确为均匀介质，且该介质的吸收率为1.0000。
接收信息定义为X射线穿过模板的过程中，被模板所吸收的信息，反映在附件2，为512×180的表格。由CT的工作原理可知，某一列的512组数据即为X射线从某一个角度发射，探测器512个单元接收到的信息。某一行的180组数据代表该行对应的单元，随着反射器转动180次而分别得到的数据。
由CT的图像重建原理可知，接收信息来源于X射线照射在被测物体上，部分能量被物体吸收而导致的射线强度衰减。这种衰减呈指数变化，如果发射器发射出的射线是单能的，并且照射在了均匀材料的物体上，则遵循Lambert-Beers吸收定律[2]：
（1）
μ为衰减系数，d为被测物体的厚度，I0为X射线穿过被测物体前的初始强度，I为X射线穿过被测物体后的射出强度。
由吸收率的含义可知，吸收率a和衰减系数μ保持一定的比例关系：
（2）
若定义P为接收信息量，且满足式（3）：
（3）
根据吸收率和接收信息的实际意义，可知接收信息P的数值大小与吸收率a以及被测物体厚度d有关。联立（1）、（2）、（3）式，可得三者之间的函数关系：
（4）
分析可知，吸收率一定的情况下，被测物体越厚，则接受信息的数值越大。
5.1.2探测器单元之间的距离
求解探测器单元之间的距离，需要将模板的真实尺寸与其投影在探测器的长度进行比较。投影在探测器上的长度与附件2的数据相关。每一列的数据为X射线从某一个角度发射，探测器的512个单元接收到的信息。由接收信息的含义可知，数值为0的单元的X射线是从发射器直接发射到探测器的，不经过模板的任意部位。而数值不为0的单元对应的射线，不论数值大小，均穿过模板。因此，只要确定出每一列的非零数据的数量，就知道有多少个单元的射线穿过模板。
由于使用Excel对每一列接收信息数据进行观察的过程比较繁琐，而且容易出现错误，因此，可以运用MATLAB导入数据，编写程序得到非零数据的数量。由于某些列的非零数据分成了两部分，所以所编写的程序将两部分拆开计数，最终得到的结果为两组数据（程序见附录2）。得到数据后进行观察，发现其中一组数据绝大多数为29，因此可以确定这组数据是射线穿过圆形模板得到的，29个单元对应的就是圆形模板的直径8mm。则探测器单元之间的距离L为： (5)
5.1.3 CT系统旋转中心在正方形托盘中的位置
由于CT系统的旋转中心一定在探测器的中垂线上，因此找任意两个探测器位置，分别做其中垂线，相交的点即为旋转中心。
由5.1.2可知，探测器单元之间的距离已知，因此，取附件2的任意列，从上向下查到非零数据点的时候，有m组数据就对应着m条射线没有穿过模板。又由5.1.1可知模板投影到探测器对应的单元个数n，且总单元数为512，由此可知模板数据下方的零数据点共有512-m-n个。
如图5.1.2所示，AF为探测器,B为正方形托盘左下角，也是坐标原点，C和E分别与椭圆相切，D点为探测器AF的中点。根据附录2的数据，可以找到一组最小，一组最大数量的数据，最小数量为109组数据，为CE段，对应上方零数据点有168个，为EF段，下方零数据点有235个，为AC段，对应的是椭圆的短轴，即AC:CE:EF三者比值为235:109:168；最大数量为289组数据，上方零数据点有69个，下方零数据点有134个，对应的是椭圆的长轴。
以正方形左下角B为原点建立直角坐标系，则旋转中心O的坐标可表示为：
（6）
同理，纵坐标可表示为：

综上，旋转中心坐标为（40.7931,56.0690），如图5.1.2所示。
A B C D E F

图5.1.2 旋转中心
5.1.4 X射线的180个方向
由5.1.3可知，最大数量的是第61组数据，最小数量的是第151组数据，并且这两组数据分别对应的是水平射线和竖直射线，即夹角为90度。
假设180次旋转，每次转动的角度相同，则：
(7)
由于发射-接收系统绕逆时针旋转，根据附件2的数据可推断出，投影长度先增加后减小再增加，因此初始位置是从正方形托盘的右下方发射，左上方接收。又因为椭圆形长轴投影对应的是第61组数据，即从初始位置旋转了60次，因此，初始位置的射线与水平方向夹60度角。
综上，X射线从与水平方向夹60度角开始旋转，每次旋转1度。
5.2 CT系统成像模型
5.2.1直接傅里叶反变换法（DFR）求衰减系数μ
直接傅里叶反变换法是由不同角度的投影的数据，恢复原始图像的方法[3]。
由于衰减系数μ在不均匀介质中会随着位置的不同而变化，因此接收信息P与μ的关系满足：
(8)
其中d为X射线穿过物体的路径，μ（x）表示d上某点x处的衰减系数。可以看出，投影过程和X射线照射过程均为线积分，因此由附件3中的接收信息P可以使用直接傅里叶反变换法方便地求解出衰减系数μ的数据。
具体过程为：首先原始图像fx,y作线积分得到投影p∅(r)，再由投影值p∅(r)求一维傅里叶变换P∅(v)，把变换所得当作二维傅里叶变换的一个切片。所有角度的切片数据放在一个极坐标系下，然后将极坐标系栅格化，即使用二维插值求出笛卡尔坐标系中相对应的值F(u,v)，再进行二维傅里叶逆变换，得到原始图像f(x,y)。具体流程如图5.2.1所示。

栅格化

二维傅里叶逆变换 一维傅里叶变换
投影

图5.2.1 流程图
分别导入附件3和附件5，求出衰减系数μ的数据表格并生成图像，发现数据量较大，因此调用MATLAB中的imresize函数运用插值法可进行一定程度的缩减，最终如图5.2.2。（程序见附录3）设定颜色越暗则衰减系数越小，反之，越亮则衰减系数越大。

图5.2.2 附件3、附件5衰减系数图
5.2.2位置和几何形状
根据图5.2.2观察可知，附件5的形状不规则，外边界近似为一个正方形，内部有许多暗色区域，分布无规律，因此很难谈论其形状以及位置信息。
附件3的几何形状近似为椭圆，因此，可用类似求解椭圆中心的方法求解该未知介质的位置问题。与附录2程序功能类似，用MATLAB编写程序将附件3的每一列数据中非零数据的数量算出，并以该数值为纵坐标，对应旋转次数为横坐标建立直角坐标系，如图5.2.3（程序见附录4）。从图5.2.3中可以清晰地看出第58组的非零数最多，第146组的非零数最少，由于该组数非常接近第61和151组，因此可近似认为这二者对应的投影长度分别为椭圆长轴和短轴投影长度。

图5.2.3投影长度—旋转次数
由于设备和旋转中心均不变，因此探测器在第61和151组的位置不变。竖直放置的探测器非零数值的数据量为296，其上方的零数值数据量有96，下方的零数值数据量为120，由5.1.3的几何关系可知，最下端到x轴的距离为10-69×829=-9.0345mm；水平放置的探测器非零数值得数据量为156，其上方的零数值数据量有215，下方的零数值数据量为141，同理可求得最左端到y轴的距离为35-235×829=29.8276mm。所以，该介质的位置，即椭圆中心的坐标为：
x=43.032+215×829-29.8276=50.9977mm (9)
y=81.662+96×829-9.0345=58.2783mm
综上，未知介质的位置，即中心坐标为（50.9977,58.2783）。
由5.2.1可得到衰减系数的数据，并可清晰地观察分析出衰减系数在不同位置的大小。分析图5.2.2可知，明亮部分衰减系数较大，说明此区域的形状为该介质的几何形状。但具体的几何形状难以通过观察分析得到，因此，可以取衰减系数值恰当，在图像中处于明和暗的交界处的点(假设为此边界衰减系数为0.15)，进行曲线拟合，得到具体的函数关系，以此来表达该介质的几何形状。
但由于该图形为封闭式图形，无法直接用函数的拟合方式进行，故可运用MATLAB编写程序将一个闭环拆分成上下两部分开环，分别进行拟合，用两个式子来完整表达形状（程序见附录5）。另外通过对图5.2.2的观察可发现，介质内部有两块区域衰减系数较小，甚至接近于外部区域，因此采用数值比较的方法将此两块区域和其他内部区域分开考虑。而这两块区域的边界也是封闭式的，同样需要拆分开进行函数关系的拟合。综上，可拟合出以下六式（图5.2.3）：

（10）
（11）
图5.2.4 外边界拟合曲线
（12）
（13）

（14）

（15）
图5.2.5 内边界拟合曲线
综上所述，该未知介质的几何形状即为（10）~（15）式联立后围住的区域，可由图5.2.5观察分析其形状特征。
5.2.3吸收率
将5.1当中模板的吸收信息按照直接傅里叶反变换法进行计算，将会得到其对应的衰减系数数据以及图像（如图5.2.6）。

图5.2.6 模板衰减系数图
由5.1可知，衰减系数和吸收率呈正比关系，则通过现在模板的衰减系数μ模和吸收率a模以及未知介质的衰减系数μ未，可以推测出未知介质的吸收率a未。具体的做法是：从模板的衰减系数数据中判断出椭圆位置，因为椭圆模板内部的吸收率均为1.0000，因此可从椭圆内部任取足够多的点求其平均值μ模，现取DF-EM列，76-109行数据，共1156个，得到μ模=0.623452。由于a模=1.0000是定值，因此，未知介质的吸收率：
a未=μ未×a模μ模=1.6040μ未 （16）
用Excel软件打开附件3、附件5对应的衰减函数图，每一个衰减系数乘一个固定的数值1.6040即可得到最终需要获得的该介质的吸收率数据表，具体数据见附件problem2.xls及problem3.xls。
5.2.4 “空点”的筛除
由于探测器的单元间隔较大，投影次数太少等因数的影响，会导致探测器探测到的吸收信息有限，从而使得利用这些信息不足以完全还原二维情形。所以才会出现吸收信息矩阵中有0值，但是变换后的系数矩阵没有0值的情况。对应的实际情况即没有物体的区域衰减系数也有非零值。所以为了精确得出二维真实情形可设计一个阈值来筛除‘’空点‘’。
利用Excel表格来确定衰减系数矩阵的最大值M、最小值m。将左右端点值分别为m、M的区间等分为N个小区间，第i个小区间记为 ni(i=0、1···N)，则第i个小区间的右端点值为ni(M-m)N+m。再以落在某个小区间内元素的个数作为该小区间的函数值。运用MATLAB绘出将横坐标区间等分成10分的图像。（如图5.2.7）

图像大致分为三段，第一段函数随着横坐标值增加，函数值迅速减小。该段表示此区间内的点绝大部分的吸收系数接近0，并且整体变化一致。故第一段应该对应“空点”即不是实体上的点。第二段前半部分曲线缓慢减小，且函数值很低，后半部分几乎不变。这段图像前半部分说明这类点数目少但吸收系数分布区间长，这类点出现的原因可从图5.2.2看出，在靠近椭圆的黑色区域有“白线”存在，说明这些区域的点较黑色区域有较高的值。但这些点是由于radon逆变换时产生的干扰，它不反应真实物体情况，所以设计的阈值应该排除掉这些点。第二段曲线后半部分可以认为真实点开始变多，从而曲线变平缓。第三段为上升段，该区域的点应为实体上的点，应该保留。综上所述，我们选择当第二段曲线的变化刚接近平缓的点0.413为阈值。在Excel中用快速分析功能的色阶分析表格可以大致得出类似于图5.2.2的图像如下，而后在其基础上令值大于阈值的格子填充为粉色。另外再作一个图，先使值大于阈值的格子填充为粉色，再用快速分析功能的色阶描述表格，使后两图相互印证。

未填充的图 先色阶分析再用粉色填充的图 先用粉色填充再用色阶分析的图
如图可看出粉色区域几乎覆盖了绿色区域，但边缘还有一些绿色，也就说明阈值有点偏高。在筛选时可能筛掉实体点，但可以近似认为阈值是合适的。用阈值对衰减系数矩阵进行筛选。
对于衰减矩阵中的值小于阈值令其为0，大于等于0则保留原值。

第三题阈值寻找方法如第二题所述，并作出N=10的曲线。

图线大致分为两端，第一段如第二题中所述对应“空点”第二段连续而缓慢上升到最高，说明此时图像干扰点不多可以忽略，此时可以只用把第一段对应点排除。综上所述以0.067为阈值。先在Excel中用快速分析功能的色阶分析表格可以大致得出类似于图5.2.2 的图像如下第一个图所示，而后在其基础上令值大于阈值的格子填充为黄色如下第二个图所示 。另外再作一个图，先使值大于阈值的格子填充为黄色，再用快速分析功能的色阶描述表格，使后两图相互印证。

未填充的图 先用色阶分析再用黄色填充的图 先用黄色填充再用色阶分析的图
从第二张看出黄色基本覆盖了绿色，第三张除了一些边缘绿色基本覆盖了黄色，这说明该阈值有点低，可能会在筛除时漏过一些“空点”但基本还适合。用与第二题同样的法则对衰减系数矩阵进行筛选

5.2.5十个位置处的吸收率
由于X射线的初始位置与竖直方向夹30度角，因此经过直接傅里叶反变换产生的图像并非实际方向，而是与竖直方向夹30度角的。因此，想要求出10个位置的吸收率，首先需要把这些点的坐标转换到Excel表格对应的坐标系——衰减率坐标系中。
已知模板的形状和位置，即椭圆长轴平行于y轴，因此，坐标系可共用一个模板，即围绕这个模板作两个坐标系。但由于两坐标系的原点不在同一位置，故需要采取一个中间坐标系进行转换。因为在衰减率坐标系中难以找到椭圆中心，因此可选取圆形的圆心作为中间变换坐标系的原点。
利用Excel快速分析中的色阶分析命令分析，可以找到圆形模板在表中所对应的位置，如图5.2.7所示。观察圆的边缘，当某一行的数据出现突变时，可以认为该行所在直线是圆的切线，记录行数。同理可以找到另一条切线并记录该切线所在行数，由此可知圆的直径大小为两行行数之差。同样，可以找到与圆相切的两列，得到直径大小。最后以二者的平均值作为圆的最终直径，两列序号的中间值为圆心的纵坐标，两行序号的中间值为圆心的横坐标。

图5.2.7 圆形模板位置
以正方形托盘的左下角为原点o，水平向右为x轴，竖直向上为y轴建立平面xoy坐标系。模板中的圆形模板在此坐标系下的圆心位置为（95,50）。再以Excel表格的左上角（1，A）为原点o2，水平向右为v轴，竖直向下为u轴建立平面uo2v坐标系。圆心在此坐标系下的坐标为（，）。
为了能实现10个位置坐标的变换，需要以圆形模板的圆心o1为原点，平行于x轴的轴为x1轴，平行于y轴的轴为y1轴建立x1o1y1坐标系。同理建立u1o1v1坐标系。坐标系的建立如图5.2.8所示。

图5.2.8 坐标
假设A点为10个点中的任一点，在xoy坐标系下的坐标为（xA，yA），在x1o1y1坐标系下的坐标为（xA1，yA1），在uo2v坐标系下的坐标为（uA，vA），在u1o1v1坐标系下的坐标为（uA1，vA1）。
已知圆形圆心o1在xoy坐标系下的坐标为（95,50），由向量相减准则可知：

o1A=oA-oo1
整理得：
xA1=xA-95yA1=yA-50（17）
因为x1o1y1和u1o1v1坐标系原点均为o1，夹角为-120度，故可用坐标变换公式求得（uA1，vA1）：

整理得：
uA1=(-12xA1-32yA1)×nvA1=(-12yA1+32yxA1)×n（18）
再运用向量相加准则：
o2A=o2o1+o1A
uA=+u1AvA=+v1A（19）
经过（17）、（18）、（19）式，可将10个点的坐标变换到u02v坐标系中，取整后查找Excel表格中对应坐标的数值，即为该点的吸收率。最后，得到坐标变换前后坐标及10个点的吸收率。如表1所示。
表一 10个点的坐标变换及吸收率

六、模型评价
首先，第一题有3个小问题，我们分别分析3个问题的精确度

第一个是求探测器单元之间的距离，我们是通过小球的长度除以投影的数量(射线通过小球投影到探测器单元的数量)即8/29.但这个数就存在误差。如图所示，只有如图一，单元格对应的射线左侧和小球相切时小球的直径才对应射线通过小球投影到探测器单元的数量，而正常情况下如图2，显然单元格射线线左侧和和小球之间存在距离，这就是是误差。而当图三所示探测器单元的右侧和小球相切时即为临界值，即为误差最大值，此时射线通过小球投影到探测器单元的数量其实为27个，探测器单元之间的距离为8/27。综上所述探测器单元之间的距离x为8/29<=x<8/27，所以它的误差为8/27-8/29=0.0204.

图一 图二

图三

第二个是求X射线的180个方向，我们是通过求射线水平射入椭圆的组数和垂直射入的组数差即为90度，即投影最大时(80)是水平射入，投影为30时是垂直切入,但是无法准确找到投影为80和30的点。只能找到79.7241和30.0690这两个点，这是误差一，而这2个点本身就是通过射线通过椭圆投影到探测器单元的数量除以探测器单元之间的距离得到的，而探测器单元之间的距离就有误差，这是误差二。而且有8组投影是79.7241,分别为58-65组，4组是30.0690，分别为149,151,153,154组。虽然我们是通过2组数据的中间值相减，但还是存在误差，这是误差三。通过误差三可以求出每次旋转的方向角a范围为0.9375<a<1.0714, 误差为1.0714-0.9375=0.1339。入射角（和水平方向的夹角）的范围为57.1875~65.3554，误差为8.1679。

第三个是求CT系统旋转中心在正方形托盘中的位置。我们是通过求射线水平射入探测器时和垂直射入时中垂线的交点即为旋转中心。这里用到探测器单元之间的距离，因为这个值的范围为8/29~8/27，存在误差。且有8组水平投影最接近80，4组垂直投影接近30，找到4组中2个临界点，得出5.1的算法的AC:CE:EF=236:109:167和234:109:169，和8组中2个临界点比例为132:289:91和136:289:87。最后可以得出最后旋转中心的横坐标的范围为40.5172~41.5185和纵坐标的范围为55.5172~60.0741，误差分别为1.0013和4.5569。

求稳定度时，我们分析附件2存在误差，即在每一列中0和非0数据之间临界值误差。因此我们把每列中0和非0相邻的0变成非0。即把每个角度射线通过图像投影到探测器单元的数量加2（上面和下面分别加1），求出每个角度投影的长度，然后得出问题一的结果。分析这个结果和我们求得的结果是否相似，以此来确定它的稳定性。
首先我们把每个角度射线通过图像投影到探测器单元的数量加2后，通过小球投影到探测器单元的数量由29变成31，所以探测器单元之间的距离为8/31。然后是求X射线的180个方向，垂直射入时的组数是不变的(因为是找投影最大值)，仍然为为58-65组，但水平射入时有变化，根据附录程序可求出第159组投影最接近30，为29.9355。以61组垂直射入时的组数，因此可求得旋转角度为90/ (159-61)=0.9184，入射角（和水平方向的夹角）为(61-1)*0.9184=55.1020。最后求CT系统旋转中心在正方形托盘中的位置，以61组为水平射入方向，159组为垂直射入方向，得出5.1的算法的AC:CE:EF分别为228:116:168和133:291:88。计算出旋转中心为(42.2258，53.3548)。
第一题求得的4个数据分别为距离8/29、旋转度为1度、入射角为60度、旋转中心(40.7931,56.0690)，因此，探测器单元之间的距离的相对误差为(8/29-8/31)/(8/29)=0.0645, 旋转角度的相对误差为为（1-0.9184）/1=0.0816，入射角的相对误差为(60-55.1020)/60=0.0816，旋转中心横坐标的相对误差为(42.2258-40.7931)/ 40.7931=0.0351纵坐标的相对误差为(53.3548-56.0690)/ 56.0690=-0.0484。因此，当初始数据有扰动时，参数依然可以标定相对稳定。

七、总结

八、参考文献
[1]郭立倩. CT系统标定与有限角度CT重建方法的研究[D].大连理工大学,2016.

[2] 卢彦斌. X射线CT成像技术与多模态层析成像技术研究[D].北京大学,2012.

[3]庄天戈. CT原理与算法[M]. 上海交通大学出版社, 1992.

九、附录
附录1：根据吸收信息求吸收率
data=xlsread('A题附件.xls','附件1');
imshow(data);
title('吸收率图像');

附录2：椭圆和圆投影所占的探测器单元数量
data=xlsread('A题附件.xls','附件2');
max=[];%定义放椭圆投影的数组
min=[];%定义放圆投影的数组
for i = 1:size(data, 2) %180列循环
num1=0; %记录2个图之一的投影
num2=0; %记录2个图之一的投影
flag=0; %标志位
for j = 1:size(data, 1) %512行循环
if data(j,i)~=0
if flag<2
flag=1;
end

if flag==1
num1 =num1+1; %如果flag为1的时候num1+1

else %如果flag不为1的时候num2+1
num2 =num2+1;

end

elseif flag==1
flag=2;
else

end
end

if num1>num2 %如果num1大于num2，num1属于max，num2属于min
max(i)=num1;
min(i)=num2;
else %反之，num2属于max，num1属于min
max(i)=num2;
min(i)=num1;
end

end
disp(max) %分别输出max和min数组
disp(min)

附录3：吸收信息转换为衰变系数
（1）data=xlsread('A题附件.xls','附件3'); %读入附件3的数据变成矩阵
show_img(data); %进入函数
（2）function [ output_args ] = show_img( x )

p_N=256; %图像默认大小256
theta_N=180; %180度平行投影
pad_N=1024; %投影后变为367*180，每列数据就是当前角度下的投影值。367<512考虑到傅里叶变换是需要基2对其，故扩展为1024*180
theta=0:(theta_N-1); %0~180度平行投影
output_args=0;
%% 根据每个投影角度下的投影数据，求其一维傅里叶变换
proj_sino=x;
if mod(length(proj_sino(:,1)),2)==1 %如果奇数个接受栅格，则需要扩展为偶数，填充0
proj_sino=[proj_sino:zeros(1,size(proj_sino,2))]; %填充一行，367*180变为368*180
end
pad_row=(pad_N-size(proj_sino,1))/2; %（1024-368）/2=328
proj_sino=padarray(proj_sino,[pad_row 0],0,'both'); % 368*180,上下填充328行，变为了1024*180
L_pad=pad_row+ceil(((p_N.*sqrt(2)+2)-p_N)/2)+1; %原始图像大小为p_N,傅里叶逆变换后为pad_N,则需要截断数据
proj_sino=ifftshift(proj_sino,1); %将实际数据移动到两端，中间的是填充的0
f_p=fft(proj_sino,[],1); %求取正弦图的一维傅里叶变换，1024*180
f_p=fftshift(f_p,1); %反向处理，中间的数据移动到两端，两端的数据移动到中间
%%极坐标化为笛卡尔坐标系
nfp=length(f_p(:,1));
omega_sino=(-(nfp-1)/2:(nfp-1)/2).*(2*pi/size(f_p,1));
%极半径范围 -pi~pi，分成1024等分
theta=theta*pi/180; %角度转化为弧度，范围0~pi,180等分
[theta_grid,omega_grid]=meshgrid(theta,omega_sino);
%网格后，omega_grid和theta_grid都是1024*180，theta_grid_x每一列都是一样的角度，omega_grid没一行都是一样的位置

omega_image=omega_sino; %根据极半径大小，建立笛卡尔坐标系
[omega_grid_x,omega_grid_y]=meshgrid(omega_image,omega_image);