CN103714578A - 针对半覆盖螺旋锥束ct的单层重排滤波反投影重建方法 - Google Patents

Abstract

本发明涉及一种针对半覆盖螺旋锥束CT的单层重排滤波反投影重建方法，首先将锥束投影重排成多排扇形束投影，每一排对应物体的一层，然后利用平行束投影的对称性质将每一排的扇形束投影重排成π范围内的平行束投影，此时平行束投影完全覆盖物体横截面，不存在数据截断，最后通过在重排前的锥束几何下计算反投影点，进行三维滤波反投影，避免了将第一次的重排误差引入重建图像。通过采用上述方法，能够实现半覆盖螺旋锥束几何下快速、高质量的三维重建，重建图像中没有明显截断伪影，重建质量不受数据截断位置的影响，能够更大程度地扩展螺旋锥束CT的成像视野。

Description

针对半覆盖螺旋锥束CT的单层重排滤波反投影重建方法

技术领域

本发明涉及一种针对半覆盖螺旋锥束CT的单层重排滤波反投影重建方法。

背景技术

螺旋锥束CT能够解决长物体的成像问题，并且具有扫描速度快、空间分辨率高和射线利用高的优点，在实际中被越来越广泛地使用。传统的螺旋锥束CT只能对在各个投影角度下横截面完全被射束覆盖的长物体进行成像，但是，在实际应用中经常会遇到应用小面阵探测器对大物体进行成像的问题，比如对发动机的无损检测，此时物体的横截面在各个投影角度下不能被射线全部覆盖，从而导致传统的螺旋锥束CT失效。

为了解决该问题，一些新的扫描方式和相应的重建方法被提出，如双螺旋扫描方式和双螺旋FDK、双螺旋Katsevich以及双螺旋BPF型重建方法，能够实现螺旋锥束CT接近2倍的横向成像视野扩展；又如多次螺旋扫描方式和相应的FBP型重建方法，能够实现更大的横向视野扩展率。但是，这两种方法在工程实施上均存在一些难度，于是一种更为简单的方法便是实现半覆盖螺旋锥束CT成像。

半覆盖螺旋锥束CT成像能够扩展接近2倍的横向视野，其可由两种方式实现，一种是载物台偏置，一种是探测器偏置，两种实现方式本质上相同，重建方法能够相互适用。例如，在载物台偏置方式下实现半覆盖螺旋锥束CT成像时，X射线源和探测器固定不动，载物台承载被扫描物体首先沿平行探测器方向偏置，使每个投影角度下X射线只覆盖被扫描物体横截面的一半，然后进行螺旋轨迹运动，同时在预定的投影角度下探测器采集投影数据。重建时，与标准螺旋锥束CT相比，半覆盖螺旋锥术CT只利用一半的投影数据重建整个超视野物体，因此能够扩展成像视野。针对半覆盖螺旋锥束CT，邹晓兵等人提出了半覆盖螺旋FDK方法，相比标准螺旋FDK，该方法使用一半的投影数据重建整个物体，但是，每个投影角度下数据存在横向截断，而该方法采用全局的斜坡滤波器，因此重建结果存在较为严重的截断伪影。对此问题，郭吉强等人又提出了改进的半覆盖螺旋FDK方法，该方法借鉴二维CT中局部滤波的思想，重建过程中采用一种新的具有局部特性的滤波器，从而能够在一定程度上抑制由于数据截断和全局滤波带来的截断伪影。但是，其重建结果中截断伪影依然存在，尤其投影数据截断位置更加接近物体中心时表现更为明显。

发明内容

为克服以上现有技术的不足，本发明要解决的技术问题是提供一种针对半覆盖螺旋锥束CT的单层重排滤波反投影重建方法，该方法能够实现半覆盖螺旋锥束几何下快速、高质量的三维重建，重建图像中没有明显截断伪影，重建质量不受数据截断位置的影响，能够更大程度地扩展螺旋锥束CT的成像视野。

本发明的技术方案是：一种针对半覆盖螺旋锥束CT的单层重排滤波反投影重建方法，该方法基于被扫描物体的螺旋锥形束投影进行重建，主要包括以下步骤：

