CN102743184B - 一种x射线锥束计算机层析成像系统的几何参数标定方法 - Google Patents

Abstract

本发明是一种X射线锥束计算机层析成像系统的几何参数标定方法，属于X射线锥形束计算机层析成像系统（简称锥束CT）领域，本发明通过相机标定技术来复现锥束CT系统中的X射线源、平板探测器和转轴之间的几何位置关系，从而对锥束CT的几何参数及其误差进行直接求解，是一种系统化的测量求解方法，可以将锥束CT抽象为基本的相机系统模型，从而同时求解系统中相互关联的多个几何参数，与现有标定方法相比，不需要逐个测量几何参数和逐个参数地反复调整，这样就大大简化了整个系统的标定测量的复杂程度，同时也提高了锥束CT几何标定的稳定性。

Description

一种X射线锥束计算机层析成像系统的几何参数标定方法

技术领域

本发明涉及一种X线锥形束计算机层析成像系统(Computerized Tomography,简称锥束CT)的标定方法，特别是指X射线源与平板探测器的相对位置固定，X射线源与平板探测器构成的成像子系统绕固定的旋转轴旋转，或者被测样本放在一个中心转台上绕固定的转轴旋转，而成像子系统固定不动，从而通过成像子系统获取各个角度的X线投影数据的小型或微型锥束CT系统几何参数的标定方法。

背景技术

计算机层析成像系统是利用已知的几何条件，通过多角度投影数据重建断层中的密度或衰减率的成像方法，在医学临床和工业检测中已经得到广泛应用。锥束CT是指利用2维的平板探测器进行X线的投影数据检测，并基于平板探测器收集的多个投影数据集进行三维重建的一种新型CT，目前其重建理论和算法都还在发展之中。几何参数和误差的准确标定对于锥束CT的精确重建至关重要。对于锥束CT的标定，已有学者提出了几种方法，这些方法可以分为两类：（1）迭代求解法：这类方法通过锥束CT成像的几何关系创建关于待求解几何参数的代价函数，利用迭代方法使代价函数达到最小点，从而得到几何参数。这种方法的问题是容易陷入局部最小点，而难以得到全局最优解，因此需要一个良好的初始值，这通常是比较难得到的。（2）解析求解法：理论上绕固定转轴旋转的空间点在平板探测器上投影的轨迹为一个椭圆，这类方法根据椭圆的参数与锥束CT几何参数之间的间接几何关系来进行求解。上述两类方法都需要假设某些参数是理想的或其误差是可以忽略的，从而只需要对其余的部分几何参数进行求解。同时，目前的标定方法都是对各个几何参数进行逐次的标定和调整，由于参数之间的关联性，往往需要逐个参数反复调整，使标定过程十分复杂。本发明通过将锥束CT系统考虑成一个计算机视觉系统，并利用计算机视觉理论进行求解，从而得到一种用于锥束CT几何参数和几何误差的系统性标定方法。

发明内容

技术问题：本发明要解决的技术问题是提出了一种基于计算机视觉理论，利用相机标定技术对锥束CT的几何参数和误差进行标定的方法。通过相机标定技术来复现锥束CT系统中的X射线源、平板探测器和转轴之间的几何位置关系，从而对锥束CT的几何参数及其误差进行直接求解。本发明是一种系统化的测量求解方法，可以将锥束CT抽象为基本的相机系统模型，从而同时求解系统中相互关联的多个几何参数，与现有标定方法相比，不需要逐个测量几何参数和逐个参数地反复调整，这样就大大简化了整个系统的标定测量的复杂程度，同时也提高了锥束CT几何标定的稳定性。

其特征在于，是在一个由平板探测器、同轴地安置有被检测物的转台、X射线源以及计算机共同构成的锥束CT系统中实现的，是一种利用相机标定技术来确定所述锥束CT系统中的所述X射线源、平板探测器和转台中的转轴之间的几何位置关系，以对所述锥束CT系统的几何参数及其误差进行直接求解的方法，步骤如下：

步骤（1），制作用于模拟所述被检测物的三维标定模块：

在所述三维标定模块上有不少于6个不全共面或不全共线的空间标志点，这里使用12个空间标志点，所述三维标定模块采用有机玻璃基体，外表面精确钻孔并镶嵌小钢珠作为空间标志点；

步骤（2），建立如下三个直角坐标系：

成像子系统坐标系，是一个所述的X射线源-平板探测器坐标系，以下称成像系，设在所述平板探测器上，用（X_C、Y_C、Z_C）表示，坐标原点在所述平板探测器的左下角，表示为(x_c0,y_c0,z_c0)，x_c0=y_c0=z_c0=0，X_C轴、Y_C轴分别对应于所述平板探测器上像素排布的方向，Z_C轴是X_C轴、Y_C轴的叉积，并垂直于探测器平面，从所述X射线源到所述平板探测器引垂线，垂足为C_C，X射线源到垂足C_C的距离为D_SD，单位为mm，在所述成像系中：

X射线源的坐标点为S_C，坐标为(x_cs,y_cs,z_cs)，z_cs=D_SD，x_cs=y_cs≠0，

垂足C_C的坐标为(x_cc,y_cc,z_cc)，z_cc=0，x_cc=x_cs，y_cc=y_cs，

在所述平板探测器的各个投影图像内，当把以mm为单位的坐标(x_ci,y_ci)转换为以像素为单位的坐标(u_ci,v_ci)时，把以像素为单位的坐标的原点设在投影图像左上角，用(u_c0,v_c0)表示，且坐标轴u向右为正方向，坐标轴v向下为正方向，则：

$u_{\mathrm{ci}}=\frac{x_{\mathrm{ci}}}{\mathrm{dx}},$ $v_{\mathrm{ci}}=\frac{D-y_{\mathrm{ci}}}{\mathrm{dy}},$ 其中：

i表示在所述投影图像内的点的序号，

D为所述平板探测器在Y_C轴上的长度，

dx，dy分别代表单个像素在X_C、Y_C方向上的尺寸，

旋转子系统坐标系，以下称旋转系，用（X_R、Y_R、Z_R）表示，主轴线为Y_R轴，与所述转台的旋转轴线重合，Z_R轴为由所述X射线源向所述转轴所引的垂线，并且指向所述X射线源，X_R轴是Y_R轴和Z_R轴的叉积，

标定物坐标系，以下称标定物系，用（X、Y、Z）表示，是一个同轴地设置在所述三维标定模块上，用于标定锥束CT参数的辅助坐标系，所述X射线源在所述标定物系中的坐标为(x_s,y_s,z_s)，所述三维标定模块上的12个空间标志点用P_i表示，坐标为(x_i,y_i,z_i)，是已知的，i=1,2,…,I，I=12，以便通过所述成像子系统对所述三维标定模块进行拍摄，从而在所述平板探测器上得到对应于12个所述空间标志点的投影位置，坐标点用P_ci表示，坐标用(x_ci,y_ci)表示；

步骤（3），把所述三维标定模块与所述转台在纵向同轴安置，在所述X射线源-平板探测器相对静止，转台转动的条件下，对所述三维标定模块成像：所述计算机启动锥束CT采集程序，获取所述平板探测器上不少于360个投影角度的投影图像，用j表示所述投影图像序号，j=1,2,L,J，J≥360，在获取到的各个所述投影图像j中抽取12个所述空间标志点的坐标点，用P_cji表示，坐标为(x_cji,y_cji)，转换为像素坐标时用(u_cji,v_cji)表示，对应的所述三维标定模块上的坐标点用P_i表示，坐标用(x_i,y_i,z_i)表示；

步骤（4），按以下步骤从所述标定物系的空间标志点P_i(x_i,y_i,z_i)和成像系中对应的像素坐标P_cji(u_cji,v_cji)中求解所述锥束CT系统的以下参数：D_SD、x_cs和y_cs，标定物系到成像系的空间转换矩阵 $T_j=\left[\begin{array}{cc}{R}_j& t_j\\ 0^T& 1\end{array}\right]$ 中的各所述投影图像的旋转矩阵R_j和位移向量t_j，其中，D_SD、x_cs和y_cs是所述锥束CT系统的内方位参数，在所述成像子系统制作完成后是不变的，但又是未知的，而对于一个设定的拍摄位置，对应于该拍摄位置的空间转换矩阵T是不变的：

步骤（4.1），按下式建立所述标定物系内已知空间坐标点P_i(x_i,y_i,z_i)与各所述投影图像j中的对应坐标点P_cji(u_cji,v_cji)之间的关系：

$z_{\mathrm{ci}}^{\′}\left[\begin{array}{c}{u}_{\mathrm{cji}}\\ v_{\mathrm{cji}}\\ 1\end{array}\right]=H_{(3\×4)j}\left[\begin{array}{c}{x}_i\\ y_i\\ z_i\\ 1\end{array}\right],$

H_(3×4)j为所述以像素为单位的坐标与标定物系坐标之间的转换关系，

$H_{(3\×4)j}=\left[\begin{array}{cccc}{D}_{\mathrm{SD}}/\mathrm{dx}& 0& x_{\mathrm{cs}}/\mathrm{dx}& {-D}_{\mathrm{SD}}\·x_{\mathrm{cs}}/\mathrm{dx}\\ 0& {-D}_{\mathrm{SD}}/\mathrm{dx}& (D-y_{\mathrm{cs}})/\mathrm{dy}& D_{\mathrm{SD}}\·y_{\mathrm{cs}}/\mathrm{dy}\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}{R}_j& t_j\\ 0^T& 1\end{array}\right]$

$=\left[\begin{array}{cccc}{h}_{j11}& h_{j12}& h_{j13}& h_{j14}\\ h_{j21}& h_{j22}& h_{j23}& h_{j24}\\ h_{j31}& h_{j32}& h_{j33}& h_{j34}\end{array}\right]$

H_(3×4)j为第j幅投影图像的变换矩阵，

R_j为第j幅投影图像的旋转矩阵，是单位正交矩阵，

t_j为第j幅投影图像的位移向量，

$\left[\begin{array}{cc}{R}_j& t_j\\ 0^T& 1\end{array}\right]$ 为第j幅投影图像的空间转换矩阵，

对于所述三维标定模块的I=12个空间标志点，得到下述对应于第j个投影图像的，包含2I个方程的联立方程组：

$\left[\begin{array}{ccccccccccc}{x}_1& y_1& z_1& 1& 0& 0& 0& 0& {-u}_{\mathrm{cjl}}{x}_1& {-u}_{\mathrm{cjl}}{y}_1& {-u}_{\mathrm{cjl}}{z}_1\\ 0& 0& 0& 0& x_1& y_1& z_1& 1& {-v}_{\mathrm{cj}1}{x}_1& {-v}_{\mathrm{cjl}}{y}_1& {-v}_{\mathrm{cjl}}{z}_1\\ & \·\·\·& \·\·\·& & & & & & & & \\ & \·\·\·& \·\·\·& & & & & & & & \\ x_I& y_I& z_I& 1& 0& 0& 0& 0& {-u}_{\mathrm{cjI}}{x}_I& {-u}_{\mathrm{cjI}}{y}_I& {-u}_{\mathrm{cjI}}{z}_I\\ 0& 0& 0& 0& x_I& y_I& z_I& 1& {-v}_{\mathrm{cjI}}{x}_I& {-v}_{\mathrm{cjI}}{y}_I& {-v}_{\mathrm{cjI}}{z}_I\end{array}\right]\left[\begin{array}{c}{h}_{j11}\\ h_{j12}\\ h_{j13}\\ h_{j14}\\ h_{j21}\\ h_{j22}\\ h_{j23}\\ h_{j24}\\ h_{j31}\\ h_{j32}\\ h_{j33}\end{array}\right]=\left[\begin{array}{c}{u}_{\mathrm{cj}1}{h}_{j34}\\ v_{\mathrm{cj}1}{h}_{j34}\\ \·\·\·\\ \·\·\·\\ u_{\mathrm{cjI}}{h}_{j34}\\ v_{\mathrm{cjI}}{h}_{j34}\end{array}\right]$

