CN104132950A - 基于原始投影信息的cl扫描装置投影旋转中心标定方法 - Google Patents

Abstract

本发明公开了一种基于原始投影信息的CL扫描装置投影旋转中心标定方法，本方法先将被扫描对象360度旋转范围内的投影信息合成为一幅椭圆投影P(x,y)，将该椭圆投影进行水平镜像得到一新的椭圆投影P(x,y)^Mirror，求取P(x,y)和P(x,y)^Mirror沿X方向的互相关函数R(τ,y)；然后将R(τ,y)沿y方向进行积分，消去y变量，得到一维函数Q^ACF(2x₀-N+τ)，求取该一维函数的最大值坐标，从而换算出投影旋转中心的位置坐标。本发明方法基于椭圆投影的对称性原理，不需要制作专门的标定模体，不需要对原始投影数据进行重排，直接利用互相关函数的最大值进行标定，标定方法具有结果唯一、抗噪性强、精度高的优点，有效保证了重建图像的精度。

Description

基于原始投影信息的CL扫描装置投影旋转中心标定方法

技术领域

本发明涉及一种适用于CL扫描装置投影旋转中心的标定，可用于医学和工业领域射线数字成像(DR-Digital Radiography)、常规锥束射线层析扫描及非常规锥束射线层析扫描成像过程中的相关测量。

背景技术

CL(Computed Tomography)是一种专门针对板状构件的层析扫描装置。该扫描装置克服了传统CT(Computed Tomography)装置中心射线必须垂直于物体旋转轴的限制，采用射线与物体旋转轴倾斜的方式，使得大部分射线穿过物体的路径长度接近于物体的厚度尺寸，从而大大降低了射线的能量，有效保证了投影图的灵敏度。该扫描装置的结构如图1所示。X射线源1发出锥束射线2，被检测样品3放置在扫描平台4上，并在扫描平台4的带动下沿X、Y、Z方向平移，以及绕旋转轴线5转动；平板探测器6安装在扫描平台4下方的弧形导轨8上。层析扫描成像开始时，被检测样品3在扫描平台4的拖动下绕旋转轴5在360度范围内旋转，平板探测器6采集被检测样品3在不同旋转位置的投影图像；基于360度旋转范围内的所有投影图像，通过CL层析扫描装置的图像重建与图像处理单元反演出被检测板状构件的断层图像。

参见图2所示，当被检测样品3被锥束射线2照射时，其透视投影会投射在平板探测器6的成像平面7上，同时，被检测样品的旋转轴线5也会被投射到成像平面7上。根据CL扫描层析重建算法要求，需要知道旋转轴线5在成像平面7上的精确位置，即旋转轴线投影9的精确位置，否则会导致图像重建精度的降低，从而影响被检测对象内部结构与质量状态信息的准确输出，对产品的质量评判与结构逆求造成不利。然而，对于实际的CL扫描层析装置，旋转轴线5与成像平面7均虚拟不可见，因此无法直接测量到旋转轴线5在成像平面7上的精确投影位置，该问题已成为决定CL层析图像质量的主要因素。

发明内容

本发明的目的是提出一种实用的CL扫描装置投影旋转中心的标定方法。本方法基于椭圆投影的对称性原理，将被检测样品360旋转范围内的投影信息合成为一幅椭圆投影P(x,y)，将该椭圆投影进行水平镜像得到一新的椭圆投影P(x,y)^Mirror，求取椭圆投影P(x,y)和P(x,y)^Mirror沿X方向的二维互相关函数R(τ,y)；将函数R(τ,y)沿y方向进行积分，消去y变量，得到一维函数Q^ACF(2x₀-N+τ)，求取该一维函数的最大值坐标，从而换算出投影旋转中心的位置坐标。

本发明是一种基于原始投影信息的CL扫描装置投影旋转中心标定方法，所述CL扫描装置包括有射线源(1)、扫描平台(4)、平板探测器(6)、弧形导轨(8)；其特征在于该标定方法包括有下列实施步骤：

步骤一：将被检测样品(3)置于扫描平台(4)上，调整扫描平台的升降位置，以保证被检测样品的投影不超出平板探测器(6)的成像平面(7)；

