CN104424625A - 一种gpu加速cbct图像重建方法和装置 - Google Patents

Abstract

本发明公开了一种GPU加速CBCT图像重建方法，包括以下步骤：读取投影数据，在GPU中采用SART算法作为逼近项，更新重建得到的体；采用自适应梯度下降的方法，使所述重建得到的体的全变差最小化；重复上述步骤，迭代N次，直至算法收敛。本发明将图像全变差最小化优化准则和SART算法相结合，不仅改善了单纯用SART算法重建的图像的质量，还可以用更少的投影数据来重建图像，从而减少成像过程中的X射线对人体的辐射，降低治病风险，同时采用GPU硬件对所提出的算法设计并行算法，有效的减少了迭代图像重建的时间。

Description

一种GPU加速CBCT图像重建方法和装置

【技术领域】

本发明涉及计算机和图像处理技术，尤其涉及一种采样不完全投影数据下的CBCT图像重建的方法和装置。

【背景技术】

基于平板探测器的三维锥束CT具有空间分辨率高、投影数据采集时间短、射线利用率高等优点，是一种很有发展空间的新型的CT设备。锥形束CT的三维图像重建算法一般分为解析类方法和迭代类方法两大类。其中解析类方法如FBP（Filter back projection，滤波反投影）算法，主要是基于傅里叶切片定理来实现，而迭代类方法如ART（Algebraic Reconstruction Technique，代数重建技术）算法和SART（Simultaneous Algebraic Reconstruction Technique，联合代数重建技术）算法，主要是基于线性方程组求解的。

与解析类方法相较，迭代类方法能和特定的成像设备及数据采集物理过程的特性相结合，并能利用某些先验知识，尤其适合于不能获得完整投影数据场合的图像重建。迭代类重建算法，在投影数较少，信噪比低的情况下也能重建出质量好的图像，然而在数学理论推导上比较复杂，计算量大，算法复杂度高，重建所需时间难以满足实际应用的需求。在算法的改进方面，美国斯坦福大学的Donoho等【D.Donoho.Compressed sensing.IEEE Transactions on InformationTheory,2006,52(4):1289-1306】从信号分解和逼近理论提出的压缩感知理论（compressed sensing,CS），证明了稀疏信号可以由远不满足香农/奈奎斯特采样定理的采样数据精确重建，很多研究者将其与已有的图像迭代重建算法相结合，用于不完全投影数据的图像重建中。Sidky EY等人【E.Y.Sidky,C.-M.Kao,X.Pan,Accurate image reconstruction from few-views and limited-angle data indivergent beam CT,J.X-Ray Sci.Technol.2006,14(2):119-139】率先将压缩传感理论引入扇束CT中，提出了基于全变分（total variation,TV)）和凸集投影(Projection onto Convex Sets,POCS)的算法，对局部平滑性很好的图像有较好的重建效果，接着提出了ASD-POCS(adaptive steepest descent-projection onto convexsets)算法【Y.Sidky Emil,Pan Xiaochuan,Image reconstruction in circularcone-beam computed tomography by constrained total variation minimization,Phys.Med.Biol.2008,53(17):4777-4807.】

近几年，随着GPU的性能的进一步提升，其作为通用处理单元的巨大潜力和强大功能在图像视频处理、模式识别、流体力学计算、生物计算等领域表现的越来越明显。三维重建在应用于临床或工程时，在重建速度、精度等方面都有严格要求，而图形处理器技术正是以追求更加真实、更为实时的3D景物模拟为目标的。三维重建中所使用的数据规模庞大，操作方式简单，非常符合流数据处理对象的特点，这使得基于GPU实现算法加速成为一种适合三维重建特点的硬件加速方法。目前国内外已经有许多学者对解析类方法进行了硬件加速的研究，并取得了很好的效果，但由于迭代类算法的结构较为复杂，设计符合图形硬件的实现方式相对较难，所以对迭代类算法的硬件加速研究较少。本发明基于CUDA技术，利用GPU硬件加速迭代类三维锥束CT重建中的SART算法，在3秒内，可以利用80幅1282的投影图像重建1283的体，而且不损失质量，相对于主流的CPU实现了100倍以上的加速效果（VIDIA8800GTX GPU），相对于Graphics-based方法实现了近3倍的加速效果[Yuqiang Lu,Weiming Wang,ShifuChen,Yongming Xie,Jing Qin,Wai-Man Pang,Pheng-Ann Heng,AcceleratingAlgebraic Reconstruction Using CUDA-Enabled GPU,in preceeding of IEEEconference on Computer Graphics,Imaging and Visualization(CGIV09).2009.]。

