CT系统参数标定

CT系统参数标定

摘要

问题一的求解：由于双球心球心投影轨迹的坐标数据是以平板探测器的坐下角为坐标原点，以水平向右方向为X轴的正方向，以垂直向上方向为Y轴的正方向，我们在原有二维坐标系的基础上，仍以探测器左下角为坐标原点，以垂直于坐标平面指向光源方向为Z轴的正方向建立坐标系。并且，我们以像素作为坐标系的长度单位。

问题二的求解：利用cholesy分解法解对称正定矩阵，对两条曲线进行椭圆曲线拟合，得到椭圆长短轴长度、椭圆中心坐标。再设出需要的参数，包括电光源与平板探测器的垂直距离、载物台中心转轴与电光源间垂直距离、玻璃台中心与载物台转轴垂直距离，通过几何关系得到未知参数间彼此的关系，联立方程解出。

问题三的求解：误差分析主要讨论了模型误差、计算误差、机械舍入误差。其中最主要的为计算误差，并进行了量化分析。

关键词：正交坐标系曲线拟合空间几何平面化误差分析

一、问题重述

CT自发明以来，被公认为自伦琴发现X射线后是放射领域的最重要的发明之一。医用CT有出现为人类的医学诊断带来福音。在实际诊断中，X射线放射源和探测器围绕人体旋转一周，就可以采集出一组图像，通过一定的重算法实现扫描切片的重建。医用CT之所以采用X射线源和探测器旋转结构，主要是为了避免人体的旋转不便。而工业CT一般采用扫描物件旋转和面探测器的结构，同样实现相应的功能，其示意图如图1所示。

图1：工业CT原理示意图

扫描物件围绕某一固定转轴旋转，每隔一定角度采集一张图像，然后根据采集的图像采用3D图像重建算法即可将原始3D物件重建出来。在实际中，由X 射线源发出X射线，经过扫描物体衰减后照射在探测器上，探测器根据接收到的光子数的计数实现光电转换，从而形成灰度图像。一般平板探测器大小为3000×2000像素，每个像素为0.127mm。

在CT系统安装过程中往往存在机械误差，而这些误差对于物件重建的准确

性往往是致关重要的，实际中就需要对安装好的CT系统进行参数标定，系统参数主要包括X射线源的位置、载物中心线的位置和探测器的位置参数等。传统的方法是采用尺子测量法，对于CT系统参数标定精度要求是不够的（通常CT重建过程中系统参数的机械误差不允许超过一个像素）。目前，一般通过实验的方法实现参数标定。所谓实验的方法标定参数就是通过扫描已知参数的体模，分析投影数据，估计系统参数的方法。针对工业CT要求解决如下问题：（1）建立合适的坐标系，正确描述CT系统的各种参数和机械误差，并建模分析这些参数的关系和可能的机械误差。

（2）通常采用轴承钢球作为CT参数标定体模。因为它具有各向投影一致，边缘清晰等优点，在里采用两个钢球实现CT的参数标定实验（如图2所示）。两个钢球置于有机玻璃管中（X射线容易透过有机玻璃，容易后期球心获取的图像处理），钢球直径为8mm±0.0008mm，两球心距离约为100mm±1mm。将有机玻璃管固定在旋转载物台上，载物台携带有机玻璃管以均匀速度旋转，与此同时X 射线源发出射线，探测器采集数据（旋转一周，等间距采集180张）。附件1给出了根据180张采集图像提取出的球心投影坐标。两个钢球球心投影的轨迹是两个椭圆（如图3所示）。试根据实验数据估计该CT系统的参数值，即给出CT系统的标定。

图2：双球体模示意图

图3：双球球心投影轨迹

（3）在参数标定过程中可能存在多种可能的机械误差，试就你的参数标定可能的误差进行分析。

二、模型假设

1.X辐射光源为点光源，忽略光源本身物理尺寸

2.旋转载物台为完全水平放置。

3.X辐射点光源，载物台中轴线所确定的平面在平板探测器的某个垂面上。

4.平板探测器表面完全平整光滑。

三、符号说明

H: 光源点与成像平板之间的垂直距离

L：光源点与圆形旋转载物台的圆心之间的水平距离

r: 旋转载物台的半径

E：两椭圆形双球球心轨迹长轴长的平均值

：上椭圆形双球球心轨迹短轴长

：下椭圆形双球球心轨迹短轴长