（1）将螺旋锥形束投影重排成多层扇形束投影，重排后每个角度下的锥形束投影数据转化为互相平行的多层扇形束投影数据，每层扇形束数据对应被扫描物体沿旋转轴方向的一层；

（2）将步骤（1）中的多层扇形束投影重排成多层平行束投影，该重排只涉及被扫描物体相同层的扇形束投影，即重排后的多层平行束投影仍对应被扫描物体的相同层；

（3）将步骤（2）中的平行束投影数据进行滤波，然后对滤波后的数据进行反投影，从而实现重建；并且，其中的反投影过程是在锥形束几何下进行的。

上述单层重排滤波反投影重建方法，其中经步骤（1）重排后的被扫描物体每层对应2π范围内的扇形束投影，且每个角度下对应的扇形束投影只覆盖被扫描物体横截面的一半。

上述单层重排滤波反投影重建方法，其中步骤（2）中的重排利用了平行束投影的对称性质，重排后的平行束投影角度范围为π，且每个角度下的投影完全覆盖被扫描物体的横截面，不存在数据截断。

上述单层重排滤波反投影重建方法，其中被扫描物体的锥形束投影数据的具体获取过程为：构建半径为R₀的圆柱体，假设被扫描物体即包含于该圆柱体中，在该圆柱体的中心位置建立固定的右手笛卡尔坐标系(x,y,z)，且被扫描物体的偏置量为D₀，成像时，被扫描物体做螺旋轨迹运动，即相当于被扫描物体固定、射线源和探测器同时做反向的螺旋轨迹运动；旋转轴与Z轴重合，射线源的旋转半径为R，射线源到探测器的距离为S，射线源到平行于探测器且过旋转轴的平面的距离为D，建立旋转的坐标系(u,v,w)，其中心为探测器中心，对于探测器上的任何一点，都可以用坐标(u,v)表示；为获得重建所需的足够投影数据，每个角度下的射线至少覆盖被扫描物体横截面的一半，即需要满足：

