CN102973291B - C型臂半精确滤波反投影断层成像方法 - Google Patents

Abstract

本发明公开了一种适合于C型臂成像设备的断层成像方法，其特点是该方法具有移不变特征，对圆弧轨迹扫描的锥束投影数据进行半精确重建。圆弧轨迹是短扫描轨迹，即π加扇角扫描轨迹；也是超短扫描轨迹，即小于短扫描的轨迹。此方法包括：权重投影数据的偏导数；沿水平与非水平方向的一维希尔伯特变换；圆弧轨迹的反投影。相对FDK类型近似重建算法，本发明提出的方法具有重建精度、离散误差小等优点。

Description

C型臂半精确滤波反投影断层成像方法

技术领域

本发明涉及一种C型臂半精确滤波反投影断层成像方法，属于生物医学成像领域，适用于C型臂X射线的断层成像。

背景技术

C型臂成像系统主要由C型臂及安装C型臂两端的X射线源和平板检测器组成。X射线源发射的X射线光信号在空间中形成锥束，C型臂可绕水平轴线在有限空间内进行最大角度范围的旋转，与此同时，平板检测器能可采集X射线衰减后的光信号。这些光信号又能转换为锥束投影数据，利用这些数据就能重建出CT断层图像。具备CT断层功能的C型臂在临床诊断和手术评估具有重要的应用价值，是一项具有广泛应用前景的重要技术。

利用C型臂进行断层重建，最实用的扫描轨迹是圆弧轨迹。目前有FDK类型和Katsevich类型两类算法可以实现相应的重建。前者是借鉴Feldmap等人在Journal of the Optical Society of America上发表的Practical cone beam algorithm论文的思想，将二维扇束重建公式向锥束投影空间的推广，属于近似锥束重建算法。后者是Guang-Hong Chen等人利用Katsevich在International Jounal of Mathematics and Mathematical Sciences 上发表的A generalscheme for constructing inversion algorithm for conebeam CT论文的思想，将精确锥束重建算法推广到圆弧轨迹，并申请美国专利cone beam filtered backprojection image reconstruction method for short trajectories，其专利号为US7,203,272 B2，属于半精度的锥束重建算法。

虽然FDK类型算法重建思想简单，速度较快，全圆轨迹的FDK算法在较小的锥角(＜±4°)下，能取得良好的重建效果，然而对于圆弧轨迹的FDK类型算法重建密度差较大的物体时在锥角±2°左右就出现了伪像。Katsevich类型重建算法的关键是计算临界面的构造因子。Guang-Hong Chen等人利用通过重建点与圆弧轨迹的端点且与轨迹平面垂直的平面将圆弧轨迹分成三段，再利用射源所在轨迹区段位置来判别构造因子的大小。因不同的重建点会将圆弧轨迹分成不同的三个区段，从而不同的投影角度以及不同的重建点都需要判别构造因子的大小，并且相应的计算关系还需要一些平方、开方以及三角函数的运算，是一个相对复杂计算判别过程。因此，尽管Katsevich类型的算法可获得较精确的重建结果，结构因子的复杂计算过程不利于算法的并行处理，从而降低了重建算法的计算效率。此外，在某些投影角度上Katsevich类型重建算法存在沿竖直或者接近竖直方向的滤波线，而Guang-Hong Chen等人提出的算法并未对投影数据沿检测器竖直方向进行重采样，影响了重建结果的精度。以上两点不利于Katsevich类型锥束重建算法在实际中的应用。

发明内容

本发明的目的是针对现有技术的不足而提供一种C型臂半精确滤波反投影断层成像方法，其特点是将Katsevich精确重建算法应用到圆弧轨迹，利用同一条投影射线上的重建点构造因子相同，将结构因子相关的计算映射到检测器平面上，获得构造因子计算规则，进而得到具有滤波反投影的重建公式。避免了现有Katsevich重建算法构造因子的复杂计算过程，提高重建算法的速度；同时也避免了FDK类型算法重建精度不高的缺点，为C型臂断层成像提供方法支撑。

本发明的目的由以下技术措施实现：

