CN105678820A - 一种基于cuda的s-bpf重建算法加速方法 - Google Patents

Abstract

本发明公开了一种基于CUDA的S-BPF重建算法加速方法，克服了现有技术中，CT成像的图像重建算法时间长的问题。该发明含有如下步骤：步骤1：从硬盘上读取若干投影并在CPU中计算用于有限Hilbert逆变换的常数C；步骤2：将若干投影从内存传入显存并在GPU中实现求导-反投影，得到DBP图像；步骤3：对步骤2得到的DBP图像进行有限Hilbert逆变换并将结果从显存传回内存。本发明解决了重建算法的GPU加速已经取得了显著的加速效果，通信时延成为了限制现有加速策略的瓶颈。实验结果表明：本发明可获得的加速比约为现有策略的2倍。

Description

一种基于CUDA的S-BPF重建算法加速方法

技术领域

该发明涉及一种CT成像的图像重建算法，特别是涉及一种基于CUDA的S-BPF重建算法加速方法。

背景技术

计算机断层成像(ComputedTomography,CT)可以在无损的条件下同时获得物体的外部和内部几何结构特征。在CT成像过程中，图像重建由于其时间复杂度高的缺点，已成为限制CT成像速度的关键因素。在医疗上，为准确地分析病情，医院通常会借助CT对病人内部组织进行成像。然而由于重建时间长，病人往往会忍受长时间的病痛后才能得到准确的治疗；在手术导航中，为了准确定位病灶位置，常使用C-armCT对病人进行成像，该应用中对重建速度具有严格的要求。因此图像重建耗时长已严重降低了CT成像的使用性。

目前，使用高性能的硬件设备对重建算法进行并行加速是解决其计算时间过长的常用方法。在众多硬件加速器(如计算机集群、DSP、FPGA、GPU等)中，GPU具有并行计算能力强、性价比高、易于开发等优势，在图像重建中得到了广泛的研究和应用。

选择性反投影滤波算法(Selective-BPF,S-BPF)是一种用于短扫描局部重建的快速重建算法。S-BPF算法主要分为投影数据的加权、求导、反投影和DBP图像的滤波四个步骤，其并行性分析如下：

步骤一：投影数据的加权。该操作是对平板探测器获取的投影进行几何修正，其计算表达式如下：

其中，R为光源到旋转轴的距离，p(β,m,n)表示在β扫描角度下平板探测器上坐标为(m,n)点的探测单元采集到的投影值对应的衰减系数，为余弦加权系数，p_w(β,m,n)为加权后的投影。由加权公式可知：1)加权与扫描角度β无关，即不同角度下的投影加权操作可以并行计算；2)加权仅与p(β,m,n)及几何参数有关，即同一角度下不同位置的投影可以并行计算。因此投影数据加权操作具有良好的并行性。

步骤二：加权后投影数据的求导。该操作是将投影沿着扫描轨迹方向求导，在圆轨迹中为沿探测器水平方向求导。该操作是对投影的一种滤波操作，在时域中求导对应在频域中乘以jω。在实际离散投影数据中，我们使用差分代替微分，其计算公式如下：

其中d_pixel为探测器的各像素尺寸，p′_w(β,m,n)为求导后的投影数据。与投影数据加权相似，不同角度下的求导操作可以并行计算，同一角度下不同位置的投影也可并行计算。

步骤三：求导后投影数据的反投影。反投影操作是使用求导后投影数据计算差分反投影(DBP)图像的过程，其计算公式如下：

其中b(x,y,z)为DBP图像，X_min,X_max,Y_min,Y_max,Z_min,Z_max分别为x,y,z的最小值和最大值，ψ和π-ψ为扫描起始扫描角度和终止扫描角度，且ψ＝arcsin(y_min/R)。m(β,x,y,z)和n(β,x,y,z)分别为空间一点(x,y,z)在扫描角度为β时的投影位置，其计算公式为：

对于单个体素的反投影过程又可以分解为以下四个步骤：

1)计算被重建物体中该体素(x,y,z)沿旋转轴旋转β度后的新坐标(x′,y′,z′)。显然，计算各体素旋转后的位置可以在不同或同一角度下并行执行。

2)计算该点在探测器上锥束投影位置的坐标索引(m′,n′)。

其中，D为光源到探测器的距离，(m₀,n₀)是中心射束投射位置的坐标索引。该过程包括几何放大以及坐标-索引的转换两项操作，计算各体素点的投影索引仅与体素旋转后的坐标和若干固定的几何参数有关，因此该过程具有良好的并行性。

3)计算投影中索引为(m′,n′)处的投影值。由于计算得到的(m′,n′)并不一定为正整数，所以需要对其进行插值得到(m′,n′)处的投影值。通常采用三线性插值，其计算公式如下：