等式两端同除以h_j34，得到：

${\mathrm{Ah}}_j^{\′}=B$

式中： $A=\left[\begin{array}{ccccccccccc}{x}_1& y_1& z_1& 1& 0& 0& 0& 0& {-u}_{\mathrm{cj}1}{x}_1& {-u}_{\mathrm{cj}1}{y}_1& {-u}_{\mathrm{cj}1}{z}_{21}\\ 0& 0& 0& 0& x_1& y_1& z_1& 1& {-v}_{\mathrm{cj}1}{x}_1& {-v}_{\mathrm{cj}1}{y}_1& {-v}_{\mathrm{cj}1}{z}_1\\ & \·\·\·& \·\·\·& & & & & & & & \\ & \·\·\·& \·\·\·& & & & & & & & \\ x_I& y_I& z_I& 1& 0& 0& 0& 0& {-u}_{\mathrm{cjI}}{x}_I& {-u}_{\mathrm{cjI}}{y}_I& {-u}_{\mathrm{cjI}}{z}_I\\ 0& 0& 0& 0& x_I& y_I& z_I& 1& {-v}_{\mathrm{cjI}}{x}_I& {-v}_{\mathrm{cjI}}{x}_I& {-v}_{\mathrm{cjI}}{z}_I\end{array}\right],$ 是2I×11的矩阵，B=[u_cj1,v_cj1,…,u_cjI,v_cjI]^T是2I×1的向量，

$h_j={\left[\begin{array}{ccccccccccc}{h}_{j11}& h_{j12}& h_{j13}& h_{j14}& h_{j21}& h_{j22}& h_{j23}& h_{j24}& h_{j31}& h_{j32}& h_{j33}\end{array}\right]}^T;$

利用最小二乘法，得到：

$h_j^{\′}={\left(A^{T}A\right)}^{-1}{A}^{T}B$

步骤（4.2），按以下步骤计算表换矩阵H_(3×4)j中的各个元素：h_j1，h_j2，h_j3，h_j14，h_j24，h_j34，其中：h_j1=[h_j11,h_j12,h_j13]^T，h_j2=[h_j21,h_j22,h_j23]^T，h_j3=[h_j31,h_j32,h_j33]^T，

步骤（4.2.1），根据步骤（4.1），H_(3×4)j表达为：

$H_{(3\×4)j}=\left[\begin{array}{cc}{h}_{j1}^T& h_{j14}\\ h_{j2}^T& h_{j24}\\ h_{j3}^T& h_{j34}\end{array}\right]\left[\begin{array}{cccc}{a}_{11}& 0& a_{13}& a_{14}\\ 0& a_{22}& a_{23}& a_{24}\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}{r}_{j1}^T& t_{\mathrm{jx}}\\ r_{j2}^T& t_{\mathrm{jy}}\\ r_{j3}^T& t_{\mathrm{jz}}\\ 0^T& 1\end{array}\right]$

其中， $a_{11}=D_{\mathrm{SD}}/\mathrm{dx},$ $a_{13}=x_{\mathrm{cs}}/\mathrm{dx},$ $a_{14}=({-D}_{\mathrm{SD}}\·x_{\mathrm{cs}})/\mathrm{dx},$ $a_{22}={-D}_{\mathrm{SD}}/\mathrm{dy},$ $a_{23}=(D-y_{\mathrm{cs}})/\mathrm{dy},$ $a_{24}=(D_{\mathrm{SD}}\·y_{\mathrm{cs}})/\mathrm{dy},$

$R_j=\left[\begin{array}{c}{r_{j1}}^T\\ {r_{j2}}^T\\ {r_{j3}}^T\end{array}\right],$ 其中，r_j1=[r_j11,r_j12,r_j13]^T，r_j2=[r_j21,r_j22,r_j23]^T，r_j3=[r_j31,r_j32,r_j33]^T

t_j=[t_jx,t_jy,t_jz]^T

从而得到：

$H_{(3\×4)j}=\left[\begin{array}{cc}{a}_{11}{r}_{j1}^T+a_{13}{r}_{j3}^T& a_{11}{t}_{\mathrm{jx}}+a_{13}{t}_{\mathrm{jz}}+a_{14}\\ a_{22}{r}_{j2}^T+a_{23}{r}_{j3}^T& a_{22}{t}_{\mathrm{jy}}+a_{23}{t}_{\mathrm{jz}}+a_{24}\\ r_{j3}^T& t_{\mathrm{jz}}\end{array}\right]$

由于R_j是单位正交矩阵，则有：

$\left|\right|r_{\mathrm{jm}}\left|\right|=1$

$r_{\mathrm{jm}}^T{r}_{\mathrm{jn}}=0$

$r_{\mathrm{jm}}\×r_{\mathrm{jm}}=0$

$\left|\right|r_{\mathrm{jm}}\×r_{\mathrm{jn}}\left|\right|=1$

其中，n=1,2,3,m=1,2,3，且n≠m，

分别对应于R_j的第m行和第n行，r_jm和r_jn相互正交，并且均为单位向量，

由上所述有：

$\left|\right|r_{j3}\left|\right|=\left|\right|h_{j3}\left|\right|=h_{j34}\left|\right|h_{j3}^{\′}\left|\right|=1$

其中 $h_{j3}^{\′}=h_{j3}/h_{j34}$

故有： $h_{j34}=1/\left|\right|h_{j3}^{\′}\left|\right|$

由上得到： $h_j=h_{j34}{h}_j^{\′}$

至此，就得到了矩阵H_(3×4)j中的各个元素，也就得到了h_j1，h_j2，h_j3，h_j14，h_j24，h_j34；

步骤（4.3），由步骤（4.2）得到H_(3×4)j中的各个元素后，进一步求解参数D_SD，x_cs，y_cs，R_j，t_j，具体步骤如下：

$a_{13}=(a_{11}{r}_{j1}^T+a_{13}{r}_{j3}^T)\·r_{j3}=h_{j1}^T{h}_{j3}$

$a_{23}=(a_{22}{r}_{j2}^T+a_{23}{r}_{j3}^t)r_{j3}=h_{j2}^T{h}_{j3}$

$a_{11}=\left|\right|(a_{11}{r}_{j1}+a_{13}{r}_{j3})\×r_{j3}\left|\right|=\left|\right|h_{j1}\×h_{j3}\left|\right|$

$a_{22}=\left|\right|(a_{22}{r}_{j2}+a_{23}{r}_{j3})\×r_{j3}\left|\right|=\left|\right|h_{j2}\×h_{j3}\left|\right|$

从上可得：

$D_{\mathrm{SD}}=a_{11}\·\mathrm{dx}={-a}_{22}\·\mathrm{dy}$

$x_{\mathrm{cs}}=a_{13}\·\mathrm{dx}$

$y_{\mathrm{cs}}=D-a_{23}\·\mathrm{dy}$

$r_{j1}=\frac{1}{a_{11}}(h_{j1}-a_{13}{h}_{j3})$

$r_{j3}=\frac{1}{a_{22}}(h_{j2}-a_{23}{h}_{j3})$

$r_{j3}=h_{j3}$

$R_j=\left[\begin{array}{c}{r}_{j1}^T\\ r_{j3}^T\\ r_{j3}^T\end{array}\right]$

从上可得：

$a_{14}=-(D_{\mathrm{SD}}\·x_{\mathrm{cs}})/\mathrm{dx}$

$a_{24}=D_{\mathrm{SD}}\·y_{\mathrm{cs}}/\mathrm{dy}$

因此，tj求解如下：

$t_{\mathrm{jx}}=\frac{1}{a_{11}}(h_{14}-a_{14}-a_{13}\·h_{34})$

$t_{\mathrm{jy}}=\frac{1}{a_{22}}(h_{24}-a_{24}-a_{23}\·h_{34})$

$t_{\mathrm{jz}}=h_{34}$

$t_j={[t_{\mathrm{jx}},t_{\mathrm{jy}},t_{\mathrm{jz}}]}^T$

由上述关系，则求解出对应各个投影位置的参数D_SD，x_cs，y_cs，R_j，t_j；

步骤（5），求解锥束CT系统中所述转轴的位置和方向

利用步骤（4.3）中的参数，得到不同投影位置时标定物系的原点在成像系中的空间位置向量：

$S_j=\left[\begin{array}{cc}{R}_j& t_j\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}0\\ 0\\ 0\\ 1\end{array}\right]=t_j$

这些向量的端点轨迹是一个空间的圆周，因此通过求解S_j的均值，得到其圆心位置C，即所述转轴通过的位置：

$c=\frac{1}{j}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{S}_j$

锥束CT系统中所述转轴的方向，则为向量S_j端点的轨迹所拟合出平面的法向量，通过求圆周上三个点所构成任意两个向量的叉积而得到，本发明使用以下公式进行求解：

$Y=\frac{1}{J}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}[(S_{j+k}-S_j)\×(S_{j+1}-S_{j+k})]$

其中，k和l分别取为1/3和2/3的投影数目，

至此，确定了锥束CT系统中所述转轴的空间向量是过C点，方向为Y的向量；

步骤（6），求解成像子系统的投影主光轴，

投影主光轴是所述X射线源到所述平板探测器的垂线，在成像系中，投影主光轴的向量为O为：

$O=C_C-S_C$

其中C_C为所述X射线源在所述平板探测器上的垂足的齐次坐标，S_C为所述X射线源的空间位置的齐次坐标，S_C=[x_cs,y_cs,D_SD,1]^T，C_C=[x_cs,y_cs,0,1]^T；

步骤（7），锥束CT几何误差的求解，

锥束CT的理想几何条件为所述投影主光轴与所述转轴相交且垂直，并且所述转轴在所述平板探测器上的投影与探测器的Y_C方向平行，这里首先判断锥束CT是否符合上述理想条件，若不满足，其误差是多大，

步骤（7.1），判断所述投影主光轴与所述转轴是否垂直，若不垂直，其夹角是多大，

设a=O^TY，

若a=0，则所述投影主光轴和所述转轴垂直，满足理想情况，

若a≠0，则所述投影主光轴和所述转轴不垂直，其夹角为α：

$\mathrm{\α}={\cos }^{-1}\frac{O^{T}Y}{{\left|\right|O\left|\right|\left|\right|Y\left|\right|}^{,}}$

步骤（7.2）判断所述投影主光轴与所述转轴是否相交，若不相交，二者之间的距离是多大，