传统的解析图像重建算法当受到噪声影响或者不完全投影数据情形下无法得到满意的结果。迭代图像重建算法如SART算法，虽然可以获得比解析法更好的重建结果，然而仍然不能满足实际的要求。而且迭代图像重建算法需要多次迭代，每次都要对大量的数据进行正投影和反投影计算，因而重建速度较慢。

【发明内容】

针对现有的迭代图像重建算法需要多次迭代，每次都要对大量的数据进行正投影和反投影计算，因而重建速度较慢，本发明实施例提供一种GPU加速CBCT图像重建方法，不仅改善了单纯用SART算法重建的图像的质量，还可以用更少的投影数据来重建图像，从而减少成像过程中的X射线对人体的辐射，降低治病风险，同时采用GPU硬件对所提出的算法设计并行算法，有效的减少了迭代图像重建的时间。

为达到上述目的，本发明解决其技术问题所采用的技术方案是，所述的一种GPU加速CBCT图像重建方法，包括以下步骤：

步骤（1），读取投影数据，在GPU中采用SART算法作为逼近项，更新重建得到的体；

步骤（2），采用自适应梯度下降的方法，使所述重建得到的体的全变差最小化；

重复步骤（1）和步骤（2），迭代N次，直至算法收敛，其中，所述N为大于一的自然数。

具体的，所述步骤（1）还包括以下步骤：

（10）：当所述重建得到的体的每个投影角度的投影数据由m*m个象素组成时，启动m*m个并行的独立线程，每一个所述独立线程对穿过所述投影角度的一条光线，进行如下计算，其中，m为自然数：

（101）：在探测器平面上计算所述投影角度的模拟投影数据；

（102）：根据真实投影数据与所述每个投影角度的模拟投影数据，计算投影误差，得到修正影像；

具体的，按如下公式计算投影误差：

${\mathrm{\ϵ}}_i=\frac{p_i-{p_i}^{\′}}{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{w}_{\mathrm{in}}}$

其中，εi为第i条光线的投影误差，i∈[1,m×m]，p_i为所述真实投影数据的值，p′_i为所述模拟投影数据的值，w_in代表了体素v_j对到达象素p′_i的射线的影响。

优选的，所述p′i，采用Ray-tracing算法。

（11）：当所述重建得到的体的大小为n*n*n个体素时，启动n*n*n个并行的独立线程，每一个所述独立线程对所述重建得到的体数据中的一个体素，进行如下计算，其中，n为自然数：

（111）：将所述投影误差依照计算得到的权重值，反投到所述重建得到的体的各个体素上，生成更新数据；

具体的，按如下公式计算反投：

其中，k为子迭代过程的序数，λ为松弛因子，在(0.0，1.0]区间内取值，v_j ^(k)为第k次迭代的第j个体素值。

（112）：更新整个所述重建得到的体；

（12）：重复步骤（10）和步骤（11），直至遍历用于所述重建的所有投影角度。

具体的，所述步骤（2）还包括以下步骤：

（22）：启动n*n*n个并行的独立线程，每一个所述独立线程计算所述重建得到的体数据的全变分图像中的一个体素的近似偏导，计算公式如下：

$\frac{\∂{\left|\right|\stackrel{\→}{v}\left|\right|}_{\mathrm{TV}}}{\∂v_{x,y,z}}\≈\frac{3\×v_{x,y,z}-v_{x-1,y,z}-v_{x,y-1,z}-v_{x,y,z-1}}{\sqrt{{(v_{x,y,z}-v_{x-1,y,z})}^2+{(v_{x,y,z}-v_{x,y-1,z})}^2+{(v_{x,y,z}-v_{x,y,z-1})}^2+\mathrm{\ξ}}}$

