DOI:10.12677/AAM.2018.74049,PDF,HTML,XML,被引量下载: 1,882浏览: 4,027科研立项经费支持
作者:李好阳,罗 政,敬端阳,吴宜均^*,左连翠：天津师范大学数学科学学院，天津
关键词:CT系统；三次样条插值；Radon变换；物体重构；CT System；Cubic Spline Interpolation；Radon Transform；Object Reconstruction

文章引用：李好阳, 罗政, 敬端阳, 吴宜均, 左连翠. 关于CT系统参数的标定及成像的数学模型[J]. 应用数学进展, 2018, 7(4): 388-400. https://doi.org/10.12677/AAM.2018.74049

1. 引言

CT (Computed Tomography)技术是20世纪50至70年代通过核物理、核医学等领域的一系列研究和实验发明的，它有效地克服了传统的X射线呈现的图像虽具有一定的分辨率，但仍不够清晰的缺点 [1] 。X射线的发射器和探测器相对位置固定不变，整个发射–接收系统绕某固定的旋转中心逆时针旋转180次。对每一个X射线方向，在具有512个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信息，如图1所示。

CT系统安装时往往存在误差，从而影响成像质量，因此需要对安装好的CT系统进行参数标定，即借助于已知结构的样品(称为模板)标定CT系统的参数，并据此对未知结构的样品进行成像。

2. 模型假设

1) 假设X射线穿过空气介质时能量不衰减；

2) 不考虑光的反射、散射、衍射等光学现象对X射线的影响；

3) 不考虑噪声对X射线传播的影响。

3. 问题重述

我们通过建立相应的数学模型和算法，解决以下问题：

1) 在正方形托盘上放置两个均匀固体介质组成的标定模板，模板的几何信息如图2所示，其中每一点对应的数值反映了该点的吸收强度，称为“吸收率”。根据这一模板及其接收信息，确定CT系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该CT系统使用的X射线的180个方向。

2) 利用上述CT系统得到某未知介质的接收信息的数据。利用(1)中得到的标定参数，确定该未知介质在正方形托盘中的位置、几何形状和吸收率等信息。另外，请具体给出图3所给的10个位置处的吸收率。

3) 利用上述CT系统得到另一个未知介质的接收信息的数据。利用(1)中得到的标定参数，给出该未知介质的相关信息。另外，请具体给出图3所给的10个位置处的吸收率。

4) 分析(1)中参数标定的精度和稳定性。在此基础上自行设计新模板、建立对应的标定模型，以改进标定精度和稳定性，并说明理由。

4. 模型的建立与求解

4.1. 问题1的求解

在正方形托盘上放置两个均匀固体介质组成的标定模板，模板的几何信息如图2所示，根据这一模板及其接收信息，确定CT系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该CT系统使用的X射线的180个方向。

此题用三次样条插值求解 [2] [3] [4] ，其相关定义如下：

设在区间 $\left[a,b\right]$上有如下分割 $\Delta$： $a=x_1\le x_2\le \cdots \le x_n=b$，如果函数 $f\left(x\right)$满足如下两个条件：

1) $f\left(x\right)$具有连续的二阶导数，即 $f\left(x\right)$， $f^{\prime }\left(x\right)$和 $f^{″}\left(x\right)$均在区间 $\left[a,b\right]$上连续；

2) $f\left(x\right)$在每一个小区间 $\left[x_j,x_{j+1}\right]$( $j=0,1,2,\cdots ,n$)是一个三次多项式，则称 $f\left(x\right)$是一个三次样条函数。

若对于一个函数 $S\left(x\right)$，上述三次样条函数 $f\left(x\right)$的函数坐标 $\left(x_i,f\left(x_i\right)\right)$满足 $f\left(x_i\right)=S\left(x_i\right)$( $i=0,1,2,\cdots ,n$)，则称 $f\left(x\right)$为 $S\left(x\right)$关于分割 $\Delta$的一个三次样条插值函数。

三次样条插值适用于分段插值生成的曲线在相邻区间的连接处也光滑的情形，我们根据EU列的数据绘制出分段插值曲线，如图4所示。

此题分段插值产生的曲线恰满足上述样条插值的特征，因而可以使用三次样条插值的方法来确定：

