CN106228584A - 锥束ct圆加直线轨迹反投影滤波重建方法 - Google Patents

Abstract

本发明涉及一种锥束CT圆加直线轨迹反投影滤波重建方法，针对圆加直线扫描轨迹的成像几何特点，选取M‑lines，通过求解穿过M‑line上重建点的R‑line的端点坐标，确定反投影积分区间，通过Hilbert变换，实现物体重建。本发明能够有效解决单圆轨迹大锥角扫描时锥角效应问题，在圆轨迹扫描的基础上，增加直线扫描轨迹，使其能够满足精确重建条件，获得完整的重建投影数据，能够较好的应用于实际扫描，有效提高大锥角重建图像质量，对于扩展锥束CT扫描应用范围具有重要的实用意义。

Description

锥束CT圆加直线轨迹反投影滤波重建方法

技术领域

本发明属于CT扫描成像技术领域，特别涉及一种锥束CT圆加直线轨迹反投影滤波重建方法。

背景技术

计算机断层成像(Computed Tomography,CT)技术能够在非接触、不破坏的情况下获得物体内部三维结构信息，目前已被广泛应用于工业无损检测和医疗影像诊断等领域。在CT系统中，重建算法的研究占有极其重要的地位。圆轨迹锥束扫描具有工程上易于实现、容易控制等优点，已成为实际CT系统应用最为广泛的扫描方式之一；然而，圆轨迹扫描并不能满足精确重建条件，其重建质量受锥角效应影响，即当锥角增大时，重建图像两端会出现严重的伪影，成像质量明显变差，进而影响CT技术检测效果。

目前，针对圆轨迹锥角效应问题的研究，依据其是否能够满足精确重建条件，大致可以分为两类。首先，依然采用圆轨迹扫描，对圆轨迹重建算法进行改进，通过一些合理的操作处理使得能够减少大锥角重建时的锥束伪影，改进方法主要有：数据重排方法、加权修正方法等，这类方法的优点在于能够保留圆轨迹扫描的优点，但是不能从根本上解决圆轨迹扫描数据缺失问题，当锥角较大时，改进效果受到限制；其次，可以在圆轨迹扫描的基础上，增加额外的扫描轨迹，使其能够满足精确重建条件，获得完整的数据，进行精确重建从而能够完全消除锥角效应，该类方法主要有圆加直线轨迹，双圆加直线轨迹，圆加圆弧轨迹等，相较于其它轨迹，直线轨迹相对简单，容易控制，且其能够更易获得缺失的数据，更适合轴向较长物体的扫描。当前对圆加直线轨迹重建算法的研究，主要集中在Katsevich算法领域，Katsevich算法最早针对螺旋轨迹提出，随后Katsevich等人将其推广到圆加直线轨迹上，Katsevich算法具有滤波反投影形式，在算法中PI线和构造因子是两个重要的概念，PI线决定了重建范围，构造因子不为零的Radon平面为重建做贡献，使用Katsevich算法进行精确图像重建时需要用到略大于TD窗的投影数据。

发明内容

为克服现有技术中的不足，本发明提供一种锥束CT圆加直线轨迹反投影滤波重建方法，有效解决现有技术中的锥角效应问题，有效提升大锥角扫描时物体重建质量。

按照本发明所提供的设计方案，一种锥束CT圆加直线轨迹反投影滤波重建方法，包含如下步骤：

步骤1、根据锥束CT系统工作原理，分别建立探测器坐标系、及针对圆加直线扫描轨迹和待扫描物体空间位置的世界坐标系；

步骤2、依据待扫描物体尺寸大小及物体重建圆柱区域，通过公式确定扫描所需的最短直线段长度，其中，H表示物体重建圆柱区域高度，r表示物体重建圆柱区域半径，R表示圆扫描轨迹半径；

步骤3、对待扫描物体进行照射投影，依据圆加直线轨迹扫描轨迹，通过M-line重建方法对物体待重建区域进行重建，得到M-line重建结果；

步骤4、对M-line重建结果进行重采样，得到最终重建物体。

上述的，步骤1中的圆加直线轨迹扫描轨迹表达式为：

其中，λ表示对应空间几何坐标变量，Z_l表示直线扫描轨迹长度；步骤2中的物体重建圆柱区域表述为：

∑＝{(x,y,z):x²+y²＜r²,-H＜z＜H}。

上述的，步骤3中具体包含如下内容：

步骤3.1、依据圆加直线轨迹扫描轨迹，选取一系列M-line，该一系列M-line覆盖整个重建物体；

