CN107194864A - 基于异构平台的ct图像三维重建加速方法及其装置 - Google Patents

Abstract

本发明涉及一种基于异构平台的CT图像三维重建加速方法及其装置，异构平台包含主机及异构OpenCL计算设备，该加速方法包含：对FDK重建算法进行算粒分解，分析各个算粒的并行计算流程；通过异构平台中的主机及异构OpenCL计算设备对各个算粒进行加速优化处理。本发明深度挖掘CT重建算法的可并行性，采用GPU+FPGA的异构计算模式，使用不同类型指令集和体系架构的计算单元组成计算***，最大程度上使算法与异构体系结构相匹配，充分利用不同加速部件的性能；同时设计适合重建算法高效运算的存储与通信方案，***支持PCI‑E/Ethernet互联，支持多块处理板通过互联总线实现多处理器高效并行处理，实现同步或异步的协同处理机制，在尽量减少的损失精度的前提下提高重建速度。

Description

基于异构平台的CT图像三维重建加速方法及其装置

技术领域

本发明属于X射线CT技术领域，特别涉及一种基于异构平台的CT图像三维重建加速方法及其装置。

背景技术

X射线计算机断层成像(Computed Tomography，CT)是一种由物体的X射线投影来反求其衰减分布的技术，涵盖核物理、数学、计算机、精密仪器等多个学科。由于CT可以在非接触、不破坏的条件下得到物体内部的高精度三维结构信息，因此自Hounsfield成功研制第一台CT以来，CT已在无损检测、医学诊断、材料分析等领域得到广泛应用。

在实际应用中，锥束CT高分辨率三维重建所需的计算资源和存储资源都非常大，随着重建规模的增加，重建的存储需求和计算量迅速增加，在很多情况下，很难满足实际应用的需求，以重建算法中通用的反投影计算为例：设待重建三维图像的各维尺寸都为N，则相应反投影的计算复杂度将高达O(N4)，重建一个分辨率为10243的三维图像，计算将循环10995亿次，在普通PC机上完成如此大的计算量是非常耗时的，难以满足实际应用的要求。因此，锥束CT重建过程加速是当前工程人员亟需解决的问题，针对CT重建算法设计相应的加速平台和加速策略具有重要的实际意义，是工业CT技术在实际应用中急需解决的一个难点问题。

发明内容

针对现有技术中的不足，本发明提供一种基于异构平台的CT图像三维重建加速方法及其装置，针对CT重建算法的特性，结合FPGA、GPU等通用加速器件性能并基于异构加速平台实现，在尽量减少损失精度的前提下提高其重建速度，性能稳定，加速效果理想。

按照本发明所提供的设计方案，一种基于异构平台的CT图像三维重建加速方法，异构平台包含主机及异构OpenCL计算设备，该加速方法包含如下内容：对FDK重建算法进行算粒分解，分析各个算粒的并行计算流程；通过异构平台中的主机及异构OpenCL计算设备对各个算粒进行加速优化处理。

上述的，所述的主机为运行主程序的CPU，OpenCL计算设备包含运行内核程序的异构容器GPU和FPGA，CPU、GPU及FPGA之间通过PCI-E总线相通信，主程序通过定义上下文来管理内核程序的运行。

优选的，对FDK重建算法进行算粒分解，包含：根据FDK算法内容，分解为：用于对投影数据进行加权的投影加权算粒，用于对加权后的投影数据进行滤波的滤波算粒，用于将滤波后的投影数据反投影到重建物体上的反投影算粒，及用于对反投影结果进行归约处理的归约算粒。

优选的，根据FDK重建算法公式：

通过对公式中积分的拆分及离散化，分为：

投影加权算粒，表示为：其中，p'(θ,u,v)表示旋转分度为θ时对投影数据进行加权后的数据，为加权系数；

滤波算粒，表示为：其中，d_f(θ,u,v)为滤波后的数据，h(u)为滤波算子的单位冲激响应，[-u_m,u_m]表示探测器每一行采集的2m个数据；

反投影算粒，表示为：其中，f(x,y,z,θ)表示重建物体在旋转分度为θ时投影点对f(x,y,z)的贡献值；

归约算粒，表示为：其中，φ_max为重建物体旋转一周时离散采得的投影分度数。

上述的，对各个算粒进行加速优化处理，包含：采用FPGA对投影加权算粒进行并行处理，通过异步传输至GPU，在传输过程中同时对滤波算粒进行处理；结合反投影时各体素点的数据并行操作，在GPU内按体素点对反投影算粒进行多线程并行反投影计算。

优选的，根据FDK重建算法中重建区域在旋转轴方向上行的各层投影数据和探测器投影数据纵轴方向上的各行数据的投影对应关系，采用分块重建策略，沿旋转轴方向将待重建区域分成若干块，重建其中一块时从外部存储器取出对应的投影数据进行重建操作。

优选的，采用FPGA对投影加权算粒进行并行处理，包含：将全局存储器分割为2块bank，通过负载分配实现随机存取存储器的访问均衡；通过常数存储器存储需要多次重复计算的中间变量。

