三源螺旋CT精确重建算法在动态重建中的应用

作者：金燕南 赵俊 刘尊钢 庄天戈    作者单位：上海交通大学 生物医学工程系 上海 200240；　基金项目: 国家自然科学基金（30570511；30770589）资助
【摘要】  目的 将三源CT精确重建算法应用到动态重建的场合，用于抑制心脏或呼吸运动在CT成像中产生的伪迹。方法 采用单源与三源CT精确重建算法对两种动态数字模型进行仿真重建。通过比较三源CT与单源CT的重建结果，评价三源精确重建算法在动态重建方面的性能。结果 三源CT的扫描方式能够有效抑制动态重建中的伪迹，改善重建质量。结论 本文的结果证实了多源CT精确重建算法在动态重建中的有效性，对于进一步改善心脏成像及小动物成像质量具有积极的实际意义。

【关键词】  三源CT;螺旋CT;精确重建算法;动态重建

 

　　Application of Exact Reconstruction Algorithm of Triplesource Helical CT in Dynamic Reconstruction.

　　JIN Yannan, ZHAO Jun, LIU Zungang, ZHUANG Tiange.

　　Space Medicine & Medical Engineering,2009,22(6):453～458

　　Abstract: Objective To reduce the artifacts produced by cardiac or respiratory motion,through the application of exact reconstruction algorithm with triplesource CT in dynamic reconstruction. Methods Simulation studies on two dynamic phantoms were performed with exact reconstruction algorithms of both singlesource and triplesource spiral CT. The performance of triplesource spiral CT in dynamic reconstruction was evaluated with comparing the reconstruction results from triplesource and singlesource configurations. Results The triplesource configuration could effectively reduce the artifact in dynamic reconstruction. Conclusion The results of this paper confirm the validity of the exact reconstruction algorithm of multisource CT, which has the potential to improve the image quality in cardiac imaging and small animal imaging.

　　Keywords: triplesource CT;helical CT;exact reconstruction algorithm;dynamic reconstruction

　　在CT成像，特别是心脏成像和小动物成像中，运动伪迹严重影响着重建图像质量。近年来，动态重建已经引起了广泛关注。在心脏成像中普遍采用门控信号加以校正［1］；Lu等［2］提出呼吸门控技术用于抑制呼吸伪迹；King等［３］采用ROI重建技术，弃用部分投影数据；Yu等［４］提出一种基于推广的HelgasonLudwig连续性条件的方法来补偿运动伪迹。西门子公司推出了双源CT并用于心脏成像，取得了很好的效果［56］。
   
　　为了进一步提高数据采集速度，抑制运动伪迹，笔者提出了一种三源螺旋CT的扫描方式并给出了相应的精确重建算法［78］。与西门子公司的双源CT采用的近似重建算法不同，三源精确重建算法从Katsevich［9］提出的FBP精确重建算法及Zou、Pan［10］提出的BPF精确重建算法中推广得到。这种算法基于跨螺旋PI线和最小检测窗的概念，其z轴方向分辨率与z轴覆盖范围无关。由于采用三源的扫描方式，获取完整投影数据的时间间隔缩短为单源情形的三分之一，更加适合动态重建的场合。
   
　　本文将在前期工作的基础上，将三源精确重建算法应用到动态重建中，并与单源精确重建算法进行比较，以评价三源精确重建算法在动态重建方面的性能。

　　方 法

　　三源CT扫描结构

　　三源螺旋CT的扫描方式如图1所示。3个X射线源对称放置，每个X射线源发出的X射线穿过被测对象后，被对侧的平板检测器接收。当X射线与检测器一起环绕被测对象作同一方向的连续旋转时，放置被测物体的检查床沿旋转轴做直线运动。这时以被测对象为参考系观察，X射线的运动轨迹为三螺旋，用数学式表示为:a1(t)=R cos(t),R sin(t),h2πt,a2(t)=R cos(t+12),R sin(t+12),h2πt,a3(t)=R cos(t+12+23),R sin(t+12+23),h2πt,(1)

　　其中R表示螺旋线的半径，h表示X射线源每旋转1圈时检查床前进的距离。

　　对于给定的投影角度t，锥形束投影变换定义为:

　　Df(aj(t),θj)=∫0∞ f(aj(t)+sθj)ds(2)
　　
　　其中θj表示从X射线源aj(t)指向空间一点x的单位向量：

　　θj=axj(t)=x-aj(t)‖x-aj(t)‖(3)

　　 三源精确重建算法