C型臂半精确滤波反投影断层成像方法包括以下步骤：检测器上的阵元接受X射线源发射的X射线衰减信号；X射线源和检测器绕水平轴线旋转，数据采集系统将阵元获取的X射线衰减信号转换为锥束投影数据；图像重建系统接收锥束投影数据，输入锥束投影数据预处理模块进行预处理；预处理后的投影数据并行输入到偏导数模块Ⅰ、偏导数模块Ⅱ和偏导数模块Ⅲ计算关于不同参数的偏导数；偏导加权模块并行接收相应的偏导数据，进行加权求和，获得投影数据的偏导数g＇(λ,u,v)；投影数据的偏导数并行输入到滤波模块Ⅰ、滤波模块Ⅱ和滤波模块Ⅲ分别计算沿三类滤波线的滤波；滤波加权反投影模块并行接收三类滤波数据，进行权重求和，产生反投影数据，再进行沿圆弧轨迹的反投影操作，获得断层图像。

所述数据采集系统是把检测器阵元接收的X射线模拟信号转换为数字信号，进而转化成锥束投影数据。

所述射源发射X射线，并绕水平轴线旋转，产生圆弧扫描轨迹，包括短扫描和超短扫描轨迹。

所述检测器是指二维平板检测器，由若干个阵元组成，每个阵元可感光X射线。

所述的沿三类滤波线是指：与圆弧轨迹相切的临界面与检测器平面的交线；过圆弧起始或终止端点的临界面与检测器平面的交线。

所述的沿三类滤波线的滤波是指：偏导数据沿与圆弧轨迹相切的临界面与检测器平面的交线进行的希尔伯特滤波；偏导数据沿过圆弧起始端点的临界面与检测器平面的交线进行的希尔伯特滤波；偏导数据沿过圆弧终止端点的临界面与检测器平面的交线进行的希尔伯特滤波。

在进行偏导数据沿过圆弧起始端点的临界面与检测器平面的交线进行的希尔伯特滤波前，需过端点在检测器上的投影点作斜率为-1的直线，将位于该直线下方的偏导数据g＇(λ,u,v)乘以-1，获得修正后偏导数据g＇_s(λ,u,v)，g＇_s(λ,u,v)沿纵向或横向坐标方向进行重采样获得横向修正后偏导数据g＇_sh(λ,u,v)或纵向修正后偏导数据为g＇_sv(λ,u,v)，再进行希尔伯特滤波，得到，其值为： $g_s^F(\mathrm{\λ},u,v)=k_1{\∫}_{-\∞}^{\∞}\frac{{\mathrm{Dg}}_{\mathrm{sh}}^{\′}(\mathrm{\λ},\tilde{u},\tilde{v})}{\sqrt{D^2+{\tilde{u}}^2+{\tilde{v}}^2}}{h}_H(u-\tilde{u})d\tilde{u}+k_2{\∫}_{-\∞}^{\∞}\frac{{\mathrm{Dg}}_{\mathrm{sv}}^{\′}(\mathrm{\λ},\tilde{u},\tilde{v})}{\sqrt{D^2+{\tilde{u}}^2+{\tilde{v}}^2}}{h}_H(v-\tilde{v})d\tilde{v}$ ，其中：k₁和k₂取值只能是0或1且k₁和k₂异或值为1。

在进行偏导数据沿过圆弧终止端点的临界面与检测器平面的交线进行的希尔伯特滤波前，需过端点在检测器上的投影点作斜率为-1的直线，将位于该直线下方的偏导数据g＇(λ,u,v)乘以-1，获得修正后偏导数据g＇_e(λ,u,v)，g＇_e(λ,u,v)沿纵向或横向坐标方向进行重采样获得横向修正后偏导数据g＇_eh(λ,u,v)或纵向修正后偏导数据为g＇_ev(λ,u,v)，再进行希尔伯特滤波，得到，其值为： $g_e^F(\mathrm{\λ},u,v)=k_1{\∫}_{-\∞}^{\∞}\frac{{\mathrm{Dg}}_{\mathrm{eh}}^{\′}(\mathrm{\λ},\tilde{u},\tilde{v})}{\sqrt{D^2+{\tilde{u}}^2+{\tilde{v}}^2}}{h}_H(u-\tilde{u})d\tilde{u}+k_2{\∫}_{-\∞}^{\∞}\frac{{\mathrm{Dg}}_{\mathrm{ev}}^{\′}(\mathrm{\λ},\tilde{u},\tilde{v})}{\sqrt{D^2+{\tilde{u}}^2+{\tilde{v}}^2}}{h}_H(v-\tilde{v})d\tilde{v}$ ，其中：k₁和k₂取值只能是0或1且k₁和k₂异或值为1。

所述的三种滤波数据的加权是指：三类滤波数据、和按照公式进行加权，获得反投影值。

与现有技术相比，本发明具有如下优点：