步骤3.2、求解R-line的端点，确定差分反投影积分区间，并求解差分反投影值，其中，R-line穿过M-line上的待重建点x；

步骤3.3、依据差分反投影值，通过Hilbert变换，得到M-line重建点的Hilbert变换值；

步骤3.4、沿着M-line方向进行Hilbert滤波，得到M-line重建结果。

上述的，步骤3.2具体包含如下内容：

步骤3.2.1、设定待重建点x的坐标表示为(x₀,y₀,z₀)，M-line对应的扫描轨迹坐标表示为直线轨迹上R-line端点坐标表示为圆轨迹上R-line端点坐标表示为根据圆加直线扫描轨迹，R-line的端点坐标分别表示为：

步骤3.2.2、求解得到变量z₁，并通过圆加直线扫描轨迹表达式，得到直线轨迹上的反投影积分区间；

步骤3.2.3、根据R-line上的两点坐标(x₀,y₀,z₀)和(R,0,z₁)，得到该R-line在空间中的直线方程，根据公式求解得到z₁；

步骤3.2.4、根据z₁及圆加直线扫描轨迹表达式，得到直线轨迹上对应的差分反投影值λ′；根据R-line在空间中的直线方程及圆轨迹上的端点坐标，求解得到圆轨迹上对应的差分反投影值λ″。

本发明的有益效果：

本发明针对圆加直线扫描轨迹的成像几何特点，选取M-lines，通过求解穿过M-line上重建点的R-line的端点坐标，确定反投影积分区间，通过Hilbert变换，实现物体重建，能够有效解决单圆轨迹大锥角扫描时锥角效应问题，在圆轨迹扫描的基础上，增加直线扫描轨迹，使其能够满足精确重建条件，获得完整的重建投影数据，通过实验验证，能够较好的应用于实际扫描，有效提高大锥角重建图像质量，对于扩展锥束CT扫描应用范围具有重要的实用意义。

附图说明：

图1为锥束CT系统原理图；

图2为本发明的流程示意图；

图3为圆加直线扫描轨迹坐标系示意图；

图4为M-line选取示意图；

图5为R-line反投影积分区间示意图；

图6为仿真实验一重建结果示意图；

图7为仿真实验一重建结果剖线图；

图8为仿真实验二重建结果示意图。

具体实施方式：

在M-line重建理论中，对于确定的扫描轨迹来说，M-line为仅与其有一个交点的直线，R-line为与其有两个交点的线段，M-line上重建点的衰减系数f在该M-line上的Hilbert变换值与M-line上点的差分反投影密切相关，其中，M-line决定Hilbert变换的方向，R-line确定反投影的积分区间范围，其具体关系表达式为：

$K^{*}\left(w\right({\mathrm{\λ}}_3,x),x)=-\frac{1}{2\mathrm{\π}}\left(b_a\right(x,{\mathrm{\λ}}_3,{\mathrm{\λ}}_1)+b_a(x,{\mathrm{\λ}}_3,{\mathrm{\λ}}_2\left)\right),---\left(1\right)$

公式(1)左边表示M-line待重建点的Hilbert变换值。如果已经得到公式左边的Hilbert变换值，对其进行有限Hilbert逆变换即可得到重建值。其中w(λ,x)表示待重建点x指向射线源点的方向向量，其具体表达式如下：

$w(\mathrm{\λ},x)=\frac{x-\stackrel{\→}{a}\left(\mathrm{\λ}\right)}{\left|\right|x-\stackrel{\→}{a}\left(\mathrm{\λ}\right)\left|\right|},---\left(2\right)$

在公式(1)等号右边，λ₁和λ₂和λ₃为扫描轨迹上的三个扫描点，使用表示扫描轨迹空间几何位置，其中，λ₃为该M-line对应的扫描角度，和的连线穿过待重建点x，即表示为穿过该重建点的R-line。公式(1)中b_a与差分反投影有以下的关系：

$b_a(x,{\mathrm{\λ}}_1,{\mathrm{\λ}}_2)=b(x,{\mathrm{\λ}}_1,{\mathrm{\λ}}_2)+\underset{q=1}{\overset{2}{\Σ}}\frac{{(-1)}^{q}g({\mathrm{\λ}}_q,w({\mathrm{\λ}}_q,x\left)\right)}{\left|\right|x-\stackrel{\→}{a}\left({\mathrm{\λ}}_q\right)\left|\right|},---\left(3\right)$

