CN103226113B - 锥束3d-ct扫描系统重建体素尺寸的自动标定方法 - Google Patents

Abstract

本发明公开了一种锥束3D-CT扫描系统重建体素尺寸自动标定的方法，该方法利用球形目标体在X轴方向不同行程距离下的DR投影的直径信息、采用最小二乘拟合法回归出不同成像位置下行程－体素尺寸关系V=b+kD，然后将V=b+kD内嵌在锥束3D-CT扫描系统的自动标定模块中。当锥束3D-CT扫描系统运行时，标定模块将自动计算出不同扫描位置对应的精确体素尺寸值，从而保证基于3D-CT图像的数据分析与特征测量的准确性。该标定方法避开了现有标定方法求取投影几何放大比难度大、精度低的难点，且无需对标准件进行CT扫描后再进行标定，具有操作简单、精度高、自动程度高的优点。

Description

锥束3D-CT扫描系统重建体素尺寸的自动标定方法

技术领域

本发明涉及一种适用于锥束3D－CT扫描系统重建体素尺寸的自动标定，可用于医学和工业领域射线数字成像（即DR投影）、二维、三维计算机断层扫描（2D-CT/3D-CT）成像过程中的相关测量。

背景技术

近些年来，随着计算机技术的飞速发展和面阵探测器的出现，锥束3D－CT（Cone-beam Three Dimensional Computed Tomography）系统日益成为NDT（Non-destructive Testing）领域内的研究热点。锥束3D－CT扫描系统的扫描原理如图1所示，锥束3D－CT扫描系统10的第一调节架11上安装有射线源1，第二调节架12上安装有面阵列探测器5，第三调节架13上安装有转台3，转台3的转轴31上放置有被检测物体4，且被检测物体4能够在转台3的带动下绕转轴31的轴线旋转。

锥束3D－CT扫描系统10在工作状态下，射线源1发出的锥束射线2对转台3上的被检测物体4进行透照，同时面阵列探测器5采集被检测物体4在不同旋转角度下的二维DR投影图像。该二维DR投影图像按照序列进行排列后被称为二维DR投影序列，利用这些二维DR投影序列进行三维图像重建，将得到被检测物体4内部结构与材质分布的信息。

3D－CT重建所得到的三维图像反映被检测物体4内部结构与材质分布信息，被检测物体4可认为是由大量的具有一定尺寸的小立方体排列堆积而成，这些小立方体被称为体素，立方体的实际物理尺寸被称为体素尺寸（Voxel Size）。体素尺寸是描述锥束3D－CT扫描系统的分辨率的重要指标，体素尺寸越小，三维图像的空间分辨率越高。然而，三维重建是在脱离了锥束3D－CT扫描系统的实际成像条件下进行的图像重建，重建所得三维图像的尺寸单位为像素，因此不可能从重建图像中直接获取到体素尺寸。而体素尺寸的准确测量，是锥束3D－CT扫描系统的分辨率标定的重要环节。同时，体素尺寸的精确测量，也决定着基于3D－CT三维重建结果的图像分析与特征测量的精度。

因此，在锥束3D－CT扫描和三维重建后，如何将像素单位的体素尺寸换算为实际的物理体素尺寸，是锥束3D－CT扫描系统标定中的重要步骤。在图1中对于锥束3D－CT扫描系统而言，当被检测物体4在射线源1与面阵列探测器5之间的不同位置（即行程距离）成像时，其对应的体素尺寸也是不同的。因此需要动态地自动标定随被检测物体4的成像位置改变的体素尺寸。

发明内容

本发明的目的是提出一种锥束3D－CT扫描系统重建体素尺寸自动标定的方法，该方法利用面阵列探测器5采集到的球形目标体（即被检测物体4）在X轴方向上的不同成像位置（即行程距离）下的二维DR投影；然后基于二维DR投影的直径信息采用最小二乘拟合法，回归出不同成像位置下行程距离与体素尺寸之间的行程－体素尺寸关系V=b+kD；最后将该行程－体素尺寸关系V=b+kD内嵌在锥束3D－CT扫描系统的自动标定模块中；锥束3D－CT扫描系统10在工作状态下，通过带有行程－体素尺寸关系V=b+kD的自动标定模块，以当前的成像位置为自变量，自动计算出当前成像位置下的体素尺寸的精确值。该标定方法避开了现有标定方法求取投影几何放大比难度大、精度低的难点，且无需对标准件进行CT扫描后再进行标定，具有操作简单、精度高、自动程度高的优点。

