基于相机模型的锥束CT重建误差校正
王梦蛟, 丁辉, 王广志
清华大学 生物医学工程系,北京 100084
通讯作者:王广志, 教授, E-mail:wgz-dea@mail.tsinghua.edu.cn

作者简介: 王梦蛟(1988—), 男(汉), 河南, 博士研究生。

摘要

该文对存在几何误差的锥束计算机层析成像(CT)的图像重建校正方法进行了研究。该文将锥束CT系统看做针孔相机,因此相机的成像模型可以用来修正三维体素与二维投影之间的映射关系。将该映射关系带入到重建算法中,就可以对存在误差的锥束CT进行有效的重建,提高重建图像的质量。计算机视觉中的相机标定技术被用来构建成像模型。该技术的可靠性保证了本文方法的精度和鲁棒性。仿真和实际系统实验表明: 本方法可以在图像重建过程中对锥束CT中的几何误差进行有效的校正。

关键词: 锥束计算机层析成像(CT)重建; 几何误差校正; 针孔相机成像模型; 相机标定技术
中图分类号:R814.42 文献标志码:A 文章编号:1000-0054(2015)01-0122-06
Misaligned cone beam computed tomography reconstruction based on a camera model
Mengjiao WANG, Hui DING, Guangzhi WANG
Department of Biomedical Engineering, Tsinghua University, Beijing 100084, China
Abstract

A correction was developed for a misaligned cone beam computed tomography (CT) system. The cone beam CT system was modeled as a pin-hole camera. A camera imaging model was used to describe the mapping information between the 3D voxels and the 2D projection. This information was integrated into the cone beam CT reconstruction to correct the misaligned cone beam CT image. A computing camera calibration technique was used to build the imaging model. This technique is reliable and accurate. Simulations and real system tests show that this method effectively corrects geometric misalignments during cone beam CT reconstruction.

Keyword: cone beam computed tomography(CT) reconstruction; geometric misalignments correction; pin-hole camera imaging model; camera calibration technique

近几年来,锥束计算机层析成像( CT)系统被越来越多地应用在医学和工业领域中。然而由于加工与组装中精度的限制,往往会使锥束 CT系统中产生几何误差,从而导致重建图像质量的下降。

之前针对锥束 CT重建误差校正的研究主要集中于对三位体素与二维投影之间映射关系的修正。 Karolczak Sun[1,2]提出首先求解出锥束 CT中的各个几何误差和参数,然后通过几何关系建立起三维体素与二维投影的映射关系。此方法需要实现对各个参数的标定,计算繁琐。并且在标定过程中通常假设一个或几个几何误差为零的理想条件,从而简化运算。这就导致构建的几何模型的不准确。 Wiesent Galigekere Navab Anne Li[3,4,5,6,7,8]提出利用投影矩阵来描述该映射关系,但其所获得的投影矩阵对锥束 CT的投影角度和标定仿体的摆放具有特异性,当二者中任意一个因素发生变化时,该投影矩阵将无法准确描述变化后的映射关系。

本文对基于投影矩阵的方法进行改进,以避免投影矩阵的求解受到投影角度变化和标定仿体摆放的影响。利用锥束 CT中固定的旋转轴线作为参考,从而定义三维体素在成像系统中的位置。在此基础上,可以得到绕转轴旋转任意角度后的体素位置,从而构建任意角度下的投影矩阵[9,10]。本文使用针孔相机成像模型实现对投影矩阵的求解。

1 标准的锥束CT图像重建方法

标准的锥束CT重建方法由Feldkamp、Davis和Kress提出的FDK算法[11],是几何条件理想的锥束CT重建中经常使用的方法。标准的FDK算法可以描述为

f FDK(x,y,z)= 02πU2(P)·{[W(P)·p F(φ,a p,b p)]·g F(a)} dφ. (1)

在上式的求解中,第一步要求解出体素 P( x, y, z)在旋转角度为 φ时所对应的二维图像中的投影 pF( φ, ap, bp),这里的 ap bp可以使用理想几何条件下的相似三角形的关系进行求解,即

ap =DSD -xsinφ+ycosφDSR+xcosφ+ysinφ, (2)

bp =z DSDDSR+xcosφ+ysinφ.(3)

其中: DSR为源到旋转中心的距离, DSD为源到探测器的距离。

而后对投影 pF( φ, ap, bp)乘上加权因子 W( P),以弥补穿过三维数据体的 X射线的差异,即

W(P)=DSDDSpF.(4)

其中 DSpF为源到二维投影位置的距离。