$-\frac{v_{x+1,y,z}-v_{x,y,z}}{\sqrt{{(v_{x+1,y,z}-v_{x,y,z})}^2+{(v_{x+1,y,z}-v_{x+1,y-1,z})}^2+{(v_{x,y,z}-v_{x,y+1,z-1})}^2+\mathrm{\ξ}}}$

$-\frac{v_{x,y+1,z}-v_{x,y,z}}{\sqrt{{(v_{x,y+1,z}-v_{x-1,y+1,z})}^2+{(v_{x,y+1,z}-v_{x,y,z})}^2+{(v_{x,y+1,z}-v_{x,y+1,z-1})}^2+\mathrm{\ξ}}}$

$-\frac{v_{x,y,z+1}-v_{x,y,z}}{\sqrt{{(v_{x,y,z+1}-v_{x-1,y,z+1})}^2+{(v_{x,y,z+1}-v_{x,y-1,z+1})}^2+{(v_{x,y,z+1}-v_{x,y,z})}^2+\mathrm{\ξ}}}$

其中，v_x,y,z为所述重建得到的体数据中的一个体素，ξ为一个避免除零的极小的正数。

优选的ξ为10的负8次方。

具体的，所述步骤（2）还包括以下步骤：

（20）：计算所述重建得到的体数据中的任一个体素的梯度，计算公式如下：

$|\▿v_{x,y,z}|=\sqrt{{(v_{x+1,y,z}-v_{x,y,z})}^2+{(v_{x,y+1,z}-v_{x,y,z})}^2+{(v_{x,y,z+1}-v_{x,y,z})}^2}$

其中，x,y,z为v_x,y,z的三维索引。

本发明的另一目的，在于提供一种GPU加速CBCT图像重建装置，包括以下部分：

图像重建模块：用于读取投影数据，在GPU中采用SART算法作为逼近项，更新重建得到的体；

除噪模块：用于采用自适应梯度下降的方法，使所述重建得到的体的全变差最小化；

迭代模块：用于迭代直至算法收敛。

具体的，所述图像重建模块还包括：

正投影单元：用于在探测器平面上计算所述投影角度的模拟投影数据，并根据真实投影数据与所述每个投影角度的模拟投影数据，计算投影误差，得到修正影像；

反投影单元：用于将所述投影误差依照预定权重值，反投到所述重建得到的体的各个体素上，生成更新数据，并更新整个所述重建得到的体。

与现有的图像重建方法和装置相比，本发明实施例具有如下优点：

1、改善了单纯用SART算法重建的图像的质量，还可以用更少的投影数据来重建图像，从而减少成像过程中的X射线对人体的辐射，降低治病风险。

2、同时采用GPU硬件对所提出的算法设计并行算法，有效的减少了迭代图像重建的时间。

【附图说明】

图1为本发明实施例提供的GPU加速CBCT图像重建方法流程图；

图2为本发明实施例提供的GPU加速CBCT图像重建方法正投影并行流程图；

图3为本发明实施例提供的GPU加速CBCT图像重建方法反投影并行流程图；

图4为本发明实施例提供的GPU加速CBCT图像重建装置结构图；

图5为本发明实施例提供的图像重建模块结构图。

【具体实施方式】

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例的目的之一，是提供一种GPU加速CBCT图像重建方法，对于采样不完全投影数据下的CBCT图像重建，结合图像全变差(TV)最小化和SART算法进行CBCT图像重建：首先采用SART算法作为逼近项用于更新重建得到的体，然后采用自适应梯度下降的方法使重建的体的全变差最小化，这两个过程交替进行，迭代一定的次数，直到算法收敛。本发明实施例不仅改善了单纯用SART算法重建的图像的质量，还可以用更少的投影数据来重建图像，从而减少成像过程中的X射线对人体的辐射，降低治病风险，同时采用GPU硬件对所提出的算法设计并行算法，有效的减少了迭代图像重建的时间。

实施例一，如图1所示，包括以下步骤：

步骤（1），读取投影数据，在GPU中采用SART算法作为逼近项，更新重建得到的体；

步骤（2），采用自适应梯度下降的方法，使所述重建得到的体的全变差最小化；

重复步骤（1）和步骤（2），迭代N次，直至算法收敛，其中，所述N为大于一的自然数。