本发明适合C型臂断层成像方法，该技术方案保持了滤波反投影的框架，计算速度较快，可以有效改善锥角伪像，提高重建精度范围；可适应短扫描以及超短扫描轨迹锥束投影数据的断层重建；投影数据轴向截断对重建结果影响不大，对投影数据的噪声不敏感。

附图说明

图1为C型臂成像设备

1.C型臂，2.X射线源，3.X射线，4.平面检测器，5.阵元, 6.手术台, 7病人, 8.基座, 9.水平轴线。

图2为C型臂重建系统框图

10.控制系统，11.采集系统，12.运动控制器，13.X射线控制器，14.手术台控制器，15.断层图像重建单元, 16计算机，17.存储设备，18.接收操作控制台, 19.观察显示器。

图3为C型臂成像系统的几何参数图

图4为几何参数在检测器上的投影示意图

图5为圆弧轨迹构造因子的确定示意图

(a)与轨迹相切的临界面

(b)过圆弧端点的临界面

图6为检测器区域的划分

(a)端点投影在检测器外侧且滤波线最大夹角小于π/2

(b)端点投影在检测器外侧且滤波线最大夹角大于π/2

(c)端点投影在检测器内

图7为本发明的半精确滤波反投影重建方法流程图

20.锥束投影数据预处理模块，21.偏导数模块，22.偏导数模块，23.偏导数模块，24.偏导数加权模块，25.滤波模块，26.滤波模块，27.滤波模块，28.滤波加权反投影模块，29.层析图像。

具体实施方式

以下通过实施例对本发明进行具体的描述，有必要在此指出的是本实施例只用于对本发明进行进一步说明，不能理解为对发明保护范围的限制，该领域的技术熟练人员可以根据上述本发明的内容作出的一些非本质的改进和调整。

下面对照附图详细描述本发明的实施例。

图1示出了C型臂成像设备，该设备包括型C臂1，C型臂1的一端配有X射线源2，用于产生X射线3，X射线3在空间中形成锥束发射；在C型臂1另一端配置二维平面检测器4，检测器由一些阵元5组成，每个阵元5可感光穿过躺在手术台6上病人7的X射线3衰减信号；C型臂1被安装在基座8上，可绕水平轴线9旋转，旋转的最大角度范围240度；随C型臂1自动旋转的同时，X射线源2形成了沿圆弧的扫描轨迹，检测器4上的阵元5可获取不同扫描位置的X射线衰减信号。

图2示出了成像系统中的控制系统10，包括：数据采集系统11，将阵元5获取的X射线2衰减信号转换为锥束投影数据；C臂运动控制器12，控制C臂运动的速度和位置；X射线控制器13，提供X射线2的能量和时间信号；手术台控制器14，控制手术台6的平移和升降运动；断层图像重建单元15，接收来自于数据采集系统11的X射线2投影数据，依据目前发明的方法对数字化后的数据进行断层成像；重建后的图像被输入到计算机内16，并进一步存储到大容量存储设备17；计算机16接收操作控制台18的输入的参数或指令，允许操作者观察显示器19重建图形和其它数据；观察者也可利用计算机16为数据采集系统11、C臂运动控制器12、X射线控制器13以及手术台控制器14提供控制指令和相应的参数。

图3示出了C型臂成像系统的几何参数图。O代表旋转中心，以其为原点建立笛卡尔xyz坐标系；投影采集过程中，检测器和射源绕z轴旋转，检测器平面始终垂直于xoy平面，检测器边长为210mm*210mm，其分辨率为1024*1024，射源与原点O的连线始终垂直于检测器中心；射源到检测器中心的距离记为D，其值为1000mm；射源起始点位于x轴的正向，其运动轨迹可描述为，其中，λ∈(λ_s,λ_e)表示射源的转角，而λ_s和λ_e分别表示圆弧起始点s和终止点e对应的转角值，两者的差值最大值为240度；重建体的密度函数记为，其中，的支撑区满足x²+y²≤R²,|z|＜H。对于锥束投影采集系统，在每个投影视角引入由轴标识的旋转坐标系，其中，和轴分别平行于检测器的水平和竖直阵列方向，轴是射源指向原点O的方向，则各坐标轴的单位矢量方向为：

