CN109741434B

CN109741434B - 一种基于体积驱动的投影方法及其系统

Info

Publication number: CN109741434B
Application number: CN201811627387.8A
Authority: CN
Inventors: 徐丹; 曾凯; 吴小页
Original assignee: Shenzhen Anke High Tech Co ltd
Current assignee: Shenzhen Anke High Tech Co ltd
Priority date: 2018-12-28
Filing date: 2018-12-28
Publication date: 2022-11-01
Anticipated expiration: 2038-12-28
Also published as: CN109741434A

Abstract

本发明公开了一种基于体积驱动的投影方法及其系统，所述方法包括步骤：对同一排像素依次计算像素在探测器上的二维投影；将二维投影扩展到三维投影。由于在二维情况下基于面积的快速正投/反投算法，并且将该算法扩展到三维情况，成为基于体积的算法。不仅大幅提高精度，而且射线与像素的相交面积时，同一排像素使用顺序计算的办法，因而大量中间计算变量可以被缓存，因而提高计算速度。

Description

一种基于体积驱动的投影方法及其系统

技术领域

本发明涉及医学成像技术领域，尤其涉及的是一种基于体积驱动的投影方法及其系统。

背景技术

一个完整的CT成像系统通常包括一个X射线源和一个探测器阵列。X射线源发出的射束被物体衰减后，到达探测器并被转化成数字信号。这些数字信号经过适当的校正算法的处理后，被送到重建系统并最终被转化成扫描图像。在目前主流的CT重建算法中，不管是解析重建还是迭代重建，正投/反投影算法是一个关键步骤。正投/反投算法直接影响到图像的生成速度，准确度和图像质量，正投/反投算法可以分为几类：像素驱动(PixelDriven或简称PD)，射线驱动(RayDriven或简称RD)，面积驱动(AreaDriven或简称AD)，距离驱动(DistanceDriven或简称DD)以及可分离足印投影(SeparableFootprint或简称SF)。现有技术中，正投/反投算法无法兼具高精度和快速的优势。

因此，现有技术还有待于改进和发展。

发明内容

本发明要解决的技术问题在于，针对现有技术的上述缺陷，提供一种基于体积驱动的投影方法及其系统，旨在解决现有技术中正投/反投算法无法兼具高精度和快速的优势的问题。

本发明解决技术问题所采用的技术方案如下：

一种基于体积驱动的投影方法，其中，包括步骤：

对同一排像素依次计算像素在探测器上的二维投影；

将二维投影扩展到三维投影。

所述的基于体积驱动的投影方法，其中，所述对同一排像素依次计算像素在探测器上的二维投影步骤具体包括：

根据上一像素在探测器上的二维投影时，没有覆盖的像素面积和像素的值的乘积之和，计算像素在探测器上的二维投影；

更新本次像素在探测器上的二维投影时，没有覆盖的像素面积和像素的值的乘积之和；

计算下一像素在探测器上的二维投影，直至完成一排像素在探测器上的二维投影。

所述的基于体积驱动的投影方法，其中，所述像素在探测器上的二维投影P_j为：

当t_U>t_L时，

当t_U≤t_L时，

t_U＝X_b(j+1)-X_U(j+1)，

t_L＝X_b(j+1)-X_L(j+1)，

其中，P_j'为本次像素在探测器上的二维投影时，投影覆盖的像素面积和像素的值的乘积之和，Δθ为探测器对点源形成的二维角，L(j)为在点源与探测器中心之间的连线、像素行的中心线中，两线相交的交点到点源的距离，B_prev为上一像素在探测器上的二维投影时，没有覆盖的像素面积和像素的值的乘积之和，B_delta为本次像素在探测器上的二维投影时，本次像素中的像素面积和像素的值的乘积之和，X_b(j+1)为本次像素右边界的坐标，X_U(j+1)为二维投影的右边界与图像上边界交点的坐标，X_L(j+1)为二维投影的右边界与图像下边界交点的坐标，t_U为X_b(j+1)与X_U(j+1)之间的距离，t_L为X_b(j+1)与X_L(j+1)之间的距离，S_U,k为，S_L,k为，k为，a为像素长度，θ_j+1为二维投影的右边界与竖直方向的夹角，μ(·)是像素的值，i为本次像素。

所述的基于体积驱动的投影方法，其中，所述本次像素在探测器上的二维投影时，没有覆盖的像素面积和像素的值的乘积之和为：

B_prev'＝B_prev+B_delta-P_j'。

述的基于体积驱动的投影方法，其中，所述将二维投影扩展到三维投影步骤采用线性插值、重叠插值或二维面积驱动中的任意一种方式进行计算。