步骤二：启动CL扫描装置的旋转控制指令，被检测样品(3)在扫描平台(4)的带动下旋转360度；扫描平台(4)步进旋转，步进角为0.5°～1°；扫描平台(4)每旋转过一个步进角，平板探测器(7)采集样品的一个投影，所述投影记为Prⁱ(x,y)，x表示X轴上的坐标变量，y表示Y轴上的坐标变量，其中i表示投影的编号，即i＝1,2,3,......T，T代表投影的总数；

步骤三：被检测样品(3)旋转360度后，平板探测器(7)将采集到的投影合成为一幅图像，即椭圆投影轨迹(11)，所述椭圆投影轨迹记为P(x,y)；P(x,y)满足的关系为 $P(x,y)=\frac{1}{T}\underset{i=1}{\overset{T}{\mathrm{\Σ}}}{\mathrm{Pr}}^i(x,y);$

步骤四：将P(x,y)沿X轴方向进行镜像变换，得到另一椭圆投影P(x,y)^Mirror，求取P(x,y)和P(x,y)^Mirror在X方向上的二维互相关函数R(τ,y)，R(τ,y)的计算式为：

$R(\mathrm{\τ},y)=\frac{1}{N}\underset{0}{\overset{N}{\∫}}P(x,y)P(2x_0-N+x+\mathrm{\τ},y)\mathrm{dx}$

步骤五：将二维互相关函数R(τ,y)沿y方向进行积分，消去y变量，得到一维函数Q^ACF(2x₀-N+τ)，函数Q^ACF(2x₀-N+τ)的值随变量τ的变化而变化，该一维函数Q^ACF(2x₀-N+τ)的计算式如为：

$Q^{\mathrm{ACF}}({2x}_0-N+\mathrm{\τ})=\frac{1}{M}\underset{0}{\overset{M}{\∫}}{R}^{\mathrm{ACF}}({2x}_0-N+\mathrm{\τ},y)\mathrm{dy}$

$\begin{array}{c}{R}^{\mathrm{ACF}}({2x}_0-N+\mathrm{\τ},y)=R(\mathrm{\τ},y)\\ =\frac{1}{N}\underset{0}{\overset{N}{\∫}}P(x,y)P({2x}_0-N+x+\mathrm{\τ},y)\mathrm{dx}\end{array}$

步骤六：统计函数Q^ACF(2x₀-N+τ)的最大值，当Q^ACF(2x₀-N+τ)取最大值时，自变量τ的对应值记为τ'，则有投影旋转中心x₀的值为

本发明标定方法的优点在于：

①该标定方法基于椭圆投影的对称性原理，不需要将原始投影数据进行重排，直接对原始投影信息进行互相关运算得到函数的最大值坐标，具有操作简单、适合工程应用的优点。

②该不需要制作专门的标定模体，直接利用被检测样品的投影信息进行标定，不仅降低了成本，也有效减缩短了系统标定时间。

③将数据叠加与互相关方法应用到测量过程中，具有结果唯一、抗噪性强、精度高的优点，有效保证了重建图像的空间分辨率。

④该发明方法可用软件自动实现，可作为板状构件层析扫描装置的系统自检模块，自动监测扫描装置投影旋转中心的变化。

附图说明

图1是CL扫描装置结构图。

图2是CL扫描几何原理图。

图3是被检测样品的一幅投影图像。

图4是合成投影图像。

图5是一维互相关函数曲线图。

图6A是利用本发明标定值的断层重建图像。

图6B是利用系统原始标定值的断层重建图像。

1.射线源	2.锥束射线	3.被检测样品
			4.扫描平台	5.旋转轴线	6.平板探测器
7.成像平面	8.弧形导轨	9.旋转轴线投影
			10.中心射线	11.椭圆投影轨迹	12.样品投影

具体实施方式

下面将结合附图和实施例对本发明做进一步的详细说明。

参见图2所示的CL扫描几何原理图，该图不同于传统的CT系统扫描方式，CL扫描方式下旋转轴线5与中心射线10的夹角小于90度(α＜90°)，而传统的CT系统旋转轴线则垂直于中心射线。

在本发明中，旋转轴线5在成像平面7上的投影为9(即旋转轴线投影)，投影平面坐标系定义为XOY，该坐标系也是投影图像的坐标系XOY。由投影几何关系可知，旋转轴线投影9垂直于X轴，因此，旋转轴线投影9在成像平面7上的位置可以用该轴线与的X轴的交点定义，交点记为点A，该交点的坐标记为A(x₀,0)。根据CL层析算法的要求，需要精确知道投影旋转中心点x₀的值。