b(x,λ₁,λ₂)表示在[λ₁,λ₂]区间上待重建点x的差分反投影，其具体的表达式与扫描轨迹相关，具体公式为：

$b(x,{\mathrm{\λ}}_1,{\mathrm{\λ}}_2)={\∫}_{{\mathrm{\λ}}_1}^{{\mathrm{\λ}}_2}\frac{D\left|\right|a^{\′}\left(\mathrm{\λ}\right)\left|\right|g_F(\mathrm{\λ},u,v)}{{\[\left(a\right(\mathrm{\λ})-x)\·{\stackrel{\→}{e}}_w\]}^2}d\mathrm{\λ}.---\left(4\right)$

其中为平板探测器的法向量。

下面结合附图和技术方案对本发明作进一步详细的说明，并通过优选的实施例详细说明本发明的实施方式，但本发明的实施方式并不限于此。

实施例一，参见图1和2所示，一种锥束CT圆加直线轨迹反投影滤波重建方法，包含如下步骤：

步骤1、根据锥束CT系统工作原理，分别建立探测器坐标系、及针对圆加直线扫描轨迹和待扫描物体空间位置的世界坐标系；

步骤2、依据待扫描物体尺寸大小及物体重建圆柱区域，通过公式确定扫描所需的最短直线段长度，其中，H表示物体重建圆柱区域高度，r表示物体重建圆柱区域半径，R表示圆扫描轨迹半径；

步骤3、对待扫描物体进行照射投影，依据圆加直线轨迹扫描轨迹，通过M-line重建方法对物体待重建区域进行重建，得到M-line重建结果；

步骤4、对M-line重建结果进行重采样，得到最终重建物体。

本发明针对圆加直线扫描轨迹的成像几何特点，选取M-lines，通过求解穿过M-line上重建点的R-line的端点坐标，确定反投影积分区间，通过Hilbert变换，实现物体重建，能够有效解决单圆轨迹大锥角扫描时锥角效应问题，在圆轨迹扫描的基础上，增加直线扫描轨迹，使其能够满足精确重建条件，获得完整的重建投影数据。

实施例二，参见图1～8所示，一种锥束CT圆加直线轨迹反投影滤波重建方法，包含如下内容：

根据锥束CT系统工作原理，分别建立探测器坐标系、及针对圆加直线扫描轨迹和待扫描物体空间位置的世界坐标系；依据待扫描物体尺寸大小及物体重建圆柱区域，确定扫描所需的最短直线段长度；对待扫描物体进行照射投影，依据圆加直线轨迹扫描轨迹，通过M-line重建方法对物体待重建区域进行重建，得到M-line重建结果；对M-line重建结果进行重采样，得到最终重建物体。

具体内容如下：

锥束CT系统的工作原理，如图1所示，圆加直线扫描轨迹描述为一个半径为R的圆C与垂直于该圆轨迹平面的直线L两部分构成，建立的坐标系如图3所示，探测器坐标系各方向向量用表示；对于世界坐标系，以圆轨迹中心O为坐标原点，以O与直线轨迹与圆轨迹的垂足连线方向为x轴正方向，以与直线轨迹平行的方向为z轴方向，满足右手坐标系，标定y轴方向，圆加直线轨迹分为三部分，具体表达式为：

其中，参数λ表示对应空间几何坐标的变量；对于待重建区域，将其看作为一个半径为r、高度为2H的圆柱体区域，中心放置于坐标系源点，重建区域具体表示为：∑＝{(x,y,z):x²+y²＜r²,-H＜z＜H},直线上半段和下半段部分具有对称关系，具体求解原理相同。

根据待重建物体尺寸大小及重建圆柱区域，求解得到对应的最短直线段长度；M-line重建中，重建与滤波均是沿着M-line方向进行，首先，选取一系列垂直于z轴的切片作为选取的M-lines，如图4所示，M-lines能够覆盖整个重建物体，对于M-line上的待重建点，满足至少一条R-line穿过，对于非圆轨迹平面上的物体部分进行重建，R-line的端点必定一个在圆轨迹上，一个在直线轨迹上，M-line对应的扫描轨迹坐标表示为直线轨迹上R-line端点坐标表示为圆轨迹上R-line端点坐标表示为差分反投影区间如图5中粗线条部分所示，设定待重建点x的坐标可表示为(x₀,y₀,z₀)，根据圆加直线轨迹几何表达式和各轨迹特点，R-line的端点坐标分别表示为： 其中，z₁和λ″是需要求解的变量，再求得z₁之后，将其带入到圆加直线扫描轨迹表达式中，即可得到直线轨迹上的反投影积分区间；由R-line上的两点坐标(x₀,y₀,z₀)和(R,0,z₁)，得到该条R-line在空间中直线方程：