　　三源的精确重建算法基于跨螺旋PI线的概念。如图2所示，考虑不同螺旋上的两个点ai(tmini)，aj(tmaxj)，i∈{1,2,3}，i≠j它们对应的投影角度满足tmaxj-tmini<2π。跨螺旋PI线定义为连接这两点的直线，该直线在螺旋所围成的圆柱以内的部分称为跨螺旋PI线线段。为了表达方便，我们把ai(tmini)记为abi，aj(tmaxj)记为aej，ai(tmini)与aj(tmaxj)之间的跨螺旋PI线表示为aej-abi。
   
　　三源精确算法用算子形式表示为:

　　f(x)=12Haekmod3+1-abk D#kD±f(x),k∈{1,2,3}(4)

　　其中算子表示投影相对于投影角度t的偏微分。反投影算子D#k定义为:

　　D#kg(x)=∑3j=1ηjk∫Ijdt‖x-aj(t)‖g(aj(t),axj(t)),k∈{1,2,3}(5)

　　ηjk表示3个X射线的投影数据在反投影上的加权系数，定义为:

　　ηjk=(-1)(J-K)mod32,k∈{1,2,3}(6)

　　Hilbert算子的表达式如下:

　　Heπ［φ(x)］=1πj∫R3sgneπ(x-y)φ(y)dy(7)

　　其中

　　sgneπ(z)=1(2π)3∫R3sgn(eπ·v)eiz·vdv(8)

　　算法实现

　　以zπ,z2,z3为坐标轴建立局部坐标系，其中zπ轴与向量eπ的方向相同，而z2,z3轴与ezπ垂直。令zπ,z2,z3在频域中对应的变量为vπ,v2,v3，则：sgneπ(z)=1(2π)3∫R3sgn(e·v)·ejz·vdv=1(2π)3∫R3sgn(vπ)·ej(zπvπ+z2v2+z3v3)dv=12π∫Rsgn(vπ)·ejzπvπdvπ
∫Rejz2v2dv2∫Rejz3v3dv3=-1jπzπδ(z2)δ(z3)(9)
   
　　令局部坐标系中点x,y的坐标为(xπ,x2,x3),(yπ,y2,y3)。将(9)式代入(7)式，得Heπ［φ(x)］=1πj∫R3sgneπ(x-y)φ(y)dy
=1πj∫R3-1jπ(xπ-yπ)δ(x2-y2)δ(x3-y3)φ(yπ)dyπdy2dy3(10)

　　根据δ函数的性质，只有当x2-y2=0,x3-y3=0(即x-y的方向与向量eπ的方向相同)时，积分号内的表达式才不为0，此时：Heπ［(x)］=1π2∫R1(xπ-yπ)(yπ)dyπ(11)

　　根据重建公式(4)，滤波可以沿三条跨螺旋PI线中的任何一条进行。不失一般性，考虑a1(t)与a2(t)之间的跨螺旋PI线ae2-ab1。令跨螺旋PI线线段两个端点对应的投影角度分别为tmax2和tmin1，用xπ12表示沿ae2-ab1的坐标轴，建立跨螺旋PI坐标系，这时，待重建物体的密度函数可以写成fπ12(xπ12,tmin1,tmax2)，反投影可以成为D#1(xπ12,tmin1,tmax2)。
   
　　BPF精确算法的实现可以分为以下4个步骤：

　　 对投影数据求导  对于平板检测器的CT装置，我们用u(j),v(j)表示aj(t)对应的检测器平面上点的二维坐标。对于一个给定的投影角度t,(u(j),v(j)）由x确定。于是，我们可以把锥形束投影表示为：

　　gj(t,u(j),v(j))=Df(aj(t),θj)(12)

　　与算子对应，首先要计算gj(t,u(j),v(j))对投影角度的微分Dgj(t,u(j),v(j))=ddtgj(t,u(j),v(j))。

　　对求导后的数据反投影  用跨螺旋PI线线段的端点对应的投影角度作为积分的上下限，式(5)的反投影算子可以改写成：

　　D#1(xπ12,tmin1,tmax2)=∑3j=1ηjk∫tjmaxtminjDgj(t,u(j),v(j)）‖x-aj(t)‖dt(13)

　　有限Hilbert变换的逆变换  从(11)式可以看到，(7)式给出的Hilbert算子实际上是沿ae2-ab1进行一维的Hilbert变换。即：

fπ12(xπ12,tmin1,tmax2)=12πH［D#1(xπ12,tmin1,tmax2)］(14)
   
　　其中H表示一维Hilbert变换，上面的变换要用到无限长的反投影数据，这显然无法实现。对公式(14)两边同时进行Hilbert变换，利用性质H［H((x))］=-(x)，得：D#1(xπ12,tmin1,tmax2)=-2∫ULdx′π12(xπ12-x′π12)fπ12(x′π12,tmin1,tmax2)(15)

　　由于重建物体的密度函数fπ12(xπ12,tmin1,tmax2)在实数集上是紧支撑的（即X0,使得xπ12>X0，fπ12(xπ12,tmin1,tmax2)=0），令L和U表示其紧支集的下界和上界，则(15)式成为求解有限Hilbert逆变换的问题，我们采用［11］中所述的方法求解。表1  动态球体模型参数设置（略）表2  仿真重建参数设置（略）
   
　　跨螺旋PI坐标转换成直角坐标  如图2所示，给定空间一点x∈ΩR3，求出该点在跨螺旋PI坐标系中对应的（xπ12,tmin1,tmax2)，由于这时求出的坐标值未必在反投影的采样网格上，因此要用到三维内插。数据重排之后得到笛卡儿坐标系中的重建灰度值f(x)。
   