进一步分析可知，当被检测样品3在扫描平台4的拖动下绕旋转轴线5转动一周后，其运动轨迹为一圆；当被检测样品3被投射到成像平面7上时，形成的投影记为样品投影12。可以看出，当被检测样品3绕轴线旋转时，样品投影12也随之在成像平面7上转动，只不过样品投影12的运动轨迹为一椭圆，该椭圆投影轨迹记为11。

进一步分析可知，椭圆投影轨迹11的长轴和短轴分别平行于投影图像坐标系XOY的X轴和Y轴，并且椭圆的短轴与旋转轴线投影9重合。因此，求取旋转轴线投影9与X轴的交点坐标A(x₀,0)等同于求样品投影12的运动轨迹的中心点B在X轴上的坐标。

由以上分析可知，样品投影12的运动轨迹为一椭圆投影轨迹11，这里记为P(x,y)。由于椭圆投影为对称图形，因此在X轴方向上该椭圆投影存在对称中心，该对称中心即为上面所述的投影旋转中心点x₀。因此存在以下关系式：

P(x,y)＝P(2x₀-x,y)        (1)

x表示X轴上的坐标变量，y表示Y轴上的坐标变量，x₀表示投影旋转中心点。

对所述的椭圆投影轨迹P(x,y)沿X轴方向进行镜像变换处理，得到X轴镜像—椭圆投影，记为P(x,y)^Mirror。

在本发明中，X轴镜像—椭圆投影P(x,y)^Mirror表达式为：

P(x,y)^Mirror＝P(N-x,y)    x∈[0 N]    (2)

x表示X轴上的坐标变量，y表示Y轴上的坐标变量，N表示椭圆轨迹投影P(x,y)在X轴方向上的长度。

在本发明中，求取椭圆投影轨迹P(x,y)和X轴镜像—椭圆投影P(x,y)^Mirror在X方向上的二维互相关函数，记为R(τ,y)，R(τ,y)的表达式为：

$R(\mathrm{\τ},y)=\frac{1}{N}\underset{0}{\overset{N}{\∫}}P(x,y)P{(x+\mathrm{\τ},y)}^{\mathrm{Mirror}}\mathrm{dx}---\left(3\right)$

τ表示互相关函数自变量，dx表示坐标变量x的差分变量。

由公式(2)得公式(3)中的P(x+τ,y)^Mirror可等效为：

P(x+τ,y)^Mirror＝P(N-x-τ,y)        (4)

联立公式(1)和(4)可得：

P(x+τ,y)^Mirror＝P(2x₀-N+x+τ,y)          (5)

将式(5)代入公式(3)得：

$R(\mathrm{\τ},y)=\frac{1}{N}\underset{0}{\overset{N}{\∫}}P(x,y)P(2x_0-N+x+\mathrm{\τ},y)\mathrm{dx}---\left(6\right)$

公式(6)亦为函数P(x,y)关于2x₀-N+τ的在X轴方向的自相关函数，该自相关函数记为R^ACF(2x₀-N+τ,y)，即：

$\begin{array}{c}{R}^{\mathrm{ACF}}({2x}_0-N+\mathrm{\τ},y)=R(\mathrm{\τ},y)\\ =\frac{1}{N}\underset{0}{\overset{N}{\∫}}P(x,y)P({2x}_0-N+x+\mathrm{\τ},y)\mathrm{dx}\end{array}---\left(7\right)$

对函数R^ACF(2x₀-N+τ,y)进行降维处理，消去y变量，得到一个一维函数Q^ACF(2x₀-N+τ)，其表达式为：

$Q^{\mathrm{ACF}}({2x}_0-N+\mathrm{\τ})=\frac{1}{M}\underset{0}{\overset{M}{\∫}}{R}^{\mathrm{ACF}}({2x}_0-N+\mathrm{\τ},y)\mathrm{dy}---\left(8\right)$

dy表示坐标变量y的差分变量，M表示椭圆轨迹投影P(x,y)在Y轴方向上的长度。

由自相关函数的特性可知，当Q^ACF(2x₀-N+τ)取最大值时，自变量τ的对应值记为τ'，此时满足以下条件：

2x₀-N+τ＝0|_τ＝τ'        (9)

由公式(9)可以得出投影旋转中心x₀的计算公式：

$x_0=\frac{N-{\mathrm{\τ}}^{\′}}{2}---\left(10\right)$

