1 引言

计算机断层成像(CT)技术能够无损地获得被检测物体的内部三维结构信息,重建物体三维图像,被广泛应用于无损检测^[1]、医疗诊断^[2]等领域。实际应用中,对集成电路、印刷电路板(PCB)等大尺寸扁平状物体进行成像时,由于受到检测空间和射线能量的限制,CT成像效果并不理想。为避免这种限制,计算机层析成像(CL)技术得到了发展,该技术对于扁平状物体的检测具有独特的优势^[3]。现有CL扫描方式主要有直线型、螺旋型和旋转型^[6]。旋转型CL技术在成像时,光源与探测器不动,旋转物体,采集被成像物体的投影,通过重建算法对采集到的投影进行重建得到物体的三维图像。实现重建算法需要已知系统的几何参数,如果几何参数存在偏差将造成重建图像出现几何伪影,因此,必须对实际CL系统进行几何参数标定。

现有CL几何参数标定方法仅能标定系统部分几何参数^[9]。通过分析CL和CT的系统结构可知,CL系统中旋转轴与探测器面不平行,存在一个倾斜角α;而在传统CT系统中旋转轴与探测器面近似平行。因此,理论上不受α角影响的CT系统几何参数标定方法都可以被CL系统所借鉴。现有CT系统几何参数标定方法中,一类通常忽略角度α,只能标定部分几何参数^[11];另一类能标定所有几何参数,但通常需要复杂的标定体模^[14]。2014年,Zhang等^[16]提出了一种迭代的基于双球定标体模的锥束CT圆轨迹几何参数标定方法。该方法所用定标体模简单,能够标定所有几何参数,但是在CL系统中该方法受转轴倾斜角影响较大,部分参数标定不准确,影响了重建图像质量,所以不能直接用于CL系统。

为了通过一次标定快速得到CL系统的所有几何参数,本文对Zhang等的方法(以下简称Zhang方法)进行了改进,提出了一种迭代的CL系统几何参数标定方法,能够一次标定系统的所有几何参数。首先分析Zhang方法不能直接用于CL系统的原因并确定敏感参数;其次,基于相对误差最小化建立包含敏感参数和受敏感参数影响的其他参数的非线性最小二乘目标函数,并用高斯-牛顿算法^[17]对目标函数进行迭代求解;最后,利用数值仿真实验对该方法的收敛速度、求解精度和稳健性进行分析,并利用实验标定的几何参数进行重建实验,验证算法的有效性。

2 锥束CL系统扫描结构及其几何参数

锥束CL系统扫描结构如图1所示。引入两个三维坐标系,XYZ是按照右手法则以探测器中心O为原点建立的世界坐标系,X'Y'Z'是按照右手法则以被扫描物体中心O'为原点建立的物体坐标系。SsX,sY,sZ为射线源焦点坐标,它在探测器上的投影点为PpX,0,pZ。Z'为被扫描物体的转轴,主射线SP与Z'轴垂直,垂足为QqX,qY,qZ。从几何结构上看,CL系统与CT系统的区别在于转轴Z'与Z轴不再平行,存在一个倾斜角α。定义Z'在探测器上的投影与Z轴的夹角为旋转轴的偏转角η,SP长度与SQ长度之比为几何放大比t。该系统需要标定的全部几何参数如下:射线源焦点坐标SsX,sY,sZ、射线源焦点投影坐标PpX,0,pZ、转轴旋转角η以及几何放大比t。

图 1. 锥束CL系统扫描结构

Fig. 1. Scanning configuration of the cone beam CL system

下载图片 查看所有图片

3 迭代的锥束CL系统几何全参数标定算法

3.1 确定CL系统敏感参数

通过对Zhang方法用于CL系统的有效性进行分析,确定CL系统几何参数中对Zhang方法敏感的参数。Zhang方法主要分为两个步骤:1)求得转轴投影方程和隐消线方程^[18],根据这两个方程可求得转轴偏转角和射线源焦点在探测器上的投影坐标;2)利用双球定标体模的投影,根据两小球球心的投影轨迹为椭圆的特性,拟合出上下两椭圆的中心坐标,即可得到过这两个中心点的直线方程,该直线方程近似为转轴投影方程,由转轴投影方程即可求得转轴偏转角。三维空间中一组平行线的投影相交于一点,该点称为隐消点,由一系列隐消点组成了隐消线^[18]。借鉴隐消线的思想,在两小球旋转轨迹中分别采样若干个点,这些点中每四个点构成一组平行线,则可以得到一系列平行线在探测器上形成的隐消线方程,隐消线方程与转轴投影方程的交点即为射线源焦点在探测器上的投影坐标。其次,通过系统几何的误差项建立包含射线源焦点坐标和系统几何放大比的目标函数,并利用高斯-牛顿迭代算法进行求解。