　　动态模型定义

　　本文采用2个动态数字模型来评价重建算法的动态特性。第一个模型由6个球体组成，其参数在表1中给出。最外层的大球（球体1、球体2）静止不动，位于中央的小球（球体7）半径随时间作正弦变化，其它4个小球（球体3～球体6）则分别沿x轴或y轴方向做正弦振动表中x0,y0,z0分别表示球心在空间的位置，r表示球体半径，ρ表示球体密度。

　　为了使仿真更加接近临床应用的情况，本文采用的第二个数字模型为动态胸部模型，该模型能够模拟人体胸部的主要结构以及心脏的运动。模型呈长方体结构，长、宽均为256像素，高140像素。模型由2种规则几何图形构成：椭球和椭圆柱（包括圆柱）。模型包括：躯干；左、右肺叶；腹腔；左、右臂；左、右肱骨；胸骨；九节脊椎骨以及一个按周期跳动的心脏。各组织采用不同的密度值以仿真实际的生理结构。
   
　　参考King等［3］在心脏动态仿真中采用的模型，心脏的半径应乘以下面的系数：A(t)=1-0.625sin(πt0.33T)  0≤t<0.33 T1 0.33T≤t    
　　即在一个周期内，心脏在前33%的时间快速收缩，随后进入静息期，半径保持不变。图3给出了动态胸部模型的横切面和透视图，动态胸部模型具体的参数设置参见［12］。

　　仿真实验参数设置

　　仿真实验中投影数据采用解析法求得。为了保证重建结果的可比性，单源与三源都采用如表2所示的参数设置。值得指出的是，在平板检测器的高度相同的情况下，三源螺旋CT的螺距是单源情形的3倍，因此成像速度也是单源情况的3倍。因此，三源CT在单位时间内的剂量增大，但由于成像时间缩短，病人接受辐射的总剂量并没有发生变化。

　　结果与讨论

　　〖HT5SS]采用解析法计算动态模型中给出的2个动态数组模型的投影值，再分别采用单源和三源精确重建算法进行重建，结果如图4、图5所示。
    图4为动态球体模型的重建结果，这里给出的是z=0对应的横截面。其中(a)对应单源重建算法的重建结果，(b)对应三源重建算法的结果，灰度显示窗口为［0,2］。从重建结果可以看出，采用三源CT的扫描方式能够有效地抑制运动伪迹。
   
　　图5给出了动态胸部模型的重建结果，灰度显示窗口为［0.5,1.5］。此时的动态伪迹主要由于心脏运动产生，仿真中设心脏每分钟跳动60次，X射线源转速为1圈/s。由于三源CT数据扫描速度更快，所以在投影采集过程心脏的运动造成的影响也比较小。图(a)与图(b)分别对应单源与三源的重建结果。比较图(a)、(b)同样可以发现，单源CT比三源CT更容易收到运动伪迹的影响。

　　图6给出了动态胸部模型重建结果的剖面线，三条剖面线的灰度值分别对应图3(a)，图5(a)和图5(b)。可以看到三源重建算法（虚线）与数字模型吻合得很好，而单源重建算法（实线）在心脏附近出现伪迹。

　　值得指出的是，被测对象某个部分的运动不仅会造成该部分及其周围组织的模糊，而且会通过Hilbert滤波的运算将运动产生的误差传递到同一条PI线（或跨螺旋PI线）上，从而影响整幅重建图像。横切面上的每一个点位于不同的PI线（或跨螺旋PI线）上，因此经过数据重排之后得到的图4、图5实际上是不同滤波方向上的运动伪迹综合作用的结果。三源的扫描方式可以减小重建某一像素用到的投影对应的数据采集时间，因此能够提高时间分辨率，提高整幅图像的重建质量。
   