本发明是一种实用的、不需要任何专用标定模体的、基于原始投影信息的CL扫描装置投影旋转中心标定技术，结合上述内容，该投影旋转中心标定方法包括有下列实施步骤：

步骤一：将被检测样品3置于扫描平台4上，调整扫描平台的升降位置，以保证被检测样品的投影不超出平板探测器6的成像平面7；

步骤二：启动CL扫描装置的旋转控制指令，被检测样品3在扫描平台4的带动下旋转360度；

扫描平台4步进旋转，步进角为0.5°～1°；

扫描平台4每旋转过一个步进角，平板探测器7则采集被检测样品3的一个投影，所述投影记为Prⁱ(x,y)，其中i表示投影的编号，即i＝1,2,3,......T，T代表投影的总数；

步骤三：被检测样品3旋转360度后，平板探测器7将采集到的投影合成为一幅图像，该图像即为图2中的椭圆投影轨迹11，所述椭圆投影轨迹记为P(x,y)。P(x,y)的计算式为 $P(x,y)=\frac{1}{T}\underset{i=1}{\overset{T}{\mathrm{\Σ}}}{\mathrm{Pr}}^i(x,y);$

步骤四：将P(x,y)沿X轴方向进行镜像变换处理，得到X轴镜像—椭圆投影P(x,y)^Mirror；

然后，求取P(x,y)和P(x,y)^Mirror在X方向上的二维互相关函数R(τ,y)， $R(\mathrm{\τ},y)=\frac{1}{N}\underset{0}{\overset{N}{\∫}}P(x,y)P{(x+\mathrm{\τ},y)}^{\mathrm{Mirror}}\mathrm{dx};$

步骤五：将二维互相关函数R(τ,y)沿y方向进行积分，消去y变量，得到一维函数Q^ACF(2x₀-N+τ)，该一维函数 $Q^{\mathrm{ACF}}({2x}_0-N+\mathrm{\τ})=\frac{1}{M}\underset{0}{\overset{M}{\∫}}{R}^{\mathrm{ACF}}({2x}_0-N+\mathrm{\τ},y)\mathrm{dy},$ 可以看出，函数Q^ACF(2x₀-N+τ)的值随变量τ的变化而变化。

x₀表示投影旋转中心点；

N表示椭圆轨迹投影P(x,y)在X轴方向上的长度；

τ表示互相关函数自变量；

M表示椭圆轨迹投影P(x,y)在Y轴方向上的长度；

y表示Y轴上的坐标变量；

dy表示坐标变量y的差分变量。

步骤六：统计函数Q^ACF(2x₀-N+τ)的最大值，当Q^ACF(2x₀-N+τ)取最大值时，自变量τ的对应值记为τ'，则有投影旋转中心x0的值为

实施案例

实验对AL-CT1303-CL型CL扫描装置的投影旋转中心重新进行标定，该装置的主要配置如下：

(1)射线源：日本产L9181S型130kV微焦点射线源，焦点尺寸5μm；

(2)面阵探测器：美国Varian公司产PaxScan1313型非晶硅面阵探测器，成像面积130×130mm²，探元尺寸0.127mm；

(3)射线源焦点到成像平面距离820mm。

实验步骤：

(1)将一被检测样品3放置在扫描平台4上，启动扫描平台旋转带动被检测样品步进旋转360度，扫描平台4的步进角为1°，共采集360(T＝360)幅投影图像。图3为被检测样品的一幅投影图像。将360幅投影图像叠加为一幅合成图像，如图4所示，该合成图像即为P(x,y)。

(2)将合成图像P(x,y)沿X轴方向进行镜像变换，得到X轴镜像—投影图像P(x,y)^Mirror，求取P(x,y)和P(x,y)^Mirror在X方向上的二维互相关函数R(τ,y)，将该二维互相关函数R(τ,y)沿y方向进行积分，消去y变量，得到关于变量τ的一维函数，该一维函数即为Q^ACF(2x₀-N+τ)，图5为Q^ACF(2x₀-N+τ)的曲线图。

(3)统计函数Q^ACF(2x₀-N+τ)的最大值，找到Q^ACF(2x₀-N+τ)为最大值时的τ值对应的值，记为τ'，由式(10)计算出投影旋转中心x₀的值。本次实验计算出的x₀值为536.540，而该AL-CT1303-CL型CL扫描装置x₀的初始标定值为536(出厂设置)。分别用x₀＝536.540和x₀＝536进行图像重建，其结果如图6所示，图6A为利用本发明标定值的重建结果，图6B为利用初始标定值的重建结果，可以看出，本发明的标定精度更高，重建图像更为清晰。