$\frac{x-R}{x_0-R}=\frac{y-0}{y_0-0}=\frac{z-z_1}{z_0-z_1}=t,$

把圆轨迹上的端点坐标带入上述的直线方程中，得到：

$z_1=\frac{t}{t-1}{z}_0,\frac{Rsin\left({\mathrm{\λ}}^{\′\′}\right)}{y_0}=t,\frac{Rcos\left({\mathrm{\λ}}^{\′\′}\right)-R}{x_0-R}=t,$

显然，z₁＞z₀，得到t＞1，进一步得到关于t的方程：并进一步得t的取值：将其带入中，得到：

将z₁的结果带入到圆加直线扫描轨迹表达式中，得到直线轨迹上对应的轨迹点λ′；对于圆轨迹上点扫描角度λ″的求解，涉及到三角函数的运算，将带入到公式得到其对应的正弦值和余弦值，然后进行正负判断，确定λ″所在象限，进而确定λ″的值；对于圆加直线扫描轨迹，从最终求得到λ′和λ″结果可以看出，对于确定的扫描轨迹来说，其仅与待重建点x的空间几何坐标相关；对于确定的待重建点x，得到的反投影积分区间λ′和λ″均具有唯一值；则对于圆加直线扫描轨迹，其重建区域(圆轨迹平面区域除外)上任意待重建点有且仅有一条R-line穿过它。将得到的λ′和λ″带入到上述的公式(1)，即得到圆加直线轨迹M-line上重建点的Hilbert变换值，具体公式表达如下：利用该公式得到重建点的Hilbert变换值，然后再沿着M-line方向进行有限Hilbert滤波，即得到重建结果；由于重建和滤波是沿着M-line进行的，得到的重建结果是M-line上采样结果，因此，最后需要对其进行重采样，得到沿着笛卡尔坐标系下的重建结果，最终得到重建物体。

下面通过具体的数字仿真实验进一步验证本发明的有效性，并将实验结果与圆轨迹反投影滤波重建方法进行比较：

如图7所示，仿真实验使用的体模为多层圆盘(multi disk)体模，实验中各参数的设置如下，圆轨迹上射线源到物体中心的距离为100mm，射线源距离探测器的距离为200mm，圆轨迹上投影数据采样间隔为1度，共采集360张投影数据；直线轨迹的长度为60mm，直线轨迹上的采样间隔为0.01，共采集600张投影数据；使用的探测器探元大小为0.2，探测器尺寸为256x256，重建图像尺寸为256x256x256；本次仿真实验的重建锥角约为15°，对得到的投影数据分别使用圆轨迹BPF算法和本文推导的圆加直线轨迹反投影滤波重建算法进行重建，图7中：(a1)、(a2)、(a2)分别对应圆盘体模在xy平面、yz平面、xz平面切片重建结果；(b1)、(b2)、(b2)分别对应使用圆轨迹反投影滤波BPF重建方法在xy平面、yz平面、xz平面切片重建结果；(c1)、(c2)、(c2)分别对应使用本发明在xy平面、yz平面、xz平面切片重建结果。从图示的重建结果中可以看出，在大锥角重建时，圆轨迹BPF算法由于采集数据不完备，其重建图像受锥角效应影响严重，具体表现为随着轴向距离的增大其重建伪影越来越严重，重建质量明显下降；而本发明所提出的圆加直线轨迹滤波反投影重建方法，属于精确重建，理论上能够完全消除锥角效应的影响，从重建结果中可以看出，所得的重建图像伪影较少，重建质量得到明显改善，能够较好的还原出物体。为更好的从数值上进行分析，图7给出了重建结果沿着z轴方向的剖线图，剖线位置位于x＝128,y＝128处，由剖线结果得到，圆轨迹BPF算法随着轴向距离的增加，其与真值误差越来越大，而本发明的剖线与与真值基本吻合，说明本发明能够有效的抑制锥角效应，提高大锥角时的成像质量。