设b=O×Y，

若

则所述投影主光轴和所述转轴平行，即二者不相交，

若

则二者之间的距离d为：

$d={(C-S_C)}^{T}b$

若d=0，则二者共面且不平行，即二者相交，

若d≠0，则二者不共面，即不相交，二者之间的距离为d；

步骤（7.3），判断所述转轴在所述平板探测器上的投影是否与所述平板探测器的Y_C方向平行，若不平行，倾斜角是多少，

由步骤（5）可知，转轴上存在两点在成像系中的坐标为C和C+Y，其在探测器上的投影为

$P=\mathrm{LC}$

$P^{\′}=L(C+Y)$

其中， $L=\left[\begin{array}{cccc}{a}_{11}& 0& a_{13}& a_{14}\\ 0& a_{22}& a_{23}& a_{24}\\ 0& 0& 1& 0\end{array}\right],$ P(u_P,v_P)，P^′(u_P′,v_P′)，

所述转轴在所述平板探测器上投影的倾斜角为β：

$\mathrm{\β}={\tan }^{-1}\frac{u_P-u_{P^{\′}}}{v_P-v_{P^{\′}}}$

若β=0，则所述转轴在所述平板探测器上的投影与探测器的Y_C方向平行；

至此，就完成了所述锥束CT系统的几何参数的标定，所述几何参数包括D_SD，x_cs，y_cs，以及所述投影主光轴与所述转轴之间的几何误差。

本发明的效果在于利用计算机视觉理论和相机标定技术，通过建立一个合理的锥束CT系统的模型，将锥束CT系统抽象为成像子系统和旋转子系统，并将系统的几何参数及其误差表达为两者间的位置关系，从而可以直接求解锥束CT的几何参数和几何误差，方法简便，易于实施。由于本发明的方法在实施中不涉及复杂的测量和迭代算法，并且相机标定技术已经发展成为成熟、可靠的理论体系，因此本发明中方法易于实施，其稳定性很好地得到保证。

附图说明

图1是锥束CT系统的几何参数和误差的示意图；

图2是“源-探测器对”的示意图；

图3是理想锥束CT模型的示意图；

图4是理想的几何条件下锥束CT采集数据的示意图；

图5是源-探测器对存在俯仰误差时锥束CT采集数据的示意图；

图6是源-探测器对存在偏转或平移误差时锥束CT采集数据的示意图；

图7是源-探测器对存在旋转误差时，锥束CT采集数据的示意图；

图8是坐标系定义的示意图；

图9是锥束CT投影过程的示意图；

图10(x_ci,y_ci)与(u_ci,v_ci)之间的关系的示意图；

图11是一种标定模块的示意图；

图12是实施的软件流程图

图13是一种典型的锥束CT系统的示意图

具体实施方式

技术方案：

1)锥束CT系统模型的建立

一个典型的锥束CT系统包括X线源，平板探测器和旋转系统三大部分。在CT重建时利用三部分之间的几何关系，根据多个投影数据来重建被检测物体的密度或X线衰减系数的空间分布。传统的标定方法是逐个考虑系统的几何参数及其误差，并进行逐个测量和调整，使其安装误差接近为零，并得到精确的成像几何参数。如图1所示，系统的坐标系是以系统的台架平面为参考系来定义的，理想的几何条件是：X线源与平板探测器中心的连线通过坐标原点，并落在台架平面为参考的XZ坐标平面内。这样，系统的几何参数和误差将包括以下几类(参见图1)：

(1)平板探测器的安装误差：即平板探测器在X,Y,Z方向偏离理想位置的位移(δx,δy,δz,)，以及平板探测器绕轴线Z的旋转角度α、绕轴线X的俯仰角度β、和绕轴线Y的偏转角度γ；

(2)投影系统参数：即X射线源到平板探测器的距离D_SD、X射线源在平板探测器上的投影坐标(x_cs,y_cs)；

(3)转轴的倾斜，即转轴与台架平面不完全垂直，所造成的偏离理想位置的角度。

由于这些误差和参数是互相关联的，因此在测量调整中需要逐次反复调整，费时费力。

为了描述锥束CT中的几何参数和误差，将锥束CT系统看作成像子系统和旋转子系统两部分的组合。其中，成像子系统包括X射线源和平板探测器，旋转子系统包括转轴。由分析可知，锥束CT的几何参数和误差完全可以由两个子系统的相对位置关系来直接表示。通过复现两个子系统之间的相对位置关系，则可以对锥束CT进行标定。实际上，如果我们从计算机视觉的角度来分析锥束CT系统，我们完全可以先脱离上述实际的几何参数及其偏差，从系统角度来抽象出数学模型和简化问题。

在一个锥束CT系统中X线源与平板探测器的相对位置是固定不动的，因此，总可以通过X线源S向平板探测器引出一条垂线，它与平板探测器垂直相交，设其交点位置为(x_cs,y_cs)。我们可以把该垂线类比为光学相机系统中的主光轴，把平板探测器考虑为光学相机系统的光学传感器，则X线源的位置相当于相机系统中的光心位置，因此称该垂线为“投影主光轴”，X线源到平板探测器的距离相当于光学相机系统中的焦距，我们把符合这样条件的X线源与平板探测器共同考虑为一个固定的系统部件，称为“源-探测器对”。这样，上述锥束CT系统中的两个几何参数，即：X射线源到平板探测器的距离D_SD、投影中心位置(x_cs,y_cs)，分别相当于相机系统中的焦距和主点。通过计算机视觉理论中的相机标定方法可以进行检测和推算。由于从一个空间点一定可以向一个平面引出一条垂线，因此满足上述条件，即X射线源到平板探测器的距离为D_SD、投影中心位置为(x_cs,y_cs)的“源-探测器对”一定存在(图2)。

另一方面，锥束CT中的旋转系统有一个固定的旋转子系统，无论是源-探测器对绕样品旋转的探测模式，还是源-探测器对不动、被探测样品自身旋转的探测模式，都可以考虑为由源-探测器对构成的成像系统从不同角度对被测样品进行拍摄。因此可以利用相机标定方法对各个拍摄位置进行检测和推算，从而得到锥束CT的旋转子系统与源-探测器对子系统之间的几何关系，进而得到系统的几何参数及其误差。

如图3所示，当把源-探测器对子系统与旋转子系统放在一起考虑时，希望的空间关系是旋转子系统上的被测物完全被X射线源与平板探测器形成的四棱锥所包围，同时源-探测器对子系统的投影主光轴与旋转子系统的旋转轴线垂直相交，且平板探测器的Y_C轴(即像素排布方向)与旋转子系统的旋转轴线在平板探测器上的投影互相平行，若完全满足这些条件，则锥束CT系统成像的几何条件是理想的，若部分偏离了这些条件，则锥束CT系统是存在几何误差的，我们希望通过标定来解析其误差究竟有多大。

具体讲，一个理想的锥束CT模型需要且只需要同时满足下列二个条件：

(1)在成像的各个投影成像位置，源-探测器对子系统的投影主光轴都与旋转子系统的旋

转轴线严格垂直相交；

(2)平板探测器的坐标轴Y_C与旋转子系统的旋转轴线在平板探测器上的投影严格平行；

同时，我们需要通过标定得到系统的以下二个参数：

(1)在理想情况下，投影主光轴准确穿过平板探测器的中心位置，即投影中心位置

(x_cs,y_cs)落在平板探测器中心。否则的话，要知道投影中心(x_cs,y_cs)的精确位置；

(2)精确掌握X射线源到平板探测器的距离D_SD，即成像子系统的焦距；

上述四个几何约束条件都可以利用计算机视觉理论中的相机标定方法进行求解，从而分析系统是否满足了必须满足的两个几何条件，若不满足，几何误差是多大，同时计算得到系统几何参数(x_cs,y_cs)和D_SD。通过对几何参数和误差的求解，达到对锥束CT系统进行标定的目的。其中(x_cs,y_cs)和D_SD相当于摄影系统中的内方位元素，而投影主光轴与旋转轴之间的关系，以及平板探测器的坐标轴YC与旋转子系统的旋转轴线在平板探测器上的投影是否严格平行，相当于摄影系统中的外方位元素。因此，通过相机标定方法得到相机的内方位元素参数，并得到相机拍摄的位置和姿态，就可以得到X线源、平板探测器与旋转子系统旋转轴线之间的相对位置关系，进而得到锥束CT的成像参数及其误差。

2)锥束CT系统的几何参数及误差与系统模型的关系

在建立了上述的源-探测器对的概念，并把系统的投影采集过程考虑为利用源-探测器对构成的相机对旋转系统进行拍摄以后，就简化了系统的几何要素，通过建立相机的模型，上述的几类系统几何误差都可以看成是源-探测器对与旋转子系统之间空间关系的变化所导致的结果，可以将系统的几何误差归纳为三项：即投影主光轴与旋转子系统的旋转轴线是否相交？在相交的情况下两个轴线是否垂直？以及平板探测器的YC轴(即探测单元排布的方向)是否与旋转子系统的旋转轴线在平板探测器上的投影平行？我们只需要关注这三个几何条件，就可以完成锥束CT的标定，这样就大大降低了对锥束CT系统分析的难度。

例如，对平板探测器的安装误差而言，平板探测器在X方向偏离理想位置的位移δx相当于投影主光轴与旋转子系统的旋转轴线不相交；而δy对系统的几何构型无影响；δz相当于使X射线源到平板探测器的距离D_SD发生变化。同时，平板探测器的旋转角度α相当于拍摄时平板探测器的Y_C轴旋转子系统的旋转轴线在平板探测器上的投影不平行，类似于拍摄时相机摆放的不在水平面内；由于规定了投影主光轴总与平板探测器垂直，因此平板探测器的俯仰角度β相当于投影主光轴与旋转子系统的旋转轴线相交但不垂直；而平板探测器的偏转角度γ相当于投影主光轴与旋转子系统的旋转轴线不相交。这样，就可以通过对等效的相机系统的内外方位元素的求解，得到锥束CT的成像参数及其误差。

对于第二类的几何参数及误差，即投影系统参数(X射线源到平板探测器的距离D_SD、投影中心(x_cs,y_cs))通过相机的标定，可以直接得到。

对于第三类的误差，即转轴的倾斜，也可以通过相机的标定验证在各个投影方位，投影主光轴与旋转子系统的旋转轴线是否总保持垂直相交，从而得到其误差。

归纳以上情况，图4是源-探测器对中的投影主光轴与旋转轴线在理想位置的锥束CT采集过程示意图。图中央的正方体代表被测物体，箭头代表转轴，外围的小菱形图案代表X射线源，最外围的平面代表平板探测器。图5是源-探测器对中的投影主光轴与旋转轴线相交但不垂直时(其误差是单一的平板探测器俯仰)，采集过程的示意图；图6是源-探测器对中的投影主光轴与旋转轴线垂直但不相交时(其误差可能是平板探测器偏转或平板探测器平移)，锥束CT采集过程示意图。图7是平板探测器存在旋转误差，即平板探测器YC轴与旋转轴线在平板探测器上的投影不平行时，锥束CT采集过程示意图。