通过将加权后的投影数据与函数 gF( a)进行卷积,可以实现滤波校正。最后通过将滤波后的投影数据乘上反投因子 U( P),进行反投影,就可以重建出三维CT数据。其中 U( P)可以表示为

U(P)=DSD(DSP·x^).(5)

其中: DSP· x^为源到体素的距离在投影主轴上的投影。这里所述的投影主轴为源到探测器的垂线。

2 存在几何误差的锥束CT重建方法

锥束 CT重建方法中重要的步骤之一为修正三维体素与二维投影之间的映射关系。在理想几何条件下,可以按照式(2)和式(3)使用相似三角形关系进行求解。但当锥束 CT中存在几何误差时,该几何关系不复存在。因此,本文使用针孔相机成像的模型来类比锥束 CT的成像过程,从而得到存在几何误差的锥束 CT中,三维体素与二维投影之间的映射关系。

2.1 锥束CT成像过程描述

在锥束CT系统中,无论是源-探测器绕被测物体旋转的模式,还是被测物体旋转的模式,其采集模式都是等价的。为了表达方便,采用被测物体旋转、源-探测器保持静止的模式进行描述。在该采集模式下,锥束CT成像的物理过程可以概括为:

1) 三维体素绕转轴进行旋转;

2) 依据投影关系,将旋转后的三维体素投影到二维图像中。

为了描述三维体素的空间位置和旋转,本文定义了世界坐标系(图1)。世界坐标系中的基准轴为锥束CT的转轴,这里表示为 ZW; 原点 OW为旋转中心,即源在旋转轴线上的垂足; XW YW可以依据Cartesian右手坐标系的准则任意定义,为了求解方便,将 XW定义为源到旋转轴线的垂线, YW依据Cartesian右手坐标系唯一确定。

图1 锥束CT成像中坐标系的定义

为了描述三维体素与二维图像之间的投影关系,仿照相机标定中的定义[12],定义了成像坐标系,其中原点 OI为射线源的位置; XI YI分别与探测器的 u轴和 v轴平行; ZI为投影主轴,与探测器垂直。

因此,锥束CT的投影成像过程在数学上,可以描述为3个步骤:

1) 在世界坐标戏中,三维体素 P3D( x, y, z)绕转轴 ZW旋转角度 φ, 得到旋转后的位置 Pφ;

2) 将 Pφ从世界坐标系中转换成像坐标系中;

3) 利用投影成像关系,得到出 Pφ在探测器上的投影位置 P2D( ap, bp)。

2.2 锥束CT中三维体素与二维投影映射关系的求解

在得到锥束 CT成像的数学描述之后,三维体素与二维投影映射关系可以表示为

sP2D =A· R· R( φ P= proj( φ P3D.(6)

其中: R( φ)、 R A分别描述了锥束投影成像的数学描述中的3个步骤。其中矩阵 R( φ)为绕 ZW轴的旋转矩阵为

R(φ)=cosφ-sinφ0sinφcosφ0001.(7)

矩阵 R描述了从世界坐标系到成像坐标系之间的转换关系。矩阵 A描述了在成像坐标系中,三维空间点与二维投影位置之间的关系。其可以表示为

A=DSDdx0u000DSDdyv000010.(8)

其中: DSD为源到探测器的垂直距离;( u0, v0)为投影中心的位置,即源在探测器上的垂足位置。

在式(6)中,矩阵 R( φ)可以由旋转角度 φ直接获得。矩阵 A R可以通过相机标定方法进行求解,在下文中将进行详细描述。

2.2.1 相机标定方法在锥束CT中的应用

相机标定方法通过三维仿体中坐标精确已知的标志点与二维投影图像中标志点的投影坐标,求解出相机的内、外方位矩阵 A Rj。内方位矩阵 A描述了相机的内部几何结构,在不同的投影位置下应当保持一致。外方位矩阵 Rj可以用来描述标定物体在成像坐标系中的位置。其中 Rj可以表示为

Rj= rjtj01.(9)

在利用相机标定方法对锥束CT系统进行标定时,首先对标志点位置精确已知的标定模体进行旋转采集,从而在各个位置求解出 A Rj, 进一步得到式(6)中的矩阵 A和矩阵 R。这部分将分别在2)和3)中进行详细描述。这里的相机标定使用文[12]中所述的方法,其要求标定模体上至少有6个三维标志点。

2.2.2 矩阵 A的求解

式(6)中的矩阵 A即为相机标定中的内方位矩阵 A。在各个投影位置计算得到内方位矩阵 A后,对其进行平均运算,从而消除随机噪声,得到用于描述锥束CT投影关系的矩阵 A