通过对工业设备连接件进行扫描试验进一步验证本发明在大锥角时重建的有效性，实验采用的射线源是Hawkeye 130KV，所用的平板探测器为Varian4030E，实验的具体参数为，射线源到旋转中心的距离为365mm，到探测器的距离为920mm，扫描的锥角大小为19.6°，在圆轨迹360°扫描范围内均匀的采集720张投影数据，采集的直线轨迹长度为13.20mm，在直线轨迹上共采集1320张投影数据，重建物体规模为880x880x1860，分别使用圆轨迹BPF方法和本发明进行重建，重建图像结果在xy平面、yz平面、xz平面切片如图8所示，图8中：(a1)、(a2)、(a2)分别对应使用圆轨迹BPF重建方法在xy平面、yz平面、xz平面切片重建图像结果，(b1)、(b2)、(b2)分别对应使用本发明在xy平面、yz平面、xz平面切片重建图像结果。由图7的重建结果可以看出，对于大锥角扫描时实际数据重建，圆轨迹BPF算法由于受到锥角效应的影响，其轴向方向重建结果具有明显的伪影，如图7中椭圆画定的区域所示，重建图像的顶端和低端成像质量变差，不能有效的重建出扫描物体该部分区域信息，严重影响了实际检测效果，这主要是因为使用圆轨迹BPF算法重建时，其轴向距离越大，缺失的数据越多；而本发明的圆加直线轨迹M-line重建算法的成像质量较好，由重建图像可以明显看出本发明重建结果几乎不受锥角效应的影响，能够对被扫描物体进行有效检测，这也证明了本发明能够很好的应用于实际扫描，能够有效的提高大锥角时的重建图像成像量，这对于扩展锥束CT扫描应用范围具有重要的实用意义。

本发明不局限于上述具体实施方式，本领域技术人员还可据此做出多种变化，但任何与本发明等同或者类似的变化都应涵盖在本发明权利要求的范围内。

Claims (4)

1.一种锥束CT圆加直线轨迹反投影滤波重建方法，其特征在于：包含如下步骤：

步骤1、根据锥束CT系统工作原理，分别建立探测器坐标系、及针对圆加直线扫描轨迹和待扫描物体空间位置的世界坐标系；

步骤2、依据待扫描物体尺寸大小及物体重建圆柱区域，通过公式确定扫描所需的最短直线段长度，其中，H表示物体重建圆柱区域高度，r表示物体重建圆柱区域半径，R表示圆扫描轨迹半径；

步骤3、对待扫描物体进行照射投影，依据圆加直线轨迹扫描轨迹，通过M-line重建方法对物体待重建区域进行重建，得到M-line重建结果；

步骤4、对M-line重建结果进行重采样，得到最终重建物体。

2.根据权利要求1所述的锥束CT圆加直线轨迹反投影滤波重建方法，其特征在于：步骤1中的圆加直线轨迹扫描轨迹表达式为：其中，λ表示对应空间几何坐标变量，Z_l表示直线扫描轨迹长度；步骤2中的物体重建圆柱区域表述为：∑＝{(x,y,z):x²+y²＜r²,-H＜z＜H}。

3.根据权利要求2所述的锥束CT圆加直线轨迹反投影滤波重建方法，其特征在于：所述步骤3中具体包含如下内容：

步骤3.1、依据圆加直线轨迹扫描轨迹，选取一系列M-line，该一系列M-line覆盖整个重建物体；

步骤3.2、求解R-line的端点，确定差分反投影积分区间，并求解差分反投影值，其中，R-line穿过M-line上的待重建点x；

步骤3.3、依据差分反投影值，通过Hilbert变换，得到M-line重建点的Hilbert变换值；

步骤3.4、沿着M-line方向进行Hilbert滤波，得到M-line重建结果。

4.根据权利要求3所述的锥束CT圆加直线轨迹反投影滤波重建方法，其特征在于：步骤3.2具体包含如下内容：

步骤3.2.1、设定待重建点x的坐标表示为(x₀,y₀,z₀)，M-line对应的扫描轨迹坐标表示为直线轨迹上R-line端点坐标表示为圆轨迹上R-line端点坐标表示为根据圆加直线扫描轨迹，R-line的端点坐标分别表示为：

步骤3.2.2、求解得到变量z₁，并通过圆加直线扫描轨迹表达式，得到直线轨迹上的反投影积分区间；

步骤3.2.3、根据R-line上的两点坐标(x₀,y₀,z₀)和(R,0,z₁)，得到该R-line在空间中的直线方程，根据公式求解得到z₁；

步骤3.2.4、根据z₁及圆加直线扫描轨迹表达式，得到直线轨迹上对应的差分反投影值λ′；根据R-line在空间中的直线方程及圆轨迹上的端点坐标，求解得到圆轨迹上对应的差分反投影值λ″。