从图4-图7可以清晰地看到我们所建立的模型的三种误差的表现形式，其误差实际上全部表现为利用源-探测器对这一成像系统对物体进行拍摄时的位置与姿态，因此可以采用相机标定的技术测定整个锥束CT系统成像的几何误差。当源-探测器对的几何位置存在混合误差时，锥束CT的采集过程也可以采用相机标定技术进行解析。

3)锥束CT系统模型的标定方法

本发明中首先定义了三个直角坐标系:(1)成像子系统坐标系(即源-探测器坐标系，简称成像系)固结在平板探测器上，用（X_C、Y_C、Z_C）表示，其中的X_C、Y_C轴分别与平板探测器阵列像素排布方向相一致，Z_C轴与探测平面垂直。假设从X线源向平板探测器引出的垂线与平板探测器相交于(x_cs,y_cs)点，X线源到平板探测器的距离为D_SD；(2)旋转子系统坐标系(简称旋转系)，用（X_R、YR_、Z_R）表示，主轴线为Y_R轴，与所述转台的旋转轴线重合，Z_R轴为由所述X射线源向所述转轴所引的垂线，并且指向所述X射线源，X_R轴是Y_R轴和Z_R轴的叉积；(3)标定物坐标系(简称标定物系)，用X,Y,Z表示，是一个与用于标定相机参数的三维物体相关联的辅助坐标系。该三维物体上固定有若干可以由X线源-探测器子系统成像的空间标志点，各个标志点在标定物坐标系内的空间位置是已知的，因此可以通过成像子系统对该三维物体进行拍摄，得到已知的各个空间标志点在平板探测器上对应的投影位置，进而利用相机标定技术解算成像子系统和旋转子系统的几何参数。坐标系的定义如图8所示

锥束CT系统几何标定的具体实现方法为：

(1)标定模块的制做

本发明所采用的相机标定方法，要求利用一个固定的三维标定模块对锥束CT系统的几何参数及其误差进行标定。所设计的三维标定模块上有6个以上可以被X线成像的空间标志点，

这些标志点应当不全共面或共线，以达到三维空间的分布，并且要求已知各个标志点的三维位置。一个具体的标定模块的例子如具体实施方式一节所描述。

(2)用相机标定方法求解锥束CT成像系统的内、外方位元素，从而得到成像系和标定物系之间的转换关系和成像子系统的结构参数(x_cs,y_cs)和D_SD

锥束CT的成像子系统可以看作光学相机模型，X射线源和平板探测器分别对应光学相机模型中的光心和CCD传感器。锥束CT采集单帧图像时，等价于相机对物体进行一次拍摄，当被测物体(或源-探测器对)旋转时，等价于相机从不同角度拍摄物体。利用计算机视觉理论中的相机标定技术，可以获得拍摄的内、外方位元素，其中内方位元素代表X射线源和平板探测器的位置关系，即源-探测器对的X射线源到平板探测器的距离D_SD和投影中心位置(x_cs,y_cs)；外方位元素代表标定物系与成像系之间的空间转换关系，而对于一个特定的拍摄位置，对应于该拍摄位置的外方位元素是不变的。

具体讲，在一个特定的拍摄位置，将标定物系作为固定坐标系XYZ，则对于标定模块上的一个三维空间标志点P_i(x_i,y_i,z_i)，设其对应的在探测器平面上的投影位置为P_ci(x_ci,y_ci)，可以建立以下坐标变换关系(参见图9)：

首先，对于源-探测器对，建立成像系统坐标系X_cY_cZ_c，则在此坐标系中X线源的位置为S_c(x_cs,y_cs,z_cs)，其中z_cs=D_SD;同时，设X线源S在标定物坐标系XYZ中的坐标位置为S(x_s,y_s,z_s)，并将成像系统坐标系与标定物坐标系之间的空间变换关系利用4×4的空间变换矩阵T表示，即：

$\left[\begin{array}{c}{x}_{\mathrm{cs}}\\ y_{\mathrm{cs}}\\ D_{\mathrm{SD}}\\ 1\end{array}\right]=T\left[\begin{array}{c}{x}_s\\ y_s\\ z_s\\ 1\end{array}\right];$

其中： $T=\left[\begin{array}{cccc}{r}_{11}& r_{12}& r_{13}& t_x\\ r_{21}& r_{22}& r_{23}& t_y\\ r_{31}& r_{32}& r_{33}& t_z\\ 0& 0& 0& 1\end{array}\right];$ 简写为： $T=\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right];---\left(2\right)$ 因此标定模块上(在标定物坐标系内)的已知点P_i(x_i,y_i,z_i)在成像系统坐标系中的坐标可以表示为：

$\left[\begin{array}{c}{x}_{\mathrm{ci}}^{\′}\\ y_{\mathrm{ci}}^{\′}\\ z_{\mathrm{ci}}^{\′}\\ 1\end{array}\right]=T\left[\begin{array}{c}{x}_i\\ y_i\\ z_i\\ 1\end{array}\right];$ 或

其中

为标定模块上的P_i点在成像坐标系中的坐标。

进一步在成像坐标系中考虑已知的P_i坐标与在平板探测器上的投影位置P_ci(x_ci,y_ci)之间的关系(参见图9)：

$\left[\begin{array}{c}{x}_{\mathrm{ci}}-x_{\mathrm{cs}}\\ y_{\mathrm{ci}}-y_{\mathrm{cs}}\\ 1\end{array}\right]=\frac{1}{z_{\mathrm{ci}}^{\′}}\left[\begin{array}{cccc}{D}_{\mathrm{SD}}& 0& 0& 0\\ 0& D_{\mathrm{SD}}& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{ci}}^{\′}-x_{\mathrm{cs}}\\ y_{\mathrm{ci}}^{\′}-y_{\mathrm{cs}}\\ z_{\mathrm{ci}}^{\′}\\ 1\end{array}\right]$

$=\frac{1}{z_{\mathrm{ci}}^{\′}}\left[\begin{array}{cccc}{D}_{\mathrm{Sd}}& 0& 0& 0\\ 0& D_{\mathrm{SD}}& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& {-x}_{\mathrm{cs}}\\ 0& 1& 0& {-y}_{\mathrm{cs}}\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{ci}}^{\′}\\ y_{\mathrm{ci}}^{\′}\\ z_{\mathrm{ci}}^{\′}\\ 1\end{array}\right]$

$=\frac{1}{z_{\mathrm{ci}}^{\′}}\left[\begin{array}{cccc}{D}_{\mathrm{SD}}& 0& 0& {-D}_{\mathrm{SD}}\·x_{\mathrm{cs}}\\ 0& D_{\mathrm{SD}}& 0& {-D}_{\mathrm{SD}}\·y_{\mathrm{cs}}\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{ci}}^{\′}\\ y_{\mathrm{ci}}^{\′}\\ z_{\mathrm{ci}}^{\′}\\ 1\end{array}\right]$

(4)通过变形可以得到：

$\left[\begin{array}{c}{x}_{\mathrm{ci}}\\ y_{\mathrm{ci}}\\ 1\end{array}\right]=\frac{1}{z_{\mathrm{ci}}^{\′}}\left[\begin{array}{ccc}1& 0& x_{\mathrm{cs}}\\ 0& 1& y_{\mathrm{cs}}\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}{D}_{\mathrm{SD}}& 0& 0& {-D}_{\mathrm{SD}}\·x_{\mathrm{cs}}\\ 0& D_{\mathrm{SD}}& 0& {-D}_{\mathrm{SD}}\·y_{\mathrm{cs}}\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{ci}}^{\′}\\ y_{\mathrm{ci}}^{\′}\\ z_{\mathrm{ci}}^{\′}\\ 1\end{array}\right]$

$=\frac{1}{z_{\mathrm{ci}}^{\′}}\left[\begin{array}{cccc}{D}_{\mathrm{SD}}& 0& x_{\mathrm{cs}}& {-D}_{\mathrm{Sd}}\·x_{\mathrm{cs}}\\ 0& D_{\mathrm{SD}}& y_{\mathrm{cs}}& {-D}_{\mathrm{SD}}\·y_{\mathrm{cs}}\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{ci}}^{\′}\\ y_{\mathrm{ci}}^{\′}\\ z_{\mathrm{ci}}^{\′}\\ 1\end{array}\right]$ (5)以mm为单位的坐标(x_ci,y_ci)与以像素为单位的坐标(u_ci,v_ci)之间的关系如图10示：(x_ci,y_ci)和(u_ci,v_ci)之间的关系为：

其中，d_x和dy分别代表像素在X_C和Y_C方向上的尺寸，D代表平板探测器在Y_C方向上的长度（mm）。因此，(x_ci,y_ci)与(u_ci,v_ci)之间的关系可重写为：

$\left[\begin{array}{c}{u}_{\mathrm{ci}}\\ v_{\mathrm{ci}}\\ 1\end{array}\right]=\left[\begin{array}{ccc}1/\mathrm{dx}& 0& 0\\ 0& -1/\mathrm{dy}& D/\mathrm{dy}\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{ci}}\\ y_{\mathrm{ci}}\\ 1\end{array}\right]---\left(6\right)$

综合公式(5)和(6)，可建立标定物坐标系内的已知点P_i(x_i,y_i,z_i)与平板探测器上探测得到的对应位置(u_ci,v_ci)之间的关系：

$z_{\mathrm{ci}}^{\′}\left[\begin{array}{c}{u}_{\mathrm{ci}}\\ v_{\mathrm{ci}}\\ 1\end{array}\right]=\left[\begin{array}{ccc}1/\mathrm{dx}& 0& 0\\ 0& -1/\mathrm{dy}& D/\mathrm{dy}\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}{D}_{\mathrm{SD}}& 0& x_{\mathrm{cs}}& {-D}_{\mathrm{SD}}\·x_{\mathrm{cs}}\\ 0& D_{\mathrm{SD}}& y_{\mathrm{cs}}& {-D}_{\mathrm{SD}}\·y_{\mathrm{cs}}\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}{x}_i\\ y_i\\ z_i\\ 1\end{array}\right]$

$=\left[\begin{array}{cccc}{D}_{\mathrm{SD}}/\mathrm{dx}& 0& x_{\mathrm{cs}}/\mathrm{dx}& {-D}_{\mathrm{SD}}\·x_{\mathrm{cs}}/\mathrm{dx}\\ 0& {-D}_{\mathrm{SD}}/\mathrm{dx}& (D-y_{\mathrm{cs}})/\mathrm{dy}& D_{\mathrm{SD}}\·y_{\mathrm{cs}}/\mathrm{dy}\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}{x}_i\\ y_i\\ z_i\\ 1\end{array}\right]$ (7)

$=\mathrm{LT}\left[\begin{array}{c}{x}_i\\ y_i\\ z_i\\ 1\end{array}\right]$

其中， $L=\left[\begin{array}{cccc}{D}_{\mathrm{SD}}/\mathrm{dx}& 0& x_{\mathrm{cs}}/\mathrm{dx}& {-D}_{\mathrm{SD}}\·x_{\mathrm{cs}}/\mathrm{dx}\\ 0& {-D}_{\mathrm{SD}}/\mathrm{dy}& (D-y_{\mathrm{cs}})/\mathrm{dy}& D_{\mathrm{SD}}\·y_{\mathrm{cs}}/\mathrm{dy}\\ 0& 0& 1& 0\end{array}\right]$ 代表由成像子系统坐标系到以像素为单位的坐标的转换关系。