经分析可知将Zhang方法直接运用于CL系统几何参数标定中存在两个问题:1)射线源焦点投影纵坐标p_Z对该方法敏感,因为CL系统存在的倾斜角α使球心投影轨迹所形成的椭圆更加扁平,隐消线求解误差增大,进而造成p_Z求解误差增大;2)Zhang方法在迭代求解射线源交点坐标SsX,sY,sZ和几何放大比t的过程中将p_Z作为真值,p_Z与后续参数的求解直接相关,p_Z的求解不准确会进一步影响SsX,sY,sZ和t的准确求解。因此,p_Z为CL系统几何参数中对Zhang方法敏感的参数,Zhang方法直接运用于CL系统几何参数标定中会使部分参数求解误差较大。为验证上述观点,在两种噪声水平下分别重复500次实验,实验条件如下:探测器像素数为3200 pixel×2304 pixel,单个像素尺寸为0.127 mm,噪声水平(n_os)分别为0.05、0.1,噪声是均值为0、标准差为0.127n_os的高斯白噪声,CL系统倾斜角α为40°,小球直径为2 mm;其他参数设置见表1和表2中的真值。实验结果如表1和表2所示。

从实验结果可以看出,0.05噪声水平下p_Z的绝对误差从0.203增加到6.984,0.1噪声水平下p_Z的绝对误差从0.304增加到6.547,p_Z误差明显增大,证明了p_Z为敏感参数。并且在迭代求解SsX,sY,sZ的过程中将p_Z作为真值,造成SsX,sY,sZ的求解误差增大。

3.2 CL系统几何参数迭代标定方法

针对Zhang方法不能直接应用于CL系统几何参数标定的问题,提出了一种迭代的几何全参数标定方法。首先,利用隐消线原理和小球投影轨迹为椭圆的几何特性计算旋转轴偏转角η和射线源焦点投影横坐标p_X。其次,构建包含敏感参数p_Z(射线源焦点投影纵坐标)和受敏感参数影响的其他参数的目标函数,并利用高斯-牛顿算法进行求解。

3.2.1 计算η和p_X

在360°范围内均匀选取120个采样点(每隔3°取一个),并根据这些采样点模拟校正体模的120个投影。体模投影为上下两个椭圆,CL系统倾斜角为40°时体模投影如图2所示。根据体模投影可以拟合出两个椭圆的中心,由此得到旋转轴投影方程,即过上下椭圆中心的直线方程为

A1X+B1Y+C1Z=0。(1)

图 2. 体模投影

Fig. 2. Projections of phantom

下载图片 查看所有图片

根据3.1节所述隐消线理论^[18]可求得隐消线方程为

A2X+B2Y+C2Z=0。(2)

由(1)式转轴投影所在直线与Z轴的夹角即可求得转轴偏转角η。射线源焦点投影为转轴投影与隐消线的交点,所以由(1)式和(2)式可求得射线源焦点投影坐标PpX,0,pZ。

3.2.2 构建目标函数

由3.1节可知射线源焦点投影坐标为PpX,0,pZ,设射线源焦点坐标为SsX,sY,sZ,几何放大比为t,SP与Z'的交点坐标为QqX,qY,qZ。根据系统的几何关系,可得到两小球旋转轨迹上每个采样点的坐标分别为( X1n, Y1n, Z1n)和( X2n, Y2n, Z2n),且到其相应旋转中心的距离分别为 r1n和 r2n,n=1,2,…,120,则上下一组相对应的采样点之间的距离可表示为

dn=(X1n-X2n)2+(Y1n-X2n)2+(Z1n-X2n)212。(3)

因为小球实际运动轨迹与理想轨迹之间存在误差,所以引入e_1n、e_2n和e_3n误差,

e1n=∑m=1nr1m-r-1r-12e2n=∑m=1nr2m-r-2r-22e3n=∑m=1ndm-d2,(4)

式中 r-1= 1n∑m=1nr1m为上小球旋转轨迹上采样点到其中心的平均距离, r-2= 1n∑m=1nr2m为下小球旋转轨迹上采样点到其旋转中心的平均距离,d为标定体模两球心距离。

将射线源焦点投影纵坐标p_Z、射线源焦点坐标 sX,sY,sZ和几何放大比t一起作为待求变量,综合(4)式给出的3种误差,建立包含p_Z、 sX,sY,sZ、t的目标函数FpZ,sX,sY,sZ,t:

FpZ,sX,sY,sZ,t=∑n=1me1n+e2n+e3n。(5)

3.2.3 求解目标函数

系统没有几何误差的情况下满足FpZ,sX,sY,sZ,t=0,所以为了使得到的几何参数尽可能准确,要使FpZ,sX,sY,sZ,t的值尽可能小。FpZ,sX,sY,sZ,t具有非线性最小二乘解的形式,因此可以采用最优化理论与算法进行求解,本文采用的求解算法为高斯-牛顿算法^[17],具体过程如下:

1) 两球球心距离d为常数,优化过程中保持不变,将通过隐消线求解得到的p_Z和存在误差的机械系统的读数作为算法初值;

2) 计算第k次迭代的雅可比矩阵J,j_ij=∂e_i/∂X_j,即计算误差向量的第i个误差关于决策向量的第j个参数的偏导;

3) 计算增益矩阵:

ΔX=-JTkJk-1JTkek;(6)

4) X_k+1=X_k+ΔX,更新参数,直到满足终止条件,终止条件设为 ΔX<ε,其中ε=10^-6。

4 实验结果与分析

为了验证迭代的锥束CL系统几何参数标定方法,通过数值仿真对本文方法的求解精度和稳健性进行定量分析。实验条件为:探测器像素数为3200 pixel×2304 pixel,单个像素尺寸为0.127 mm,噪声水平n_os分别为0.05、0.1,噪声是均值为0、标准差为0.127n_os的高斯白噪声,α=40°,小球直径为2 mm,在两种噪声水平下分别重复实验500次,其他参数的设置见表3和表4中的真值。图3为目标函数F的收敛曲线,图4为F满足收敛终止条件后最终的相对误差,表3和表4为不同噪声水平下几何参数的标定结果。

图 3. F的收敛曲线

Fig. 3. Convergence curve of F

下载图片 查看所有图片

图 4. 最终相对误差

Fig. 4. Final relative error

下载图片 查看所有图片

图3表明本文方法的收敛速度较快,在迭代40多次时目标函数趋向收敛。图4表明最终的相对误差较小,在±0.05%范围内。表3和表4表明,与Zhang方法相比,在0.05噪声水平和0.1噪声水平下,本文方法对射线源焦点投影纵坐标p_Z的标定精度有明显提升,p_Z的绝对误差分别从6.984降低到0.013、从6.547降低到0.039。受益于p_Z标定精度的提升,本文方法对射线源焦点坐标 sX,sY,sZ和几何放大比t的标定精度也得到了提高。表3表明,在0.05噪声水平下,本文方法对 sX,sY,sZ以及t的标定精度较Zhang方法均有提高,s_X的绝对误差从0.164降低到0.055,s_Y的绝对误差从4.165降低到1.483,s_Z的绝对误差从5.771降低到1.739,t的绝对误差从0.032降到0.011。表4表明,在0.1噪声水平下,本文方法对 sX,sY,sZ以及t的标定精度较Zhang方法均有提高,s_X的绝对误差从0.529降低到0.089,s_Y的绝对误差从2.619降低到0.921,s_Z的绝对误差从7.069降低到1.071,t的绝对误差从0.041降到0.017。在两种噪声水平下,本文方法求解精度变化不大,因此本文方法具有较强的稳健性。

为进一步验证算法的有效性,利用噪声水平较高时(表4)的Zhang方法校正结果及本文方法校正结果分别对三维Shepp-Logan体模和PCB体模进行仿真重建实验。三维Shepp-Logan体模重建矩阵规模为256×256×256,PCB体模重建矩阵规模为360×360×360,重建结果分别如图5和图6所示。从图5和图6重建图像可以看出,Zhang方法校正结果不准确,造成重建图像中存在明显的几何伪影,图像质量差。而本文方法校正结果精确,重建图像中没有几何伪影,验证了本文方法的有效性,但CL特性使得扫描投影数据会有部分缺失,造成重建图像中出现数据缺失引起的伪影和层间混叠。

图 5. Shepp-Logan体模重建结果。(a) XY面切片图像;(b) XZ面切片图像;(c) YZ面切片图像

Fig. 5. Reconstruction results of Shepp-Logan phantom. (a) Images of XY slice; (b) images of XZ slice; (c) images of YZ slice

下载图片 查看所有图片

图 6. PCB体模第二层重建结果。(a)理想图像;(b) Zhang方法校正后图像;(c)本文方法校正后图像

Fig. 6. Reconstruction results of the PCB phantom. (a) Ideal image; (b) image corrected by Zhang's method; (c) image corrected by the proposed method

下载图片 查看所有图片

5 结论

基于锥束CT系统中经典的几何参数标定方法,提出了一种锥束CL系统几何全参数标定方法,能够通过一次标定得到系统的所有几何参数。仿真结果表明,本文方法收敛速度快、稳健性较强、计算精度高。相比Zhang方法,本文方法对敏感参数以及受敏感参数影响的几何参数的标定精度均有明显提升。利用本文方法校正系统几何参数后对Shepp-Logan体模和PCB体模进行重建,重建图像中没有几何伪影,验证了本文方法的有效性。但CL系统扫描结构特性使得采集到的投影数据有一部分缺失,造成重建图像中存在伪影和层间混叠,如何改善因数据缺失引起的伪影和层间混叠是下一步需要深入研究的工作。