1) 找出探测器接受信息中所有数据中的最大值所对应的列。

2) 确定X射线经过小圆介质在探测器上的垂直投影所对应的探测器单元为46~74，以万分之一为步长，进行三次样条插值。

3) 找出三次样条插值函数所对应的最大值和迭代的次数。

4) 得出小圆的圆心在探测器上的垂直投影所对应的探测器编号x

x = 46 + (迭代次数 − 1) ´ 步长

5) 确定X光经过椭圆介质在探测器上的垂直投影所对应的探测器单元为169~277，以万分之一为步长，进行三次样条插值。

6) 找出三次样条插值函数所对应的最大值和迭代的次数。

7) 得出椭圆的中心在探测器上的垂直投影所对应的探测器编号y

y = 169 + (迭代次数 − 1) ´ 步长

8) 已知小圆的圆心和椭圆的中心之间的距离是45 mm，计算探测器两相邻探测单元的距离d，为 $d=\frac{45}{x-y}$，如图5所示。

9) 计算CT系统旋转中心的横坐标 $x_c$

$x_c=50-\left(256.5-x\right)\times \frac{45}{x-y}$

10) 顺时针旋转90˚，确定X射线经过椭圆介质在探测器上的垂直投影所对应的探测器单元编号为90~378，以万分之一为步长，进行三次样条插值。

11) 找出三次样条插值函数所对应的最大值和迭代的次数，以得到椭圆中心在探测器上的垂直投影所对应的探测器单元编号z

z = 90 + (迭代次数 − 1) ´ 步长

12) 计算CT系统旋转中心的纵坐标 $y_c$

$y_c=50+\left(256.5-z\right)\times \frac{45}{x-y}$

我们用MATLAB求得 $x=60.2822$， $y=222.8767$， $d=0.2768$， $x_c=40.6943$， $z=235.0094$， $y_c=55.9478$。

由模板示意图和X射线的能量衰减规律并通过计算可知，椭圆的长轴为80 mm，而椭圆的短轴长与圆的直径的和为38 mm。通过比较，当X射线的方向与椭圆的长轴方向平行时，探测器的接收信息值最大。在第一个模板的接收数据中，找出所有接收信息中的最大值为141.7794，它位于第233行EU列(即151列)。

CT系统是均匀旋转的，所以我们可以认为该CT系统每一次旋转的角度为1˚。如图6，我们定义CT系统使用的X射线的初始方向为θ。第一个模板的接收数据中数据的最大值位于EU列(即第151列)，该列数据可视为X射线从初始位置逆时针旋转150次(即150˚)得到的。因而将图7中X射线的方向顺时针旋转150˚后就是我们所要得到的该CT系统的初始位置。

如图6，记椭圆长轴的下半轴为0˚，则X射线初始位置的 $\theta ={30}^{\circ }$，该CT系统使用的X射线的180个方向 ${\theta }_n\left(0<n\le 180\right)$为

${\theta }_n={\left(30+n\right)}^{\circ },0<n\le 180$

4.2. 问题2的求解

由物理学中X射线的能量衰减规律知：X射线在穿过均匀介质时其能量的衰减率与能量本身的关系式为：

$I=I_0{\text{e}}^{-{\mu }_{0}d}$

在穿过非均匀介质时X射线的能量衰减率与能量本身的关系式为：

$I\left(x\right)=I_0{\text{e}}^{-{\displaystyle {\int }_L\mu \left(x\right)\text{d}x}}$

4.2.1. 增益处理的探究

标定模板是由两个均匀固体介质组成的，X射线穿过它的长度d与探测器接收到的原数据是一一对应的，但是，探测器接收到的原数据与实际接收到的数据并不相等，由题目知实际的接收信息是探测器接收到的原数据经过增益等处理后得到的，我们通过函数拟合的方法给出这种增益处理的量化表达。

1) 计算X射线被标定模板所截长度 $d_i$

问题1中的标定模板是均匀固体介质，且由第1个模板的接收数据知在有介质的位置X射线吸收率 ${\mu }_0=1$，在无介质的位置X射线吸收率 ${\mu }_0=0$。相对于椭圆而言，计算X射线被圆所截的长度的误差更小且更容易计算。因此我们计算X射线被圆所截的长度，即圆的弦长 $d_i$。