优选的，在GPU内按体素点对反投影算粒进行多线程并行反投影计算，包含：采用基于体素驱动方式，按重建体数据输出对GPU进行任务划分；将计算中与体素无关的变量进行分离和合并，并在反投影前计算出来存储于GPU的常数存储器中，反投影计算时，直接读取常数存储器中的变量参与计算；优化内核程序中一次反投影的张数。

更进一步地，将计算中与体素无关的变量进行分离和合并，包含内容如下：体数据中任意一点(x，y，z)在投影角度为θ时投射到探测器上的点(u，v)，投影点(u,v)的计算为：

u＝(x-vCenter)×cos(θ)+(y-vCenter)sin(θ)+pCenter

dis＝(u-pCenter)×a

v＝(z-(s₀+θ×h)-γ×h/a)×w+γ×h/a+pCenter

，分离和合并变量后为:

u＝x×A[0]+y×A[1]+A[2]+pCenter

dis＝(u-pCenter)×a

v＝(z-A[4]-γ×A[5])×w+γ×A[5]+pCenter

，其中，vCenter表示体数据中心，pCenter为投影数据中心，α为体素大小，θ为投影角度，r为旋转射线源旋转半径，h为螺距，γ为射束在中心平面上投影与中心射束的夹角。

一种基于异构平台的CT图像三维重建加速装置，异构平台采用PCI-Express作为传输数据信号和控制信号的互连总线，并以Ethernet作为与外部进行联网控制和数据传输的附加总线；异构平台的架构包含提供功能模块的应用层，提供面向应用层各功能模块基于不同处理器的构件库及重建算法所需接口规范的构件层，及面向构件层和应用层提供服务的支撑层，支撑层包含执行主程序的CPU及执行内核程序的多个OpenCL计算设备，CPU和OpenCL计算设备相通信连接，所述的OpenCL计算设备包含GPU、FPGA；基于所述的异构平台的架构的CT图像三维重建加速装置包含如下内容：

算粒分解模块，用于根据FDK重建算法内容将该算法以算粒形式进行拆分和离散化，分解为用于对投影数据进行数据校正的投影加权算粒，用于对加权后的投影数据进行滤波的滤波算粒，用于将滤波后的投影数据反投影到重建物体上的反投影算粒，及用于对反投影结果进行归约处理的归约算粒；

加速处理模块，用于对通过附加总线传输至异构平台计算节点上的投影数据，根据用户设置和重建性能评估，通过互连总线进行协调处理，分别在GPU和FPGA加速部件上完成算粒加速，重建数据实时存储并反馈至用户。

本发明的有益效果：

本发明深度挖掘CT重建算法的可并行性，采用GPU+FPGA的异构计算模式，使用不同类型指令集和体系架构的计算单元组成计算***，最大程度上与异构体系结构相匹配，充分利用不同加速部件的性能，同时设计了重建算法高效运算的存储与通信方案，支持PCI-E/Ethernet互联，支持多块处理板通过互联总线实现多处理器高效并行处理，实现同步或异步的协同处理机制，在尽量减少损失精度的前提下提高其重建速度，满足用户使用需求。

附图说明：

图1为本发明的方法流程示意图；

图2为FDK算法的几何关系示意图；

图3为FDK算法的总体算粒流程图；

图4为异构平台框架模型；

图5为异构平台软件框图；

图6为本发明的装置示意图；

图7为石油岩心重建结果三维剖面图。

具体实施方式：

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

如图1所示，在一个实施例中，提供一种基于异构平台的CT图像三维重建加速方法。本实施例中异构平台包含主机及异构OpenCL计算设备，所述的主机为运行主程序的CPU，OpenCL计算设备包含运行内核程序的异构容器GPU和FPGA，CPU、GPU及FPGA之间通过PCI-E总线相通信，主程序通过定义上下文来管理内核程序的运行。

FPGA采用逻辑单元阵列LCA(Logic Cell Array)，内部包括：可配置逻辑模块CLB(Configurable Logic Block)、输出输入模块IOB(Input Oouput Block)和内部连线(Interconnect)三个部分。现场可编程门阵列(FPGA)是可编程器件，与传统逻辑电路和门阵列(如PAL,GAL,及CPLD器件)相比，FPGA具有不同的结构，FPGA利用小型查找表(16×1RAM)来实现组合逻辑，每个查找表连接到一个D触发器的输入端，触发器再来驱动其他逻辑电路或驱动I/O，构成既可实现组合逻辑功能又可实现时序逻辑功能的基本逻辑单元模块，这些模块间利用金属连线互相连接或连接到I/O模块。FPGA的逻辑是通过向内部静态存储单元加载编程数据来实现的，存储在存储器单元中的值决定了逻辑单元的逻辑功能以及各模块之间或模块与I/O间的联接方式，并最终决定了FPGA所能实现的功能，FPGA允许无限次的编程。GPU体系结构由一个高度线程化的多核流处理器(Streaming Multiprocess，SM)阵列组成。两个SM形成一个构建块，在基于CUDA技术的GPU的每一代之间，每个构建块中SM的数量可能不同。此外，每个SM又包含了多个流处理器(Streaming Processor，SP)，它们之间共享控制逻辑和指令缓存。每个GPU都带有若干千兆字节(GB)的图形双数据速率(Graphics Double Data Rate，GDDR)DRAM，称为全局存储器(Global Memory)，GPU中这些GDDR DRAM完全不同于CPU体系中安装在主板上的***DRAM，它们主要用于图形处理的帧缓冲区存储器。在图形应用程序中，它们用来保存视频图像和用于3D渲染的纹理信息；而对于计算，它们可以作为带宽芯片外存储器。尽管比典型***存储器的延迟要长，大规模并行应用程序通常高带宽来弥补时延。