上述方程中D，dx，dy为已知量，D_SD，x_cs，y_cs以及T是带求解的参数。其中D_SD，x_cs，y_cs是成像子系统的内方位参数，在系统制作完成后是不变的，但是是未知的。对于一个特定的拍摄位置，成像子系统的空间位置姿态都不改变，因此空间变换矩阵T是不变的。因此若有足够多的标志点，就可以根据(7)式，由已知的标志点坐标(x_i,y_i,z_i)和对应的像素坐标(u_ci,v_ci)求解系统的参数，包括D_SD，x_cs，y_cs和空间变换阵 $\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]$ 中的R阵和位移向量t。

方程(7)就简化为：

$z_{\mathrm{ci}}^{\′}\left[\begin{array}{c}{u}_{\mathrm{ci}}\\ v_{\mathrm{ci}}\\ 1\end{array}\right]=H_{3\×4}\left[\begin{array}{c}{x}_i\\ y_i\\ z_i\\ 1\end{array}\right]---\left(8\right)$ 其中:

$H_{3\×4}=\mathrm{LT}$

$=\left[\begin{array}{cccc}{D}_{\mathrm{Sd}}/\mathrm{dx}& 0& x_{\mathrm{cs}}/\mathrm{dx}& {-D}_{\mathrm{SD}}\·x_{\mathrm{cs}}/\mathrm{dx}\\ 0& {-D}_{\mathrm{SD}}/\mathrm{dy}& (D-y_{\mathrm{cs}})/\mathrm{dy}& D_{\mathrm{SD}}\·y_{\mathrm{cs}}/\mathrm{dy}\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]$ (9)

$=\left[\begin{array}{cccc}{h}_{11}& h_{12}& h_{13}& h_{14}\\ h_{21}& h_{22}& h_{23}& h_{24}\\ h_{31}& h_{32}& h_{33}& h_{34}\end{array}\right]$

从(8)式中可以消去zci，得到方程组:

$\left\{\begin{array}{c}{h}_{11}{x}_i+h_{12}{y}_i+h_{13}{z}_i+h_{14}-h_{31}{u}_{\mathrm{ci}}{x}_i-h_{32}{u}_{\mathrm{ci}}{y}_i-h_{33}{u}_{\mathrm{ci}}{z}_i=h_{34}{u}_{\mathrm{ci}}\\ h_{21}{x}_i+h_{22}{y}_i+h_{23}{z}_i+h_{24}-h_{31}{v}_{\mathrm{ci}}{x}_i-h_{32}{v}_{\mathrm{ci}}{y}_i-h_{33}{v}_{\mathrm{ci}}{z}_i=h_{34}{v}_{\mathrm{ci}}\end{array}\right.---\left(10\right)$

对于标定模块上的I个标志点，可以得到包含2I个方程的联立方程组：

$\left[\begin{array}{ccccccccccc}{x}_1& y_1& z_1& 1& 0& 0& 0& 0& {-u}_{\mathrm{cjl}}{x}_1& {-u}_{\mathrm{cjl}}{y}_1& {-u}_{\mathrm{cjl}}{z}_1\\ 0& 0& 0& 0& x_1& y_1& z_1& 1& {-v}_{\mathrm{cj}1}{x}_1& {-v}_{\mathrm{cjl}}{y}_1& {-v}_{\mathrm{cjl}}{z}_1\\ & \·\·\·& \·\·\·& & & & & & & & \\ & \·\·\·& \·\·\·& & & & & & & & \\ x_I& y_I& z_I& 1& 0& 0& 0& 0& {-u}_{\mathrm{cjI}}{x}_I& {-u}_{\mathrm{cjI}}{y}_I& {-u}_{\mathrm{cjI}}{z}_I\\ 0& 0& 0& 0& x_I& y_I& z_I& 1& {-v}_{\mathrm{cjI}}{x}_I& {-v}_{\mathrm{cjI}}{y}_I& {-v}_{\mathrm{cjI}}{z}_I\end{array}\right]\left[\begin{array}{c}{h}_{j11}\\ h_{j12}\\ h_{j13}\\ h_{j14}\\ h_{j21}\\ h_{j22}\\ h_{j23}\\ h_{j24}\\ h_{j31}\\ h_{j32}\\ h_{j33}\end{array}\right]=\left[\begin{array}{c}{u}_{\mathrm{cj}1}{h}_{j34}\\ v_{\mathrm{cj}1}{h}_{j34}\\ \·\·\·\\ \·\·\·\\ u_{\mathrm{cjI}}{h}_{j34}\\ v_{\mathrm{cjI}}{h}_{j34}\end{array}\right]---\left(11\right)$

因此可以通过标定模块上已知空间位置的多个标志点，以及其对应的在平板探测器上的投影位置，求解出成像系统每一次特定位置拍摄的所有内、外方位元素，从而得到成像子系统的参数。

由公式（11）可知，每个空间标志点可以列出两个方程。在求解公式（11）所述的线性方程组时，通常设h₃₄=1。公式（11）中共有11个未知数，故至少需要6个空间标志点才可以求解系统的几何参数。

在依据公式（11）求解出H_3×4后，则可以根据公式（9）求解出特定位置时的空间变换阵T和D_SD，x_cs，y_cs。当从多个位置进行类似的拍摄时，每一个位置都可以得到一组空间变换阵T和D_SD，x_cs，y_cs三个变量，可以通过对各个位置计算得到的D_SD，x_cs，y_cs进行均值计算，减少测量标定中的随机误差，从而得到精化的成像系统几何参数D_SD，x_cs，y_cs。

另一方面，通过对各个位置得到的多组空间变换阵T的分析，可以解析出锥束CT系统的最核心的几何参数和误差，从而了解系统是否满足理想的锥束CT所需要满足的两个几何约束条件(即：①在成像的各个投影成像位置，成像子系统的投影主光轴都与旋转子系统的旋转轴线严格垂直相交；②平板探测器的坐标轴Y_C与旋转子系统的旋转轴线在平板探测器上的投影严格平行)，下面具体分析。

(3)锥束CT系统转轴的位置方向的求解

实际上，从第1）和第2）节的分析可知，无论锥束CT系统是否满足理想的两个几何约束条件，标定物体上的空间标志点的空间轨迹都是一个空间圆周运动。因此从各次投影得到的标定物体上空间标志点的空间轨迹就可以直接求解出锥束CT系统转轴位置。具体讲，在每个投影位置，在成像子系统坐标系中标定物系的原点的空间位置可以由公式（2），通过空间变换阵T得到：

$S_j=\left[\begin{array}{cc}{R}_j& t_j\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}0\\ 0\\ 0\\ 1\end{array}\right]=t_i$

S_j的轨迹是一个空间的圆周，通过求均值可以得到其圆心位置，即锥束CT系统转轴的位置；通过计算向量端点轨迹所拟合出平面的法向量，可以得到锥束CT系统转轴的方向。

(4)成像子系统投影主光轴的求解

成像子系统投影主光轴是X线源到平板探测器的垂线。因此，成像子系统投影主光轴的向量为：

$O=C_C-S_C$

其中C_C为所述X线源在所述平板探测器上的垂足的齐次坐标，S_C为所述X射线源的空间位置的齐次坐标，S_C=[x_cs,y_cs,D_SD,1]^T，C_C=[x_cs,y_cs,0,1]^T；

(5)成像子系统投影主光轴与旋转子系统旋转轴线几何关系的求解

在得到旋转子系统旋转轴线的空间向量和成像子系统投影主光轴的空间向量后，就可以求解二者之间的几何关系，并验证两个轴线是否垂直相交。利用成像子系统投影主光轴向量与旋转子系统旋转轴线向量之间的点积可以分析两个轴线是否相互垂直；通过计算两个向量之间的距离，可以知道两者是否共面(相交)，从而分析系统是否存在偏离理想位置的位移δx或存在偏转角度γ等几何误差。

(6)平板探测器像素方向与旋转子系统的旋转轴线在平板探测器上的投影之间几何关系

的求解

由于已经解析出X线源的位置，同时旋转子系统旋转轴线是已知的(并且是固定的)，因此可以以X线源为投影点，将转轴向平板探测器进行投影，从而分析平板探测器像素方向YC与旋转子系统的旋转轴线在平板探测器上的投影是否满足平行的几何关系，进而分析系统是否存在几何误差。

(7)锥束CT系统几何参数的获取和误差矫正

在得到了锥束CT系统转轴的位置方向和成像子系统的投影主光轴与旋转子系统的旋转轴线之间关系以后，可以计算出平板探测器的安装误差、投影系统参数、转轴的倾斜等必要的几何参数和所对应的误差，从而对锥束CT系统的机械安装进行调整，以便减小或消除误差。同时，也可以将这些不够理想的几何参数作为CT重建的几何条件，优化重建算法，从而得到更高质量的重建图像。

本发明利用计算机视觉理论和相机标定技术，通过建立一个合理的锥束CT系统的模型，将锥束CT系统抽象为成像子系统和旋转子系统，并将系统的几何参数及其误差表达为两者间的位置关系，从而可以直接求解锥束CT的几何参数和几何误差，方法简便，易于实施。由于本发明的方法在实施中不涉及复杂的测量和迭代算法，并且相机标定技术已经发展成为成熟、可靠的理论体系，因此本发明中方法易于实施，其稳定性很好地得到保证。

以下结合本发明的技术方案和附图详细说明本发明的具体实施方式。其中的矩阵和向量运算均由Matlab Main Toolbox（Matlab主工具箱）进行实现。

第1步，制做标定模块。所设计的三维标定模块上至少有6个以上可以被X线成像的空间标志点，这些标志点应当不全共面或共线，以达到三维空间的分布，并且要求已知各个标志点的三维位置。作为一个具体的例子，发明人设计的一个标定模块包含有12个空间标志点，模块采用有机玻璃作为基体，在其外表面精确钻孔并镶嵌小钢珠作为空间标志点。

第2步，利用锥束CT系统对标定模块进行成像。将模块放在锥束CT成像系统的视场中央附近，启动锥束CT的数据采集程序，获取在锥束CT整个旋转范围内不少于360个不同投影角度的投影图像，这些投影应当比较均匀地覆盖整个投影范围。

第3步，在各个获取的投影图像中抽取各个标志点的二维投影位置坐标，对第j个投影图像中的第i个标志点的投影位置坐标，记为(u_cji,v_cji)，其对应的标定模块上的空间标志点坐标是(x_i,y_i,z_i)。

第4步，利用公式（11）求解每个投影位置的H阵的诸元素。对投影位置j的具体计算公式为：