本实施例通过将SART算法和TV最小化算法相结合，提出的算法的伪代码如下。

算法伪代码：

ng nd←4，λ←0.25，γ_max←0.95,γ_red←0.95

$\stackrel{\→}{v}\←0$

Iteration_num←10//主循环迭代次数

for(i=0;i<Iteration_num;i++)do

${\stackrel{\&RightArrow;}{v}}_0\&LeftArrow;\stackrel{\&RightArrow;}{v}$

for(j=0;j<Proj_num;j++)do//SART循环

1)正投：在探测器平面上计算该投影角度的所有模拟投影数据；

2)修正：根据实际的“投影数据”与相应角度的模拟投影数据计算出“误差数据”，得到修正影像；

3)反投：将误差数据按照一定的“权重”反投到该目标重建体的各个体素上形成“更新数据”，公式如下，λ为松弛因子；

4)更新：更新整个目标重建体；end for

for(j=0;j<Vol_num;j++)do//保证数据非负

if v_j＜0then v_j←0end if

end for

$\mathrm{dp}\&LeftArrow;\stackrel{\&RightArrow;}{v}-{\stackrel{\&RightArrow;}{v}}_0$

if first iteration then dtvg←0end if

${\stackrel{\&RightArrow;}{v}}_0\&LeftArrow;\stackrel{\&RightArrow;}{v}$

for(i=0;i<nd;i++)do//TV下降循环

$\stackrel{\&RightArrow;}{d}v\&LeftArrow;{\&dtri;}_{\stackrel{\&RightArrow;}{v}}{\left|\right|\stackrel{\&RightArrow;}{v}\left|\right|}_{\mathrm{TV}}$

$\hat{d}v\&LeftArrow;\stackrel{\&RightArrow;}{d}v/\left|\stackrel{\&RightArrow;}{d}v\right|$

$\stackrel{\&RightArrow;}{v}\&LeftArrow;\stackrel{\&RightArrow;}{v}-\mathrm{dtvg}*\hat{d}v$

end for

$\mathrm{dg}\&LeftArrow;\stackrel{\&RightArrow;}{v}-{\stackrel{\&RightArrow;}{v}}_0$

if dg>γ_max then dtvg←dtvg*γ_red end if

end for

return

因此，三维CT系统可以描述为一个线性系统，即

Wv＝p    （1）

其中v=[v1,v2,v3]T是一个待求的体素值向量，它包括N=n3个元素，这些体素构成了一个n*n*n大小的目标重建体；p=[p1,p2,p3]T是投影数据向量，如果共有M个投影角度，每个投影角度的投影数据由L(m*m)个象素组成，则R=ML，对于每个投影数据值都是在投影角度为时获得的，这里的是索引i的一个子集；W是一个R×N大小的权重矩阵，它的每一个元素w_ij代表了体素v_j对到达象素（对应于一个探测器）p_i的射线的影响。

在SART算法中，每一次迭代过程主要包括四步：正投、修正、反投、更新。为了设计并行的算法，本发明将前两步合并，称之为正投影，目的是计算修正影像；将后两步合并，称之为反投影，目的是更新体数据。这样的划分是为了程序更好的并行执行，并能够更好的符合GPU的硬件架构，以实现更快的重建速度。

实施例二，如图2所示，在SART循环部分的GPU加速方面，在正投影过程中，设计为光线驱动的方法，即一个线程计算一根光线。这样，对于m*m的投影图片，本发明就可以得到m*m个线程。

一种GPU加速CBCT图像重建方法，包括以下步骤：

步骤（1），读取投影数据，在GPU中采用SART算法作为逼近项，更新重建得到的体；

步骤（2），采用自适应梯度下降的方法，使所述重建得到的体的全变差最小化；

重复步骤（1）和步骤（2），迭代N次，直至算法收敛，其中，所述N为大于一的自然数。

具体的，所述步骤（1）还包括以下步骤：

（10）：当所述重建得到的体的每个投影角度的投影数据由m*m个象素组成时，启动m*m个并行的独立线程，每一个所述独立线程对穿过所述投影角度的一条光线，进行如下计算，其中，m为自然数：

（101）：在探测器平面上计算所述投影角度的模拟投影数据；