d: 旋转载物台上两钢球之间的垂直距离

D：上椭圆形轨迹的下短轴端点与下椭圆形轨迹的上短轴端点之间的距离

：上椭圆形轨迹的下短轴端点与光源点在平板上的垂直投影之间的距离

：下椭圆形轨迹的上短轴端点与光源点在平板上的垂直投影之间的距离

四、模型的建立与求解

问题一：

1.椭圆曲线拟合原理

用matlab对已知数据进行描点，观察图形，近似为椭圆。坐标系利用题中给出的坐标系，暂不考虑坐标系变换。

图4：matlab描点作图

我们已知的是上方球心坐标与下方球心坐标各180组数据，需要由此拟合出相应两个椭圆的方程以获取我们需要的参数。

从数学的角度看，考虑到实际的观测误差，本问题可以归结为：已知M个二维数据，找一个恰当的椭圆

使这些观察点最接近曲线。

选择平方偏差函数以度量曲线拟合的可靠性。

要做的是，找到参数组，使得

公式(1)的两边同时乘以一个常数并不影响问题的实质，不妨假定(1)中的A

取值为1。其次，具体考虑CT系统，因为载物盘上为一圆盘在转动，所以椭圆的长轴相对于水平，即我们所取的平面是平行的。因此取项前的系数2B为0。

综上，目标函数(2)是一个关于参数的二次多元函数。问题归结为，

求参数，使得(3)式成立。

由多元函数极值定理，要求满足一下四个方程组成的方程组：

引进列向量和未知参数列向量，则前面的方程组可以写成：

引进列向量和矩阵，这样就有

相应的矩阵表示为

容易知道对称正定线性方程组(4)的解是存在的。如果矩阵U是满秩的，则方程(4)的解是唯一的，在此我们假设这个条件是满足的。

下介绍解对称正定线性方程组所用到的cholesy分解法。讨论一般的对称线性方程组为

这里为系数矩阵，其为对称正定矩阵，x为未知量，b为右端项。

将S分解为

其中，L为下对角阵，

比较(6)式两端对应元素可得

于是得到求L各项的递推公式

这样，求解方程组(5)就化为以下三角形方程组

其中

至此完成椭圆曲线拟合，具体以C语言实现，在/code中。

2.光源点坐标、载物台位置及半径的确定

首先，我们需要建立一个坐标系，由于双球心球心投影轨迹的坐标数据是以平板探测器的坐下角为坐标原点，以水平向右方向为X轴的正方向，以垂直向上方向为Y轴的正方向，我们在原有二维坐标系的基础上，仍以探测器左下角为坐标原点，以垂直于坐标平面指向光源方向为Z轴的正方向建立坐标系。并且，我们以像素作为坐标系的长度单位。

之后，建立如图所示的解析图，其中图5为光源所在的与平板平面垂直的剖面图，而图6为旋转平台所在平面的切面图。

A

由于图中各三角形之间存在相似关系，可得到一系列的比例关系。可以得到：

对上述两式进行变换，可以得到如下关系

又因为

可以将与用如下关系式表示：

同样，我们可以发现与之间存在如下关系：

于是，可以将与以如下方式表示：

之后，我们可以继续考查图中的比例关系：

由EHD~HAI可得

①

由ADF~ABG可得

②由PST~PQO可得

③

由①、②、③可得

令，则L=(+1)r

另一方面，由于椭圆的标准方程为

其中，点(a,b)为椭圆的中心，A为椭圆的长半轴长，B为椭圆的短半轴长（当A>B时），通过转化，可以得到如下形式：

即为：

由上下两个椭圆的拟合方程可以发现：

对于上面一个椭圆，有

对于下面一个椭圆，有

解得：

对于上面一个椭圆，其中心坐标以及长短半轴长如下：

对于下面一个椭圆，其中心坐标以及长短半轴长如下：

设光源点在平台上的垂直映射点的坐标为(X,Y)。

根据假设，光源点与旋转台圆心的连线在平台的垂面上，我们取两拟合椭圆中心的X坐标的平均值，可得

又因为

取两者的平均数作为光源点的映射点的Y坐标，得

同时，我们可以确定如下的关系：

从而得到E，D，，的值分别为：

由于两球心之间的距离约为99mm~101mm，现取d=100mm，转化为像素作为单位，得

又这些已知数据，利用先前所求得的H，r，L以及光源映射点的坐标为：