$\left[\begin{array}{ccccccccccc}{x}_1& y_1& z_1& 1& 0& 0& 0& 0& {-u}_{\mathrm{cjl}}{x}_1& {-u}_{\mathrm{cjl}}{y}_1& {-u}_{\mathrm{cjl}}{z}_1\\ 0& 0& 0& 0& x_1& y_1& z_1& 1& {-v}_{\mathrm{cj}1}{x}_1& {-v}_{\mathrm{cjl}}{y}_1& {-v}_{\mathrm{cjl}}{z}_1\\ & \·\·\·& \·\·\·& & & & & & & & \\ & \·\·\·& \·\·\·& & & & & & & & \\ x_I& y_I& z_I& 1& 0& 0& 0& 0& {-u}_{\mathrm{cjI}}{x}_I& {-u}_{\mathrm{cjI}}{y}_I& {-u}_{\mathrm{cjI}}{z}_I\\ 0& 0& 0& 0& x_I& y_I& z_I& 1& {-v}_{\mathrm{cjI}}{x}_I& {-v}_{\mathrm{cjI}}{y}_I& {-v}_{\mathrm{cjI}}{z}_I\end{array}\right]\left[\begin{array}{c}{h}_{j11}\\ h_{j12}\\ h_{j13}\\ h_{j14}\\ h_{j21}\\ h_{j22}\\ h_{j23}\\ h_{j24}\\ h_{j31}\\ h_{j32}\\ h_{j33}\end{array}\right]=\left[\begin{array}{c}{u}_{\mathrm{cj}1}{h}_{j34}\\ v_{\mathrm{cj}1}{h}_{j34}\\ \·\·\·\\ \·\·\·\\ u_{\mathrm{cjI}}{h}_{j34}\\ v_{\mathrm{cjI}}{h}_{j34}\end{array}\right]$

等式两端同除以hj34，得到：

${\mathrm{Ah}}_j^{\′}=B$

式中： $A=\left[\begin{array}{ccccccccccc}{x}_1& y_1& z_1& 1& 0& 0& 0& 0& {-u}_{\mathrm{cj}1}{x}_1& {-u}_{\mathrm{cj}1}{y}_1& {-u}_{\mathrm{cj}1}{z}_{21}\\ 0& 0& 0& 0& x_1& y_1& z_1& 1& {-v}_{\mathrm{cj}1}{x}_1& {-v}_{\mathrm{cj}1}{y}_1& {-v}_{\mathrm{cj}1}{z}_1\\ & \·\·\·& \·\·\·& & & & & & & & \\ & \·\·\·& \·\·\·& & & & & & & & \\ x_I& y_I& z_I& 1& 0& 0& 0& 0& {-u}_{\mathrm{cjI}}{x}_I& {-u}_{\mathrm{cjI}}{y}_I& {-u}_{\mathrm{cjI}}{z}_I\\ 0& 0& 0& 0& x_I& y_I& z_I& 1& {-v}_{\mathrm{cjI}}{x}_I& {-v}_{\mathrm{cjI}}{x}_I& {-v}_{\mathrm{cjI}}{z}_I\end{array}\right],$ 是2I×11的矩阵，B=[u_cj1,v_cj1,…,u_cjI,v_cjI]^T是2I×1的向量，

$h_j={\left[\begin{array}{ccccccccccc}{h}_{j11}& h_{j12}& h_{j13}& h_{j14}& h_{j21}& h_{j22}& h_{j23}& h_{j24}& h_{j31}& h_{j32}& h_{j33}\end{array}\right]}^T;$

利用最小二乘法，得到：

$h_j^{\′}={\left(A^{T}A\right)}^{-1}{A}^{T}B$

按以下步骤计算表换矩阵H_(3×4)j中的各个元素：h_j1，h_j2，h_j3，h_j14，h_j24，h_j34，其中：h_j1=[h_j11,h_j12,h_j13]^T，h_j2=[h_j21,h_j22,h_j23]^T，h_j3=[h_j31,h_j32,h_j33]^T，

根据公式（9），H_(3×4)j表达为：

$H_{(3\×4)j}=\left[\begin{array}{cc}{h}_{j1}^T& h_{j14}\\ h_{j2}^T& h_{j24}\\ h_{j3}^T& h_{j34}\end{array}\right]\left[\begin{array}{cccc}{a}_{11}& 0& a_{13}& a_{14}\\ 0& a_{22}& a_{23}& a_{24}\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}{r}_{j1}^T& t_{\mathrm{jx}}\\ r_{j2}^T& t_{\mathrm{jy}}\\ r_{j3}^T& t_{\mathrm{jz}}\\ 0^T& 1\end{array}\right]$

其中， $a_{11}=D_{\mathrm{SD}}/\mathrm{dx},$ $a_{13}=x_{\mathrm{cs}}/\mathrm{dx},$ $a_{14}=({-D}_{\mathrm{SD}}\·x_{\mathrm{cs}})/\mathrm{dx},$ $a_{22}={-D}_{\mathrm{SD}}/\mathrm{dy},$ $a_{23}=(D-y_{\mathrm{cs}})/\mathrm{dy},$ $a_{24}=(D_{\mathrm{SD}}\·y_{\mathrm{cs}})/\mathrm{dy},$

$R_j=\left[\begin{array}{c}{r_{j1}}^T\\ {r_{j2}}^T\\ {r_{j3}}^T\end{array}\right],$ 其中，r_j1=[r_j11,r_j12,r_j13]^T，r_j2=[r_j21,r_j22,r_j23]^T，r_j3=[r_j31,r_j32,r_j33]T

t_j=[t_jx,t_jy,t_jz]^T

从而得到：

$H_{(3\×4)j}=\left[\begin{array}{cc}{a}_{11}{r}_{j1}^T+a_{13}{r}_{j3}^T& a_{11}{t}_{\mathrm{jx}}+a_{13}{t}_{\mathrm{jz}}+a_{14}\\ a_{22}{r}_{j2}^T+a_{23}{r}_{j3}^T& a_{22}{t}_{\mathrm{jy}}+a_{23}{t}_{\mathrm{jz}}+a_{24}\\ r_{j3}^T& t_{\mathrm{jz}}\end{array}\right]$

由于R_j是单位正交矩阵，则有：

$\left|\right|r_{\mathrm{jm}}\left|\right|=1$

$r_{\mathrm{jm}}^T{r}_{\mathrm{jn}}=0$

$r_{\mathrm{jm}}\×r_{\mathrm{jm}}=0$

$\left|\right|r_{\mathrm{jm}}\×r_{\mathrm{jn}}\left|\right|=1$

其中，n=1,2,3,m=1,2,3，且n≠m，

分别对应于R_j的第m行和第n行，r_jm和r_jn相互正交，并且均为单位向量，

由上所述有：

$\left|\right|r_{j3}\left|\right|=\left|\right|h_{j3}\left|\right|=h_{j34}\left|\right|h_{j3}^{\′}\left|\right|=1$

其中 $h_{j3}^{\′}=h_{j3}/h_{j34}$

故有： $h_{j34}=1/\left|\right|h_{j3}^{\′}\left|\right|$

由上得到： $h_j=h_{j34}{h}_j^{\′}$

至此，就得到了矩阵H_(3×4)j中的各个元素，也就得到了h_j1，h_j2，h_j3，h_j14，h_j24，h_j34；

第5步，得到H_(3×4)j中的各个元素后，进一步求解参数D_SD，x_cs，y_cs，R_j，t_j，具体步骤如下：

$a_{13}=(a_{11}{r}_{j1}^T+a_{13}{r}_{j3}^T)\·r_{j3}=h_{j1}^T{h}_{j3}$

$a_{23}=(a_{22}{r}_{j2}^T+a_{23}{r}_{j3}^t)r_{j3}=h_{j2}^T{h}_{j3}$

$a_{11}=\left|\right|(a_{11}{r}_{j1}+a_{13}{r}_{j3})\×r_{j3}\left|\right|=\left|\right|h_{j1}\×h_{j3}\left|\right|$

$a_{22}=\left|\right|(a_{22}{r}_{j2}+a_{23}{r}_{j3})\×r_{j3}\left|\right|=\left|\right|h_{j2}\×h_{j3}\left|\right|$

从上可得：

$D_{\mathrm{SD}}=a_{11}\·\mathrm{dx}={-a}_{22}\·\mathrm{dy}$

$x_{\mathrm{cs}}=a_{13}\·\mathrm{dx}$

$y_{\mathrm{cs}}=D-a_{23}\·\mathrm{dy}$

$r_{j1}=\frac{1}{a_{11}}(h_{j1}-a_{13}{h}_{j3})$

$r_{j3}=\frac{1}{a_{22}}(h_{j2}-a_{23}{h}_{j3})$

$r_{j3}=h_{j3}$

$R_j=\left[\begin{array}{c}{r}_{j1}^T\\ r_{j3}^T\\ r_{j3}^T\end{array}\right]$

从上可得：

$a_{14}=-(D_{\mathrm{SD}}\·x_{\mathrm{cs}})/\mathrm{dx}$

$a_{24}=D_{\mathrm{SD}}\·y_{\mathrm{cs}}/\mathrm{dy}$

因此，t_j求解如下：

$t_{\mathrm{jx}}=\frac{1}{a_{11}}(h_{14}-a_{14}-a_{13}\·h_{34}),$

$t_{\mathrm{jy}}=\frac{1}{a_{22}}(h_{24}-a_{24}-a_{23}\·h_{34})$

$t_{\mathrm{jz}}=h_{34}$

$t_j={[t_{\mathrm{jx}},t_{\mathrm{jy}},t_{\mathrm{jz}}]}^T$

由上述关系，则求解出对应各个投影位置的参数D_SD，x_cs，y_cs，R_j，t_j；

第6步，求解锥束CT系统旋转子系统转轴的位置和方向

利用第5步中的参数，得到不同投影位置时标定物系的原点在成像系中的空间位置向量S_j：

$S_j=\left[\begin{array}{cc}{R}_j& t_j\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}0\\ 0\\ 0\\ 1\end{array}\right]=t_j$

这些向量的端点轨迹是一个空间的圆周，因此通过求解S_j的均值，得到其圆心位置C，即所述转轴通过的位置C：

$c=\frac{1}{j}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{S}_j$

上式中J代表一共采集了J个投影，J应不少于360，

锥束CT系统中所述转轴的方向，则为向量S_j端点的轨迹所拟合出平面的法向量，通过求所述圆周上三个点所构成任意两个向量的叉积而得到，本发明使用以下公式进行求解：

$Y=\frac{1}{J}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}[(S_{j+k}-S_j)\×(S_{j+1}-S_{j+k})]$

其中，k和l分别取为1/3和2/3的投影数目，

至此，确定了锥束CT系统中所述转轴的空间向量是过C点，方向为Y的向量；

第7步，求解成像子系统的投影主光轴，

投影主光轴是所述X射线源到所述平板探测器的垂线，在成像系中，投影主光轴的向量为O为：

$O=C_C-S_C$

其中C_C为所述X线源在所述平板探测器上的垂足的齐次坐标，S_C为所述X射线源的空间位置的齐次坐标，S_C=[x_cs,y_cs,D_SD,1]^T，C_C=[x_cs,y_cs,0,1]^T；

第8步，锥束CT几何误差的求解，

锥束CT的理想几何条件为所述投影主光轴与所述转轴相交且垂直，并且所述转轴在所述平板探测器上的投影与平板探测器的YC方向平行，这里依据第6步和第7步的结果，首先判断锥束CT是否符合上述理想条件，若不满足，其误差是多大，