$\left\{\begin{array}{c}{\stackrel{\→}{e}}_u\left(\mathrm{\λ}\right)=(-\sin \mathrm{\λ},\cos \mathrm{\λ},0)\\ {\stackrel{\→}{e}}_v=\left(\mathrm{0,0,1}\right)\\ {\stackrel{\→}{e}}_w\left(\mathrm{\λ}\right)=(-\cos \mathrm{\λ},-\sin \mathrm{\λ},0)\end{array}\right.\mathrm{\λ}\∈({\mathrm{\λ}}_s,{\mathrm{\λ}}_e)$

由此，射源和重建点确定的单位矢量方向可表达为：

则过射源沿单位矢量方向的锥束投影可表示为：

$g(\mathrm{\λ},\stackrel{\→}{\mathrm{\β}})={\∫}_0^{\∞}f(\stackrel{\→}{\mathrm{\α}}\left(\mathrm{\λ}\right)+t\stackrel{\→}{\mathrm{\β}})\mathrm{dt},\stackrel{\→}{\mathrm{\β}}\∈S^2$

上式中的S²表示三维空间的单位球。

过源点且包含单位矢量的平面记为，其单位法向矢量记为，而θ表示平面在垂直于矢量局部坐标系中的极角。

Katsevich在A general scheme for constructing inversionalgorithm for cone beam CT论文中推导出通用的具有移不变滤波反投影结构的精确锥束重建算法，其重建公式可表示为：

$f\left(\stackrel{\→}{x}\right)=-\frac{1}{4{\mathrm{\π}}^2}{\∫}_{I\left(\stackrel{\→}{x}\right)}\mathrm{d\λ}\underset{m}{\mathrm{\Σ}}\frac{c_m(\mathrm{\λ},\stackrel{\→}{x})}{|\stackrel{\→}{x}-\stackrel{\→}{\mathrm{\α}}\left(\mathrm{\λ}\right)|}{\∫}_0^{2\mathrm{\π}}\frac{\mathrm{d\γ}}{\sin \mathrm{\γ}}\frac{\∂}{\∂q}g(q,\mathrm{\Θ}(\mathrm{\λ},\stackrel{\→}{x},\mathrm{\γ})){|}_{q=\mathrm{\λ}}$

式中，Θ表示平面上任意投影射线矢量方向，γ则为矢量方向Θ与之间的夹角；称之为构造因子，其表达式为：

$c_m(\mathrm{\λ},\stackrel{\→}{x})=\mathrm{sgn}(\stackrel{\→}{n}\left(\mathrm{\θ}\right)\·{\stackrel{\→}{\mathrm{\α}}}^{\′}\left(\mathrm{\λ}\right))w(\mathrm{\λ},\stackrel{\→}{x},\stackrel{\→}{n}\left(\mathrm{\θ}\right)){|}_{\mathrm{\θ}={\mathrm{\θ}}_m-\mathrm{\ϵ}}^{\mathrm{\θ}={\mathrm{\θ}}_m+\mathrm{\ϵ}}$

构造因子中，表示平面的权重值。由于过重建点的平面与轨迹的交点可能不止一个，则权重值反映了平面对重建点的冗余贡献，需要满足下面关系：

$\underset{i}{\mathrm{\Σ}}w(\mathrm{\λ},\stackrel{\→}{x},\stackrel{\→}{n}\left(\mathrm{\θ}\right))=1$

在本发明中选择等权重策略，即满足下式：

$w(\mathrm{\λ},\stackrel{\→}{x},\stackrel{\→}{n}\left(\mathrm{\θ}\right))=\frac{1}{N}$

其中，N代表了平面与扫描轨迹的交点个数。

从构造因子公式中，可以看出平面绕矢量分别顺时针和逆时针作无限小ε角度旋转，两者之间的符号权重差值即为值。Katsevich指出重建过程仅发生构造因子不为零的，而临界面是与轨迹相切或者通过轨迹端点的平面。对于圆弧轨迹上存在三类临界面：与轨迹相切以及分别通过圆弧起始和终止端点的平面。为此需判别这三类临界面上构造因子的大小来确定投影数据的滤波方向。

图4示出了当源点为时，锥束扫描几何参数在检测器平面上的投影，上标符号^分别表示相应参数在检测器上的投影，而∏_t、∏_s和∏_e分别标识了过且与轨迹相切的平面、过和圆弧起始端点s的平面和过和圆弧终止端点e的平面与检测器的交线。显然，平面法线矢量在检测器上的投影与平面同检测器的交线相垂直。此外，根据构造因子的含义可知，同一条投影射线上的重建点具有相同的构造因子，并且这些重建点在检测器上的投影地址相同，因此，可以将构造因子存储在检测器对应的阵元上。

图5（a）示出了与圆弧轨迹相切的临界面同检测器平面的交线是沿水平方向的。当该平面绕投影射线分别沿顺时针与逆时针旋转后，与轨迹的交点个数均为2，因此，；而的符号值发生改变，由负变为正，由此可计算出该临界面上的结构因子值为1，并且相应的滤波线是沿水平方向。

图5（b）示出了过圆弧端点的临界面相关参数在检测器上的投影。从图中可看出，重建点在检测器上的投影在检测器轴上坐标正负的不同，临界面平面穿过圆弧端点前后与轨迹的交点个数变化不同；此外，投影点位置的不同，临界面上的的符号不同。因此，计算过圆弧端点的临界面上构造因子相对要复杂些。

过圆弧端点的临界面与检测器的交线决定了投影数据的滤波方向，亦称为滤波线。端点的投影一定在轴上，根据它沿轴的坐标的位置，将滤波线在检测器上的分布划分为三种情况，如图6所示。图6(a)示出了端点投影在检测器外，最外侧两条滤波线的夹角小于等于π/2；图6 (b)也示出了端点投影仍然在检测器外，但最外侧两条滤波线夹角小于等于π且大于π/2；图6 (c) 示出了端点投影在检测器内，滤波线通过端点投影并沿整周分布。

检测器上的投影数据需要沿这些滤波线进行重采样，而且重采样后的投影数据被权重及滤波后，还需要反采样到检测器阵列单元上以方便反投影运算。为了减少采样过程中数据的缺失，对于斜率绝对值小于等于1的滤波线，投影数据应保持检测器水平阵列线的索引不变，而沿着轴进行重采样；对于斜率绝对值大于1的滤波线，投影数据应保持检测器垂直阵列线的索引不变，而沿着轴进行重采样。

研究上述三种情况构造因子的取值，得出过圆弧端点处构造因子的绝对值均为0.5，而正负取值仅与过圆弧端点在检测器上的投影的斜率为-1的直线相关。对于圆弧的起始端点，则该直线的上方对应的构造因子为0.5，而在该直线下方对应的构造因子为-0.5；而对于圆弧的终止端点，则在直线的上方对应的构造因子为-0.5，而对应在该直线下方的构造因子为0.5。

确定了构造因子不为零的Radon平面以及相应的数值后，还需要推导适合平面检测器上的重建公式。在Katsevich重建公式中，对轨迹参数λ的偏导数相关部分为：

$g^{\′}(\mathrm{\λ},\mathrm{\Θ}(\mathrm{\λ},\stackrel{\→}{x},\mathrm{\γ}))=\frac{\∂}{\∂q}g(q,\mathrm{\Θ}(\mathrm{\λ},\stackrel{\→}{x},\mathrm{\γ})){|}_{q=\mathrm{\λ}}$

上式中，当采用平板检测器时，投影数据g(λ,Θ)又可表达为g(λ,u,v)，其中，u,v分别表示在检测器平面上沿和轴的坐标。通过和，利用偏导数链式法则，可推导出投影数据的偏导数为：

$g^{\′}(\mathrm{\λ},u,v)=(\frac{\∂}{\∂\mathrm{\λ}}+\frac{D^2+u^2}{D}\frac{\∂}{\∂u}+\frac{\mathrm{uv}}{D}\frac{\∂}{\∂v})g(\mathrm{\λ},u,v)$

构造因子不为零的Radon平面，其上相关的运算才会对重建点有所贡献。实际上，它与检测器平面的交线确定了投影数据的滤波方向。Katsevich重建公式中与滤波运算相关的表达式是：

${\∫}_0^{2\mathrm{\π}}\frac{\∂}{\∂q}g(q,\mathrm{\Θ}(\mathrm{\λ},\stackrel{\→}{x},\mathrm{\γ})){|}_{q=\mathrm{\λ}}\frac{\mathrm{d\γ}}{\sin \mathrm{\γ}}$

对于斜率绝对值小于1的滤波线，可以推导出：

$\frac{\mathrm{d\γ}}{\sin \mathrm{\γ}}\≈\frac{\mathrm{du}}{(\tilde{u}-u)}\frac{\sqrt{D^2+u^2+v^2}}{\sqrt{D^2+{\tilde{u}}^2+{\tilde{v}}^2}}$

利用Hilbert变换，进一步推导出投影数据的滤波公式：

${\∫}_{-\∞}^{\∞}\frac{\mathrm{\π}\sqrt{D^2+u^2+v^2}}{\sqrt{D^2+{\tilde{u}}^2+{\tilde{v}}^2}}{g}^{\′}(\mathrm{\λ},\tilde{u},\tilde{v})h_H(\tilde{u}-u)d\tilde{u}$

其中，h_H(·)表示Hilbert变换滤波核。记g^F(λ,u,v)为投影数据偏导数的滤波函数，其表达式为：

$g^F(\mathrm{\λ},u,v)={\∫}_{-\∞}^{\∞}\frac{{\mathrm{Dg}}^{\′}(\mathrm{\λ},\tilde{u},\tilde{v})}{\sqrt{D^2+{\tilde{u}}^2+{\tilde{v}}^2}}{h}_H(u-\tilde{u})d\tilde{u}$

在检测器平面上，由于通过圆弧端点的滤波线存在斜率绝对值大于1的情况。投影数据需要按轴进行重排，并也需要推导相应的公式。依据上述相同的推导方式，同样获得斜率绝对值大于1的滤波线对应的滤波函数公式g^F(λ,u,v)：

$g^F(\mathrm{\λ},u,v)={\∫}_{-\∞}^{\∞}\frac{{\mathrm{Dg}}^{\′}(\mathrm{\λ},\tilde{u},\tilde{v})}{\sqrt{D^2+{\tilde{u}}^2+{\tilde{v}}^2}}{h}_H(v-\tilde{v})d\tilde{v}$

由此，Katsevich类型适合圆弧轨迹的重建公式：

$f\left(\stackrel{\→}{x}\right)=\frac{1}{4\mathrm{\π}}\underset{j=0}{\overset{2}{\mathrm{\Σ}}}{c}_j(\mathrm{\λ},\stackrel{\→}{x}){\∫}_{I_j\left(\stackrel{\→}{x}\right)}\frac{g^F(\mathrm{\λ},u,v)}{D+\stackrel{\→}{x}\·{\stackrel{\→}{e}}_w}\mathrm{d\λ}$

其中，c₀=1表示投影数据沿水平方向进行滤波时构造因子的取值；c₁和c₂分别表示投影数据沿过圆弧起始和终止端点在检测器上的投影点和的直线方向进行滤波时构造因子的取值，依据在检测器上相应的分布规律取0.5或-0.5。

根据上述的分析，重建过程主要可划分为投影数据的偏导数、滤波以及反投影三个步骤。由此，图7示出了本发明适合C型臂系统的半精确滤波反投影方法的流程图。数据采集系统11将检测器4上的阵元5采集到的X射线3衰减信号转换为锥束投影数据并输入到断层图像重建单元15。

在图像重建单元15中，锥束投影数据经过锥束投影预处理模块20进行诸如偏置校正、增益校正、坏像素校正及log运算等预处理。获取适合断层重建的投影数据，并按照检测器阵元布置的二维数组形式进行存储。

预处理后的数据并行进入偏导数模块Ⅰ21、偏导数模块Ⅱ22和偏导数模块Ⅲ23进行偏导数的运算。偏导数模块Ⅰ21进行关于轨迹参数λ的偏导数， 其运算公式为：

$\frac{\∂}{\∂\mathrm{\λ}}g(\mathrm{\λ},u,v)=\frac{g(\mathrm{\λ}+\mathrm{\Δ\λ},u,v)-g(\mathrm{\λ}-\mathrm{\Δ\λ},u,v)}{\mathrm{\Δ\λ}}$

偏导数模块Ⅱ22进行关于检测器局部坐标的横向坐标参数u的偏导数，其运算公式为：

$\frac{\∂}{\∂u}g(\mathrm{\λ},u,v)=\frac{g(\mathrm{\λ},u+\mathrm{\Δu},v)-g(\mathrm{\λ},u-\mathrm{\Δu},v)}{\mathrm{\Δu}}$

偏导数模块Ⅲ23进行关于检测器局部坐标的横向坐标参数v的偏导数，其运算公式为：

$\frac{\∂}{\∂v}g(\mathrm{\λ},u,v)=\frac{g(\mathrm{\λ},u,v+\mathrm{\Δv})-g(\mathrm{\λ},u,v-\mathrm{\Δv})}{\mathrm{\Δv}}$

偏导数加权模块24接收以上三个偏导数据进行加权求和，获得投影数据的偏导数，相应的公式为：

$g^{\′}(\mathrm{\λ},u,v)=(\frac{\∂}{\∂\mathrm{\λ}}+\frac{D^2+u^2}{D}\frac{\∂}{\∂u}+\frac{\mathrm{uv}}{D}\frac{\∂}{\∂v})g(\mathrm{\λ},u,v)$

投影数据的偏导数并行进入滤波模块Ⅰ25，滤波模块Ⅱ26和滤波模块Ⅲ27，进行相关滤波运算。滤波模块Ⅰ25是计算与轨迹相切的临界面上相关数据的滤波，滤波线是沿检测器横向阵列方向，滤波核是希尔伯特滤波核，相应的公式为：

$g_t^F(\mathrm{\λ},u,v)={\∫}_{-\∞}^{\∞}\frac{{\mathrm{Dg}}^{\′}(\mathrm{\λ},u,v)}{\sqrt{D^2+{\tilde{u}}^2+{\tilde{v}}^2}}{h}_H(u-\tilde{u})d\tilde{u}$

滤波模块Ⅱ26是计算过圆弧起始端点的临界面上相关数据的滤波，滤波线是沿在检测器平面内过圆弧起始端点在检测器平面上投影点的任意直线，滤波核是希尔伯特滤波核。在进行滤波前，过投影点作斜率为-1的直线，将位于该直线下方的投影数据的偏导数乘以-1，记修正后偏导数据为g＇_s(λ,u,v)。此外，由于可能存在沿检测器纵向采样的滤波线，如图6（b）和（c）所示，g＇_s(λ,u,v)进一步划分横向修正后偏导数据g＇_sh(λ,u,v)和纵向修正后偏导数据为g＇_sv(λ,u,v)。相应的滤波公式为：

$g_s^F(\mathrm{\λ},u,v)=k_1{\∫}_{-\∞}^{\∞}\frac{{\mathrm{Dg}}_{\mathrm{sh}}^{\′}(\mathrm{\λ},\tilde{u},\tilde{v})}{\sqrt{D^2+{\tilde{u}}^2+{\tilde{v}}^2}}{h}_H(u-\tilde{u})d\tilde{u}+k_2{\∫}_{-\∞}^{\∞}\frac{{\mathrm{Dg}}_{\mathrm{sv}}^{\′}(\mathrm{\λ},\tilde{u},\tilde{v})}{\sqrt{D^2+{\tilde{u}}^2+{\tilde{v}}^2}}{h}_H(v-\tilde{v})d\tilde{v}$

k₁和k₂取值只能是0或1。当偏导数据沿检测器横向方向采样，k₁取值为1，k₂取值为0；当偏导数据沿检测器纵向方向采样，k₁取值为0，k₂取值为1。

滤波模块Ⅲ27是计算过圆弧终止端点的临界面上相关数据的滤波，滤波线是沿在检测器平面内过圆弧终止端点在检测器平面上投影点的任意直线，滤波核是希尔伯特滤波核。在进行滤波前，过投影点作斜率为-1的直线，将位于该直线下方的投影数据的偏导数乘以-1，记修正后偏导数据为g＇_e(λ,u,v)。此外，由于可能存在沿检测器纵向采样的滤波线，如图6（b）和（c）所示，g＇_e(λ,u,v)进一步划分横向修正后偏导数据g＇_eh(λ,u,v)和纵向修正后偏导数据为g＇_ev(λ,u,v)。相应的滤波公式为：

$g_e^F(\mathrm{\λ},u,v)=k_1{\∫}_{-\∞}^{\∞}\frac{{\mathrm{Dg}}_{\mathrm{eh}}^{\′}(\mathrm{\λ},\tilde{u},\tilde{v})}{\sqrt{D^2+{\tilde{u}}^2+{\tilde{v}}^2}}{h}_H(u-\tilde{u})d\tilde{u}+k_2{\∫}_{-\∞}^{\∞}\frac{{\mathrm{Dg}}_{\mathrm{ev}}^{\′}(\mathrm{\λ},\tilde{u},\tilde{v})}{\sqrt{D^2+{\tilde{u}}^2+{\tilde{v}}^2}}{h}_H(v-\tilde{v})d\tilde{v}$

k₁和k₂取值只能是0或1。当偏导数据沿检测器横向方向采样，k₁取值为1，k₂取值为0；当偏导数据沿检测器纵向方向采样，k₁取值为0，k₂取值为1。

获得的滤波数据，并行输入到滤波加权反投影模块28进行加权求和，求和公式为：

$g_t^F(\mathrm{\λ},u,v)+0.5g_s^F(\mathrm{\λ},u,v)-0.5g_e^F(\mathrm{\λ},u,v)$

获得在检测坐标u,v处的反投影值，再按照u,v与重建点满足的投影映射关系：

$u=D\frac{\stackrel{\→}{x}\·{\stackrel{\→}{e}}_u}{D+\stackrel{\→}{x}\·{\stackrel{\→}{e}}_w},v=D\frac{\stackrel{\→}{x}\·{\stackrel{\→}{e}}_v}{D+\stackrel{\→}{x}\·{\stackrel{\→}{e}}_w}$

进行反投影运算。

在X射线源2每个采样点进行上述同样的步骤，即可所获得断层图像29，这些图像经过计算机16存到大容量存储器17中并输入到显示器19进行显示和相关操作。

Claims (7)

1.一种C型臂半精确滤波反投影断层成像方法，其特征在于成像方法的步骤包括：检测器(4)上的阵元(5)接受X射线源(2)发射的X射线(3)衰减信号；X射线源(2)和检测器(4)绕水平轴线(9)旋转，数据采集系统(11)将阵元(5)获取的X射线(3)衰减信号转换为锥束投影数据；图像重建系统(15)接收锥束投影数据，输入锥束投影数据预处理模块(20)进行预处理；预处理后的投影数据并行输入到偏导数模块Ⅰ(21)、偏导数模块Ⅱ(22)和偏导数模块Ⅲ(23)计算关于不同参数的偏导数；偏导加权模块(24)并行接收相应的偏导数据，进行加权求和，获得投影数据的偏导数g'(λ,u,v)；投影数据的偏导数并行输入到滤波模块Ⅰ(25)、滤波模块Ⅱ(26)和滤波模块Ⅲ(27)分别计算沿三类滤波线的滤波；滤波加权反投影模块(28)并行接收三类滤波数据，进行权重求和，产生反投影数据，再进行沿圆弧轨迹的反投影操作，获得断层图像(29)； 

其中，在进行偏导数据沿过圆弧起始端点的临界面与检测器平面的交线进行的希尔伯特滤波前，需过端点在检测器上的投影点作斜率为-1的直线，将位于该直线下方的偏导数据g'(λ,u,v)乘以-1，获得修正后偏导数据g'_s(λ,u,v)，g'_s(λ,u,v)沿纵向或横向坐标方向进行重采样获得横向修正后偏导数据g'_sh(λ,u,v)或纵向修正后偏导数据为g'_sv(λ,u,v)，再进行希尔伯特滤波，得到其值为： 

其中：k₁和k₂取值只能是0或1且k₁和k₂异或值为1； 

在于在进行偏导数据沿过圆弧终止端点的临界面与检测器平面的交线进行的希尔伯特滤波前，需过端点在检测器上的投影点作斜率为-1的直线，将位于该直线下方的偏导数据g'(λ,u,v)乘以-1，获得修正后偏导数据g'_e(λ,u,v)，g'_e(λ,u,v)沿纵向或横向坐标方向进行重采样获得横向修正后偏导数据g'_eh(λ,u,v)或纵向修正后偏导数据为g'_ev(λ,u,v)，再进行希尔伯特滤波，得到其值为： 

其中：k₁和k₂取值只能是0或1且k₁和k₂异或值为1。 

2.根据权利要求1所述C型臂半精确滤波反投影断层成像方法，其特征在于所述数据采集系统(11)是把检测器阵元接收的X射线模拟信号转换为数字信号，进而转化成锥束投影数据。 

3.根据权利要求1所述C型臂半精确滤波反投影断层成像方法，其特征在于所述X射线射源发射X射线(3)，并绕水平轴线旋转，产生圆弧扫描轨迹，包括短扫描和超短扫描轨迹。 

4.根据权利要求1所述C型臂半精确滤波反投影断层成像方法，其特征在于所述检测器(4)是指二维平板检测器，由若干个阵元组成，每个阵元可感光X射线。 

5.根据权利要求1所述的C型臂半精确滤波反投影断层成像方法，其特征在于所述的沿三类滤波线是指：与圆弧轨迹相切的临界面与检测器平面的交线；过圆弧起始或终止端点的临界面与检测器平面的交线。 

6.根据权利要求1所述的C型臂半精确滤波反投影断层成像方法，其特征在于所述的沿三类滤波线的滤波是指：偏导数据沿与圆弧轨迹相切的临界面与检测器平面的交线进行的希尔伯特滤波偏导数据沿过圆弧起始端点的临界面与检测器平面的交线进行的希尔伯特滤波偏导数据沿过圆弧终止端点的临界面与检测器平面的交线进行的希尔伯特滤波

7.根据权利要求1所述的C型臂半精确滤波反投影断层成像方法，其特征在于所述的三类滤波数据的加权是指：三类滤波数据和 按照公式 进行加权，获得反投影值。