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

摘要

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

主权项

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_cx=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向下为正方向，则：uci=xcidx,]]>vci=D-ycidy,]]>其中：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,…,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，标定物系到成像系的空间转换矩阵Tj=Rjtj0T1]]>中的各所述投影图像的旋转矩阵R_j和位移向量t_j，其中，D_SD、x和y_cs是所述锥束CT系统的内方位参数，在所述成像子系统制作完成后是不变的，但又是未知的，而对于一个特定的拍摄位置，对应于该拍摄位置的空间转换矩阵T是不变的；步骤（4.1），按下式建立所述标定物系内已知空间坐标点P_i(x_i,y_i,z_i)与各所述投影图像j中的对应坐标点P_cji(u_cji，v_cji)之间的关系：zci′ucjivcji1=H(3×4)jxiyizi1,]]>H_(3×4)j为所述以像素为单位的坐标与标定物系坐标之间的转换关系，H(3×4)j=DSD/dx0xcs/dx-DSD·xcs/dx0-DSD/dy(D-ycs)/dyDSD·ycs/dy0010Rjtj01]]>=hj11hj12hj13hj14hj21hj22hj23hj24hj31hj32hj33hj34]]>H_(3×4)j为第j幅投影图像的所述变换矩阵，R_j为第j幅投影图像的旋转矩阵，是单位正交矩阵，t_j为第j幅投影图像的位移向量，Rjt0T1]]>为第j幅投影图像的空间转换矩阵，对于所述三维标定模块的I=12个空间标志点，得到下述对应于第j个投影图像的，包含2I个方程的联立方程组：x1y1z110000-ucj1x1-ucj1y1-ucj1z10000x1y1z11-vcj1x1-vcj1y1-vcj1z1............xIyIzI10000-ucjIxI-ucjIyI-ucjIzI0000xIyIzI1-vcjIxI-vcjIyI-vcjIzIhj11hj12hj13hj14hj21hj22hj23hj24hj31hj32hj33=ucj1hj34vcj1hj34......ucjIhj34vcjIhj34]]>等式两端同除以h_j34，得到：Ah′_j=B式中：A=x1y1z110000-ucj1x1-ucj1y1-ucj1z10000x1y1z11-vcj1x1-vcj1y1-vcj1z1............xIyIzI10000-ucjIxI-ucjIyI-ucjIzI0000xIyIzI1-vcjIxI-vcjIyI-vcjIzI,]]>是2I×11的矩阵，B=[u_cj1,v_cj1,…,u_cjI,v_cjI]^T是2I×1的向量，h′_j＝h_j/h_j34，h_j=[h_j11 h_j12 h_j13 h_j14 h_j21 h_j22 h_j23 h_j24 h_j31 h_j32 h_j33]^T；利用最小二乘法，得到：h′_j=(A^TA)^-1A^TB步骤（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=hj1Thj14hj2Thj24hj3Thj34=a110a13a140a22a23a240010rj1Ttjxrj2Ttjyrj3Ttjz0T1]]>其中，a₁₁=D_SD/dx，a₁₃=x_cs/dx，a₁₄＝(-D_SD·x_cs)/dx，a₂₂=-D_SD/dy，a₂₃＝(D-y_cs)/dy，a₂₄＝(D_SD·y_cs)/dy，Rj=rj1Trj2Trj3T,]]>其中，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]^Tt_j=[t_jx,t_jy,t_jz]^T从而得到：H(3×4)j=a11rj1T+a13rj3Ta11tjx+a13tjz+a14a22rj2T+a23rj3Ta22tjy+a23tjz+a24rj3Ttjz]]>由于R_j是单位正交矩阵，则有：||r_jm||＝1rjmTrjn=0]]>r_jm×r_jm=0||r_jm×r_jn||＝1其中，n=1,2,3,m=1，2,3，且n≠m，分别对应于R_j的第m行和第n行，r_jm和r_jn相互正交，并且均为单位向量，由上所述有：||r_h3||＝||h_j3||=h_j34||h′_j3=1其中h′_j3＝h_j3/h_j34故有：h_j34＝1/||h′_j3‖由上得到：h_j=h_j34h′_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，具体步骤如下：a13=(a11rj1T+a13rj3T).rj3=hj1Thj3]]>a23=(a22rj2T+a23rj3T).rj3=hj2Thj3]]>a₁₁=‖(a₁₁r_j1+a₁₃r_j3)×r_j3‖=‖h_j1×h_j3‖a₂₂=‖(a₂₂r_j2+a₂₃r_j3)×r_j3‖=‖h_j2×h_j3‖从上可得：D_SD=a₁₁·dx=-a₂₂·dyx_cs=a₁₃·dxy_cs=D-a₂₃·dyrj1=1a11(hj1-a13hj3)]]>rj2=1a22(hj2-a23hj3)]]>r_j3＝h_j3Rj=rj1Trj2Trj3T]]>从上可得：a₁₄＝-(D_SD·x_cs)/dxa₁₄＝-(D_SD·x_cs)/dxa₂₄＝D_SD·y_cs/dya₂₄＝D_SD·y_cs/dy因此，t_i求解如下：tjx=1a11(h14-a14-a13.h34)]]>tjy=1a22(h24-a24-a23.h34)]]>t_jz=h₃₄t_j=[t_jx，t_jy，t_jz]^T由上述关系，则求解出对应各个投影位置的参数D_SD，x_cs，y_cs，R_j，t_j；步骤（5），求解锥束CT系统中所述转轴的位置和方向利用步骤（4.3）中的参数，得到不同投影位置时标定物系的原点在成像系中的空间位置向量：Sj=Rjtj0T10001=tj]]>这些向量的端点轨迹是一个空间的圆周，因此通过求解S_j的均值，得到其圆心位置C，即所述转轴通过的位置：C=1JΣj=1JSj]]>锥束CT系统中所述转轴的方向，则为向量S_j端点的轨迹所拟合出平面的法向量，通过求圆周上三个点所构成任意两个向量的叉积而得到，本发明使用以下公式进行求解：Y=1JΣj=1J[(Sj+k-Sj)×(Sj+1-Sj+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，则所述投影主光轴和所述转轴不垂直，其夹角为α：α=cos-1OTY||O||||Y||,]]>步骤（7.2）判断所述投影主光轴与所述转轴是否相交，若不相交，二者之间的距离是多大，设b=O×Y，若||b‖=0，则所述投影主光轴和所述转轴平行，即二者不相交，若||b‖≠0，则二者之间的距离d为：d=(C-S_C)^Tb若d=0，则二者共面且不平行，即二者相交，若d≠0，则二者不共面，即不相交，二者之间的距离为d；步骤（7.3），判断所述转轴在所述平板探测器上的投影是否与所述平板探测器的Y_C方向平行，若不平行，倾斜角是多少，由步骤（5）可知，转轴上存在两点在成像系中的坐标为C和C+Y，其在探测器上的投影为P=LCP′=L(C+Y)其中，L=a110a13a140a22a23a240010,]]>P(u_P,v_P)，P′(u_P′,v_P′)，所述转轴在所述平板探测器上投影的倾斜角为β：β=tan-1uP-uP′vP-vP′]]>若β=0，则所述转轴在所述平板探测器上的投影与探测器的Y_C方向平行；至此，就完成了所述锥束CT系统的几何参数的标定，所述几何参数包括D_SD，x_cs，y_cs，以及所述投影主光轴与所述转轴之间的几何误差。