8.1，判断所述投影主光轴与所述转轴是否垂直，若不垂直，其夹角是多大，

设a=O^TY，

若a=0，则所述投影主光轴和所述转轴垂直，满足理想情况，

若a≠0，则所述投影主光轴和所述转轴不垂直，其夹角为α：

$\mathrm{\α}={\cos }^{-1}\frac{O^{T}Y}{{\left|\right|O\left|\right|\left|\right|Y\left|\right|}^{,}}$

8.2，判断所述投影主光轴与所述转轴是否相交，若不相交，二者之间的距离是多大，

设b=O×Y，

若

则所述投影主光轴和所述转轴平行，即二者不相交，

若

则二者之间的距离d为：

$d={(C-S_C)}^{T}b$

若d=0，则二者共面且不平行，即二者相交，

若d≠0，则二者不共面，即不相交，二者之间的距离为d；

8.3，判断所述转轴在所述平板探测器上的投影是否与所述平板探测器的Y_C方向平行，若不平行，倾斜角是多少，

由第6步可知，转轴上存在两点在成像系中的坐标为C和C+Y，其在平板探测器上的投影为

$P=\mathrm{LC}$

$P^{\′}=L(C+L)$

其中， $L=\left[\begin{array}{cccc}{a}_{11}& 0& a_{13}& a_{14}\\ 0& a_{22}& a_{23}& a_{24}\\ 0& 0& 1& 0\end{array}\right],$ P(u_P,v_P)，P^′(u_P′,v_P′)，

所述转轴在所述平板探测器上投影的倾斜角为β：

$\mathrm{\β}={\tan }^{-1}\frac{u_P-u_{P^{\′}}}{v_P-v_{P^{\′}}}$

若β=0，则所述转轴在所述平板探测器上的投影与探测器的Y_C方向平行；

至此，就完成了所述锥束CT系统的几何参数的标定，所述几何参数包括D_SD，x_cs，y_cs，以及所述投影主光轴与所述转轴之间的几何误差。

Claims (1)

1.一种X射线锥束计算机层析成像系统的几何参数标定方法，其特征在于，是在一个由平板探测器、同轴地安置有被检测物的转台、X射线源以及计算机共同构成的锥束CT系统中实现的，是一种利用相机标定技术来确定所述锥束CT系统中的所述X射线源、平板探测器和转台中的转轴之间的几何位置关系，以对所述锥束CT系统的几何参数及其误差进行直接求解的方法，步骤如下：

步骤（1），制作用于模拟所述被检测物的三维标定模块：

在所述三维标定模块上有不少于6个不全共面或不全共线的空间标志点，这里使用12个空间标志点，所述三维标定模块采用有机玻璃基体，外表面精确钻孔并镶嵌小钢珠作为空间标志点；

步骤（2），建立如下三个直角坐标系：

成像子系统坐标系，是一个所述的X射线源-平板探测器坐标系，以下称成像系，设在所述平板探测器上，用（X_C、Y_C、Z_C）表示，坐标原点在所述平板探测器的左下角，表示为(x_c0,y_c0,z_c0)，x_c0=y_c0=z_c0=0，X_C轴、Y_C轴分别对应于所述平板探测器上像素排布的方向，Z_C轴是X_C轴、Y_C轴的叉积，并垂直于探测器平面，从所述X射线源到所述平板探测器引垂线，垂足为C_C，X射线源到垂足C_C的距离为D_SD，单位为mm，在所述成像系中：

X射线源的坐标点为S_C，坐标为(x_cs,y_cs,z_cs)，z_cs=D_SD，x_cs=y_cs≠0，

垂足C_C的坐标为(x_cc,y_cc,z_cc)，z_cc=0，x_cc=x_cs，y_cc=y_cs，

在所述平板探测器的各个投影图像内，当把以mm为单位的坐标(x_ci,y_ci)转换为以像素为单位的坐标(u_ci,v_ci)时，把以像素为单位的坐标的原点设在投影图像左上角，用(u_c0,v_c0)表示，且坐标轴u向右为正方向，坐标轴v向下为正方向，则：

$u_{\mathrm{ci}}=\frac{x_{\mathrm{ci}}}{\mathrm{dx}},$ $v_{\mathrm{ci}}=\frac{D-y_{\mathrm{ci}}}{\mathrm{dy}},$ 其中：

i表示在所述投影图像内的点的序号，

D为所述平板探测器在Y_C轴上的长度，

dx，dy分别代表单个像素在X_C、Y_C方向上的尺寸，

旋转子系统坐标系，以下称旋转系，用（X_R、Y_R、Z_R）表示，主轴线为Y_R轴，与所述转台的旋转轴线重合，Z_R轴为由所述X射线源向所述转轴所引的垂线，并且指向所述X射线源，X_R轴是Y_R轴和Z_R轴的叉积，

标定物坐标系，以下称标定物系，用（X、Y、Z）表示，是一个同轴地设置在所述三维标定模块上，用于标定锥束CT参数的辅助坐标系，所述X射线源在所述标定物系中的坐标为(x_s,y_s,z_s)，所述三维标定模块上的12个空间标志点用Pi表示，坐标为(x_i,y_i,z_i)，是已知的，i=1,2,…,I，I=12，以便通过所述成像子系统对所述三维标定模块进行拍摄，从而在所述平板探测器上得到对应于12个所述空间标志点的投影位置，坐标点用P_ci表示，坐标用(x_ci,y_ci)表示；

步骤（3），把所述三维标定模块与所述转台在纵向同轴安置，在所述X射线源-平板探测器相对静止，转台转动的条件下，对所述三维标定模块成像：所述计算机启动锥束CT采集程序，获取所述平板探测器上不少于360个投影角度的投影图像，用j表示所述投影图像序号，j=1,2,…,J，J≥360，在获取到的各个所述投影图像j中抽取12个所述空间标志点的坐标点，用P_cji表示，坐标为(x_cji,y_cji)，转换为像素坐标时用(u_cji,v_cji)表示，对应的所述三维标定模块上的坐标点用P_i表示，坐标用(x_i,y_i,z_i)表示；

步骤（4），按以下步骤从所述标定物系的空间标志点P_i(x_i,y_i,z_i)和成像系中对应的像素坐标P_cji(u_cji,v_cji)中求解所述锥束CT系统的以下参数：D_SD、x_cs和y_cs，标定物系到成像系的空间转换矩阵 $T_j=\left[\begin{array}{cc}{R}_j& t_j\\ 0^T& 1\end{array}\right]$ 中的各所述投影图像的旋转矩阵R_j和位移向量t_j，其中，D_SD、x_cs和y_cs是所述锥束CT系统的内方位参数，在所述成像子系统制作完成后是不变的，但又是未知的，而对于一个设定的拍摄位置，对应于该拍摄位置的空间转换矩阵T是不变的；

步骤（4.1），按下式建立所述标定物系内已知空间坐标点P_i(x_i,y_i,z_i)与各所述投影图像j中的对应坐标点P_cji(u_cji,v_cji)之间的关系：

$z_{\mathrm{ci}}^{\′}\left[\begin{array}{c}{u}_{\mathrm{cji}}\\ v_{\mathrm{cji}}\\ 1\end{array}\right]=H_{(3\×4)j}\left[\begin{array}{c}{x}_i\\ y_i\\ z_i\\ 1\end{array}\right],$

H_(3×4)j为所述以像素为单位的坐标与标定物系坐标之间的转换关系，

$H_{(3\×4)j}=\left[\begin{array}{cccc}{D}_{\mathrm{SD}}/\mathrm{dx}& 0& x_{\mathrm{cs}}/\mathrm{dx}& {-D}_{\mathrm{SD}}\·x_{\mathrm{cs}}/\mathrm{dx}\\ 0& {-D}_{\mathrm{SD}}/\mathrm{dx}& (D-y_{\mathrm{cs}})/\mathrm{dy}& D_{\mathrm{SD}}\·y_{\mathrm{cs}}/\mathrm{dy}\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}{R}_j& t_j\\ 0^T& 1\end{array}\right]$ $=\left[\begin{array}{cccc}{h}_{j11}& h_{j12}& h_{j13}& h_{j14}\\ h_{j21}& h_{j22}& h_{j23}& h_{j24}\\ h_{j31}& h_{j32}& h_{j33}& h_{j34}\end{array}\right]$

H_(3×4)j为第j幅投影图像的变换矩阵，

R_j为第j幅投影图像的旋转矩阵，是单位正交矩阵，

t_j为第j幅投影图像的位移向量，

$\left[\begin{array}{cc}{R}_j& t_j\\ 0^T& 1\end{array}\right]$ 为第j幅投影图像的空间转换矩阵，

对于所述三维标定模块的I=12个空间标志点，得到下述对应于第j个投影图像的，包含2I个方程的联立方程组：

$\left[\begin{array}{ccccccccccc}{x}_1& y_1& z_1& 1& 0& 0& 0& 0& {-u}_{\mathrm{cjl}}{x}_1& {-u}_{\mathrm{cjl}}{y}_1& {-u}_{\mathrm{cjl}}{z}_1\\ 0& 0& 0& 0& x_1& y_1& z_1& 1& {-v}_{\mathrm{cj}1}{x}_1& {-v}_{\mathrm{cjl}}{y}_1& {-v}_{\mathrm{cjl}}{z}_1\\ & \·\·\·& \·\·\·& & & & & & & & \\ & \·\·\·& \·\·\·& & & & & & & & \\ x_I& y_I& z_I& 1& 0& 0& 0& 0& {-u}_{\mathrm{cjI}}{x}_I& {-u}_{\mathrm{cjI}}{y}_I& {-u}_{\mathrm{cjI}}{z}_I\\ 0& 0& 0& 0& x_I& y_I& z_I& 1& {-v}_{\mathrm{cjI}}{x}_I& {-v}_{\mathrm{cjI}}{y}_I& {-v}_{\mathrm{cjI}}{z}_I\end{array}\right]\left[\begin{array}{c}{h}_{j11}\\ h_{j12}\\ h_{j13}\\ h_{j14}\\ h_{j21}\\ h_{j22}\\ h_{j23}\\ h_{j24}\\ h_{j31}\\ h_{j32}\\ h_{j33}\end{array}\right]=\left[\begin{array}{c}{u}_{\mathrm{cj}1}{h}_{j34}\\ v_{\mathrm{cj}1}{h}_{j34}\\ \·\·\·\\ \·\·\·\\ u_{\mathrm{cjI}}{h}_{j34}\\ v_{\mathrm{cjI}}{h}_{j34}\end{array}\right]$

等式两端同除以h_j34，得到：

${\mathrm{Ah}}_j^{\′}=B$