p′_w(β,m′,n′)＝(n-[n])((m-[m])p′_w(β,[m]+1,[n]+1)+([m]+1-m)p′_w(β,[m],[n]+1))+

([n]+1-n)((m-[m])p′_w(β,[m]+1,[n])+([m]+1-m)p′_w(β,[m],[n]))

\*MERGEFORMAT(7)

其中，[a]表示不大于a的最大整数。显然，各点的三线性插值操作是可并行的。

4)计算体素(x,y,z)在扫描角度β下的反投加权系数并对体素进行累加。由公式可知，反投加权系数为

该加权系数与步骤3)得到的投影值进行相乘后累加到体素(x,y,z)上。该过程中权值的计算在不同角度或同一角度都是可并行的。该过程中的累加操作在同角度下是可并行的，但是对某一体素同时反投不同角度的投影时，存在读写冲突问题，因此该操作不能在不同角度下实现并行。

步骤四：DBP图像的滤波。该步骤分别在DBP图像上沿各PI线上进行有限Hilbert逆变换。为了计算方便，我们将PI线选择成平行于y轴的线段。该操作可以分为以下5个步骤：

1)对PI线上的DBP图像进行加权

2)对b_w(x,y,z)沿y方向进行FFT

fb_w(x,ω,z)＝fft[b_w(x,y,z)]\*MERGEFORMAT(10)

3)在频域中对fb_w进行希Hilbert滤波

fb′_w(x,ω,z)＝-j×fb_w(x,ω,z)\*MERGEFORMAT(11)

4)对fb′_w沿y方向进行IFFT

b′_w(x,y,z)＝ifft[b_w(x,ω,z)]\*MERGEFORMAT(12)

5)对b′_w进行加权

通过并行性分析，很容易知道以上5个步骤均可实现步骤内部并行执行，但步骤之间必须严格按照以上流程执行。而各PI线之间相互独立，因此各PI线之间的执行是可并行的。

发明内容

本发明克服了现有技术中，CT成像的图像重建算法时间长的问题，提供一种效果明显的基于CUDA的S-BPF重建算法加速方法。

含有如下步骤：

步骤1：从硬盘上读取若干投影到内存，再传入显存中并在CPU中计算用于有限Hilbert逆变换的常数C；

步骤2：GPU使用已传入的若干投影做求导-反投影操作，得到DBP图像；

步骤3：GPU对步骤2计算得到的DBP图像做有限Hilbert逆变换，同时将计算好的结果传回内存，供应用处理或使用。

所述步骤2中GPU执行求导-反投影的特征包含以下步骤：

2.1、对投影数据进行加权和求导；

2.2、将结果从全局存储器拷贝到纹理存储器中；

2.3、执行反投影操作。

所述步骤3中的有限Hilbert逆变换，需要将数据装载到一个特定的数据类型cufftComplex中，具体的执行流程为：

3.1、对PI线上的DBP数据加权并赋值到cufftComplex中；

3.2、使用CUFFT实现快速傅里叶变换；

3.3、在频域进行希尔伯特滤波；

3.4、使用CUFFT库做逆傅里叶变换；

3.5、对其进行加权并赋值DBP图像；

3.6、将结果拷贝到CPU中。

所述步骤1的通过使用两倍投影大小的内存空间解除读取投影操作和计算常数C操作之间的耦合，再利用多线程技术实现两者的并行执行。

所述步骤2通过使用步骤1中两倍投影大小的显存空间解除内存-显存的投影数据传输和GPU中求导-反投影操作之间的耦合，再利用CUDA中的异步技术实现两者的并行执行。

所述GPU执行求导-反投影方法通过使用两倍于步骤1中投影大小的显存空间来解除步骤2.1、步骤2.2、步骤2.3中的耦合，再利用CUDA中的多流技术和异步技术实现GPU中数据传输和计算的并行执行。

所述GPU执行求导-反投影方法中，为了完全隐藏步骤2中数据传输时延，如CPU-GPU数据传输时延和GPU-GPU数据传输时延，根据GPU的性能和图像重建的规模对反投影进行了分块处理。

所述步骤3中的有限Hilbert逆变换方法中，为节约显存，在一定图像重建规模下，根据GPU计算和数据传输性能，对有限Hilbert逆变换进行分块处理，实现了完全隐藏数据传输时延。

所述步骤3中的有限Hilbert逆变换方法中，利用两倍分块大小的显存空间保证两组有限Hilbert逆变换-回传操作的独立性，再利用CUDA中的多流技术和异步技术实现一组有限Hilbert逆变换与另一组的结果回传操作的并行执行。