本发明的一种锥束3D－CT扫描系统重建体素尺寸自动标定的方法，所述锥束3D－CT扫描系统至少包括有射线源（1）、转台（3）、被检测物体（4）、面阵列探测器（5），其标定方法包括有下列实施步骤：

步骤一：初始行程距离下的体素尺寸标定

步骤101：将一球形目标体安装在转台（3）的转轴（31）上，将转台（3）移动至X轴方向的起始成像位置P₀，记录下行程距离D₁、计算出首次体素尺寸V₁；

步骤102：启动锥束3D－CT扫描系统，获取球形目标体的二维DR数字的投影图像；

步骤103：利用图像处理和轮廓追踪方法获取该起始成像位置P₀下的图像轮廓点；

步骤104：根据记录的在起始成像位置P₀下的所有轮廓点坐标值，利用非线性最小二乘拟合法求取出球形目标体的投影直径的像素单位长度W₁′；

步骤二：最大行程距离下的体素尺寸标定

步骤201：将转台（3）移动至X轴方向的行程最大位置点P_max，记录下最长行程距离D_max、计算出最后次体素尺寸V_max；

步骤202：启动锥束3D－CT扫描系统，获取在行程最大位置点P_max下的球形目标体的二维DR数字的投影图像；

步骤203：利用图像处理和轮廓追踪方法获取该行程最大位置点P_max下的图像轮廓点；

步骤204：根据记录的在行程最大位置点P_max下的所有轮廓点坐标值，利用非线性最小二乘拟合法求取出球形目标体的投影直径的像素单位长度W_max′；

步骤三：中间行程距离下的体素尺寸标定

步骤301：将转台（3）移动至X轴方向上的介于P₀和P_max之间的中间成像位置点P_i，记录下中间行程距离D_i、计算出中间次体素尺寸V_i；

步骤302：启动锥束3D－CT扫描系统，获取在介于P₀和P_max之间的中间成像位置点P_i下的球形目标体的二维DR数字的投影图像；

步骤303：利用图像处理和轮廓追踪方法获取该中间成像位置点P_i下的图像轮廓点；

步骤304：根据记录的在中间成像位置点P_i下的所有轮廓点坐标值，利用非线性最小二乘拟合法求取出球形目标体的投影直径的像素单位长度W_i′；

步骤四：任意行程下的体素尺寸获取

将转台（3）在X轴方向上不同的成像位置计算出的体素尺寸值和对应的行程距离值组成一方程组为：

$\left[\begin{array}{cc}{D}_1& 1\\ D_2& 1\\ ...& 1\\ D_i& 1\\ ...& 1\\ D_{\max }& 1\end{array}\right]\×\left[\begin{array}{c}k\\ b\end{array}\right]=\left[\begin{array}{c}{V}_1\\ V_2\\ ...\\ V_i\\ ...\\ V_{\max }\end{array}\right];$

从体素尺寸值和行程距离值组成一方程组中即可得到零次项系数b、一次项系数k的值，从而得到不同成像位置下行程－体素尺寸关系V＝b+kD。

本发明体素尺寸自动标定方法的优点在于：

①本发明应用行程－体素尺寸关系V=b+kD标定方法，避开了现有标定方法求取投影几何放大比难度大、精度低的难点，且无需对标准件进行CT扫描后再进行标定，具有操作简单、精度高、自动程度高的优点。以扫描系统控制编码器或光栅尺精确输出的X轴方向的行程距离为自变量，回归出对应成像位置下的体素尺寸精确值。

②本发明方法实施简单，不需要对专门的标准件进行CT扫描后实施标定，仅利用球形目标体的若干个DR投影即可。

③利用二维DR投影的直径信息约束下的线性最小二乘拟合得出行程－体素尺寸关系，使得通过行程－体素尺寸关系计算出的体素尺寸精度高，重复性好。