2.2.3 矩阵 R的求解

矩阵 R代表了从世界坐标系到成像坐标系的转换关系。通过求解世界坐标系中的原点及各个坐标轴在成像坐标系中对应的单位向量,即可得到转换矩阵 R

首先求解世界坐标系的基准轴,即转轴 ZW。锥束CT采集过程中,标定仿体绕转轴旋转,在空间形成圆形轨迹。其圆心决定了转轴的位置,圆形平面的法向量决定了转轴的方向。在第 j个投影,得到外方位矩阵 Rj后,则可以求出仿体的原点在该投影下的位置即为式(9)中的 tj。若在圆形轨迹上均匀采集2 n个投影则圆形轨迹的圆心可以表示为

(10)

圆形平面的法向量,即 ZW在成像坐标系中的向量可表示为

(11)

其中: N = 2 n, k l分别取为1 /3和2 /3的投影数目。

依据世界坐标系的定义,原点 OW为源 OI在转轴 ZW上的垂足位置,有

OWI= CI- (OI-CI),ZWI)ZWI2· ZWI. (12)

X W轴在成像坐标系中的向量可以表示为

XWI= OI- OWI. (13)

Y W轴在成像坐标系中的向量可以表示为

YWI= ZWI× XWI. (14)

因此,世界坐标系原点和各个轴在成像坐标系中的单位向量可以组合为

VWI=[ XWI/‖ XWI‖, YWI/‖ YWI‖, ZWI/‖ ZWI‖, OWI]. (15)

世界坐标系到成像坐标系之间的转换矩阵可以表示为

R=VWI. (16)

至此,就求解出了式(6)中的矩阵 A R R( φ), 从而求解出三维体素与二维投影之间的映射关系。

2.3 加权因子 W( P)和反投影因子 U( P)的修正

在系统存在几何误差的情况下,三维空间体素的投影关系可以由图1表示。在图1的基础上,加权因子 W( P)可以表示为

W(P)=DSDOI-P2D. (17)

其中: D SD可以依据式(8)从矩阵 A中得到, P2D为空间体素的二维投影位置。

反投影因子 U( P)可以表示为

U(P)=DSDOI-P3D'. (18)

其中: P3D'为空间体素 P3D在投影主轴 ZI上的投影,其可以表示为

P3D'=( R· R( φ P3D, ZI). (19)

3 实验结果
3.1 仿真实验结果

在仿真实验中,使用镶嵌了螺旋分布标志点的仿体进行旋转采集,从而基于相机标定技术求解出用于描述三维体素和二维投影之间映射关系的矩阵 A和矩阵 R、 加权因子 W( P)、 反投影因子 U( P)。这里使用Shepp-Logan头模型来验证校正方法的有效性。

在本文中,分别对不同的锥束CT几何误差进行仿真,包括: 探测器俯仰5°、探测器偏转5°和探测器旋转5°, 以及( DSR+10mm)、( DSD+10mm)和投影中心 PR偏移10个像素( PR-(10,10))。本文中分别对未校正的重建结果、无噪声时的重建校正结果和存在噪声时的重建校正进行了仿真。误差主要来自于相机标定过程中标志点二维投影的拾取误差。这里利用方差为0.5像素的Gauss噪声对误差进行仿真。各种情况下的重建结果对比如图2所示。

图2 不同几何误差下的锥束CT重建校正结果

图2中可以看出,本文的方法对不同的几何误差对重建造成的影响均可以进行有效的校正。

此外,本文对不同重建结果中的系统调制传递函数(modulation transfer function, MTF)进行了分析,如图3所示。从图3中可以看出,使用本文方法进行几何误差校正后的MTF得到了有效的恢复。在不同的几何误差中,探测器的旋转误差和投影中心的偏移对重建图像的MTF影响最大,而探测器的俯仰误差对重建图像的MTF影响最小。

图3 不同几何误差下的MTF对比

不同几何误差情况下的空间分辨率如表1所示。空间分辨率为调制传递函数中50%的截止频率所对应的空间频率,以线对数(lp/mm)为单位。

图2图3表1中,可以看出,本文提出的方法在不同的几何误差的情况下,有效地提高了重建图像的质量和系统的空间分辨率。

表1 不同几何误差校正前后的空间分辨率(线对数)对比(lp·mm-1)
3.2 实际系统校正结果

本文中所使用的实际系统包括日本滨松公司(C7921-02, Hamamatsu, Japan)生产的CMOS 平板探测器, 牛津仪器(XTF5011, Oxford Instruments, USA)生产的固定阳极X射线源、转台系统。标定仿体使用有机玻璃圆柱体作为基底, 16个0.5mm的钢珠呈螺旋状分布在圆柱体的表面。