在某个方向上，通过确定圆心在探测器上的垂直投影的位置编号x(由三次样条插值法得到)来判断穿过圆的X射线对应的探测器单元的位置 $x_i$，并且利用圆的平分弦定理计算探测器单元所在直线被圆所截的弦长 $d_i$，其计算公式如下：

$d_i=2\sqrt{r^2-{\left[\left(x_i-x\right)\times d\right]}^2}$

其中r为圆形介质的半径，本题中 $r=2$，d为相邻两个探测器之间的距离，本题中 $d=0.2768$。在第一个模板的接收数据中随机找一列如EU列，分别计算每个 $x_{1i}\left(46\le i\le 74\right)$对应的 $d_{1i}$。再通过计算每个方向上穿过圆的所有X射线对应的探测器接收点相距最远两个的接收信息的差值，筛选出其中最小的列为第DX列，并计算该列中 $x_{2i}\left(46\le i\le 74\right)$对应的 $d_{2i}$。

2) 寻找 $d_i$与探测器的接收信息 $y_i$之间的关系

将第一个模板的接收数据中EU列和DX列中满足 $46\le i\le 74$的 $y_{1i}$和 $y_{2i}$分别与(1)中得到的 $d_{1i}$， $d_{2i}$进行函数拟合，分别如图8、图9所示：

由图可得， $d_i$与探测器的接收信息 $y_i$存在明显的线性关系，且函数图像的截距很小，对函数的影响忽略不计。因为此时介质的吸收率为1，从而穿过的均匀介质长度等于探测器上得到的原数据。从而，X射线穿过均匀介质后探测器上得到的原数据与经过增益等处理后得到的接收信息成近似正比关系。

radon变换及其逆变换具有线性 [5] ，从而对接收信息进行radon变换的逆变换得到的值与对未处理的原数据进行radon变换的变换后得到的值作比，其比值为一个常数，该常数为拟合图像的斜率，近似于2。

3) 计算问题2中X射线穿过某未知介质后探测器上得到的原数据

利用(2)中的结论，X射线穿过介质后探测器上得到的原数据与经过增益等处理后得到的接收信息成近似正比关系，我们可以推算出问题二中的接收信息所对应的原数据。

4.2.2. 对未知介质的重构

radon变换和其逆变换被广泛应用于医学成像、光学和全系干涉量度学等领域中进行图像分析和信号重构，其原理如下 [6] ：

1) 对椭圆介质进行radon变换得到一个与介质中心和旋转角度有关的函数

二维的radon变换表示为如下的线积分：

$p_{\varphi }\left(\rho \right)={\displaystyle {\int }_{L_{\varphi ,\rho }}f\left(x,y\right)\text{d}s},\rho \in \left[-R,R\right],\varphi \in \left[0,\text{π}\right]$(2.1)

其中

$L_{\varphi ,\rho }=\left\{\left(x,y\right)|x\cos \varphi +y\sin \varphi =\rho \right\}$(2.2)

是在平面上，将直线L进行参数化后得到的表达式。ρ是原点到直线L的距离，f是直线L与x轴正半轴的夹角。

$f\left(x,y\right)$表示均匀椭圆介质的吸收率，其满足如下条件：