④将行程－体素尺寸关系内嵌在锥束3D－CT扫描系统的自动标定模块中，通过软件编程自动实现，降低了对锥束3D－CT扫描系统的改进成本，同时可作为锥束3D－CT扫描系统的辅助校验模块，自动监测转台的成像位置变化。

附图说明

图1是锥束3D－CT扫描系统的工作原理示意图。

图1A是面阵列探测器的探测单元布局示意图。

图1B是面阵探测器的单个探测单元物理尺寸示意图。

图2A是球形目标体的二维DR投影线条示意图。

图2B是球形目标体的二维DR投影图。

图2C是球形目标体的二维DR轮廓线条示意图。。

图2D是球形目标体的二维DR轮廓图。

图3是经本发明标定方法标定后的锥束3D－CT扫描系统的行程－体素尺寸关系图。

10.锥束3D－CT扫描系统	11.第一调节架
		12.第二调节架	13.第三调节架
1.射线源	2.锥束射线
		3.转台	31.转轴
4.被检测物体	5.面阵列探测器
		51.探测单元

具体实施方式

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

参见图1、图1A所示，在锥束3D－CT扫描系统中，当选取出面阵列探测器5的型号后，则确定出了面阵列探测器5的探测单元51的物理尺寸，记为U，U为一个常数。如图1A所示，探测器5可认为是由许多个尺寸相同的正方形探测单元排列而成，探测单元51为一个正方形，其边长即为探测单元的物理尺寸，即等同于常数U，图1B为探测单元的物理尺寸示意图。探测单元51的长是以锥束3D－CT扫描系统坐标系O-XYZ中的Z轴方向布局，探测单元51的宽是以锥束3D－CT扫描系统坐标系O-XYZ中的Y轴方向布局。在锥束3D－CT扫描系统中，射线源1出射的锥束射线2的射线源焦点记为O₁，面阵列探测器5的中心点记为O₅，射线源焦点O₁到探测器中心点O₅之间的距离记为D_FOD，D_FOD简称为系统焦距。

在锥束3D－CT扫描系统中，转台3的初始成像位置记为P₀，初始成像位置P₀下的行程距离记为D₀，所述D₀简称为初始行程距离。

在锥束3D－CT扫描系统中，转台3的最大成像位置记为P_MAX，最大成像位置P_MAX下的行程距离记为D_MAX，所述D_MAX简称为最大行程距离。

在锥束3D－CT扫描系统中，除了初始成像位置P₀和最大成像位置P_MAX以外，转台3需要放置在射线源1与面阵列探测器5之间的任意一位置P_i，该位置也称为中间任意成像位置P_i，i表示初始成像位置P₀与最大成像位置P_MAX之间的任意一成像位置点；在不同中间成像位置P_i下，转台3沿X轴方向的行程距离记为D_i，则被检测物体4的成像距离记为L_i。

所述成像距离L_i是指转台3的转轴31的中心线与射线源焦点O₁之间的距离。

所述行程距离D_i是指转台3在电机驱动下，相对于转台的初始行程D₀移动后的距离，该距离值由转台3的电机运动控制编码器精确输出，或者由X轴方向的光栅尺精确输出。

在锥束3D－CT扫描系统中，当面阵列探测器5型号确定后，其探测单元51的物理尺寸U就随之确定，因此任意中间成像位置P_i下的体素尺寸V_i的标定也就是成像几何放大比M_i的标定，即：

$M_i=\frac{D_{\mathrm{FOD}}}{L_i}---\left(1\right)$

M_i表示任意一中间成像位置P_i下的成像几何放大比。

D_FOD表示系统焦距，一般取值为1500mm～700mm，D_FOD值在锥束3D－CT扫描系统装配完成后即为一个常数，一般锥束3D－CT扫描系统出厂时会给出其精确值。

在本发明中，所述成像距离L_i在锥束3D－CT扫描系统中无法直接测量得到，因此体素尺寸V_i的标定的难点就是如何快速精确求得式(1)中L_i的取值。由图1可知：

L_i＝D₀+D_i                              (2)