　　赵俊等［78］在数学上证明了三源精确重建的可行性，同时给出了SheppLogan模型的仿真重建结果。然而，由于SheppLogan模型为静态模型，无法评估重建算法的时间分辨率。本文采用的动态球体模型可以直观地检验算法在不同方向上的时间分辨率，动态胸部模型可以仿真心脏搏动造成的动态伪迹，由此检验三源精确重建算法的性能。动态胸部模型的尺寸，材料及心动周期都能很好地模拟人体胸部的生理特征［13］，三源CT的扫描参数（即几何结构与扫描速度）也与目前的医用CT相符。因此本文给出的仿真重建结果对于三源精确重建算法的临床应用具有参考价值。
    三源精确重建算法采用沿着跨螺旋PI线的一维滤波，与迭代重建算法［13］相比，算法重建速度更快，效率更高。值得注意的是，三源精确重建算法在提高成像速度的同时保持了相同的X射线辐射剂量［78］，这对于心脏的三维重建［14］有着重要的实际意义。

　　结 论

　　本文采用三源CT的扫描方式对动态数字模型进行仿真重建，并与单源精确重建算法的结果进行比较。仿真结果表明，三源CT能够有效抑制动态重建中的伪迹，改善重建质量。本文的结果证实了多源CT精确重建算法在动态重建中的有效性，对于进一步改善心脏成像及小动物成像质量具有积极的实际意义。

【参考文献】
  　　［1］Kachelrie M, Ulzheimer S, Kalender WA. ECGcorrelated image reconstruction from subsecond multislice spiral CT scans of the heart. ［J］. Medical Physics, 2000, 27:18811902.

　　［2］Lu W, Parikh PJ, Hubenschmidt JP, et al. Reduction of motion blurring artifacts using respiratory gated CT in sinogram space: A quantitative evaluation. ［J］. Medical Physics, 2005, 32:32953304.

　　［3］King M, Pan X, Yu L, Giger M. Regionofinterest reconstruction of motioncontaminated data using a weighted backprojection filtration algorithm. ［J］. Medical Physics, 2006,33:12221238.

　　［4］Yu H, Wei Y, Hsieh J, et al. Data consistency based translational motion artifact reduction in fanbeam CT. ［J］. IEEE Trans Med Imaging, 2006,25(6):792803.

　　［5］Flohr TG, McCollough CH, Bruder H, et al. First performance evaluation of a dualsource CT (DSCT) system［J］. European Radiology, 2006,16(2):25668.

　　［6］Johnson TRC, Nikolaou K, Wintersperger BJ, et al. Dualsource CT cardiac imaging: initial experience. ［J］. European Radiology, 2006,16(7):140915.

　　［7］Zhao J, Jiang M, Zhuang T, et al. Minimum detection window and interhelix PIline with triplesource helical conebeam scanning ［J］. Proceedings of SPIE. 2004, 5535:599610.

　　［8］Zhao J, Jiang M, Zhuang T, et al. An exact reconstruction algorithm for triplesource helical conebeam CT［J］. Journal of XRay Science and Technology, 2006, 14(3):191206.

　　［9］Katsevich A. An improved exact filtered backprojection algorithm for spiral computed tomography ［J］. Advances in Applied Mathematics, 2004, 32(4):68197.

　　［10］Zou Y, Pan X. Exact image reconstruction on PIlines from minimum data in helical conebeam CT. ［J］. Physics in Medicine and Biology, 2004, 49(6):94159.

　　［11］Noo F, Clackdoyle R, Pack JD. A twostep Hilbert transform method for 2 D image reconstruction ［J］. Physics in Medicine and Biology, 2004, 49(17):39033923.

　　［12］Zhao J, Liu Z, Zhuang T. Method for evaluating the performance of image reconstruction algorithm using dynamic chest phantom ［P］. Chinese Patent Application. Application No.200610027986.7, Open No.CN1877638, 2006.

　　［13］Gao X, Xia S, Wang Y. Multiobjective vector optimization algorithm for medical image reconstruction from incomplete projections［J］. Space Medicine & Medical Engineering, 2004, 17(4): 270274

　　［14］WANG Gang, RAO Nini. Research on registration algorithm for threedimension cardiac electrical maps［J］.Space Medicine & Medical Engineering, 2008, 21(6):523528.