$f\left(x,y\right)=\{\begin{array}{l}c,在介质内\\ 0,在介质外\end{array}$

$f\left(x,y\right)=\{\begin{array}{l}c,在介质内\\ 0,在介质外\end{array}$(2.3)

当椭圆介质的中心记为坐标原点 $\left(x_1=y_1=0\right)$，旋转角度在初始位置时，将(2.3)和(2.2)式代入(2.1)式可得

$p_{\varphi }\left(\rho \right)=\{\begin{array}{l}\frac{2cAB}{a^2\left(\varphi \right)}\sqrt{a^2\left(\varphi \right)-{\rho }^2},\left|\rho \right|\le a\left(\varphi \right)\\ 0,其他\end{array}$

当且仅当 $a^2\left(\varphi \right)=A^2{\cos }^2\varphi +B^2{\sin }^2\varphi$时等号成立。(其中A、B分别为椭圆介质的半长轴与半短轴的长度。)

如果中心为 $\left(x_1,y_1\right)$且旋转α度，同样地将将(2.3)式和(2.2)式代入(2.1)式可得

$p_{\varphi }\left(\rho \right)=p_{\varphi -\alpha }\left(\rho -s\cos cr-\varphi \right),s=\sqrt{x_1^2+y_1^2},r={\tan }^{-1}\left(\frac{y_1}{x_1}\right)$

引入狄拉克函数 $\delta \left(x,y\right)$，(2.1)式等价于

$p_{\varphi }\left(\rho \right)={\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }f\left(x,y\right)\delta \left(x\cos \varphi +y\sin \varphi -\rho \right)\text{d}x\text{d}y}}$

将其与 $I_{\varphi ,\rho }=I_0{\text{e}}^{-{\displaystyle {\int }_{I_{\varphi ,\rho }}\mu \left(x,y\right)\text{d}s}}$联立可解得:

$f\left(x,y\right)=\mu \left(x,y\right)$, $p_{\varphi }\left(\rho \right)=-\log \frac{I_{\varphi ,\rho }}{I_0}$,

即：

$p_{\varphi }\left(\rho \right)={\displaystyle {\int }_{L_{\varphi ,\rho }}\mu \left(x,y\right)\text{d}s},\rho \in \left[-R,R\right],\varphi \in \left[0^{\circ },\text{π}\right]$

2) 将上述radon变换得到的函数进行一维傅里叶变换，得到二维傅里叶变换的一个切片。

由傅立叶切片定理知：某图像函数 $\mu \left(x,y\right)$在某一方向上的投影函数 $p_{\varphi }\left(\rho \right)$的关于ρ的一维傅立叶变换，给出 $\mu \left(x,y\right)$的二维傅立叶变换 $F\left(u,v\right)$的一个切片，该切片与u轴相交成f角，且通过坐标原点。具体表达式为：

${\hat{p}}_{\varphi }\left(\omega \right)={\displaystyle {\int }_{-\infty }^{\infty }{p}_{\varphi }\left(\rho \right){\text{e}}^{-i2\text{π}\omega \rho }\text{d}\rho }={F\left(u,v\right)|}_{u=\omega \cos \varphi ,v=\omega \sin \varphi }$

3) 在全部180个方向上进行上述变换，拼接起来构成一个立体图形。

对平面上180个不同的f对应的 $\left(u,v\right)=\left(\omega \cos \varphi ,\omega \sin \varphi \right)$在二维傅立叶空间上，关于u轴进行旋转形成的截面进行拼接重构得到一个立体图形。

即 ${\hat{p}}_{\varphi }\left(\omega \right)=F\left(\omega \cos \varphi ,\omega \sin \varphi \right)$, $\varphi \in \left(0,\text{π}\right)$, $\omega \in \left(-\infty ,+\infty \right)$

4) 通过radon变换的逆变换，并添加色阶，得到色阶图。

radon变换的逆变换还原出介质的几何形状，通过添加色阶，得到一个色阶图。其中，radon变换的逆变换公式为：

$\mu \left(x,y\right)=\frac{1}{4{\text{π}}^2}{\displaystyle {\int }_0^{\text{2π}}\text{d}\varphi }{\displaystyle {\int }_{-\infty }^{+\infty }\frac{\frac{\partial }{\partial \rho }{p}_{\varphi }\left(\rho \right)}{\rho -\left(x\cos \varphi +y\sin \varphi \right)}\text{d}\rho }$

在实际的操作过程中，我们对探测器接收的数据进行radon变换的逆变换，发现得出的图形并不与问题一中的初始方向一致，如果我们强制将其角度旋转30˚，会造成数据的丢失(圆形介质区域的接收信息丢失)。但是，第一个模板和第二个模板的接收数据的矩阵边缘部分均为0。因此，我们可以通过添加零元素来扩大第二个模板的接收数据的矩阵的维数，再对矩阵进行切割，以此来避免数据的丢失。具体求解过程如下：

1) 将第二个模板的接收数据每一列增加100个零，共计180列，得到一个新的矩阵。

2) 对新的矩阵进行radon变换的逆变换，顺时针旋转30˚，得到的图形即为我们所要的初始方向。