L_i表示成像距离，一般地0＜L_i＜D_FOD。

D₀表示射线源焦点O₁与转台初始位置P₀之间的行程距离，简称为初始行程距离；所述转台初始位置是指在锥束3D－CT扫描系统进入工作状态时，转台3在X轴方向上运动的起点位置，该位置在3D－CT扫描系统设计时确定，当3D－CT扫描系统出厂后该位置即被固定，不可改变。初始行程距离即为转台3位于初始位置时其转轴31的中心线与射线源焦点O₁之间的行程距离。X轴方向为射线束的照射方向。

D_i表示任意一成像位置P_i处的行程距离。

在本发明中，转台3在X轴方向上的行程的起点位置记为P₀，也称为初始成像位置P₀。在锥束3D－CT扫描系统工作中，转台3沿X轴方向运动的起始位置P₀需要根据锥束3D－CT扫描系统的设计指标预先设定，因此当锥束3D－CT扫描系统装配完成后，初始行程距离D₀变成常数量。但该常数的量值却无法精确测量得到。

因此，成像位置P_i处所对应的体素尺寸V_i的关系为：

$V_i=\frac{U}{M_i}---\left(3\right)$

联立式(1)、式(2)和式(3)得到：

$V_i=U\×\frac{D_0}{D_{\mathrm{FOD}}}+\frac{U}{D_{\mathrm{FOD}}}\×D_i---\left(4\right)$

设 显然，零次项系数b、一次项系数k为常数，那么式(4)变化为：

V_i＝b+kD_i                              (5)

由此得出体素尺寸V_i与行程距离D_i之间的通用关系式：

V＝b+kD                                (6)

V表示体素尺寸变量；V_i表示体素尺寸变量V的某一个取值。

D表示行程距离变量（即转台3的行程），D_i表示行程距离变量D的某一个取值。

公式(6)为一次方程，如果得到b、k的值，就可以利用该式得出转台3沿X轴方向移动到任意成像位置时的体素尺寸值。根据线性拟合方法可知，如果知道若干个体素尺寸变量V的离散值为[V₁,V₂,......V_i]，以及相对应的行程距离变量D的离散值为[D₁,D₂,......D_i]，即可利用最小二乘拟合方法得到b、k的值，得到方程(6)的表达式，从而解决锥束3D－CT扫描系统体素尺寸的自动标定问题。

在本发明中，体素尺寸V_i与行程距离D_i之间的通用关系也简称为行程－体素尺寸关系V=b+kD。

在本发明中，为了事先对锥束3D－CT扫描系统进行体素尺寸的自动标定，采用球形目标体代替被检测物体4，球形目标体由金属材质制成。

本发明是一种适用于锥束3D－CT扫描系统重建体素尺寸的自动标定，该标定方法包括有下列实施步骤：

步骤一：初始行程距离下的体素尺寸标定

步骤101：如图1所示，将一球形目标体安装在转台3的转轴31上，所述球形目标体的直径可精确测量得到，其值为常数，记为W。将转台3移动至X轴方向的初始成像位置P₀；该初始成像位置P₀对应的转台3的行程距离为D₁，显然D₁＝0mm；该初始成像位置P₀下的体素尺寸记为V₁（简称为首次体素尺寸V₁）。

步骤102：启动锥束3D－CT扫描系统，获取球形目标体的二维DR数字透视投影图像，该投影图像为圆形投影，如图2A、图2B。图2A为用线条表示的二维DR投影图像，图2B为实际成像得到的二维DR投影图像，由于图2B有背底色，所以用图2A来进行参考说明。

步骤103：利用图像处理和轮廓追踪方法（锥束3D－CT扫描系统中已有的图像轮廓获取方法）获取该初始成像位置P₀下的圆形透视图像的轮廓点（简称为图像轮廓点），如图2C和图2D所示，由于图2D有背底色，所以用图2C来进行参考说明。每个图像轮廓点的坐标值记为(x_j,y_j)，j表示图像轮廓点的个数，j＝1，2，……，n，n为轮廓点总数目，x_j为投影图像坐标系o-xy下的在x轴上的像素点位置坐标值，y_j为投影图像坐标系o-xy下的在y轴上的像素点位置坐标值。例如，第一个图像轮廓点（j＝1）的坐标值记为(x₁,y₁)，第二个图像轮廓点（j＝2）的坐标值记为(x₂,y₂)，……，第j个图像轮廓点的坐标值记为(x_j,y_j)，为了方便说明，(x_j,y_j)也表示任意一个图像轮廓点的坐标。