海量数据传输是重建加速的难点和重点。PCI-E总线技术是一种高性能通用I/O互连总线，单通道传输速率可达到2.5Gbit/s。因此，本发明采用PCI-E接口作为主要数据传输总线。GPU在通信接口方面已有商业化的PCI-E驱动，因此，FPGA和PCI-E总线间的高速数据传输设计是平台设计的关键。在FPGA上实现PCI-E接口的传输一般有2种方案：一是采用专用的桥接芯片，另一种是采用可实现PCI-E物理接口的特定FPGA。由于目前桥接芯片支持的PCI-E通道数和传输速率较低，故本项目采用后一种方案。数据传输可分为2种模式：普通模式和DMA(Direct Memory Access，直接内存访问)传输模式。其中普通模式主要实现主机和设备之间的指令通信；DMA模式主要针对大块数据传输，传输过程中不需要经过CPU，数据直接从设备传送到内存中，数据传输的速率快，可以充分利用PCI-E数据带宽。因此，FPGA数据采用DMA传输模式。

下面具体描述本实施例的CT图像重建的实现过程：

S1、对FDK重建算法进行算粒分解，分析各个算粒的并行计算流程：

CT图像重建算法主要包括两类：解析类算法和迭代类算法。相对于迭代算法，解析算法数学形式简单，重建速度快，实现方便，是实际CT***中应用的主流算法。在各种基于滤波背投影的解析重建算法当中，FDK算法计算效率高，实现起来容易，而且当锥角较小时，能够取得较好的重建效果，在实际中应用最为广泛。FDK算法几何关系示意图如图2所示，图中S表示射线源，f(x，y，z)表示重建物体的某一点，成为体素，P(θ，u，v)表示在物体旋转分度为θ时探测器采集到f(x，y，z)点的投影数据。根据FDK算法的原理，其主要包括3个步骤：投影加权、一维滤波和反投影。通过对算法的分析可知：每个投影点的加权是相互独立的，每个分度下投影数据是按行进行滤波，滤波过程完全相同，且体素点在每个分度下的反投影过程只与该点映射到探测器坐标系上的坐标有关，具有独立性。也就是说，FDK算法的三个步骤都具有良好的并行性，适合用并行处理的方式进行算法加速。

在一个实施中，如图3所示，根据FDK算法内容，分解为：用于对投影数据进行加权的投影加权算粒，用于对加权后的投影数据进行滤波的滤波算粒，用于将滤波后的投影数据反投影到重建物体上的反投影算粒，及用于对反投影结果进行归约处理的归约算粒。通过分析算粒的存储量、计算量、通信量，使得在异构平台上高效实现。

在本发明的其中一些实施例中，FDK重建算法公式：

通过对公式中积分的拆分及离散化，分为：

投影加权算粒，表示为：其中，p'(θ,u,v)表示旋转分度为θ时对投影数据进行加权后的数据，为加权系数；

滤波算粒，表示为：其中，d_f(θ,u,v)为滤波后的数据，h(u)为滤波算子的单位冲激响应，[-u_m,u_m]表示探测器每一行采集的2m个数据；

反投影算粒，表示为：其中，f(x,y,z,θ)表示重建物体在旋转分度为θ时投影点对f(x,y,z)的贡献值；

归约算粒，表示为：其中，φ_max为重建物体旋转一周时离散采得的投影分度数。

加权系数只与投影射线在探测器平面的相对位置和射线源到原点的距离有关，因此每个投影点下的投影数据的加权具有独立性；根据投影数据的频率滤波公式，每个分度的投影数据在滤波前都需进行傅里叶变换，而每个分度下投影数据的傅里叶变换只与投影数据在面阵探测器上的横坐标Y有关，而频域滤波函数都是相同的，因此投影数据的频域滤波在每一行下具有独立性；体像素点在每个分度下的反投影计算只与该像素点映射到探测器坐标系上的坐标有关，因此反投影计算在每个体素点下是独立的。

综上所述，FDK算法的三个步骤都具有良好的并行性，适合用并行处理的方式进行算法加速。在细粒度的解剖下，FDK算法的总体算粒流程可如图3所示