3) 对2)所得的矩阵进行切割，并使切割之后的矩阵成为一个 $256\times 256$的矩阵，根据第1个模板和第2个模板的接收数据之间的关系，对矩阵的数据进行增益处理，得到一个新的矩阵，该矩阵中的数据即为我们所求的吸收率。

4) 由于标定模板是 $100\times 100$的，所以我们需要对给定的十个坐标数据扩大2.56倍，并进行相应的坐标变换，使得矩阵的(0，0)位置与我们之前建立坐标系的原点重合，进而求得给定的十个坐标数据以及所给的10个位置的吸收率。

基于上述四步，我们求得的结果如下：

a) 该未知介质的几何形状及在正方形托盘中的位置如图10所示。

b) 该未知介质在图3所给的10个位置处的吸收率分别是(见表1)。

4.3. 问题3的求解

与问题2的解题思路基本一致，不同的是两个介质的接收信息。

1) 该未知介质的几何形状及在正方形托盘中的位置如图11所示(见表2)。

2) 该未知介质在图3所给的10个位置处的吸收率分别是(见表2)。

4.4. 问题4的求解

CT的扫描、重建、显示等参数对图像质量、机器负荷、辐射剂量的影响存在着相互依赖又相互制约的关系 [7] [8] 。以下几个因素会对参数标定的精度和稳定性产生影响：

1) 模板具有多条对称轴，以确保能找出多组容易比较的数据进行分析。

2) 模板所占的面积应该足够多，以确保能使得数据经过重构处理后的图像更趋近于光滑，从而减小边缘的误差。

3) 模板应该具有一定的弧度，以确保能使得数据有起伏，具有区分度。

为了更好的改进参数标定的精度和稳定性，我们设计了这个模板，以正方形托盘左下角为原点建立直角坐标系，模版以 $\left(50,50\right)$为中心的三个圆(半径分别为3，4，5 mm)呈正三角形排放。记X射线与x轴的夹角为θ。探测器从与x轴垂直时开始向逆时针方向旋转。当探测器每旋转60˚时，其检测到的数据会出现容易比较的情况，呈现一定的变化规律。

在问题一中，方向的确定仅仅依靠在接收信息中找出最大值所在的角度，以此角度当作垂直正方形托盘过长轴的角度，从而逆推回探测器的初始方向。而当容易提取的数值所对应的角度多了，可以通过多组数据来确定初始方向，这样可以减少误差，得到更精确的初始角度，使数据更接近准确值。

圆可以使角度变化时所对应的探测器的接收信息呈现出连续的趋势，这样可以在radon逆变换后所得的边缘处的数值偏差降低。多个圆可以增大所占面积，从而得到充分多的数据。而呈正三角形的排放可以使得它具有多条对称轴。如图12所示。这些都改进了标定的精度和稳定性。

资助信息

国家级大学生创新创业训练计划项目资助(项目号：201710065038)。

致谢

感谢天津师范大学教务处的老师们对我们的大力支持。

NOTES

^*通讯作者。

参考文献

[1]	姜启源, 谢金星, 叶俊. 数学模型[M]. 第四版. 北京: 高等教育出版社, 2011: 187-191.
[2]	三次样条差值[EB/OL].https://wenku.baidu.com/view/c8171f53ad02de80d4d84064.html, 2017-09-17.
[3]	许小勇, 钟太勇. 三次样条差值函数的构造与Matlab实现[J]. 自动测量与控制, 2006, 25(11): 76-78.
[4]	朱立勋. 三次样条插值的收敛性及一类三次广义样条插值的误差估计[D]: [硕士学位论文]. 长春: 吉林大学数学研究所, 2006.
[5]	李静. 二维Radon变换在图像重建中的重要性质及定理[J]. 河南教育学院学报, 2014, 24(4): 21-23.
[6]	Training School: Scientific Computing for X-Ray Computed Tomography (CT) [EB/OL].http://www2.compute.dtu.dk/~pcha/HDtomo/SCforCT.html, 2017-09-17.
[7]	王甜甜. CT图像质量评估及测试体模研究[D]: [硕士学位论文]. 广州: 南方医科大学, 2012.
[8]	刘景鑫, 杨海山. 优化CT系统参数设置的可行性探讨[J]. 中华放射学杂志, 1996, 30(1): 62-64.