步骤104：根据记录的在初始成像位置P₀下的所有轮廓点坐标值，利用非线性最小二乘拟合法求取出球形目标体的投影直径的像素单位长度W₁′。

在本发明中，非线性最小二乘拟合法回归出球形目标体的投影直径方法如下：

设圆方程为(x_j-x₀)²+(y_j-y₀)²＝r²，其中：(x₀,y₀)为圆形透视图像的圆心坐标；r为半径；(x_j,y_j)为圆形透视图像的圆轮廓点坐标，也就是通过对球形目标体的投影进行图像处理与轮廓追踪后得到的轮廓点坐标。建立投影图像轮廓误差函数E，其表达式为：

$E=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{[{(x_j-x_0)}^2+{(y_j-y_0)}^2-r^2]}^2---\left(7\right)$

n表示求和元素，即轮廓点总数目；j表示求和指标；

令： $m=x_0^2+y_0^2-r^2$ ，则式(7)变换为：

$E=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{(x_j^2-{2x}_j{x}_0+y_j^2-{2y}_j{y}_0+m)}^2---\left(8\right)$

基于误差函数E对x₀求偏导，记为，其中表示求偏导运算。

基于误差函数E对y₀求偏导，记为，其中表示求偏导运算。

基于误差函数E对m求偏导，记为，其中表示求偏导运算。

令： $\frac{\∂E}{{\∂x}_0}=0,$ $\frac{\∂E}{{\∂y}_0}=0,$ $\frac{\∂E}{\∂m}=0,$ 即：