与现有技术相比，本发明基于CUDA的S-BPF重建算法加速方法具有以下优点：目前重建算法的GPU加速已经取得了显著的加速效果，通信时延成为了限制现有加速策略的瓶颈。本发明基于选择性反投影滤波算法的并行性分析及CPU-GPU异构平台特点，利用多线程技术实现CPU的运算与硬盘-内存的数据传输并行执行，利用CUDA异步技术实现CPU的运算与内存-显存的数据传输并行执行，利用流处理技术和异步技术实现GPU的运算与内存-显存、显存-显存的输出传输并行执行。本发明使用多重事件重叠的方式，在保证计算效率不变的前提下，最大可能地隐藏数据在硬盘、内存和显存之间的通信时延和CPU的运算时间，进而提高重建效率。实验结果表明：本发明可获得的加速比约为现有策略的2倍。在以上步骤中，使用多线程、多流、异步、两倍空间解耦合等技术和方法，最大程度隐藏数据在硬盘、内存、显存之间的传输时间。

附图说明

图1是本发明基于CUDA的S-BPF重建算法加速方法中CPU与GPU的采取并行策略后的时序图；

图2是本发明基于CUDA的S-BPF重建算法加速方法中GPU中针对DBP获取通信时延的隐藏策略示意图；

图3是本发明基于CUDA的S-BPF重建算法加速方法中GPU中针对DBP图像滤波的通信时延的隐藏策略结构示意图；

图4从左至右分别为本发明基于CUDA的S-BPF重建算法加速方法中体模、采用本文策略前以及采用本文策略后的重建结果对比图。

具体实施方式

下面结合附图和具体实施方式对本发明基于CUDA的S-BPF重建算法加速方法作进一步说明：如图所示，本发明为缩短图像重建的整体时间，针对S-BPF算法提出了深度隐藏通信时延优化策略的加速方法。本发明针对S-BPF的两大步骤：获取DBP图像以及有限Hilbert逆变换，设计如下加速策略。

在获取DBP图像过程中，一次反投影必须按以下6个流程顺序执行：1、CPU从硬盘中读取若干张投影数据；2、计算用于滤波的常数C；3、将投影数据从CPU中传入GPU中；4、对投影数据进行加权和求导；5、将结果拷贝到cudaArrary中，为线性滤波纹理做准备；6、对各体素分别计算投影位置，并累加相应的投影值。重复执行1-6直到反投完相应角度范围内的投影。为了方便描述，我们使用P1(Cm)、P2(Cc)、P3(Gm)、P4(Gc)、P5(Gm)、P6(Gc)分别表示以上6个操作，其中C表示该操作在CPU中执行，G表示该操作在GPU中执行；c表示该操作为计算性操作，m表示该操作为数据传输性操作。利用CPU与GPU可以同时进行运算的原理，可以将初始并行策略设计如图1所示。

其中p1-p2为CPU端执行，p3-p6为GPU端执行，在数据不冲突的情况下便可设计如图1并行策略。在一次反投影中，两次计算和两次数据传输分别在GPU上交替执行：P3(Gm)、P4(Gc)、P5(Gm)、P6(Gc)。通常两次计算的执行时间大于两次数据传输的执行时间，并且CUDA模型中多流技术可实现数据传输和计算的并行执行。因此，合理的优化策略可以使数据传输的执行时间完全被隐藏到计算的执行中。在一次反投影中，数据传输工作在大部分时间处于空闲状态。而下一次的反投影中，P4必须在P3完成后执行，P6必须在P5完成后执行。因此，P3可在上一次反投影中数据传输的空闲时间内提前执行，从而下一次反投影中P4可无需等待P3的完成。P6和P5也按照相同的方法执行。通常，P3的执行时间大于P4的执行时间。以上策略中，在两次反投影中有一次P3没有被完全隐藏。为了完全隐藏P3，必须用一个比P3更耗时的操作与其同时执行，因此我们将P6划分成若干个子操作，让部分P6的子操作与P3同时执行，如图2所示。

对于有限Hilbert逆变换，由于有Hilbert逆变换中存在傅里叶变换和逆傅里叶变换，为了提高其执行效率，我们CUDA中自带的CUFFFT库对其进行实现。该操作需要实现将数据装载到一个特定的数据类型cufftComplex中，我们得到较优的执行流程为：1.对PI线上的DBP数据加权并赋值到cufftComplex中；2.使用CUFFT实现快速傅里叶变换；3.在频域进行希尔伯特滤波；4.使用CUFFT库做逆傅里叶变换；5.对其进行加权并赋值DBP图像；6.将结果拷贝到CPU中。为了方便描述，我们使用F1(c)、F2(c)、F3(c)、F4(c)、F5(c)、F6(m)分别表示以上6个流程，其中c和m分别表示计算操作和数据传输操作。(F2和F4中较长的CUFFT的长度有助于提高重建图像质量，但会产生额外的时间开销。通常CUFFT的长度采用2的整数次方且大于等于两倍PI线长度。)