（102）：根据真实投影数据与所述每个投影角度的模拟投影数据，计算投影误差，得到修正影像；

具体的，按如下公式计算投影误差：

${\mathrm{\&epsiv;}}_i=\frac{p_i-{p_i}^{\&prime;}}{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}_{\mathrm{in}}}---\left(2\right)$

其中，εi为第i条光线的投影误差，i∈[1,m×m]，p_i为所述真实投影数据的值，p′_i为所述模拟投影数据的值，w_in代表了体素v_j对到达象素p′_i的射线的影响。

优选的，所述p′_i，采用Ray-tracing算法。

对于正投过程中线积分的计算，也就是此处的p′_i，采用Ray-tracing方法。该方法是基于采样的点模型，插值计算采样点。根据精度的要求，本发明可以方便的控制采样的频率，能够取得更为精细的效果。

如图3所示，在反投影过程中，设计为体素驱动的方法，即一个线程计算体数据中的一个体数。这样，对于n*n*n的三维物体，本发明就可以得到n*n*n个线程，即：

步骤（11）：当所述重建得到的体的大小为n*n*n个体素时，启动n*n*n个并行的独立线程，每一个所述独立线程对所述重建得到的体数据的全变分图像中的一个体素，进行如下计算，其中，n为自然数：

（111）：将所述投影误差依照计算得到的权重值，反投到所述重建得到的体的各个体素上，生成更新数据；

具体的，按如下公式计算反投：

其中，k为子迭代过程的序数，λ为松弛因子，在(0.0，1.0]区间内取值，v_j ^(k)为第k次迭代的第j个体素值。

（112）：更新整个所述重建得到的体；

（12）：重复步骤（10）和步骤（11），直至遍历用于所述重建的所有投影角度。

在TV最小化过程中，对一个给定的三维数据体v中的每个体素v_x,y,z，梯度可按以下式计算：

$|\&dtri;v_{x,y,z}|=\sqrt{{(v_{x+1,y,z}-v_{x,y,z})}^2+{(v_{x,y+1,z}-v_{x,y,z})}^2+{(v_{x,y,z+1}-v_{x,y,z})}^2}---\left(4\right)$

其中x,y,z为该体素的三维索引。

公式（4）只对非边界体素有效。v的梯度也可以看作一个三维体数据，该体数据中的每个体素的近似偏导，由于在TV下降循环中，最耗时的过程是偏导的计算，本发明采用体素驱动的方法，一个线程计算一个体素的偏导，而且线程之间完全独立，这样就可以有效的减小运算时间，即：

具体的，所述步骤（2）还包括以下步骤：启动n*n*n个并行的独立线程，每一个所述独立线程计算所述重建得到的体数据的全变分图像中的一个体素的近似偏导，计算公式如下：

$\frac{\&PartialD;{\left|\right|\stackrel{\&RightArrow;}{v}\left|\right|}_{\mathrm{TV}}}{\&PartialD;v_{x,y,z}}\&ap;\frac{3\&times;v_{x,y,z}-v_{x-1,y,z}-v_{x,y-1,z}-v_{x,y,z-1}}{\sqrt{{(v_{x,y,z}-v_{x-1,y,z})}^2+{(v_{x,y,z}-v_{x,y-1,z})}^2+{(v_{x,y,z}-v_{x,y,z-1})}^2+\mathrm{\&xi;}}}$

$-\frac{v_{x+1,y,z}-v_{x,y,z}}{\sqrt{{(v_{x+1,y,z}-v_{x,y,z})}^2+{(v_{x+1,y,z}-v_{x+1,y-1,z})}^2+{(v_{x,y,z}-v_{x,y+1,z-1})}^2+\mathrm{\&xi;}}}---\left(5\right)$

$-\frac{v_{x,y+1,z}-v_{x,y,z}}{\sqrt{{(v_{x,y+1,z}-v_{x-1,y+1,z})}^2+{(v_{x,y+1,z}-v_{x,y,z})}^2+{(v_{x,y+1,z}-v_{x,y+1,z-1})}^2+\mathrm{\&xi;}}}$