一种基于体积驱动的投影系统，其中，包括：处理器，以及与所述处理器连接的存储器，

所述存储器存储有基于体积驱动的投影程序，所述基于体积驱动的投影程序被所述处理器执行时实现以下步骤：

构建投影的二维模型；

对同一排像素依次计算像素在探测器上的二维投影；

将二维投影扩展到三维投影。

所述基于体积驱动的投影系统，其中，所述基于体积驱动的投影程序被所述处理器执行时，还实现以下步骤：

根据上一像素在探测器上的二维投影时，没有覆盖的像素面积和像素的值的乘积之和，计算像素在探测器上的投影；

所述基于体积驱动的投影系统，其中，所述像素在探测器上的二维投影P_j为：

当t_U>t_L时，

当t_U≤t_L时，

t_U＝X_b(j+1)-X_U(j+1)，

t_L＝X_b(j+1)-X_L(j+1)，

所述基于体积驱动的投影系统，其中，所述本次像素在探测器上的二维投影时，没有覆盖的像素面积和像素的值的乘积之和为：

B_prev'＝B_prev+B_delta-P_j'。

所述基于体积驱动的投影系统，其中，所述将二维投影扩展到三维投影步骤采用线性插值、重叠插值或二维面积驱动中的任意一种方式进行计算。

有益效果：由于在二维情况下基于面积的快速正投/反投算法，并且将该算法扩展到三维情况，成为基于体积的算法。不仅大幅提高精度，而且射线与像素的相交面积时，同一排像素使用顺序计算的办法，因而大量中间计算变量可以被缓存，因而提高计算速度。。

附图说明

图1是本发明中基于体积驱动的投影方法较佳实施例的流程图。

图2是本发明中基于体积驱动的投影方法的第一原理示意图。

图3是本发明中基于体积驱动的投影方法的第二原理示意图。

图4是本发明中基于体积驱动的投影方法的第三原理示意图。

图5是本发明中基于体积驱动的投影系统的功能原理框图。

具体实施方式

为使本发明的目的、技术方案及优点更加清楚、明确，以下参照附图并举实施例对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

请同时参阅图1-图4，本发明提供了一种基于体积驱动的投影方法的一些实施例。

本发明描述了一种CT系统中的基于体积驱动的正投/反投模型。该模型假设X射线是一个点源A，重建图像中的像素是相邻的均匀长方体，探测器阵列中的探测器也是相邻的均匀长方形。该模型可以扩展到具有面积的X射线源。在此假设的基础上，体积驱动的正投模型认为，一个像素在一个探测器单元上面的投影与这个像素落在点源A和探测器表面构成的锥体内的体积成正比。

具体地，如图2所示，在一个CT系统的重建算法中，X射线源往往被认为是一个点源A，探测器往往也被认为是一个点，但是在我们的模型中，我们认为每个探测器是一个均匀的长方形。重建图像中的每个像素我们用下标i来区分，探测器阵列中的每个探测器我们用下标j来区分。

像素i到探测器j的投影我们可以认为是在整个探测器表面的积分平均。这个积分我们在立体角中表示就是：

其中，P_ij是像素i在探测器j上的投影，μ(i)是像素i的值，Δθ_j是探测器j对点源A形成的立体角。因为图像的像素相对于像素到点源A的距离通常都非常小，可以认为在整个像素内，每一点到点源A的距离都可以用这个像素中心到点源A的距离L_i来近似。那么我们可以得到：

其中，V_ij是像素i落在探测器j和点源A形成的五面锥体中的体积。也就是说，像素i在探测器单元j上面的投影与V_ij成正比。

在二维情况下，上述方程可以简化为：

其中，S_ij为像素i落在探测器j和点源A形成的三角形中的面积，Δθ_j是探测器j对点源A形成的二维角。

P_ij是单个像素i在探测器j上的投影。为了计算整幅图像在探测器j上的投影，需要将所有在探测器j覆盖范围内的像素i进行求和：

为了简化计算，我们将像素i到点源的距离L_i用一个近似值来代替，并将其移出求和以获得加速。这个近似值定义为L(j)。我们先找到点源与探测器j中心的连线与当前像素行的中心线的交点，该交点到点源的距离则为L(j)(如图3所示)。这样，P_j的计算可以简化为：

反投的情况相对简单。探测器j对像素i的贡献B_ij可以通过简单的归一化获得：

其中p(j)是探测器j上的值，T_i是像素i的体积。二维情况就是：

其中T_i是像素i的面积。

如图1所示，一种基于体积驱动的投影方法包括步骤：

S100、对同一排像素依次计算像素在探测器上的二维投影；