S2、通过异构平台中的主机及异构OpenCL计算设备对各个算粒进行加速优化处理：

异构平台框架模型如图4所示。根据CT重建对数据传输速率的要求，为了在各异构处理器上高效分配和调度重建任务，异构平台采用PCI-Express作为互连总线，传输数据信号和控制信号。同时，异构平台作为总体平台的一个计算节点，以Ethernet作为附加总线，与外部进行联网控制或数据传输。根据CT重建算法的粒度、负载分布和并行化特性，本发明采用GPU和FPGA作为主要的加速部件。CT采集的投影数据通过Ethernet总线传输进入加速平台节点，根据用户设置和重建性能评估，***通过PCI-E总线进行协调处理，分别在GPU和FPGA加速部件上完成相应的算粒加速，重建数据实时存储并反馈回用户。

在一个实施例中，采用FPGA对投影加权算粒进行并行处理，通过异步传输至GPU，在传输过程中同时对滤波算粒进行处理；结合反投影时各体素点的数据并行操作，在GPU内按体素点对反投影算粒进行多线程并行反投影计算。

进一步地，CT图像的投影数据需要对几何关系进行数据校正；而根据FDK算法原理，其可分为3个步骤：加权、滤波和反投影。因此，需要合理分配三个步骤所消耗的资源，以达到加速性能最优。另一方面，FPGA器件具有可重构性，灵活性强，且适合进行并行流水处理，GPU器件有大量的高性能流处理器，适合数据并行操作模式，在设计时应充分发挥各加速器件的特长。将数据加权校正，滤波和反投影3个算粒分别进行优化实现。数据加权校正是对投影数据进行预处理，由于该算粒只是对投影的行数据进行流水操作，适合在FPGA内部采用流水线的方式进行处理。因此，采用FPGA对数据加权校正进行并行处理，并采用了异步传输的方式，在传输过程中同时流水进行加权和滤波计算，减少了处理延迟，发挥了FPGA流水处理的特点。而加权滤波以及反投影涉及大量的数据运算和存储问题，结合反投影时各体素点的数据并行操作特点，在实现时我们在GPU内按体素点进行多线程并行反投影计算，提升了资源利用率和重建速度。

更进一步，根据FDK重建算法中重建区域在旋转轴方向上行的各层投影数据和探测器投影数据纵轴方向上的各行数据的投影对应关系，采用分块重建策略，沿旋转轴方向将待重建区域分成若干块，重建其中一块时从外部存储器取出对应的投影数据进行重建操作。

当待重建规模很大时，GPU/FPGA板卡上的存储规模可能难以满足一次重建，因此，需采用分块重建的策略。根据FDK算法的公式，重建区域在旋转轴方向上的各层数据和探测器投影数据按纵轴方向上的各行数据有严格的投影对应关系。因此，可沿旋转轴方向将待重建区域分成若干块，重建其中一块时仅需从外部存储器取出对应的投影数据进行重建操作即可。

更进一步，综合考虑FPGA的硬件特点，采用FPGA对投影加权算粒进行并行处理，将全局存储器分割为2块bank，通过负载分配实现随机存取存储器的访问均衡；通过常数存储器存储需要多次重复计算的中间变量。将全局存储器分割为2块bank，通过负载分配以实现2块DDR2SDRAM的访问平衡，从而提升存储器的访问带宽；利用常数存储器存储运算过程中需要多次重复计算的中间变量，节省计算资源；优化内核函数中一次反投影的张数，提升对投影数据存储的访问带宽，同时减少对重建数据存储的访问，通过调整以达到对全局存储器访问的最优。

更进一步，在GPU内按体素点对反投影算粒进行多线程并行反投影计算，采用基于体素驱动方式，按重建体数据输出对GPU进行任务划分；将计算中与体素无关的变量进行分离和合并，并在反投影前计算出来存储于GPU的常数存储器中，反投影计算时，直接读取常数存储器中的变量参与计算；优化内核程序中一次反投影的张数。

GPU的任务划分可以按照输入划分或者按输出划分。对于反投影算法来说，输入为投影数据，输出为重建体数据，GPU的两种任务分配方式实质上反映了两种不同的反投影实现方式：基于射线驱动和基于体素驱动。基于射线驱动的方式按投影数据进行任务划分，一个或几个线程完成一条射线上所有体素的反投，反投时首先计算出当前射线穿过的所有体素，然后将这些体素的值赋为当前射线投影值，由于不同射线对应的线程或同一射线对应的不同线程可能同时给一个体素赋值，即这种方式存在“写竞争”问题；基于体素驱动的方式按体数据划分任务，一个线程完成一个或几个体素的反投，反投时首先计算出当前体素的正投影位置，然后取该位置的投影值赋给当前体素。基于体素驱动的反投影方式不存在“写竞争”问题，不需要设计额外的归约步骤，因此这种任务分配方式更适合于GPU加速。