本文对标定仿体在0°~360°的旋转角度中,以30°为间隔,均匀地采集标定仿体的二维投影图像,从而标定出重建所需的投影矩阵和参数。而后本文分别对线仿体和小鼠进行成像,具体结果如图4示。

图4 实际系统中线模体与小鼠的成像结果对比

线模体的半峰全宽(FWHM)在未校正之前为1.03, 使用本文方法校正后为0.14。使用本文方法对线仿体半峰全宽的提高幅度为86.41%。结合图4的结果,可以看到本文的方法对实际系统中的几何误差也得到了良好的校正效果。

4 结论

本文对存在几何误差的锥束CT的图像重建校正校正方法进行了研究。本文中,锥束CT成像系统被看做是针孔相机。相机成像模型被用来求解三维体素与二维投影之间的映射关系,通过将这种映射关系带入到FDK重建算法中,就可以有效地对几何误差造成校正,从而提高重建图像的质量。本文使用成熟的相机标定方法求解成像模型,保证了方法的精度和可靠性。

The authors have declared that no competing interests exist.

参考文献
[1] Karolczak M, Schaller S, Engelke K, et al. Implementation of a cone-beam reconstruction algorithm for the single-circle source orbit with embedded misalignment correction using homogeneous coordinates[J]. Medical Physics, 2001, 28(10): 2050-2069. [本文引用:1] [JCR: 3.012]
[2] Sun Y, Hou Y, Hu J. Reduction of artifacts induced by misaligned geometry in cone-beam CT[J]. IEEE Transactions on Biomedical Engineering, 2007, 54(8): 1461-1471. [本文引用:1] [JCR: 2.233]
[3] Wiesent K, Barth K, Navab N, et al. Enhanced 3-D-recon-struction algorithm for C-arm systems suitable for interventional procedures[J]. IEEE Transactions on Medical Imaging, 2000, 19(5): 391-403. [本文引用:1] [JCR: 3.799]
[4] Galigekere R R, Wiesent K, Holdsworth D W. Cone-beam reprojection using projection-matrices[J]. IEEE Transactions on Medical Imaging, 2003, 22(10): 1202-1214. [本文引用:1] [JCR: 3.799]
[5] Navab N, Bani-Hashemi A, Nadar M S, et al. 3D reconstruction from projection matrices in a C-arm based 3D-angiography system [C]//Medical Image Computing and Computer-Assisted Intervention—MICCAI'98. Boston, USA: Medical Image Computing and Computer Assisted Intervention Society, 1998: 119-129 [本文引用:1]
[6] Navab A, Bani-Hashemi A, Mitshke M, et al. Dynamic geometric calibration for 3-D cerebral angiography [C]//Proceedings of SPIE Medical Imaging. Bellingham, USA: SPIE Medical Imaging, 1996: 361-370 [本文引用:1]
[7] Anne R, Picard C, Ponchut C, et al. Geometrical calibration of X-ray imaging chains for three-dimensional reconstruction[J]. Computed Medical Imaging and Graphics, 1993, 17(4): 295-300 [本文引用:1] [JCR: 1.496]
[8] Li X, Zhang D, Liu B. A generic geometric calibration method for tomographic imaging systems with flat-panel detectors—A detailed implementation guide[J]. Medical Physics, 2010, 37(7): 3844. [本文引用:1] [JCR: 3.012]
[9] Wang M, Ding H, Wang G. An improved FDK algorithm using camera calibration technique for misaligned CBCT reconstruction [C]//34th Annual International Conference of the IEEE Engineering in Medicine and Biology Society. San Diego, USA: IEEE Engineering in Medicine and Biology Society, 2012: 5991-5994. [本文引用:1]
[10] Wang M, Ding H, Wang G. A novel method using camera calibration technique for complete geometric calibration of cone-beam CT [C]//IFMBE Proceedings of World Congress on Medical Physics and Biomedical Engineering. Beijing, China: International Federation for Medical and Biological Engineering, 2013, 842-845 [本文引用:1]
[11] Feldkamp L A, Davis L C, Kress J W. Practical cone-beam algorithm[J]. Journal of Optics Society of American A, 1984, 1(6): 612-619. [本文引用:1]
[12] Abdel-Aziz Y I, Karara H M. Direct linear transformation into object space coordinates in close-range photogrammetry [C]//Proceedings of the Symposium on Close-Range Photo-grammetry. Chicago, USA: University of Illinois at Urbana-Champaign, 1971: 1-18. [本文引用:2]