式中： $A=\left[\begin{array}{ccccccccccc}{x}_1& y_1& z_1& 1& 0& 0& 0& 0& {-u}_{\mathrm{cj}1}{x}_1& {-u}_{\mathrm{cj}1}{y}_1& {-u}_{\mathrm{cj}1}{z}_{21}\\ 0& 0& 0& 0& x_1& y_1& z_1& 1& {-v}_{\mathrm{cj}1}{x}_1& {-v}_{\mathrm{cj}1}{y}_1& {-v}_{\mathrm{cj}1}{z}_1\\ & \·\·\·& \·\·\·& & & & & & & & \\ & \·\·\·& \·\·\·& & & & & & & & \\ x_I& y_I& z_I& 1& 0& 0& 0& 0& {-u}_{\mathrm{cjI}}{x}_I& {-u}_{\mathrm{cjI}}{y}_I& {-u}_{\mathrm{cjI}}{z}_I\\ 0& 0& 0& 0& x_I& y_I& z_I& 1& {-v}_{\mathrm{cjI}}{x}_I& {-v}_{\mathrm{cjI}}{x}_I& {-v}_{\mathrm{cjI}}{z}_I\end{array}\right],$ 是2I×11的矩阵，B=[_ucj1,v_cj1,…,u_cjI,v_cjI]^T是2I×1的向量， $h_j^{\′}=h_j/h_{j34}$

$h_j={\left[\begin{array}{ccccccccccc}{h}_{j11}& h_{j12}& h_{j13}& h_{j14}& h_{j21}& h_{j22}& h_{j23}& h_{j24}& h_{j31}& h_{j32}& h_{j33}\end{array}\right]}^T;$

利用最小二乘法，得到：

$h_j^{\′}={\left(A^{T}A\right)}^{-1}{A}^{T}B$

步骤（4.2），按以下步骤计算表换矩阵H_(3×4)j中的各个元素：h_j1，h_j2，h_j3，h_j14，h_j24，

h_j34，其中：h_j1=[h_j11，h_j12，h_j13]^T，h_j2=[h_j21，h_j22，h_j23]^Th_j3=[h_j31，h_j32，h_j33]^T，

步骤（4.2.1），根据步骤（4.1），H_(3×4)j表达为：

$H_{(3\×4)j}=\left[\begin{array}{cc}{h}_{j1}^T& h_{j14}\\ h_{j2}^T& h_{j24}\\ h_{j3}^T& h_{j34}\end{array}\right]\left[\begin{array}{cccc}{a}_{11}& 0& a_{13}& a_{14}\\ 0& a_{22}& a_{23}& a_{24}\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}{r}_{j1}^T& t_{\mathrm{jx}}\\ r_{j2}^T& t_{\mathrm{jy}}\\ r_{j3}^T& t_{\mathrm{jz}}\\ 0^T& 1\end{array}\right]$

其中， $a_{11}=D_{\mathrm{SD}}/\mathrm{dx},$ $a_{13}=x_{\mathrm{cs}}/\mathrm{dx},$ $a_{14}=({-D}_{\mathrm{SD}}\·x_{\mathrm{cs}})/\mathrm{dx},$ $a_{22}={-D}_{\mathrm{SD}}/\mathrm{dy},$ $a_{23}=(D-y_{\mathrm{cs}})/\mathrm{dy},$ $a_{24}=(D_{\mathrm{SD}}\·y_{\mathrm{cs}})/\mathrm{dy},$

$R_j=\left[\begin{array}{c}{r_{j1}}^T\\ {r_{j2}}^T\\ {r_{j3}}^T\end{array}\right],$ 其中，r_j1=[r_j11,r_j12,r_j13]^T，r_j2=[r_j21,r_j22,r_j23]^T，r_j3=[r_j31,r_j32,r_j33]^T

t_j=[t_jx,t_jy,t_jz]^T

从而得到：

$H_{(3\×4)j}=\left[\begin{array}{cc}{a}_{11}{r}_{j1}^T+a_{13}{r}_{j3}^T& a_{11}{t}_{\mathrm{jx}}+a_{13}{t}_{\mathrm{jz}}+a_{14}\\ a_{22}{r}_{j2}^T+a_{23}{r}_{j3}^T& a_{22}{t}_{\mathrm{jy}}+a_{23}{t}_{\mathrm{jz}}+a_{24}\\ r_{j3}^T& t_{\mathrm{jz}}\end{array}\right]$

由于R_j是单位正交矩阵，则有：

$\left|\right|r_{\mathrm{jm}}\left|\right|=1$

$r_{\mathrm{jm}}^T{r}_{\mathrm{jn}}=0$

$r_{\mathrm{jm}}\×r_{\mathrm{jm}}=0$

$\left|\right|r_{\mathrm{jm}}\×r_{\mathrm{jn}}\left|\right|=1$

其中，n=1,2,3,m=1,2,3，且n≠m，分别对应于R_j的第m行和第n行，r_jm和r_jn相互正交，并且均为单位向量，

由上所述有：

$\left|\right|r_{j3}\left|\right|=\left|\right|h_{j3}\left|\right|=h_{j34}\left|\right|h_{j3}^{\′}\left|\right|=1$

其中 $h_{j3}^{\′}=h_{j3}/h_{j34}$

故有： $h_{j34}=1/\left|\right|h_{j3}^{\′}\left|\right|$

由上得到： $h_j=h_{j34}{h}_j^{\′}$

至此，就得到了矩阵H_(3×4)j中的各个元素，也就得到了h_j1，h_j2，h_j3，h_j14，h_j24，h_j34；

步骤（4.3），由步骤（4.2）得到H_(3×4)j中的各个元素后，进一步求解参数D_SD，x_cs，y_cs，R_j，t_j，具体步骤如下：

$a_{13}=(a_{11}{r}_{j1}^T+a_{13}{r}_{j3}^T)\·r_{j3}=h_{j1}^T{h}_{j3}$

$a_{23}=(a_{22}{r}_{j2}^T+a_{23}{r}_{j3}^t)r_{j3}=h_{j2}^T{h}_{j3}$

$a_{11}=\left|\right|(a_{11}{r}_{j1}+a_{13}{r}_{j3})\×r_{j3}\left|\right|=\left|\right|h_{j1}\×h_{j3}\left|\right|$

$a_{22}=\left|\right|(a_{22}{r}_{j2}+a_{23}{r}_{j3})\×r_{j3}\left|\right|=\left|\right|h_{j2}\×h_{j3}\left|\right|$

从上可得：

$D_{\mathrm{SD}}=a_{11}\·\mathrm{dx}={-a}_{22}\·\mathrm{dy}$

$x_{\mathrm{cs}}=a_{13}\·\mathrm{dx}$

$y_{\mathrm{cs}}=D-a_{23}\·\mathrm{dy}$

$r_{j1}=\frac{1}{a_{11}}(h_{j1}-a_{13}{h}_{j3})$

$r_{j3}=\frac{1}{a_{22}}(h_{j2}-a_{23}{h}_{j3})$

$r_{j3}=h_{j3}$

$R_j=\left[\begin{array}{c}{r}_{j1}^T\\ r_{j3}^T\\ r_{j3}^T\end{array}\right]$

从上可得：

$a_{14}=-(D_{\mathrm{SD}}\·x_{\mathrm{cs}})/\mathrm{dx}$ $a_{24}=D_{\mathrm{SD}}\·y_{\mathrm{cs}}/\mathrm{dy}$

因此，t_j求解如下：

$t_{\mathrm{jx}}=\frac{1}{a_{11}}(h_{14}-a_{14}-a_{13}\·h_{34})$

$t_{\mathrm{jy}}=\frac{1}{a_{22}}(h_{24}-a_{24}-a_{23}\·h_{34})$

$t_{\mathrm{jz}}=h_{34}$

$t_j={[t_{\mathrm{jx}},t_{\mathrm{jy}},t_{\mathrm{jz}}]}^T$

由上述关系，则求解出对应各个投影位置的参数D_SD，x_cs，y_cs，R_j，t_j；

步骤（5），求解锥束CT系统中所述转轴的位置和方向

利用步骤（4.3）中的参数，得到不同投影位置时标定物系的原点在成像系中的空间位置向量：

$S_j=\left[\begin{array}{cc}{R}_j& t_j\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}0\\ 0\\ 0\\ 1\end{array}\right]=t_j$

这些向量的端点轨迹是一个空间的圆周，因此通过求解S_j的均值，得到其圆心位置C，即所述转轴通过的位置：

$c=\frac{1}{j}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{S}_j$

锥束CT系统中所述转轴的方向，则为向量S_j端点的轨迹所拟合出平面的法向量，通过求圆周上三个点所构成任意两个向量的叉积而得到，本发明使用以下公式进行求解：

$Y=\frac{1}{J}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}[(S_{j+k}-S_j)\×(S_{j+1}-S_{j+k})]$

其中，k和l分别取为1/3和2/3的投影数目，

至此，确定了锥束CT系统中所述转轴的空间向量是过C点，方向为Y的向量；

步骤（6），求解成像子系统的投影主光轴，

投影主光轴是所述X射线源到所述平板探测器的垂线，在成像系中，投影主光轴的向量为O为：

$O=C_C-S_C$

其中C_C为所述X射线源在所述平板探测器上的垂足的齐次坐标，S_C为所述X射线源的空间位置的齐次坐标，S_C=[x_cs,y_cs,D_SD,1]^T，C_C=[x_cs,y_cs,0,1]^T；

步骤（7），锥束CT几何误差的求解，

锥束CT的理想几何条件为所述投影主光轴与所述转轴相交且垂直，并且所述转轴在所述平板探测器上的投影与探测器的Y_C方向平行，这里首先判断锥束CT是否符合上述理想条件，若不满足，其误差是多大，

步骤（7.1），判断所述投影主光轴与所述转轴是否垂直，若不垂直，其夹角是多大，

设a=O^TY，

若a=0，则所述投影主光轴和所述转轴垂直，满足理想情况，

若a≠0，则所述投影主光轴和所述转轴不垂直，其夹角为α：

$\mathrm{\α}={\cos }^{-1}\frac{O^{T}Y}{{\left|\right|O\left|\right|\left|\right|Y\left|\right|}^{,}}$

步骤（7.2）判断所述投影主光轴与所述转轴是否相交，若不相交，二者之间的距离是多大，

设b=O×Y，

若

则所述投影主光轴和所述转轴平行，即二者不相交，

若

则二者之间的距离d为：

$d={(C-S_C)}^{T}b$

若d=0，则二者共面且不平行，即二者相交，

若d≠0，则二者不共面，即不相交，二者之间的距离为d；

步骤（7.3），判断所述转轴在所述平板探测器上的投影是否与所述平板探测器的Y_C方向平行，若不平行，倾斜角是多少，

由步骤（5）可知，转轴上存在两点在成像系中的坐标为C和C+Y，其在探测器上的投影为

$P=\mathrm{LC}$

$P^{\′}=L(C+Y)$

其中， $L=\left[\begin{array}{cccc}{a}_{11}& 0& a_{13}& a_{14}\\ 0& a_{22}& a_{23}& a_{24}\\ 0& 0& 1& 0\end{array}\right],$ P(u_P,v_P)，P′(u_P′,v_P′)，

所述转轴在所述平板探测器上投影的倾斜角为β：

$\mathrm{\β}={\tan }^{-1}\frac{u_P-u_{P^{\′}}}{v_P-v_{P^{\′}}}$

若β=0，则所述转轴在所述平板探测器上的投影与探测器的Y_C方向平行；

至此，就完成了所述锥束CT系统的几何参数的标定，所述几何参数包括D_SD，x_cs，y_cs，以及所述投影主光轴与所述转轴之间的几何误差。

Images