$-\frac{v_{x,y,z+1}-v_{x,y,z}}{\sqrt{{(v_{x,y,z+1}-v_{x-1,y,z+1})}^2+{(v_{x,y,z+1}-v_{x,y-1,z+1})}^2+{(v_{x,y,z+1}-v_{x,y,z})}^2+\mathrm{\&xi;}}}$

其中ξ为一个避免除零的极小的正数，这儿取10的负8次方。

相对于Sidky EY等提出的ASD-POCS算法，本发明实施例在数据的逼近项上采用的是SART算法，与ASD-POCS算法里使用的ART算法相比，SART算法中目标重建体是逐一角度更新的，当一个投影角度内的所有探测器误差数据都计算完成后，才进行该角度的反投影过程。逐角度更新显然是更加合理的，因为在某一角度，一个体素实际上可能被很多光线穿越，而不是一根光线。而且，对于本发明的并行算法设计来说，逐角度更新将提供更多更有效的可并行线程。更重要的是，本发明通过对SART算法和TV最小化的过程进行了并行算法设计，实现了完整的SART-TV算法的GPU加速，相比单纯的SART算法，可以用更少的投影数获得较好的CBCT重建结果。

实施例三，相应的，本发明实施例还提供一种GPU加速CBCT图像重建装置，如图4所示，包括以下部分：

图像重建模块：用于读取投影数据，在GPU中采用SART算法作为逼近项，更新重建得到的体；

除噪模块：用于采用自适应梯度下降的方法，使所述重建得到的体的全变差最小化；

迭代模块：用于迭代直至算法收敛。

如图5所示，具体的，所述图像重建模块还包括：

正投影单元：用于在探测器平面上计算所述投影角度的模拟投影数据，并根据真实投影数据与所述每个投影角度的模拟投影数据，计算投影误差，得到修正影像；

反投影单元：用于将所述投影误差依照预定权重值，反投到所述重建得到的体的各个体素上，生成更新数据，并更新整个所述重建得到的体。

本发明实施例将图像全变差最小化优化准则和SART算法相结合，不仅改善了单纯用SART算法重建的图像的质量，还可以用更少的投影数据来重建图像，从而减少成像过程中的X射线对人体的辐射，降低治病风险，同时采用GPU硬件对所提出的算法设计并行算法，有效的减少了迭代图像重建的时间。

以上所述本发明的具体实施方式，并不构成对本发明保护范围的限定。任何根据本发明的技术构思所作出的各种其他相应的改变与变形，均应包含在本发明权利要求的保护范围内。

Claims (10)

1.一种GPU加速CBCT图像重建方法，其特征在于，包括以下步骤：

（1）：读取投影数据，在GPU中采用SART算法作为逼近项，更新重建得到的体；

（2）：采用自适应梯度下降的方法，使所述重建得到的体的全变差最小化；

（3）：重复步骤（1）和步骤（2），迭代N次，直至算法收敛，其中，所述N为大于一的自然数。

2.如权利要求1所述的GPU加速CBCT图像重建方法，其特征在于，所述步骤（1）还包括以下步骤：

（10）：当所述重建得到的体的每个投影角度的投影数据由m*m个象素组成时，启动m*m个并行的独立线程，每一个所述独立线程对穿过所述投影角度的一条光线，进行如下计算，其中，m为自然数：

（101）：在探测器平面上计算所述投影角度的模拟投影数据；

（102）：根据真实投影数据与所述每个投影角度的模拟投影数据，计算投影误差，得到修正影像；

（11）：当所述重建得到的体的大小为n*n*n个体素时，启动n*n*n个并行的独立线程，每一个所述独立线程对所述重建得到的体数据中的一个体素，进行如下计算，其中，n为自然数：

（111）：将所述投影误差依照计算得到的权重值，反投到所述重建得到的体的各个体素上，生成更新数据；

（112）：更新整个所述重建得到的体；

（12）：重复步骤（10）和步骤（11），直至遍历用于所述重建的所有投影角度。

3.如权利要求1或2所述的GPU加速CBCT图像重建方法，其特征在于，所述步骤（2）还包括以下步骤：

（21）：启动n*n*n个并行的独立线程，每一个所述独立线程计算所述重建得到的体数据的全变分图像中的一个体素的近似偏导，计算公式如下：