由此，可以知道光源点在该坐标系中的坐标为(1424.109253, 815.0699336, 7764.118028)，光源点距离探测平板的垂直距离约为986.042mm。

旋转载物台地半径为667.072664像素(约合84.718mm)，旋转载物台地旋转中心轴距离探测平面的垂直距离为H-L=3884.64567像素（约合493.35mm），旋转中心轴大约位于光源点与平板探测器的中间位置。

五、误差分析

讨论该方法的精度即是考虑该方法标定下带来的误差，误差来源于很多方面，大致归结为以下几类：

1.模型误差

忽略光源物理尺寸以及假设载物平台完全水平是我们为了简化问题而进行的假设。在这种假设下，小球旋转过程中在探测平板上出现的图像将是较为严格的椭圆，并且以平板的边缘横向和纵向建立平面坐标系，计算拟合出的椭圆曲线方程将没有形式上的交叉项。而实际在CT系统的安装过程中肯定不能保证载物平台的完全水平。当载物平台不是完全水平状态时，对每个小球而言，平板探测器上得到的投影将不再是单个椭圆，这将会对曲线拟合的过程带来更大的难度。同时，载物平台不水平将会影响后续具体计算参数时建立的几何平面模型。

2.计算误差

此误差有计算软件的精度引起，但此误差较小，可以忽略不计。

3.方法误差

此误差为参数标定中的主要误差，有我们计算过程中选定的方法，建立的计算模型引起，下面将对此进行分析。

椭圆曲线拟合引起的误差。

在使用给定数据拟合小球投影曲线时，我们采用了最小二乘法。拟合的理想

结果使所有的离散数据点尽可能接近你和得到的曲线，但显然仍然有很多落于曲

线两侧。，假设落在平面上的点为，待拟合在曲线上的点的坐标，则拟合误差为：

计算拟合得到的椭圆的参数时引起的误差。

在我们的计算模型中，需要用到椭圆的长轴长度。而我们假设点光源与旋转载物台的中轴线所在的面为探测器平板的垂面，那么对于两个小球而言，他们的投影轨迹所在的椭圆曲线的长轴应该是相等的。但是由于两个拟合得到的椭圆曲线的长度并不一致，故我们将两个曲线的长轴取平均数得到平均的长轴长度作为计算用长度。这样的计算过程肯定存在误差，但取平均数能减少数据误差。

计算误差。

由离散点拟合得到的椭圆曲线并不能显示地得到长短轴信息，需要转换成为标准的椭圆曲线形式才能得到需要的数据，由于在曲线方程的转换过程中求得长短轴数据时会丢失精度，故也会产生计算误差。

在使用平面几何计算点光源的空间坐标以及平台的中轴线位置坐标过程中，由于探测平板上给出的单位为像素，而其他数据单位为毫米，故需要进行单位统一化转换，在单位统一化转换过程中会丢失数据精度，从而产生计算误差。

在计算各参数过程中，两钢球之间距离我们选择中间值100mm作为已知数据，实际CT系统参数标定过程中，钢球之间的距离并不是一个十分精确的值，由于这个值是测量产生，故会产生测量误差。我们在计算过程中使用这个值将会是结果得到的参数产生误差。下面给出由于两钢球之间距离引起最终参数的误差公式,其中误差传递公式如下:

根据此公式，可知，已知参数d对于计算的到参数H无影响，而对于计算得到的参数r和L的影响如下：

其中根据已知条件d的绝对误差限是，为1mm，代入即可得到由于

d的绝对误差限引起的r和L的绝对误差限各为：

ε(r) = 0.8471822833mm

ε(L) = 4.926929895mm

六、模型的评价与改进

优点：

1.空间几何平面化，降低了问题的难度。

2.曲线拟合时，利用cholesy分解求解对称正定矩阵方程，大大简化了计算量，

并有较好的可靠性。

缺点：

1.假设太多，考虑的因素还不够多，模型要求太苛刻。

2.曲线拟合只采用了一种方法。

3.误差控制还不够好。

七、参考文献

[1] 石振东，误差理论与曲线拟合，哈尔滨，哈尔冰工程大学出版社，2010.1

[2] 张朝宗，工业CT技术和原理，北京，科学出版社，2009.6

[3] 刘士庆，基于给定轮廓线的散乱数据点的曲线拟合方法及其应用，[D]，2009.