本发明属于图像处理领域里的计算机层析成像(computedtomography，ct)技术领域，涉及利用具有特征纹理的标定模板对锥束ct系统进行探测器角度偏差测量的方法。
背景技术：
计算机层析成像(computedtomography，ct)系统的几何角度偏差校正问题是ct研究领域所要解决的基本问题之一。工业x射线ct系统作为一种检测设备，可对待测物体内部的尺寸、缺陷、特征进行无损检测。在工业ct系统中，射线源为x射线源，发出的光线为锥束。锥束ct系统在理想成像关系下要求探测器无角度偏差。但是，在工业ct系统的安装过程中，探测器不可避免的会产生角度偏差，而这种角度偏差会对重建图像的准确性造成影响，使图像中出现伪影，影响图像重建的质量。因此，在系统运行之前对探测器进行角度偏差的测量和校正至关重要。
为了准确的测量出锥束ct系统中探测器的角度偏差，研究人员发明了很多测量方法。目前广泛采用的测量方法是基于特定形状模板的检测方法，该方法中需要设计特定形状的模板，根据模板的先验知识，采集单一角度或多个角度的投影，根据投影的信息进行求解。目前基于特定形状模板的探测器角度偏差测量方法存在以下问题：仅对部分参数进行准确求解；模板复杂度较高，操作较困难。
技术实现要素：
本发明提供了一种基于特征纹理模板的锥束ct系统探测器角度偏差测量方法，该方法使用的标定模板中同时包含旋转对称纹理及平移对称纹理，将探测器的角度偏差与两种纹理的对称性之间建立联系，求解探测器的全部角度偏差。该方法模板制作容易，测量步骤简单，解决了模板制作困难及测量操作步骤复杂的不足。
本发明的技术方案：
基于特征纹理模板的锥束ct系统探测器角度偏差测量方法，步骤如下：
s1将标定模板置于置物台上，使用x射线照射标定模板，得到标定模板的投影图像，投影图像为非理想偏转图像i1；
s2初始化参数，探测器角度偏差初始值θx＝0,θy＝0,θz＝0，归一化位移参数初始值
s3将步骤s1中的投影图像i1和步骤s2中的初始化参数作为条件，通过低秩矩阵恢复方法，计算出θx和θy；
s4使用θx和θy将非理想偏转图像i1变换为仅包含面内转角偏差的非理想偏转图像i2，变换的方式为式(7)
s5对非理想偏转图像i2进行二维离散傅里叶变换，得到其相位谱；
s6对非理想偏转图像i2以垂直图像平面的轴进行旋转，每次旋转1度，且每次旋转获取旋转后图像的傅里叶变换幅度谱，共旋转90度；对90张旋转后的图像幅度谱进行加权平均，即式(9)
加权平均后的投影图像幅度谱中旋转对称纹理特征保持不变，平移对称纹理特征大大减少。使用该幅度谱与步骤s5中i2的傅里叶变换相位谱进行二维离散傅里叶逆变换，得到逆变换后的投影图像；
s7设定阈值，使用该阈值对s6逆变换后得到的投影图像进行处理，消除图像中少量的平移对称纹理，得到仅包含旋转对称纹理的投影图像ir；
s8从非理想偏转图像i2中减去s7中获得的投影图像ir，得到仅包含平移对称纹理的投影图像it；
s9使用低秩矩阵恢复方法对平移对称纹理图像it进行计算，得到θz，并通过ct系统中的旋钮对探测器角度进行校正，校正顺序为先绕xc轴旋转-θx，再绕yc轴旋转-θy，最后绕zw轴旋转-θz。
本发明的效果和益处是利用具有特征纹理的标定模板，采集单角度投影数据求解得出锥束ct系统中探测器的角度偏差，模板制作容易，操作简单，解决了模板制作复杂及测量操作复杂的不足。
附图说明
图1为锥束ct几何成像系统示意图。s为射线源，图中ow-xwywzw为世界坐标系，uov为图像坐标系，xcocyc为探测器成像平面坐标系。世界坐标系用于描述空间中点的位置。图像坐标系是以像素为单位，u轴指向探测器上的列数增大的方向，v轴指向探测器上的行数增大的方向。探测器成像平面坐标系表示与图像坐标系平行，原点oc在探测器中心。s'为射线源在图像坐标系上的投影点。
图2为锥束ct系统中，探测器的角度偏差。图2(a)为探测器围绕zw轴的角度偏差，即探测器面内偏转角；图2(b)为探测器围绕xc轴的角度偏差，即探测器面外俯仰角；图2(c)为探测器围绕yc轴的角度偏差，即探测器面外偏转角。
图3为本发明使用的标定模板示意图。
图4为不同情况下的投影图像。图4(a)表示探测器三个角度偏差均不为0时的非理想偏转图像，图4(b)表示仅包含面内偏转角偏差的非理想偏转图像，图4(c)表示理想正对图像。
图5为投影图像的fe变换过程。图5(a)表示理想正对图像的fe变换，图5(b)表示非理想偏转图像的fe变换。
图6为标定模板纹理组成示意图。
具体实施方式
以下结合技术方案和附图详细叙述本发明的具体实施方式。
基于特征纹理模板的锥束ct系统探测器角度偏差测量方法，步骤如下：
(1)基本概念
标定模板：用于测量锥束ct系统探测器角度偏差的物体；
模板纹理：模板中的图案为棋盘格纹理与同心圆环纹理的叠加纹理，二者中心重合；
平移对称纹理：模板中的棋盘格纹理；
旋转对称纹理：模板中的同心圆环纹理；
投影图像：实际场景中标定模板经过x射线照射在探测器中所成的像；
平移对称纹理图像：包含平移对称纹理的投影图像；
旋转对称纹理图像：包含旋转对称纹理的投影图像；
非理想偏转图像：探测器角度偏差不为零时标定模板的投影图像；
理想正对图像：探测器无角度偏差时标定模板的投影图像；
(2)技术方案的具体原理：
锥束ct系统中探测器的角度偏差不同，投影图像也不同，本发明根据标定模板的投影图像来计算探测器的角度偏差。本发明的基本原理为：模板纹理由平移对称纹理和旋转对称纹理组合而成，两种纹理的中心重合，且平移对称纹理的尺寸与旋转对称纹理的尺寸相对固定。拍摄一副标定模板的投影图像，由于投影图像中包含旋转对称纹理图像，可以建立旋转对称纹理的对称性与经过fe变换后的投影图像的秩之间的关系，采用低秩矩阵恢复原理，求解出探测器部分角度偏差。同时，投影图像中还包含平移对称纹理图像，可以建立平移对称纹理的对称性与投影图像秩之间的关系，采用低秩矩阵恢复原理，求解出探测器其余角度偏差。下面将对技术方案原理进行详细介绍，首先介绍锥束ct几何成像系统的定义及理想正对图像与非理想偏转图像之间的关系，然后介绍探测器角度偏差求解。
1)锥束ct几何成像系统定义及投影图像之间的关系
锥束ct几何成像系统，包括x射线源、置物台、探测器和标定模板；s为x射线源焦点，在标定模板平面上建立世界坐标系ow-xwywzw，原点ow为标定模板平面的中心点，yw轴为旋转轴，xwowyw平面为标定模板平面；在探测器平面上建立探测器成像坐标系xcoyc，其中，原点oc为中心射线在探测器上的投影位置，xc轴平行于探测器平面水平方向，yc轴平行于探测器平面垂直方向；在探测器平面上建立图像坐标系uov，原点o位于探测器平面左上角，u轴平行于xc轴，v轴平行于yc轴；s点到xwowyw平面的距离为r,s点在图像坐标系uov内的坐标为(u0,v0)；探测器通过旋钮进行绕xc轴、yc轴、zw轴角度的调节，探测器绕xc轴旋转角度为面外俯仰角θx，绕yc轴旋转角度为面外偏转角θy，绕zw轴旋转角度为面内偏转角θz；
对标定模板进行x射线照射，标定模板上某点(xw,yw,zw)在探测器上得到的对应投影点坐标为(u,v)，用齐次坐标表示两者的映射关系为：
其中，s为一个常数因子，n为内参数矩阵，d表示射线源到探测器的距离，(u0,v0)表示射线源焦点在图像坐标系中的坐标，px和py分别表示探测器u轴方向和v轴方向的像元尺寸，本发明中令px＝py；t是平移矩阵，r是一个3×3的旋转矩阵；
定义探测器的旋转顺序为：先绕zw轴顺时针旋转，再绕yc轴顺时针旋转，最后绕xc轴顺时针旋转，用探测器面外俯仰角θx、面外偏转角θy及面内偏转角θz表示旋转矩阵r为：
标定模板的投影图像是由平移对称纹理图像(棋盘格)和旋转对称纹理图像(同心圆环)组合而成，两种纹理的中心重合，且棋盘格纹理中方格的长宽尺寸与同心圆环纹理的半径尺寸相对固定；使用x射线照射标定模板，得到非理想偏转图像i1，如图4(a)所示。假设探测器角度偏差θx、θy为0，得到的投影图像为仅包含θz偏差的非理想偏转图像i2，如图4(b)所示。无角度偏差时，投影得到的图像为理想正对图像i3，如图4(c)所示。世界坐标系中同一点p(xw,yw,zw)在三种投影图像中的投影点分别为点p1(u1,v1)、p2(u2,v2)与p3(u3,v3)，则有
其中，sa为射线源到旋转轴距离；则根据式(3)、(4)、(5)对三种投影图像进行转换；
2)探测器角度偏差求解
当标定模板的投影图像为理想正对图像，其面外俯仰角θx与面外偏转角θy均为0，此时同心圆环纹理经过fe变换后为平移对称纹理，而平移对称纹理为低秩纹理；当投影图像为非理想偏转图像，其θx与θy不为0，经过fe变换后的纹理不是平移对称纹理，不具有低秩性；
使用的特征纹理模板的fe变换如图5所示。其中(a)为理想正对图像的fe变换结果，(b)表示非理想偏转图像的fe变换结果。对特征纹理模板投影图像进行fe变换，令变换后的图像大小为m×n。由图5(a)可知，当投影图像为理想正对图像i3时，i3可分为8个图像块，每个图像块均可通过一个图像块经过旋转、翻转、平移得到，则其fe变换图像i3f具有周期性，且周期t＝8。根据图像的列秩的定义，对于i3f的第n/8+1列至第n列像素，均可由第1列至第n/8列像素经过线性相关得到，因此rank(i3f)＝n/8。而对于图5(b)，投影图像i2为非理想偏转图像，其fe变换图像i2f不具有周期性，rank(i2f)>n/8。所以，当投影图像为理想正对图像时，其秩最小，可用低秩矩阵恢复原理式(6)对非理想偏转图像进行计算，得到探测器的面外俯仰角θx和面外偏转角θy。
其中，i表示投影图像，i0表示探测器无角度偏差时的投影图像，e为噪声矩阵，f{}表示fe变换，τ表示探测器角度偏差。
由式(3)(4)可得探测器三个角度偏差均不为0时点p1与探测器角度偏差中θx＝0,θy＝0,θz≠0时点p2之间转换关系为
在得到探测器的θx和θy之后，通过式(7)将三个角度偏差均不为0的投影图像转换为θx＝0,θy＝0,θz≠0的投影图像。此时的投影图像中仅包含探测器的面内偏转角θz偏差，可用于θz的测量。
本发明中使用的特征纹理模板投影图像是由两种具有不同特征的纹理图像组合而成的，理想正对图像可以表示如下：
i3＝it+ir(8)
其中it表示平移对称纹理图像，即棋盘格纹理图像。ir表示旋转对称纹理图像，即同心圆环纹理图像。i3表示组合后的混合纹理图像，如图6所示。
对于平移对称纹理图像it，若以垂直于图像平面的轴进行旋转，纹理会发生变化，其傅里叶变换幅度谱也会发生变化，即其傅里叶变换幅度谱的变化反应了其纹理特征的变化。若取权重为1，对经过不同角度旋转后的平移对称纹理图像傅里叶变换幅度谱进行加权平均，得到的幅度谱中包含的平移对称纹理特征会大大减少。而对于旋转对称纹理图像ir，若进行同样的旋转，纹理不会发生变化，其傅里叶变换幅度谱也不会发生明显变化。因此，对经过不同角度旋转后的旋转对称纹理图像傅里叶变换幅度谱进行加权平均，得到的幅度谱仍具有丰富的旋转对称纹理特征。
若对混合纹理图像i3绕垂直于图像的轴进行多次旋转，并对旋转后图像的傅里叶变换幅度谱进行加权平均，由式(8)可得
其中
实施例
射线源为commet射线源，使用paxscan1313dx探测器采集投影图像，分辨率为1024×1024，像元尺寸为0.127mm，接口为千兆网口。射线源到探测器估计距离为1000mm。标定模板上的纹理为棋盘格纹理与同心圆环纹理的叠加，模板大小为100mm×100mm，棋盘格为10×10，每个小方格长宽均为8mm，同心圆环共5个圆环，以棋盘格中心为圆心，最小半径为8mm，其它半径为以最小半径为基础，间隔8mm递增。模板上每条线宽为1mm。
步骤1、获得标定模板投影图像，确定初始值。角度参数初始值θx＝0,θy＝0,θz＝0，归一化位移参数初始值
步骤2、对投影图像进行fe变换，使用矩阵低秩恢复方法计算出面外俯仰角θx和面外偏转角θy。
步骤3、使用θx和θy将投影图像变换为仅包含θz偏差的非理想偏转图像。
步骤4、对非理想偏转图像进行二维离散傅里叶变换，得到其相位谱。
步骤5、对非理想偏转图像进行旋转，每次旋转1度，共旋转90次。每次旋转计算其傅里叶变换幅度谱并记录。对90张幅度谱进行加权平均，使用得到的幅度谱与步骤4中的相位谱进行二维离散傅里叶逆变换，得到包含丰富旋转对称纹理特征与少量平移对称纹理特征的图像。
步骤6、设定阈值，对步骤5中得到的图像进行处理，消除图像中的少量平移对称纹理，得到旋转对称纹理图像。
步骤7、从仅包含θz偏差的非理想偏转图像中减去得到的旋转对称纹理图像，分离出平移对称纹理图像。
步骤8、使用低秩矩阵恢复方法对分离出的平移对称纹理进行计算，得到面内转角θz。
根据求得的探测器角度偏差的值θx、θy、θz调节探测器，通过ct系统中的旋钮对探测器角度进行校正，校正顺序为先绕xc轴旋转-θx，再绕yc轴旋转-θy，最后绕zw轴旋转-θz，此时探测器的角度偏差校正完毕。