Claims (5)

1.一种基于原始投影信息的CL扫描装置投影旋转中心标定方法，所述CL扫描装置包括有射线源(1)、扫描平台(4)、平板探测器(6)、弧形导轨(8)；其特征在于该标定方法包括有下列实施步骤：

步骤一：将被检测样品(3)置于扫描平台(4)上，调整扫描平台的升降位置，以保证被检测样品的投影不超出平板探测器(6)的成像平面(7)；

步骤二：启动CL扫描装置的旋转控制指令，被检测样品(3)在扫描平台(4)的带动下旋转360度；扫描平台(4)步进旋转，步进角为0.5°～1°；扫描平台(4)每旋转过一个步进角，平板探测器(7)采集样品的一个投影，所述投影记为Prⁱ(x,y)，x表示X轴上的坐标变量，y表示Y轴上的坐标变量，其中i表示投影的编号，即i＝1,2,3,......T，T代表投影的总数；

步骤三：被检测样品(3)旋转360度后，平板探测器(7)将采集到的投影合成为一幅图像，即椭圆投影轨迹(11)，所述椭圆投影轨迹记为P(x,y)；P(x,y)满足的关系为 $P(x,y)=\frac{1}{T}\underset{i=1}{\overset{T}{\mathrm{\Σ}}}{\mathrm{Pr}}^i(x,y);$

步骤四：将P(x,y)沿X轴方向进行镜像变换，得到另一椭圆投影P(x,y)^Mirror，求取P(x,y)和P(x,y)^Mirror在X方向上的二维互相关函数R(τ,y)，R(τ,y)的计算式为：

$R(\mathrm{\τ},y)=\frac{1}{N}\underset{0}{\overset{N}{\∫}}P(x,y)P(2x_0-N+x+\mathrm{\τ},y)\mathrm{dx}$

x表示X轴上的坐标变量；

y表示Y轴上的坐标变量；

τ表示互相关函数自变量；

x₀表示投影旋转中心点；

dx表示坐标变量x的差分变量；

N表示椭圆轨迹投影P(x,y)在X轴方向上的长度；

步骤五：将二维互相关函数R(τ,y)沿y方向进行积分，消去y变量，得到一维函数Q^ACF(2x₀-N+τ)，函数Q^ACF(2x₀-N+τ)的值随变量τ的变化而变化，该一维函数Q^ACF(2x₀-N+τ)的计算式如为：

$Q^{\mathrm{ACF}}({2x}_0-N+\mathrm{\τ})=\frac{1}{M}\underset{0}{\overset{M}{\∫}}{R}^{\mathrm{ACF}}({2x}_0-N+\mathrm{\τ},y)\mathrm{dy}$

$\begin{array}{c}{R}^{\mathrm{ACF}}({2x}_0-N+\mathrm{\τ},y)=R(\mathrm{\τ},y)\\ =\frac{1}{N}\underset{0}{\overset{N}{\∫}}P(x,y)P({2x}_0-N+x+\mathrm{\τ},y)\mathrm{dx}\end{array}$

步骤六：统计函数Q^ACF(2x₀-N+τ)的最大值，当Q^ACF(2x₀-N+τ)取最大值时，自变量τ的对应值记为τ'，则有投影旋转中心x0的值为

2.根据权利要求1所述的基于原始投影信息的CL扫描装置投影旋转中心标定方法，其特征在于：该标定方法不需要专门的标定模体，直接利用被扫描样品的投影信息即可实现精确标定。

3.根据权利要求1所述的基于原始投影信息的CL扫描装置投影旋转中心标定方法，其特征在于：被检测样品(3)在扫描平台(4)的带动下实现360度旋转，扫描平台(4)步进旋转，步进角为0.5°～1°。

4.根据权利要求1所述的基于原始投影信息的CL扫描装置投影旋转中心标定方法，其特征在于：被检测样品(3)的投影不能超出平板探测器(6)的成像平面(7)。

5.根据权利要求1所述的基于原始投影信息的CL扫描装置投影旋转中心标定方法，其特征在于：旋转轴线(5)与中心射线(10)的夹角小于90度。