步骤S100包括：

步骤S110、根据上一像素在探测器上的二维投影时，没有覆盖的像素面积和像素的值的乘积之和，计算像素在探测器上的二维投影。

该算法对二维图像进行逐行处理，并且每一行的处理都互相独立。因此该算法很容易实现并行。下面的描述中，我们假设图像只有一行。如图3所示，像素依次排列，i-3、i-2、i-1、i、i+1。

计算探测器j与点源A形成二维投影具有两个边界(分别为二维投影左边界和二维投影右边界)。与二维投影右边界相交的最右边的像素是新增像素(即本次像素i)。每个像素都是边长为a的正方形。二维投影左边界与这行图像的上或下边界的交点的坐标分别是X_U(j)、X_L(j)，二维投影右边界与这行图像的上或下边界的交点的坐标分别是X_U(j+1)、X_L(j+1)。本次像素i的右边界是X_b(j)，本次像素i的右边界是X_b(j+1)。二维投影左边界与竖直方向的夹角是θ_j，二维投影右边界与竖直方向的夹角是θ_j+1。

X_b(j+1)与X_U(j+1)之间的距离定义为：

t_U＝X_b(j+1)-X_U(j+1)，

X_b(j+1)与X_L(j+1)之间的距离定义为：

t_L＝X_b(j+1)-X_L(j+1)。

步骤S110具体包括：

步骤S111、计算X_U(j+1)和X_L(j+1)。

步骤S112、根据本次像素i的右边界X_b(j+1)更新B_delta。本次像素在探测器上的二维投影时，本次像素中的像素面积和像素的值的乘积之和为B_delta。在上一次像素计算完成或初始状态时，B_delta被置零，即B_delta＝0，在本次像素计算时，根据本次像素i的右边界X_b(j+1)更新B_delta。

步骤S113、计算t_U和t_L。

步骤S114、比较t_U和t_L，计算P_j'，Pj'为本次像素在探测器上的二维投影时，投影覆盖的像素面积和像素的值的乘积之和。

当t_U>t_L时，

当t_U≤t_L时，

其中，

B_prev为上一像素在探测器上的二维投影时，没有覆盖的像素面积和像素的值的乘积之和，μ(·)是像素的值，S_U,k为，S_L,k为，k为。

步骤S115、计算P_j，P_j为像素在探测器上的二维投影。

Δθ是探测器对点源A形成的二维角，L(j)是像素中心到点源A的距离。

步骤S120、更新B_prev'，B_prev'为本次像素在探测器上的二维投影时，没有覆盖的像素面积和像素的值的乘积之和。

B_prev'＝B_prev+B_delta-P_j'，

更新B_prev'用于计算下一像素在探测器上的二维投影。

步骤S130、计算下一像素在探测器上的二维投影，直至完成一排像素在探测器上的二维投影。

具体步骤参见步骤S110和步骤S120。通过把图像旋转90度的方法，我们总是可以把点源A控制在-45°和45°之间。在实际实现上述算法时，可以假设在边界j+1最多与M个像素点相交，M取决于系统的扇形角大小。例如CT系统的扇形角不能超过53°时，任意边界最多与4个像素点相交。这个限制刚刚能满足目前所有商业CT的要求。在这种情况下，上述算法可以通过循环展开(loopunrolling)的方法进行进一步加速，循环展开是代码编译器中的一个选项。以这种方式编译的可执行程序，实际执行时并不进行循环，而是把循环中的代码顺序执行。比如在本例中，代码中可以把四个像素点的计算写在一个循环中，但是编译器将这个循环编译成四段顺序执行的程序。在某些特殊情况下，这种编译方式能带来效率提升。这种算法也非常适合利用图形加速卡(GPU)来进行加速。

在上述算法中，所有几何长度都可以预先对像素长度a进行归一化，从而减少运算量。1/(2tanθ)，Δθ，和L(j)都可以预先计算，并且它们都不依赖于z，因此从二维到三维的扩展中这些变量可以被重复使用。

步骤S200、将二维投影扩展到三维投影。

在三维情况下，图像由许多层组成，并且探测器由许多排组成。我们先对每一层图像单独进行二维处理，并假设探测器的厚度与图像的厚度一致。如图4所示，对于某一列的探测器以及某一行的图像，由上述的算法，我们可以计算出一系列平行的投影。这些投影从一系列“虚拟点源A'”出发(这些“虚拟点源A'”与实际点源A到旋转中心具有相同的距离)，到达一系列“虚拟探测器”(这些“虚拟探测器”具有与图像相同的厚度)。这些投影在Z方向上的位置是