本实施例中，采用基于体素的驱动方式，即按输出对任务划分。在线程分配时需要考虑一个很重要的原则，即SM的占用率不能太低，SM占用率指的是每个SM上的活动Warp数量与允许的最大活动Warp数量之比。因为GPU是通过线程间的切换来隐藏长延时操作(访问全局存储器、线程间同步等)的，当一个Warp中的线程进行长延时操作时，另外一个活动Warp中的线程会自动进入运算状态，这样便可以隐藏一部分延迟。但这并不代表SM占用率越高越好，SM占用率越高则每个线程占用的GPU资源越少，即每个线程完成的计算量越少，而每个SM上的最大活动Warp数量是一定的，因此就可能出现即使SM中活动的Warp数量达到最大，由于Warp中每个线程计算量太少，从而使所有活动Warp中线程同时进入长延时操作，不能充分地隐藏延迟。因此，在SM占用率和每个线程完成的计算量间选择一个平衡点，才能使GPU的性能发挥到最佳。经过多次实验，可采取如下的线程分配方式，假设体数据的规模为N3，则Block的大小恒定为(16,16,2)，Grid的尺寸随N的改变而变化为(N/16,N/16,1)。

反投影需要计算体数据中任意一点(x，y，z)在投影角度为θ时投射到探测器上的点(u，v)，在这个运算过程中需要多次计算只与投影角度和几何关系有关的三角函数和其他中间变量。对于每个体素，这些相同的变量都被计算一次，假设体数据规模为N³，则这些值会被重复计算N³次，造成***资源的极大浪费。针对这个问题，在一个实施例中，将运算中与体素(x,y,z)无关的变量进行分离和合并，并在反投影前计算出来存于GPU的常数存储器中，反投影运算时，直接读取常数存储器中的变量参与计算。

进一步地，将计算中与体素无关的变量进行分离和合并，包含内容如下：体数据中任意一点(x，y，z)在投影角度为θ时投射到探测器上的点(u，v)，投影点(u,v)的计算为：

u＝(x-vCenter)×cos(θ)+(y-vCenter)sin(θ)+pCenter

dis＝(u-pCenter)×a

v＝(z-(s₀+θ×h)-γ×h/a)×w+γ×h/a+pCenter

，分离和合并变量后为:

u＝x×A[0]+y×A[1]+A[2]+pCenter

dis＝(u-pCenter)×a

v＝(z-A[4]-γ×A[5])×w+γ×A[5]+pCenter

，其中，vCenter表示体数据中心，pCenter为投影数据中心，α为体素大小，θ为投影角度，r为旋转射线源旋转半径，h为螺距，γ为射束在中心平面上投影与中心射束的夹角。

分离和合并变量后每一个角度的反投影运算可以提取出6个与体素(x，y，z)无关的中间变量，假设投影数为360，则整个反投影共有2160个变量(单精度浮点型)需要在反投影前计算出来并存于GPU常数存储器中。常数存储器是GPU特有的只读存储空间，能够被缓存，并且来自同一half-Warp的线程访问常数存储器中的同一数据时，如果发生缓存命中，只需要一个周期就可以获得数据。一般来说，常数存储器空间较小，如Tesla C1060中只有64KB，但完全能够满足本实施例中的需要。反投影时，只需要读取常数存储器中的变量值参与常规的乘加运算便可得到最终的反投影参数。因此，这种加速策略不仅可以避免用GPU去计算开销极大的三角函数，而且避免了GPU的重复计算，在提升反投影运算效率方面具有较好效果。

重建体数据通常存放于GPU全局存储器内，全局存储器占据了显存绝大部分，可以用来存放大规模的数据，但全局存储器没有缓存加速，虽然合并访问方式可以极大地提高访问速度，但通常仍然会有几百个周期的访问延迟。研究表明，GPU用于高性能计算的瓶颈不是计算消耗而是访存消耗。因此，如何缩减访问全局存储器的时间是GPU加速的关键。本实施例中，设计了在一个Kernel同时反投m幅投影图像的加速策略。

通常情况下，一个Kernel中只对1幅投影图像进行运算，对于360幅投影数据，整个反投影过程需要读/写全局存储器360×N³次。实施例中，一个Kernel完成m个角度投影图像的反投影，每个Kernel需要计算m个角度投影图像的反投影参数，但只读写全局存储器N3次，即可以将读/写全局存储器的次数变为原来的1/m。在减少全局存储器读写次数的同时，算法会增加Kernel中每个线程的计算负担，若一味的增大m，势必会降低整个GPU中活动Block和活动Warp的数量，活动Block和活动Warp数量减少又会反过来影响GPU隐藏长延时操作(访问全局存储器)的优势。所以，寻找一个适中的m显得尤为重要。通过多次试验发现，当一个Kernel中同时完成的反投影图像为3幅时，双方达到一个平衡，加速效果最为理想。