为隐藏图3中F6的时间开销，采用多流技术，实现上一次滤波中的数据传输操作(F6)和下一次滤波中的计算操作(F1～F5)同时执行。为了保证结果的正确性，采用cudaEvent确保同一次滤波中的F6在F5完成后才执行。

结果：实验采用长宽高均为12.8mm的标准Shepp-Logan头模作为测试体模。系统参数设定为：光源到旋转中心的距离为750mm；光源到探测器的距离为1500mm；探测器大小为512×512像素，其像素大小为0.1mm。仿真实验中使用Liu的方法在182度扫描角度范围内均匀采集182张投影。测试平台：惠普Z820工作站。开发环境：VisualStudio2008，CUDA6.5runtimeAPI。使用的GPU为基于Tesla框架的K20c。重建规模分别为512³和1024³体素，各体素大小分别为0.05mm和0.025mm。图4中从左至右分别为体模、采用本文策略前以及采用本文策略后的重建结果。重建时间如表1所示，各时间消耗包含了数据传输、CPU的计算和GPU的计算。

表1选择性反投影滤波重建算法并行加速实验结果

Claims (9)

1.一种基于CUDA的S-BPF重建算法加速方法，其特征是含有如下步骤：

步骤1：从硬盘上读取若干投影并在CPU中计算用于有限Hilbert逆变换的常数C；

步骤2：将若干投影从内存传入显存并在GPU中实现求导-反投影，得到DBP图像；

步骤3：对步骤2得到的DBP图像进行有限Hilbert逆变换并将结果从显存传回内存。

2.根据权利要求1所述的基于CUDA的S-BPF重建算法加速方法，其特征在于：所述步骤2中GPU执行求导-反投影的特征包含以下步骤：

2.1、对投影数据进行加权和求导；

2.2、将结果从全局存储器拷贝到纹理存储器中；

2.3、执行反投影操作。

3.根据权利要求1所述的基于CUDA的S-BPF重建算法加速方法，其特征在于：所述步骤3中的有限Hilbert逆变换，需要将数据装载到一个特定的数据类型cufftComplex中，具体执行流程为：

3.1、对PI线上的DBP数据加权并赋值到cufftComplex中；

3.2、使用CUFFT实现快速傅里叶变换；

3.3、在频域进行希尔伯特滤波；

3.4、使用CUFFT库做逆傅里叶变换；

3.5、对其进行加权并赋值DBP图像；

3.6、将结果拷贝到CPU中。

4.根据权利要求1所述的基于CUDA的S-BPF重建算法加速方法，其特征在于：所述步骤1的通过使用两倍投影大小的内存空间解除读取投影操作和计算常数C操作之间的耦合，再利用多线程技术实现两者的并行执行。

5.根据权利要求1所述的基于CUDA的S-BPF重建算法加速方法，其特征在于：所述步骤2通过使用步骤1中两倍投影大小的显存空间解除内存-显存的投影数据传输和GPU中求导-反投影操作之间的耦合，再利用CUDA中的异步技术实现两者的并行执行。

6.根据权利要求2所述的基于CUDA的S-BPF重建算法加速方法，其特征在于：所述GPU执行求导-反投影方法通过使用两倍于步骤1中投影大小的显存空间来解除步骤2.1、步骤2.2、步骤2.3中的耦合，再利用CUDA中的多流技术和异步技术实现GPU中数据传输和计算的并行执行。

7.根据权利要求1所述的基于CUDA的S-BPF重建算法加速方法，其特征在于：所述GPU执行求导-反投影方法中，为了完全隐藏步骤2中数据传输时延，如CPU-GPU数据传输时延和GPU-GPU数据传输时延，根据GPU的性能和图像重建的规模对反投影进行了分块处理。

8.根据权利要求1所述的基于CUDA的S-BPF重建算法加速方法，其特征在于：所述步骤3中的有限Hilbert逆变换方法中，为节约显存，在一定图像重建规模下，根据GPU计算和数据传输性能，对有限Hilbert逆变换进行分块处理，实现了完全隐藏数据传输时延。

9.根据权利要求1所述的基于CUDA的S-BPF重建算法加速方法，其特征在于：所述步骤3中的有限Hilbert逆变换方法中，利用两倍分块大小的显存空间保证两组有限Hilbert逆变换-回传操作的独立性，再利用CUDA中的多流技术和异步技术实现一组有限Hilbert逆变换与另一组的结果回传操作的并行执行。