其中N_Z是图像的层数。需要指出的是，因为在Z方向可以放在算法最里面的循环中，上述的算法对于计算这样的投影效率是非常高的。从二维到三维的扩展就是要把这一系列二维投影在Z方向上重排，从而得到从实际点源A到实际探测器的投影值。这个重排的过程有几种实现方式：

第一、线性插值。

最简单的重排方法就是线性插值。从真实点源A到真实探测器中心的连线与图像的当前行位置相交于

其中T_Nk是探测器的排数。每个探测器与点源A形成的锥角是

假设对于探测器k，T_k位于Z_m1与Z_m2之间，线性插值可以被描述为：

第二、重叠插值。

重叠插值比线性插值更接近理论模型。它在计算探测器k上的投影值时，使用每个图像层与当前探测器的重合部分来作为权重进行计算。假设探测器k在图像当前行上的长度是D_k，图像上从m₁到m₂层与探测器k有重叠，那么重叠插值可以被描述为：

其中，R_k,m是图像层m与探测器k的重叠长度。

第三、二维面积驱动。

如果想要得到更精确的投影值，在Z方向还可以进行一次基于面积驱动的投影。从二维算法中得到的投影值我们可以认为是一幅新的二维图像。重新应用上述二维面积驱动算法，我们可以得到三维情况的投影值。这个投影值非常接近于基于体积驱动的正投模型的理论值。但是这种算法需要对每一列探测器及每一行图像都做面积驱动的算法，因此计算量比较大。这种算法的优势在于如果CT系统的锥角比较大时，因为实际投影与虚拟投影之间大夹角的影响，这种算法的精度会比线性插值或重叠插值要高。

值得说明的是，本发明基于体积驱动的投影方法的精度与SF算法相当。并且提供一种快速计算的算法，其计算速度与DD算法相当。

另外，这种快速算法可以被在计算机的中央处理器(CPU)或者图形处理器(GPU)上实现。这种算法容易在图形处理器(GPU)上实现，利用其大量并行计算的能力，能获得比DD更快的执行速度。射线与像素的相交面积时，不同于传统的每条射线单独计算的方式，同一排像素使用顺序计算的办法，因而大量中间计算变量可以被缓存，因而提高计算速度。

本发明中阐述了一种二维情况下基于面积的快速正投/反投算法，并且将该算法扩展到三维情况，成为基于体积的算法。本发明的主要优势在于在对性能损失不大的情况下，将CT重建中的正投/反投算法的精度大幅提高，为提高图像质量和降低伪影提供了基础。并且本发明中的算法可以方便地在GPU中实现，带来更快的执行速度。

其他领域中，包括但是不局限于，工业CT成像，辐射探伤，其他医用成像设备(PET，SPECT，乳腺断层成像等等)，也可以用到本发明。

本发明还提供了一种基于体积驱动的投影系统的较佳实施例：

如图5所示，本发明实施例所述一种基于体积驱动的投影系统，包括：处理器10，以及与所述处理器10连接的存储器20，

所述存储器20存储有基于体积驱动的投影程序，所述基于体积驱动的投影程序被所述处理器10执行时实现以下步骤：

构建投影的二维模型；

对同一排像素依次计算像素在探测器上的二维投影；

将二维投影扩展到三维投影，具体如上所述。

所述基于体积驱动的投影程序被所述处理器10执行时，还实现以下步骤：

计算下一像素在探测器上的二维投影，直至完成一排像素在探测器上的二维投影，具体如上所述。

所述像素在探测器上的二维投影P_j为：

当t_U>t_L时，

当t_U≤t_L时，

t_U＝X_b(j+1)-X_U(j+1)，

t_L＝X_b(j+1)-X_L(j+1)，

其中，P_j'为本次像素在探测器上的二维投影时，投影覆盖的像素面积和像素的值的乘积之和，Δθ为探测器对点源A形成的二维角，L(j)为在点源与探测器中心之间的连线、像素行的中心线中，两线相交的交点到点源的距离，B_prev为上一像素在探测器上的二维投影时，没有覆盖的像素面积和像素的值的乘积之和，B_delta为本次像素在探测器上的二维投影时，本次像素中的像素面积和像素的值的乘积之和，X_b(j+1)为本次像素右边界的坐标，X_U(j+1)为二维投影的右边界与图像上边界交点的坐标，X_L(j+1)为二维投影的右边界与图像下边界交点的坐标，t_U为X_b(j+1)与X_U(j+1)之间的距离，t_L为X_b(j+1)与X_L(j+1)之间的距离，S_U,k为，S_L,k为，k为，a为像素的边长，θ_j+1为二维投影右边界与竖直方向的夹角，μ(·)是像素的值，i为本次像素，具体如上所述。