$\left\{\begin{array}{c}{D}_0<D\tan \left({\mathrm{\&gamma;}}_m\right)\\ (D\tan \left({\mathrm{\&gamma;}}_m\right)+D_0)\cos \left({\mathrm{\&gamma;}}_m\right)>R_0\end{array}\right.$

其中γ_m表示最大的扇角的一半。

上述单层重排滤波反投影重建方法，其中步骤（1）中的由锥形束投影到扇形束投影的具体重排过程为：每层扇形束数据对应被扫描物体沿Z轴方向的一层，重排后的扇形束几何与锥形束几何相匹配，即虚拟的扇形束源点同样位于距旋转轴距离为R的位置，线阵探测器与面探测器的行匹配同时对应被扫描物体沿Z轴方向的一层；重排前的锥形束投影数据为b(λ,u,v)，其中λ表示投影角度，(u,v)表示面探测器上的点，重排后的扇形束投影数据为p_z(θ,u)，其中z表示物体沿Z轴方向的层数，θ表示扇形投影角度，且θ=λ%2π，u表示线阵探测器探元位置，则该步骤下的重排过程可由下面公式表示：

$p_z(\mathrm{\&theta;},u)=\frac{\sqrt{u^2+S^2}}{\sqrt{u^2+v^2+S^2}}b(\mathrm{\&lambda;},u,v)$

其中

Δz表示锥形束源点与虚拟扇形束源点的距离。

上述单层重排滤波反投影重建方法，其中步骤（2）中的由扇形束投影到平行束投影的具体重排过程为：

$g_z(\mathrm{\&phi;},t)=\left\{\begin{array}{c}{p}_z(\mathrm{\&phi;}-\mathrm{\&gamma;},D\tan \left(\mathrm{\&gamma;}\left(t\right)\right)),t<=0\\ p_z(\mathrm{\&pi;}+\mathrm{\&pi;}-\mathrm{\&gamma;},D\tan \left(\mathrm{\&gamma;}(-t)\right)),t>0\end{array}\right.$

其中，g_z(φ,t)表示重排后的平行束投影，φ表示平行束投影角度，t表示旋转中心到每条射线的距离，并且φ∈[0,π]，t∈[-R₀,R₀]，γ表示带符号的扇角。

上述单层重排滤波反投影重建方法，其中在步骤（3）中，采用标准的斜坡滤波器分别对被扫描物体每一层的每个投影角度的数据进行滤波，即

${\tilde{g}}_z(\mathrm{\&phi;},t)=g_z(\mathrm{\&phi;},t)*h\left(t\right)$

其中，

表示滤波后的数据，h(t)为斜坡滤波器；并且，所述反投影的具体过程为

其中，f(x,y,z)表示三维重建结果，

表示平行束投影角度，

表示滤波后的平行束投影数据，

和

表示点(x,y,z)在进行步骤（1）重排前的锥形束几何下计算的投影位置。

本发明的有益效果是：通过采用上述方法，在实现半覆盖螺旋锥束CT成像的过程中，不仅可以扩展螺旋锥束CT的横向视野，使能够检测横截面尺寸大于成像视野的长物体，并且对截断伪影抑制能力强，重建图像中没有明显截断伪影，重建质量不受数据截断位置的影响，从而能够快速、高质量地重建三维超视野物体。

附图说明

下面结合附图对本发明的具体实施方式作进一步详细的说明。

图1是载物台偏置下的半覆盖螺旋锥束CT成像示意图；

图2是半覆盖螺旋锥束CT几何示意图；

图3是锥形束投影到扇形束投影的重排示意图，其中（a）为锥形束投影，（b）为扇形束投影，（c）为重排几何关系示意图；

图4是扇形束投影到平行束投影的重排几何关系示意图；

图5是三组仿真实验的重建结果示意图，其中（a）、（b）、（c）分别对应第一组、第二组和第三组仿真实验，并且每一组中的第一列、第二列和第三列分别为重建结果的x-y、y-z和x-z中心切片；

图6是三组仿真实验的重建结果中x-y中心切片的中心纵向剖线图，其中（a）、（b）、（c）分别对应第一组、第二组和第三组仿真实验；

图7是实际实验下的增压器叶轮半覆盖螺旋扫描重建结果示意图。

具体实施方式

如图1所示，本发明在载物台偏置方式下，建立一种针对半覆盖螺旋锥束CT的单层重排滤波反投影重建方法。该方法是基于被扫描物体的螺旋锥形束投影进行重建的，首先对半覆盖螺旋锥束CT几何做必要的说明。如图2所示，构建半径为R₀的圆柱体，假设被扫描物体即包含于该圆柱体中，在该圆柱体的中心位置建立固定的右手笛卡尔坐标系(x,y,z)，为了实现半覆盖扫描必须实现被扫描物体的偏置，偏置量用D₀表示。成像时，被扫描物体做螺旋轨迹运动，即相当于被扫描物体固定、射线源和探测器同时做反向的螺旋轨迹运动，旋转轴与Z轴重合，射线源的旋转半径为R。定义射线源到探测器的距离为S，射线源到平行于探测器且过旋转轴的平面的距离为D，建立旋转的坐标系(u,v,w)，其中心为探测器中心，对于探测器上的任何一点，都可以用坐标(u,v)表示。

为获得重建所需的足够投影数据，每个角度下的射线至少覆盖被扫描物体横截面的一半，即偏置量需要满足：

$\left\{\begin{array}{c}{D}_0<D\tan \left({\mathrm{\&gamma;}}_m\right)\\ (D\tan \left({\mathrm{\&gamma;}}_m\right)+D_0)\cos \left({\mathrm{\&gamma;}}_m\right)>R_0\end{array}\right.$

其中γ_m表示最大的扇角的一半。

在获得足够的投影数据后，本发明通过以下3个步骤实现快速、高质量的三维重建：1、锥形束投影到扇形束投影的重排，重排后每个角度下的锥形束投影数据转化为互相平行的多层扇形束投影数据，每层扇形束数据对应被扫描物体沿旋转轴方向的一层；2、扇形束投影到平行束投影的重排，该重排只涉及被扫描物体相同层的扇形束投影，即重排后的多层平行束投影仍对应被扫描物体的相同层；3、三维滤波反投影，即首先对平行束投影数据进行滤波，然后对滤波后的数据进行反投影，从而实现重建。并且，其中的反投影过程是在锥形束几何下进行的。具体实施过程分别如下。

（1）锥形束投影到扇形束投影的重排。如图3所示，将每个角度下的锥形束投影数据转化为互相平行的多层扇形束投影数据，每层扇形束数据对应被扫描物体沿Z轴方向的一层，重排后的扇形束几何与锥形束几何相匹配，即虚拟的扇形束源点同样位于距旋转轴距离为R的位置，线阵探测器与面探测器的行匹配同时对应被扫描物体沿Z轴方向的一层；重排前的锥形束投影数据为b(λ,u,v)，其中λ表示投影角度，并且对于该投影角度λ，可通过

和

分别表示3个轴的方向向量，其中平行探测器的法向量方向，

沿探测器的行方向，沿探测器的列方向；(u,v)表示面探测器上的点，重排后的扇形束投影数据为p_z(θ,u)，其中z表示物体沿Z轴方向的层数，θ表示扇形投影角度，且θ=λ%2π(％表示求余运算），u表示线阵探测器探元位置。重排后的扇形束投影只覆盖物体横截面的一半，同样存在数据截断。

进一步地，该步骤下的重排过程只涉及垂直方向，每一条扇形束射线由一条倾斜的锥形束射线来估计，倾斜的锥形束射线从距离扇形束源点最近的位于扇形束源点正上方或正下方的锥束源点发出，同时穿过扇形束射线与被扫描物体交线的中心M。对于全角度重建，被扫描物体每层对应2π范围内的扇形束投影，即需要2π范围内的锥形束投影来进行估计。所以，锥形束源点到其所生成的虚拟扇形束源点的最大距离为螺距的一半。

该步骤下的重排过程可由下面公式表示：

$p_z(\mathrm{\&theta;},u)=\frac{\sqrt{u^2+S^2}}{\sqrt{u^2+v^2+S^2}}b(\mathrm{\&lambda;},u,v)$

其中

Δz表示锥形束源点与虚拟扇形束源点的距离。

（2）扇形束投影到平行束投影的重排。经过步骤（1）后得到了多层扇形束投影p_z(θ,u)，物体的每一层z均对应2π范围内的扇束投影，但是每个角度下扇形束只覆盖横截面的一半，同样存在数据截断。为了消除数据截断，本发明进一步将扇形束投影重排成平行束投影。重排只涉及被扫描物体相同层z的扇束投影，即重排后的平行束投影对应物体的第z层，定义重排后的平行束投影为g_z(φ,t)，其中φ表示平行束投影角度，t表示旋转中心到每条射线的距离，如图4所示。

重排利用了平行束的对称性质，重排后的平行束投影角度范围为π，但每个角度下投影完全覆盖物体的横截面，不存在数据截断，即在g_z(φ,t)中，φ∈[0,π]，t∈[-R₀,R₀]。具体的重排过程为：

由图4可知

β=α-γ

u=d×tanγ

其中，γ表示带符号的扇角，可以由下面公式计算

$\mathrm{\&gamma;}\left(t\right)=\arccos \left(\frac{t}{R}\right)-\arccos \left(\frac{D_0}{R}\right)$

β同样表示平行束投影角度，β和φ的关系可表示如下

$\mathrm{\&beta;}=\left\{\begin{array}{c}\mathrm{\&phi;},\mathrm{\&beta;}\&Element;[0,\mathrm{\&pi;}]\\ \mathrm{\&phi;}+\mathrm{\&pi;},\mathrm{\&beta;}\&Element;(\mathrm{\&pi;},2\mathrm{\&pi;}]\end{array}\right.$

同时，平行束的对称性可由下面公式表示

g_z(φ,t)=g_z(π+φ,-t)

因此，重排过程可表示为

$g_z(\mathrm{\&phi;},t)=\left\{\begin{array}{c}{p}_z(\mathrm{\&phi;}-\mathrm{\&gamma;},D\tan \left(\mathrm{\&gamma;}\left(t\right)\right)),t<=0\\ p_z(\mathrm{\&pi;}+\mathrm{\&pi;}-\mathrm{\&gamma;},D\tan \left(\mathrm{\&gamma;}(-t)\right)),t>0\end{array}\right.$

（3）三维滤波反投影。经过步骤（2）的重排后，平行束投影完全覆盖物体的横截面，不存在数据截断。在此基础上，本发明进一步采用简单高效的滤波形式，即采用标准的斜坡滤波器分别对被扫描物体每一层的每个投影角度的数据进行滤波，即

${\tilde{g}}_z(\mathrm{\&phi;},t)=g_z(\mathrm{\&phi;},t)*h\left(t\right)$

其中，

表示滤波后的数据，h(t)为斜坡滤波器；然后，对滤波的数据进行反投影即可完成重建，由于步骤（1）中的重排是近似的，且重排误差随锥角的增大而显著变大，因此，本发明的反投影过程是在重排前的锥形束几何而不是平行束几何下进行。该反投影的具体过程为

其中，f(x,y,z)表示三维重建结果，

表示平行束投影角度，表示滤波后的平行束投影数据，

和

表示点(x,y,z)在重排前的锥形束几何而不是最终平行束几何下计算的投影位置。

进一步地，本发明给出了以下仿真实验和实际试验结果，以验证上述方法的有效性。

表1第一组仿真实验参数

仿真实验以三维Shepp-Logan体模为被扫描物体。仿真实验共有三组，三组实验的参数中除了探测器像素矩阵和物体偏移量两个参数不同外，其他参数的数值均相同。第一组仿真实验的参数如表1所示，其中物体的横向尺寸（像素数）为256，探测器的横向尺寸（像素数）为160，因此每个角度下射线只能覆盖被扫描物体的部分，覆盖率为0.625；第二组仿真实验参数中，探测器像素矩阵为200×150，物体偏移量为-2.7mm，对应的射线覆盖率为0.585；第三组仿真实验的探测器像素矩阵为200×140，物体偏移量为-2.9mm，对应的射线覆盖率为0.547。因此，可以看出随着三组仿真实验中探测器的横向尺寸相比物体越来越小，产生数据截断的位置越来越靠近物体中心，截断伪影对有效数据的影响逐步增大。三组仿真实验的重建结果如图5所示，重建结果的x-y中心切片的中心纵向剖线如图6所示。

而在实际试验中，采用的是微焦斑锥束CT系统（射线源YXLON225.48，探测器Varian4030E），其扫描参数为：射线源到旋转轴距离358mm，射线源到探测器距离1438mm，螺距36mm，每个角度下的射线覆盖率为0.521。被成像的物体是增压器叶轮，其重建结果切片如图7所示。

上面结合附图对本发明优选的具体实施方式和实施例作了详细说明，但是本发明并不限于上述实施方式和实施例，在本领域技术人员所具备的知识范围内，还可以在不脱离本发明构思的前提下作出各种变化。

Claims (7)

1.一种针对半覆盖螺旋锥束CT的单层重排滤波反投影重建方法，该方法基于被扫描物体的螺旋锥形束投影进行重建，其特征在于包括以下步骤：

（1）将螺旋锥形束投影重排成多层扇形束投影，重排后每个角度下的锥形束投影数据转化为互相平行的多层扇形束投影数据，每层扇形束数据对应被扫描物体沿旋转轴方向的一层；

（2）将步骤（1）中的多层扇形束投影重排成多层平行束投影，该重排只涉及被扫描物体相同层的扇形束投影，即重排后的多层平行束投影仍对应被扫描物体的相同层；

（3）将步骤（2）中的平行束投影数据进行滤波，然后对滤波后的数据进行反投影，从而实现重建；并且，其中的反投影过程是在锥形束几何下进行的。

2.如权利要求1所述的单层重排滤波反投影重建方法，其特征在于：经步骤（1）重排后的被扫描物体每层对应2π范围内的扇形束投影，且每个角度下对应的扇形束投影只覆盖被扫描物体横截面的一半。

3.如权利要求1所述的单层重排滤波反投影重建方法，其特征在于：所述步骤（2）中的重排利用了平行束投影的对称性质，重排后的平行束投影角度范围为π，且每个角度下的投影完全覆盖被扫描物体的横截面，不存在数据截断。

4.如权利要求1至3中任一项所述的单层重排滤波反投影重建方法，其特征在于：被扫描物体的锥形束投影数据的具体获取过程为：构建半径为R₀的圆柱体，假设被扫描物体即包含于该圆柱体中，在该圆柱体的中心位置建立固定的右手笛卡尔坐标系(x,y,z)，且被扫描物体的偏置量为D₀，成像时，被扫描物体做螺旋轨迹运动，即相当于被扫描物体固定、射线源和探测器同时做反向的螺旋轨迹运动；旋转轴与Z轴重合，射线源的旋转半径为R，射线源到探测器的距离为S，射线源到平行于探测器且过旋转轴的平面的距离为D，建立旋转的坐标系(u,v,w)，其中心为探测器中心，对于探测器上的任何一点，都可以用坐标(u,v)表示；为获得重建所需的足够投影数据，每个角度下的射线至少覆盖被扫描物体横截面的一半，即需要满足：

$\left\{\begin{array}{c}{D}_0<D\tan \left({\mathrm{\&gamma;}}_m\right)\\ (D\tan \left({\mathrm{\&gamma;}}_m\right)+D_0)\cos \left({\mathrm{\&gamma;}}_m\right)>R_0\end{array}\right.$

其中γ_m表示最大的扇角的一半。

5.如权利要求4所述的单层重排滤波反投影重建方法，其特征在于：所述步骤（1）中由锥形束投影到扇形束投影的具体重排过程为：每层扇形束数据对应被扫描物体沿Z轴方向的一层，重排后的扇形束几何与锥形束几何相匹配，即虚拟的扇形束源点同样位于距旋转轴距离为R的位置，线阵探测器与面探测器的行匹配同时对应被扫描物体沿Z轴方向的一层；重排前的锥形束投影数据为b(λ,u,v)，其中λ表示投影角度，(u,v)表示面探测器上的点，重排后的扇形束投影数据为p_z(θ,u)，其中z表示物体沿Z轴方向的层数，θ表示扇形投影角度，且θ=λ%2π，u表示线阵探测器探元位置，则该步骤下的重排过程可由下面公式表示：

$p_z(\mathrm{\&theta;},u)=\frac{\sqrt{u^2+S^2}}{\sqrt{u^2+v^2+S^2}}b(\mathrm{\&lambda;},u,v)$

其中

Δz表示锥形束源点与虚拟扇形束源点的距离。

6.如权利要求5所述的单层重排滤波反投影重建方法，其特征在于：所述步骤（2）中由扇形束投影到平行束投影的具体重排过程为：

$g_z(\mathrm{\&phi;},t)=\left\{\begin{array}{c}{p}_z(\mathrm{\&phi;}-\mathrm{\&gamma;},D\tan \left(\mathrm{\&gamma;}\left(t\right)\right)),t<=0\\ p_z(\mathrm{\&pi;}+\mathrm{\&pi;}-\mathrm{\&gamma;},D\tan \left(\mathrm{\&gamma;}(-t)\right)),t>0\end{array}\right.$

其中，g_z(φ,t)表示重排后的平行束投影，φ表示平行束投影角度，t表示旋转中心到每条射线的距离，并且φ∈[0,π]，t∈[-R₀,R₀]，γ表示带符号的扇角。

7.如权利要求6所述的单层重排滤波反投影重建方法，其特征在于：在所述步骤（3）中，采用标准的斜坡滤波器分别对被扫描物体每一层的每个投影角度的数据进行滤波，即

${\tilde{g}}_z(\mathrm{\&phi;},t)=g_z(\mathrm{\&phi;},t)*h\left(t\right)$

其中，

表示滤波后的数据，h(t)为斜坡滤波器；并且，所述反投影的具体过程为

其中，f(x,y,z)表示三维重建结果，

表示平行束投影角度，

表示滤波后的平行束投影数据，

和

表示点(x,y,z)在进行步骤（1）重排前的锥形束几何下计算的投影位置。

Images