在另一个实施例中，如图6所示，提供一种基于异构平台的CT图像三维重建加速装置，异构平台采用PCI-Express作为传输数据信号和控制信号的互连总线，并以Ethernet作为与外部进行联网控制和数据传输的附加总线；异构平台的架构包含提供功能模块的应用层，提供面向应用层各功能模块基于不同处理器的构件库及重建算法所需接口规范的构件层，及面向构件层和应用层提供服务的支撑层，支撑层包含执行主程序的CPU及执行内核程序的多个OpenCL计算设备，CPU和OpenCL计算设备相通信连接，所述的OpenCL计算设备包含GPU、FPGA；基于所述的异构平台的架构的CT图像三维重建加速装置包含如下内容：

算粒分解模块201，用于根据FDK重建算法内容将该算法以算粒形式进行拆分和离散化，分解为用于对投影数据进行数据校正的投影加权算粒，用于对加权后的投影数据进行滤波的滤波算粒，用于将滤波后的投影数据反投影到重建物体上的反投影算粒，及用于对反投影结果进行归约处理的归约算粒；

加速处理模块202，用于对通过附加总线传输至异构平台计算节点上的投影数据，根据用户设置和重建性能评估，通过互连总线进行协调处理，分别在GPU和FPGA加速部件上完成算粒加速，重建数据实时存储并反馈至用户。

需要说明的，本发明的其中一些实施例的CT图像三维重建加速装置的具体实现过程与CT图像三维重建加速方法部分相同，具体可参见方法部分实施例，这里不再赘述。

异构平台的软件设计框图如图5所示。软件架构共分为三层，其中应用层主要为用户软件所具备的功能模块，分别为投影数据加权校正、投影数据滤波、三维反投影重建和重建图像可视化4个模块；构件层为面向应用的基于不同处理器的构件库，主要为重建算法各算粒的函数代码库与相应的接口规范；最后一层为CPU、GPU、FPGA等容器组成的支撑层，为构件和应用软件提供服务和支撑。

在构件层设计时，拟采用OpenCL框架进行开发。OpenCL全称为Open ComputingLanguage，即开放式计算语言，最早是由Apple公司提出，是一种全新的计算应用程序接口(API)。OpenCL的主要作用是为通用计算领域提供一个跨平台的统一标准语言，它所支持的异构平台可由多核CPU、GPU、DSP、Cell/B.E.processor或其他类型的处理器组成。OpenCL为并行程序的开发提供了跨厂商的非专有软件解决方案，使得程序具备较好的可移植性；同时，跨平台的异构框架利于同时发挥***中各种设备的性能潜力。OpenCL平台模型由主机(Host)及与其相连的一个或多个OpenCL计算设备(Compute Device)构成。其中，每个计算设备由一个或多个计算单元(Computer Unit)，每个计算单元又可以进一步分割成一个或多个处理单元(Processing Element)，各种计算操作都是在处理单元中完成的。主机端管理着整个平台上的所有计算资源。应用程序会从主机端向各个OpenCL设备的处理单元发送计算命令。在一个计算单元内的所有处理单元会执行完全相同的一套指令流程。异构平台的并行编程实现中，将CPU作为主机，并将GPU和FPGA作为OpenCL设备。OpenCL的编程模型可以分为两部分，一部分是在CPU上执行的主程序(Host program)，另一部分是在Device上执行的核函数(Kernel)。主程序通过定义上下文(Context)并管理内核程序在Device上的执行。在OpenCL编程中，在主机端创建一个核函数之前必须先为该核函数创建一个索引空间(Index Space)，该索引空间可以是一维、二维或者三维的，核函数会在该索引空间的每个节点(Work Item)上执行。每个工作节点相应维度上的索引都被定义为节点在该维度上的全局ID(Global ID)。所有工作节点都将执行相同的核函数程序，各工作节点即等价于不同的执行线程，通过大量节点间的并行程序执行从而达到计算加速的目的。OpenCL在对索引空间提供了全局索引之外也提供了较小粒度的工作组(work group)空间。每个工作组在都有一个唯一的工作组ID(Work Group ID)，工作组内部的节点相对该工作组的位置索引被称作局部ID(Local ID)。核函数的设计既可以选择工作节点间的并行模式，也可以选择两层并行模式，即索引空间中的工作组间并行和工作组中内部节点间并行。这既可充分利用GPU和FPGA底层硬件的计算资源，也增加了程序设计的灵活性。

在多处理异构和OpenCL编程环境下，对CPU主机端的编程采用标准C/C++语言，对GPU和FPGA协处理器的编程采用基于OpenCL规范的描述语言。OpenCL标准规范的编程语言的抽象等级远高于VHDL和Verilog等硬件描述语言。传统的编程方式需要对FPGA底层硬件单元按照时序周期进行编程描述，对于复杂的算法执行，需要设计状态机控制数据通路，同时需要处理各级接口约束和时序同步问题，编程难度大、耗时长，且程序维护和升级复杂，非常不利于实际产品的快速应用。而采用OpenCL编程方式，无需关注底层时序级的硬件设计，可根据反投影算法设计高级语言描述的类C代码，OpenCL编译器则可自动实现由OpenCL代码转化为硬件描述语言和执行程序的步骤。