$\frac{\&PartialD;{\left|\right|\stackrel{\&RightArrow;}{v}\left|\right|}_{\mathrm{TV}}}{\&PartialD;v_{x,y,z}}\&ap;\frac{3\&times;v_{x,y,z}-v_{x-1,y,z}-v_{x,y-1,z}-v_{x,y,z-1}}{\sqrt{{(v_{x,y,z}-v_{x-1,y,z})}^2+{(v_{x,y,z}-v_{x,y-1,z})}^2+{(v_{x,y,z}-v_{x,y,z-1})}^2+\mathrm{\&xi;}}}$

$-\frac{v_{x+1,y,z}-v_{x,y,z}}{\sqrt{{(v_{x+1,y,z}-v_{x,y,z})}^2+{(v_{x+1,y,z}-v_{x+1,y-1,z})}^2+{(v_{x,y,z}-v_{x,y+1,z-1})}^2+\mathrm{\&xi;}}}$

$-\frac{v_{x,y+1,z}-v_{x,y,z}}{\sqrt{{(v_{x,y+1,z}-v_{x-1,y+1,z})}^2+{(v_{x,y+1,z}-v_{x,y,z})}^2+{(v_{x,y+1,z}-v_{x,y+1,z-1})}^2+\mathrm{\&xi;}}}$

$-\frac{v_{x,y,z+1}-v_{x,y,z}}{\sqrt{{(v_{x,y,z+1}-v_{x-1,y,z+1})}^2+{(v_{x,y,z+1}-v_{x,y-1,z+1})}^2+{(v_{x,y,z+1}-v_{x,y,z})}^2+\mathrm{\&xi;}}}$

其中，v_x,y,z为所述重建得到的体数据中的一个体素，ξ为一个避免除零的极小的正数。

4.如权利要求3所述的GPU加速CBCT图像重建方法，其特征在于，所述步骤（2）还包括以下步骤：

（20）：计算所述重建得到的体数据中的任一个体素的梯度，计算公式如下：

$|\&dtri;v_{x,y,z}|=\sqrt{{(v_{x+1,y,z}-v_{x,y,z})}^2+{(v_{x,y+1,z}-v_{x,y,z})}^2+{(v_{x,y,z+1}-v_{x,y,z})}^2}$

其中，x,y,z为v_x,y,z的三维索引。

5.如权利要求4所述的GPU加速CBCT图像重建方法，其特征在于，ξ为10的负8次方。

6.如权利要求2所述的GPU加速CBCT图像重建方法，其特征在于，所述步骤（102）中，按如下公式计算投影误差：

${\mathrm{\&epsiv;}}_i=\frac{p_i-{p_i}^{\&prime;}}{\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}_{\mathrm{in}}}$

其中，ε_i为第i条光线的投影误差，i∈[1,m×m]，p_i为所述真实投影数据的值，p′_i为所述模拟投影数据的值，w_in代表体素v_j对到达象素p′_i的射线的影响。

7.如权利要求6所述的GPU加速CBCT图像重建方法，其特征在于，所述p′_i，采用Ray-tracing算法。

8.如权利要求2所述的GPU加速CBCT图像重建方法，其特征在于，所述步骤（111）中，按如下公式计算反投：

其中，k为子迭代过程的序数，λ为松弛因子，在(0.0，1.0]区间内取值，v_j ^(k)为第k次迭代的第j个体素值。

9.一种GPU加速CBCT图像重建装置，其特征在于，包括以下部分：

图像重建模块：用于读取投影数据，在GPU中采用SART算法作为逼近项，更新重建得到的体；

除噪模块：用于采用自适应梯度下降的方法，使所述重建得到的体的全变差最小化；

迭代模块：用于迭代直至算法收敛。

10.如权利要求9所述的一种GPU加速CBCT图像重建装置，其特征在于，所述图像重建模块，还包括以下部分：

正投影单元：用于在探测器平面上计算所述投影角度的模拟投影数据，并根据真实投影数据与所述每个投影角度的模拟投影数据，计算投影误差，得到修正影像；

反投影单元：用于将所述投影误差依照计算得到的权重值，反投到所述重建得到的体的各个体素上，生成更新数据，并更新整个所述重建得到的体。