所述本次像素在探测器上的二维投影时，没有覆盖的像素面积和像素的值的乘积之和为：

B_prev'＝B_prev+B_delta-P_j'，具体如上所述。

所述将二维投影扩展到三维投影步骤采用线性插值、重叠插值或二维面积驱动中的任意一种方式进行计算，具体如上所述。

综上所述，本发明所提供的一种基于体积驱动的投影方法及其系统，所述方法包括步骤：对同一排像素依次计算像素在探测器上的二维投影；将二维投影扩展到三维投影。由于在二维情况下基于面积的快速正投/反投算法，并且将该算法扩展到三维情况，成为基于体积的算法。不仅大幅提高精度，而且射线与像素的相交面积时，同一排像素使用顺序计算的办法，因而大量中间计算变量可以被缓存，因而提高计算速度。

应当理解的是，本发明的应用不限于上述的举例，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims

1.一种基于体积驱动的投影方法，其特征在于，包括步骤：

对同一排像素依次计算像素在探测器上的二维投影；

将二维投影扩展到三维投影；

所述对同一排像素依次计算像素在探测器上的二维投影步骤具体包括：

计算下一像素在探测器上的二维投影，直至完成一排像素在探测器上的二维投影；

所述像素在探测器上的二维投影P_j为：

当t_U>t_L时，

当t_U≤t_L时，

t_U＝X_b(j+1)-X_U(j+1)，

t_L＝X_b(j+1)-X_L(j+1)，

其中，P_j'为本次像素在探测器上的二维投影时，投影覆盖的像素面积和像素的值的乘积之和，Δθ为探测器对点源形成的二维角，L(j)为在点源与探测器中心之间的连线、像素行的中心线中，两线相交的交点到点源的距离，B_prev为上一像素在探测器上的二维投影时，没有覆盖的像素面积和像素的值的乘积之和，B_delta为本次像素在探测器上的二维投影时，本次像素中的像素面积和像素的值的乘积之和，X_b(j+1)为本次像素右边界的坐标，X_U(j+1)为二维投影的右边界与图像上边界交点的坐标，X_L(j+1)为二维投影的右边界与图像下边界交点的坐标，t_U为X_b(j+1)与X_U(j+1)之间的距离，t_L为X_b(j+1)与X_L(j+1)之间的距离，S_U,k、S_L,k均为用于计算P_j'的中间变量，k为整数，a为像素的边长，θ_j+1为二维投影右边界与竖直方向的夹角，μ(·)是像素的值，i为本次像素；

B_prev′＝B_prev+B_delta-P′_j。

2.根据权利要求1所述的基于体积驱动的投影方法，其特征在于，所述将二维投影扩展到三维投影步骤采用线性插值、重叠插值或二维面积驱动中的任意一种方式进行计算。

3.一种基于体积驱动的投影系统，其特征在于，包括：处理器，以及与所述处理器连接的存储器，

构建投影的二维模型；

对同一排像素依次计算像素在探测器上的二维投影；

将二维投影扩展到三维投影；

所述基于体积驱动的投影程序被所述处理器执行时，还实现以下步骤：

所述像素在探测器上的二维投影P_j为：

当t_U>t_L时，

当t_U≤t_L时，

t_U＝X_b(j+1)-X_U(j+1)，

t_LX_b(j+1)-X_L(j+1)，

其中，P_j'为本次像素在探测器上的二维投影时，投影覆盖的像素面积和像素的值的乘积之和，Δθ为探测器对点源形成的二维角，L(j)为在点源与探测器中心之间的连线、像素行的中心线中，两线相交的交点到点源的距离，B_prev为上一像素在探测器上的二维投影时，没有覆盖的像素面积和像素的值的乘积之和，B_delta为本次像素在探测器上的二维投影时，本次像素中的像素面积和像素的值的乘积之和，X_b(j+1)为本次像素右边界的坐标，X_U(j+1)为二维投影的右边界与图像上边界交点的坐标，X_L(j+1)为二维投影的右边界与图像下边界交点的坐标，t_U为X_b(j+1)与X_U(j+1)之间的距离，t_L为X_b(j+1)与X_L(j+1)之间的距离，S_U,k为，S_L,k为，k为，a为像素的边长，θ_j+1为二维投影右边界与竖直方向的夹角，μ(·)是像素的值，i为本次像素；

B_prev′＝B_prev+B_delta-P′_j。

4.根据权利要求3所述基于体积驱动的投影系统，其特征在于，所述将二维投影扩展到三维投影步骤采用线性插值、重叠插值或二维面积驱动中的任意一种方式进行计算。