为进一步验证本发明的有效性，下面通过具体实验作进一步解释说明：

如图7所示，选择石油岩心作为测试对象，对被测物体采集全角度下360张投影进行重建，重建规模为500³，重建结果中间三维剖面图7所示：a)、b)、c)分别为重建的三维图像中每一个中心平面上的剖面图，即对于世界坐标系上的x＝0、y＝0、z＝0三个平面上的切面图。

针对不同重建规模分别采用“CPU”、“GPU+FPGA”两种方式进行重建测试。每组测试重复10次，平均重建时间测试结果如下：

表1重建时间测试结果

从结果中我们可以看出，采用GPU+FPGA的加速方式可以大幅度提升传统CPU计算的速度，使得用户可以快速得到三维重建结果。同时在重建结果图像中我们可以看出所得重建结果可以清晰显示三维内部信息，符合三维重建实际应用。

以上各个实施例在具体说明中仅只针对相应步骤的实现方式进行了阐述，然后在逻辑不相矛盾的情况下，上述各个实施例是可以相互组合而形成新的技术方案的，而该新的技术方案依然在本具体实施方式的公开范围内。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例方法可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件，但很多情况下前者是更佳的实施方式。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品承载在一个非易失性计算机可读存储载体(如ROM、磁碟、光盘、服务器云空间)中，包括若干指令用以使得一台终端设备(可以是手机，计算机，服务器，或者网络设备等)执行本发明各个实施例所述的方法。

以上所述实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

Claims (10)

1.一种基于异构平台的CT图像三维重建加速方法，其特征在于，异构平台包含主机及异构OpenCL计算设备，该加速方法包含如下内容：对FDK重建算法进行算粒分解，分析各个算粒的并行计算流程；通过异构平台中的主机及异构OpenCL计算设备对各个算粒进行加速优化处理。

2.根据权利要求1所述的基于异构平台的CT图像三维重建加速方法，其特征在于，所述的主机为运行主程序的CPU，OpenCL计算设备包含运行内核程序的异构容器GPU和FPGA，CPU、GPU及FPGA之间通过PCI-E总线相通信，主程序通过定义上下文来管理内核程序的运行。

3.根据权利要求2所述的基于异构平台的CT图像三维重建加速方法，其特征在于，对FDK重建算法进行算粒分解，包含：根据FDK算法内容，分解为：用于对投影数据进行加权的投影加权算粒，用于对加权后的投影数据进行滤波的滤波算粒，用于将滤波后的投影数据反投影到重建物体上的反投影算粒，及用于对反投影结果进行归约处理的归约算粒。

4.根据权利要求3所述的基于异构平台的CT图像三维重建加速方法，其特征在于，根据FDK重建算法公式：

通过对公式中积分的拆分及离散化，分为：

投影加权算粒，表示为：其中，p'(θ,u,v)表示旋转分度为θ时对投影数据进行加权后的数据，为加权系数；

滤波算粒，表示为：其中，d_f(θ,u,v)为滤波后的数据，h(u)为滤波算子的单位冲激响应，[-u_m,u_m]表示探测器每一行采集的2m个数据；

反投影算粒，表示为：其中，f(x,y,z,θ)表示重建物体在旋转分度为θ时投影点对f(x,y,z)的贡献值；

归约算粒，表示为：其中，φ_max为重建物体旋转一周时离散采得的投影分度数。

5.根据权利要求3所述的基于异构平台的CT图像三维重建加速方法，其特征在于，对各个算粒进行加速优化处理，包含：采用FPGA对投影加权算粒进行并行处理，通过异步传输至GPU，在传输过程中同时对滤波算粒进行处理；结合反投影时各体素点的数据并行操作，在GPU内按体素点对反投影算粒进行多线程并行反投影计算。

6.根据权利要求5所述的基于异构平台的CT图像三维重建加速方法，其特征在于，根据FDK重建算法中重建区域在旋转轴方向上行的各层投影数据和探测器投影数据纵轴方向上的各行数据的投影对应关系，采用分块重建策略，沿旋转轴方向将待重建区域分成若干块，重建其中一块时从外部存储器取出对应的投影数据进行重建操作。

7.根据权利要求5所述的基于异构平台的CT图像三维重建加速方法，其特征在于，采用FPGA对投影加权算粒进行并行处理，包含：将全局存储器分割为2块bank，通过负载分配实现随机存取存储器的访问均衡；通过常数存储器存储需要多次重复计算的中间变量。

8.根据权利要求5所述的基于异构平台的CT图像三维重建加速方法，其特征在于，在GPU内按体素点对反投影算粒进行多线程并行反投影计算，包含：采用基于体素驱动方式，按重建体数据输出对GPU进行任务划分；将计算中与体素无关的变量进行分离和合并，并在反投影前计算出来存储于GPU的常数存储器中，反投影计算时，直接读取常数存储器中的变量参与计算；优化内核程序中一次反投影的张数。