$\left\{\begin{array}{c}\frac{\∂E}{{\∂x}_0}=2\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(x_j^2-{2x}_j{x}_0+y_j^2-{2y}_j{y}_0+m)\×\left({-2x}_j\right)=0\\ \frac{\∂E}{{\∂y}_0}=2\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(x_j^2-{2x}_j{x}_0+y_j^2-{2y}_j{y}_0+m)\×\left({-2y}_j\right)=0\\ \frac{\∂E}{\∂m}=2\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(x_j^2-{2x}_j{x}_0+y_j^2-{2y}_j{y}_0+m)=0\end{array}\right.----\left(9\right)$

解析式(9)，即可得到球形目标体的投影的半径值r，从而得到其直径值W₁′＝2×r，该值的单位为像素。那么初始行程距离下的体素尺寸值记为V₁，可以用如下关系式求得：

$V_1=\frac{W}{{W_1}^{\′}}---\left(10\right)$

W₁′表示初始行程距离下的球形目标体的投影直径的像素单位长度。

W表示球形目标体的实际直径值。

步骤二：最大行程距离下的体素尺寸标定

步骤201：如图1所示，将转台3移动至X轴方向的行程最大位置点，该行程最大位置点的成像位置记为P_max（简称为最大成像位置P_max），该最大成像位置P_max对应的转台3的行程距离记为D_max（简称为最长行程距离D_max）；该最大成像位置P_max下的体素尺寸记为V_max（简称为最后次体素尺寸V_max）。

步骤202：启动锥束3D－CT扫描系统，获取在最大成像位置P_max下的球形目标体的二维DR数字透视投影图像，该投影图像为圆形投影，如图2A、图2B。图2A为用线条表示的二维DR投影图像，图2B为实际测量得到的二维DR投影图像，由于图2B有背底色，所以用图2A来进行参考说明。

步骤203：利用图像处理和轮廓追踪方法（锥束3D－CT扫描系统中已有的图像轮廓获取方法）获取该最大成像位置P_max下的圆形透视图像的轮廓点（简称为图像轮廓点），如图2C和图2D所示，由于图2D有背底色，所以用图2C来进行参考说明。每个图像轮廓点的坐标值记为(x_j,y_j)，j表示图像轮廓点的个数，j＝1，2，……，n，n为轮廓点总数目，x_j为投影图像坐标系o-xy下的在x轴上的像素点位置坐标值，y_j为投影图像坐标系o-xy下的在y轴上的像素点位置坐标值。例如，第一个图像轮廓点（j＝1）的坐标值记为(x₁,y₁)，第二个图像轮廓点（j＝2）的坐标值记为(x₂,y₂)，……，第j个图像轮廓点的坐标值记为(x_j,y_j)，为了方便说明，(x_j,y_j)也表示任意一个图像轮廓点的坐标。

步骤204：根据记录的在最大成像位置P_max下的所有轮廓点坐标值，利用非线性最小二乘拟合法求取出球形目标体的投影直径的像素单位长度W_max′。

在本发明中，步骤202、步骤203、步骤204与步骤102、步骤103和步骤104是相同的获取方式，即非线性最小二乘拟合法回归出球形目标体的投影直径方法是相同的。因此，最长行程距离D_max下的最后次体素尺寸V_max，可以用如下关系式求得：

$V_{\max }=\frac{W}{{W_{\max }}^{\′}}---\left(11\right)$

W_max′表示最长行程距离D_max下的球形目标体的投影直径的像素单位长度。

步骤三：中间行程距离下的体素尺寸标定

步骤301：如图1所示，将转台3移动至X轴方向上的介于P₀和P_max之间的任一中间位置P_i（简称为中间成像位置P_i），该位置P_i对应的转台3的行程距离记为D_i（简称为中间行程距离D_i）。该中间成像位置P_i下的体素尺寸记为V_i（简称为中间次体素尺寸V_i）。

步骤302：启动锥束3D－CT扫描系统，获取在中间成像位置P_i下的球形目标体的二维DR数字透视投影图像，该投影图像为圆形投影，如图2A、图2B。图2A为用线条表示的二维DR投影图像，图2B为实际测量得到的二维DR投影图像，由于图2B有背底色，所以用图2A来进行参考说明。

步骤303：利用图像处理和轮廓追踪方法（锥束3D－CT扫描系统中已有的图像轮廓获取方法）获取该中间成像位置P_i下的圆形透视图像的轮廓点（简称为图像轮廓点），如图2C和图2D所示，由于图2D有背底色，所以用图2C来进行参考说明。每个图像轮廓点的坐标值记为(x_j,y_j)，j表示图像轮廓点的个数，j＝1，2，……，n，n为轮廓点总数目，x_j为投影图像坐标系o-xy下的在x轴上的像素点位置坐标值，y_j为投影图像坐标系o-xy下的在y轴上的像素点位置坐标值。例如，第一个图像轮廓点（j＝1）的坐标值记为(x₁,y₁)，第二个图像轮廓点（j＝2）的坐标值记为(x₂,y₂)，……，第j个图像轮廓点的坐标值记为(x_j,y_j)，为了方便说明，(x_j,y_j)也表示任意一个图像轮廓点的坐标。

步骤304：根据记录的在中间成像位置P_i下的所有轮廓点坐标值，利用非线性最小二乘拟合法求取出球形目标体的投影直径的像素单位长度W_i′。

在本发明中，步骤302、步骤303、步骤304与步骤102、步骤103和步骤104是相同的获取方式，即非线性最小二乘拟合法回归出球形目标体的投影直径方法是相同的。因此，中间行程距离D_i下的中间次体素尺寸V_i，可以用如下关系式求得：

$V_i=\frac{W}{{W_i}^{\′}}---\left(12\right)$

W_i′表示中间行程距离D_i下的球形目标体的投影直径的像素单位长度。

重复步骤三的处理，获取介于P₀和P_max之间的多个不同位置下的球形目标体的二维DR数字透视投影图像，从而得到这些成像位置下的体素尺寸值。

步骤四：任意行程下的体素尺寸获取

将转台3在X轴方向上不同的成像位置计算出的体素尺寸值和行程距离值组成一方程组为：

$\left[\begin{array}{cc}{D}_1& 1\\ D_2& 1\\ ...& 1\\ D_i& 1\\ ...& 1\\ D_{\max }& 1\end{array}\right]\×\left[\begin{array}{c}k\\ b\end{array}\right]=\left[\begin{array}{c}{V}_1\\ V_2\\ ...\\ V_i\\ ...\\ V_{\max }\end{array}\right]---\left(13\right)$

求解式（13），即可得到零次项系数b、一次项系数k的值，从而得到式(6)的表达式V＝b+kD。在得到式(6)的表达式后，只要输入转台3在X轴方向的任意行程距离值，即可得到该位置下的体素尺寸值，从而实现锥束3D－CT扫描系统体素尺寸的自动标定。

实施例

利用本发明的标定方法对某锥束3D－CT扫描系统的体素尺寸进行自动标定。验证采用的球形目标体为直径为19mm的钢珠。探测单元物理尺寸为0.2mm。图2A、图2B、图2C和图2D为该钢珠的在某个成像位置的DR图像以及图像处理后提取出的轮廓图，其对应D_i＝330mm。表1为验证中转台不同行程距离下的测量结果。绘制表1的D_i-V_i关系图，如图3所示，可以看出，转台行程距离与体素尺寸的关系为线性一次函数关系，证明上述推论的正确性。

将利用表1数据拟合出的线性方程输入锥束3D－CT扫描系统的自动标定模块。当系统正常运行时，标定模块将根据当前转台3在X轴方向的行程距离，自动计算出该扫描位置对应的精确体素尺寸值，从而保证了基于3D－CT图像的数据分析与特征测量的准确性。

表1D_i-V_i测试结果表

本发明提出的一种锥束3D－CT扫描系统重建体素尺寸自动标定的方法，该方法利用面阵列探测器采集到的球形目标体（即被检测物体）在X轴方向上的不同成像位置（即行程距离）下的二维DR投影；然后依据二维DR投影的直径信息采用最小二乘拟合法，回归出不同成像位置下行程距离与体素尺寸之间的行程－体素尺寸关系V=b+kD；最后将该行程－体素尺寸关系V=b+kD内嵌在锥束3D－CT扫描系统的自动标定模块中；锥束3D－CT扫描系统10在工作状态下，通过带有行程－体素尺寸关系V=b+kD的自动标定模块，以当前的成像位置为自变量，自动计算出当前成像位置下的体素尺寸的精确值。该标定方法避开了现有标定方法求取投影几何放大比难度大、精度低的难点，且无需对标准件进行CT扫描后再进行标定，具有操作简单、精度高、自动程度高的优点。

Claims (1)

1.一种锥束3D-CT扫描系统重建体素尺寸自动标定的方法，所述锥束3D-CT扫描系统至少包括有射线源(1)、转台(3)、被检测物体(4)、面阵列探测器(5)；

被检测物体(4)为球形目标体，且球形目标体由金属材质制成；

球形目标体安装在转台(3)的转轴(31)上；

面阵列探测器(5)为平板探测器，探测器成像面由离散的正方形探测单元整齐排列而成；

在利用锥束3D-CT扫描系统重建体素尺寸的自动标定过程中，包括有：

步骤一：初始行程距离下的体素尺寸标定

步骤101：将转台(3)移动至X轴方向的初始成像位置P₀，记录下初始行程距离D₁、记录下首次体素尺寸V₁；

步骤102：启动锥束3D-CT扫描系统，获取球形目标体的二维DR数字的投影图像；

步骤103：利用图像处理和轮廓追踪方法获取该初始成像位置P₀下的图像轮廓点；

步骤104：根据记录的在初始成像位置P₀下的所有轮廓点坐标值，利用非线性最小二乘拟合法求取出球形目标体的投影直径的像素单位长度W₁′；

步骤二：最大行程距离下的体素尺寸标定

步骤201：将转台(3)移动至X轴方向的最大成像位置P_max，记录下最长行程距离D_max、记录下最后次体素尺寸V_max；

步骤202：启动锥束3D-CT扫描系统，获取在最大成像位置P_max下的球形目标体的二维DR数字的投影图像；

步骤203：利用图像处理和轮廓追踪方法获取该最大成像位置P_max下的图像轮廓点；

步骤204：根据记录的在最大成像位置P_max下的所有轮廓点坐标值，利用非线性最小二乘拟合法求取出球形目标体的投影直径的像素单位长度W_max′；

步骤三：中间行程距离下的体素尺寸标定

步骤301：将转台(3)移动至X轴方向上的介于P₀和P_max之间的中间成像位置P_i，记录下中间行程距离D_i、记录下中间次体素尺寸V_i，i表示初始成像位置P₀与最大成像位置P_MAX之间的任意一成像位置点；

步骤302：启动锥束3D-CT扫描系统，获取在介于P₀和P_max之间的中间成像位置P_i下的球形目标体的二维DR数字的投影图像；

步骤303：利用图像处理和轮廓追踪方法获取该中间成像位置P_i下的图像轮廓点；

步骤304：根据记录的在中间成像位置P_i下的所有轮廓点坐标值，利用非线性最小二乘拟合法求取出球形目标体的投影直径的像素单位长度W_i′；

其特征在于：在利用锥束3D-CT扫描系统重建体素尺寸的自动标定过程中，还包括有：

步骤四：任意行程下的体素尺寸获取

将转台(3)在X轴方向上不同的成像位置计算出的体素尺寸值和对应的行程距离值组成一方程组为：

$\left[\begin{array}{cc}{D}_1& 1\\ D_2& 1\\ \mathrm{...}& 1\\ D_i& 1\\ \mathrm{...}& 1\\ D_{\max }& 1\end{array}\right]\×\left[\begin{array}{c}k\\ b\end{array}\right]=\left[\begin{array}{c}{V}_1\\ V_2\\ \mathrm{...}\\ V_i\\ \mathrm{...}\\ V_{\max }\end{array}\right];$

从体素尺寸值和行程距离值组成一方程组中即可得到零次项系数b、一次项系数k的值，从而得到不同成像位置下行程－体素尺寸关系V＝b+kD；

所述非线性最小二乘拟合法回归出球形目标体的投影直径方法如下：

设圆方程为(x_j-x₀)²+(y_j-y₀)²＝r²，其中：(x₀,y₀)为圆形透视图像的圆心坐标；r为半径；(x_j,y_j)为圆形透视图像的圆轮廓点坐标，也就是通过对球形目标体的投影进行图像处理与轮廓追踪后得到的轮廓点坐标；建立投影图像轮廓误差函数E，其表达式为n表示求和元素，即轮廓点总数目；j表示求和指标；

令： $m=x_0^2+y_0^2-r^2,$ 则投影图像轮廓误差函数 $E=\underset{j=1}{\overset{n}{\Σ}}{\[{\left(x_i-x_0\right)}^2+{\left(y_j-y_0\right)}^2-r^2\]}^2$ 变换为 $E=\underset{j=1}{\overset{n}{\Σ}}{(x_j^2-2x_j{x}_0+y_j^2-2y_j{y}_0+m)}^2;$

基于误差函数E对x₀求偏导，记为其中表示求偏导运算；

基于误差函数E对y₀求偏导，记为其中表示求偏导运算；

基于误差函数E对m求偏导，记为其中表示求偏导运算；

令： $\frac{\∂E}{\∂x_0}=0,\frac{\∂E}{\∂y_0}=0,\frac{\∂E}{\∂m}=0,$ 即：

$\left\{\begin{array}{c}\frac{\∂E}{\∂x_0}=2\underset{j=1}{\overset{n}{\Σ}}\left(x_j^2-2x_j{x}_0+y_j^2-2y_j{y}_0+m\right)\×\left(-2x_j\right)=0\\ \frac{\∂E}{\∂y_0}=2\underset{j=1}{\overset{n}{\Σ}}\left(x_j^2-2x_j{x}_0+y_j^2-2y_j{y}_0+m\right)\×\left(-2y_j\right)=0\\ \frac{\∂E}{\∂m}=2\underset{j=1}{\overset{n}{\Σ}}\left(x_j^2-2x_j{x}_0+y_j^2-2y_j{y}_0+m\right)=0\end{array}\right.$

解析上式得到球形目标体的投影的半径值r，从而得到其直径值W₁′＝2×r，该值的单位为像素。