9.根据权利要求8所述的基于异构平台的CT图像三维重建加速方法，其特征在于，将计算中与体素无关的变量进行分离和合并，包含内容如下：体数据中任意一点(x，y，z)在投影角度为θ时投射到探测器上的点(u，v)，投影点(u,v)的计算为：

u＝(x-vCenter)×cos(θ)+(y-vCenter)sin(θ)+pCenter

dis＝(u-pCenter)×a

<mrow> <mi>w</mi> <mo>=</mo> <mfrac> <msqrt> <mrow> <msup> <mi>r</mi> <mn>2</mn> </msup> <mo>-</mo> <msup> <mi>dis</mi> <mn>2</mn> </msup> </mrow> </msqrt> <mrow> <msqrt> <mrow> <msup> <mi>r</mi> <mn>2</mn> </msup> <mo>-</mo> <msup> <mi>dis</mi> <mn>2</mn> </msup> </mrow> </msqrt> <mo>+</mo> <mi>a</mi> <mo>&amp;times;</mo> <mrow> <mo>(</mo> <mo>-</mo> <mo>(</mo> <mrow> <mi>x</mi> <mo>-</mo> <mi>v</mi> <mi>C</mi> <mi>e</mi> <mi>n</mi> <mi>t</mi> <mi>e</mi> <mi>r</mi> </mrow> <mo>)</mo> <mo>&amp;times;</mo> <mi>s</mi> <mi>i</mi> <mi>n</mi> <mo>(</mo> <mi>&amp;theta;</mi> <mo>)</mo> <mo>+</mo> <mo>(</mo> <mrow> <mi>y</mi> <mo>-</mo> <mi>v</mi> <mi>C</mi> <mi>e</mi> <mi>n</mi> <mi>t</mi> <mi>e</mi> <mi>r</mi> </mrow> <mo>)</mo> <mi>c</mi> <mi>o</mi> <mi>s</mi> <mo>(</mo> <mi>&amp;theta;</mi> <mo>)</mo> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow> ，

v＝(z-(s₀+θ×h)-γ×h/a)×w+γ×h/a+pCenter

分离和合并变量后为:

u＝x×A[0]+y×A[1]+A[2]+pCenter

dis＝(u-pCenter)×a

<mrow> <mi>w</mi> <mo>=</mo> <mfrac> <msqrt> <mrow> <msup> <mi>r</mi> <mn>2</mn> </msup> <mo>-</mo> <msup> <mi>dis</mi> <mn>2</mn> </msup> </mrow> </msqrt> <mrow> <msqrt> <mrow> <msup> <mi>r</mi> <mn>2</mn> </msup> <mo>-</mo> <msup> <mi>dis</mi> <mn>2</mn> </msup> </mrow> </msqrt> <mo>+</mo> <mi>a</mi> <mo>&amp;times;</mo> <mrow> <mo>(</mo> <mo>-</mo> <mi>x</mi> <mo>&amp;times;</mo> <mi>A</mi> <mo>&amp;lsqb;</mo> <mn>1</mn> <mo>&amp;rsqb;</mo> <mo>+</mo> <mi>y</mi> <mo>&amp;times;</mo> <mi>A</mi> <mo>&amp;lsqb;</mo> <mn>0</mn> <mo>&amp;rsqb;</mo> <mo>+</mo> <mi>A</mi> <mo>&amp;lsqb;</mo> <mn>3</mn> <mo>&amp;rsqb;</mo> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow> ，

v＝(z-A[4]-γ×A[5])×w+γ×A[5]+pCenter

其中，vCenter表示体数据中心，pCenter为投影数据中心，α为体素大小，θ为投影角度，r为旋转射线源旋转半径，h为螺距，γ为射束在中心平面上投影与中心射束的夹角。

10.一种基于异构平台的CT图像三维重建加速装置，其特征在于，异构平台采用PCI-Express作为传输数据信号和控制信号的互连总线，并以Ethernet作为与外部进行联网控制和数据传输的附加总线；异构平台的架构包含提供功能模块的应用层，提供面向应用层各功能模块基于不同处理器的构件库及重建算法所需接口规范的构件层，及面向构件层和应用层提供服务的支撑层，支撑层包含执行主程序的CPU及执行内核程序的多个OpenCL计算设备，CPU和OpenCL计算设备相通信连接，所述的OpenCL计算设备包含GPU、FPGA；基于所述的异构平台的架构的CT图像三维重建加速装置包含如下内容：

算粒分解模块，用于根据FDK重建算法内容将该算法以算粒形式进行拆分和离散化，分解为用于对投影数据进行数据校正的投影加权算粒，用于对加权后的投影数据进行滤波的滤波算粒，用于将滤波后的投影数据反投影到重建物体上的反投影算粒，及用于对反投影结果进行归约处理的归约算粒；

加速处理模块，用于对通过附加总线传输至异构平台计算节点上的投影数据，根据用户设置和重建性能评估，通过互连总线进行协调处理，分别在GPU和FPGA加速部件上完成算粒加速，重建数据实时存储并反馈至用户。