CN105374006A - 基于遗传算法的ct图像重建反投影加速方法 - Google Patents

Abstract

本<b>发明</b>公开了一种<b>基于遗传算法的</b><b>CT</b>图像<b>重建反投影加速方法</b>，克服了现有技术中，利用CUDA加速CT重建反投影的实现方式存在访存带宽受限、自适应程度低等问题。该发明其在CUDA编程框架下，基于CUDA的反投影并行计算策略如下：（1）一次反投影核函数运算处理多个投影；（2）合理设计线程网格（grid）和线程块（block）的维度和顺序来保证全局存储器的合并访问；（3）使用常数存储器保存反投影中的常数；（4）使用纹理存储器实现双线性插值，减少算术指令的运算时间；（5）约简反投影运算中的相同因子。本发明不需要人为来识别，加速策略的自适应程度高。

Description

基于遗传算法的CT图像重建反投影加速方法

技术领域

该发明涉及一种CT图像重建的方法，特别是涉及一种基于遗传算法的CT图像重建反投影加速方法。

背景技术

计算机断层成像(ComputedTomography,CT)以其无损、精确、三维可视化等优点在工业无损检测和医疗诊断领域获得了广泛的应用。然而CT设备在提高成像质量和重建分辨率的同时，相应需求的计算资源和存储资源都非常大，重建速度已成为CT技术走向实用的一个瓶颈，如何对CT重建过程加速是目前该领域研究的一个难点问题。

目前，实际CT系统中应用的主流算法主要集中于以滤波反投影为基础的解析类重建算法上，如FBP、FDK等算法。在这类算法中，反投影过程所占计算消耗的比例最高，对反投影计算的加速是实际应用中重建加速的关键点。

在半导体工艺和并行计算技术进步的推动下，特别是以GPU为代表的众核处理器的推出，使得图像重建速度已不再是一个瓶颈。在众多高性能计算方案中，基于CUDA的图像重建算法加速效果十分显著，利用CUDA实现对CT重建的反投影过程进行加速是该领域的重要研究方向。

经过现有技术的文献检索发现，利用CUDA加速CT重建反投影的文章采取的实现方式存在访存带宽受限、自适应程度低等不足。2011年，Papenhausen在WorkshoponHighPerformanceImageReconstruction上发表《GPU-acceleratedback-projectionrevisited:squeezingperformancebycarefultuning》，通过常数存储器、纹理存储器和一次反投影多个投影等更精细的优化策略，使得GPU在图像重建中的优势得到更好的发挥。2014年，Zhang等人在AppliedMechanicsandMaterials上发表《ComparisonofParallelComputingMethodsforFastCone-BeamReconstructionwithSimilarOptimizationStrategies》，充分利用了GPU受限于访存带宽而不是计算带宽的特点，通过一次反投影更多投影的方式，对于5123的重建图像仅需要2.1秒。但上述方法的性能瓶颈是都需要人为来识别，加速策略的自适应程度低。

发明内容

本发明克服了现有技术中，利用CUDA加速CT重建反投影的实现方式存在访存带宽受限、自适应程度低等问题，提供一种自适应程度高的基于遗传算法的CT图像重建反投影加速方法。

本发明的技术解决方案是，提供一种具有以下步骤的基于遗传算法的CT图像重建反投影加速方法：其在CUDA编程框架下，基于CUDA的反投影并行计算策略如下：

(1)一次反投影核函数运算处理多个投影，减少重建图像体素值更新时对全局存储器的访问次数；

(2)合理设计线程网格(grid)和线程块(block)的维度和顺序来保证全局存储器的合并访问，降低重建图像体素值更新时对全局存储器的访问延迟；

(3)使用常数存储器保存反投影中的常数，比如三角函数值，减少算术指令的运算时间；

(4)使用纹理存储器实现双线性插值，减少算术指令的运算时间；

(5)约简反投影运算中的相同因子，如相同的除法，减少算术指令的运算时间。

所述定义策略(1)中一次反投影处理I个投影，策略(2)中block的维度为(Bx,By)、一个线程更新K个体素，建立以(I,K,Bx,By)为参数、以反投影计算时间T(I,K,B_x,B_y)为代价函数的最优化模型

${(I,K,B_x,B_y)}^{opt}=\underset{(I,K,Bx,By)}{\arg }\min T(I,K,B_x,B_y),$

$s.t.\left\{\begin{array}{c}I,K,B_x,B_y\&Element;N\\ I\&Element;\&lsqb;1,min(N_p,N_{layer},N_m)\&rsqb;,K\&Element;\&lsqb;1,N_z\&rsqb;\\ B_x\&Element;\&lsqb;1,min(N_x,N_t)\&rsqb;,B_y\&Element;\&lsqb;1,min(N_y,N_t)\&rsqb;\\ B_x*B_y\&Element;\&lsqb;1,N_t\&rsqb;,\mathrm{mod}(B_x*B_y,32)=0\end{array}\right..$

公式中N代表自然数空间，(N_x,N_y,N_z)为重建图像沿三个坐标轴方向的体素数量，N_p为CT采集的投影数量，N_layer为GPU硬件对于多层纹理层数的上限，N_m为GPU全局存储器容量存储投影数量的上限，N_t为GPU硬件对于block中所包含线程数量的上限，约束mod(B_x*B_y,32)＝0是为了保证block中的线程数量为warp(包含32个线程)的整数倍。

所述最优化模型采用SGA求解步骤如下：

步骤1：随机生成M个个体(I,K,Bx,By)作为初始种群，并对每个个体进行编码。由于个体(I,K,Bx,By)全部是整数，因此直接采用每个参数的二进制表示作为基因编码，并将其串联在一起组成基因序列；

步骤2：编译每个个体对应的反投影核函数为动态链接库，参数最优化搜索程序调用该动态链接库并记录反投影时间，从而得到每个个体的适应度；

步骤3：当终止条件未达到时，执行下述操作：

子步骤3.1：基于轮盘赌算法执行选择与复制算子；

子步骤3.2：以概率Pc执行交叉算子；

子步骤3.3：以概率Pm执行变异算子；

子步骤3.4：验证新产生的个体是否满足模型的约束条件。如果满足，则执行下一步；否则重新执行步骤3.2和3.3，直至满足约束条件；

子步骤3.5：对每个个体的基因编码进行解码得到(I,K,Bx,By)，然后重复步骤2，计算新一代种群中每个个体的适应度；

步骤4：当终止条件达到时，输出适应度最小的个体，并对其进行解码得到最优的配置参数。

本方法其在CUDA编程框架下，基于CUDA的反投影并行计算分为三个优化部分：

优化1：在第一次反投影核函数运算之前，首先完成I个投影的数据上传，然后在核函数运算的同时上传下次需要处理的投影数据；

优化2：在接下来的反投影中，GPU在对上次上传的投影数据进行反投影的同时，上传下次需要处理的投影数据；

优化3：对于最后一次反投影，将重建图像分成均等的两个部分，采用两个核函数分别完成对应的反投影，当第一个核函数完成反投影之后即开始下载数据到内存，同时第二个核函数计算另一半重建图像，待反投影结束后下载数据到内存。

与现有技术相比，本发明基于遗传算法的CT图像重建反投影加速方法具有以下优点：本发明反投影过程所占计算消耗的比例低，通过常数存储器、纹理存储器和一次反投影多个投影等更精细的优化策略，使得GPU在图像重建中的优势得到更好的发挥；不需要人为来识别，加速策略的自适应程度高。提出的流水处理方法优点是时序控制简单，容易实现。

附图说明

图1是本发明基于遗传算法的CT图像重建反投影加速方法的加速流程图；

图2是本发明基于遗传算法的CT图像重建反投影加速方法采用SGA来求解模型的算法流程图；

图3是本发明基于遗传算法的CT图像重建反投影加速方法中反投影流水处理方法的示意图；

图4是本发明基于遗传算法的CT图像重建反投影加速方法中RabbitCT的(a)xy切片、(b)xz切片、(c)yz切片、(d)三维渲染图像的重建图像。

具体实施方式

下面结合附图和具体实施方式对本发明基于遗传算法的CT图像重建反投影加速方法作进一步说明：一种基于遗传算法的CT图像重建反投影加速方法，在CUDA编程框架下，基于CUDA的反投影并行计算策略包括以下几点：(1)一次反投影核函数运算处理多个投影，减少重建图像体素值更新时对全局存储器的访问次数；(2)合理设计线程网格(grid)和线程块(block)的维度和顺序来保证全局存储器的合并访问，降低重建图像体素值更新时对全局存储器的访问延迟；(3)使用常数存储器保存反投影中的常数(比如三角函数值等)，减少算术指令的运算时间；(4)使用纹理存储器实现双线性插值，减少算术指令的运算时间；(5)约简反投影运算中的相同因子(如相同的除法)，减少算术指令的运算时间。

对于不同的GPU平台，上述并行计算策略中的(3)-(5)都是不变的，但策略(1)和(2)需要调整其中的配置参数才能使反投影并行计算达到最优。目前，调整策略(1)和(2)中的配置参数完全由人工尝试，若有一种方法能够使得策略(1)和(2)自动地随GPU硬件的变化而改变，则图像重建算法在GPU上的最优并行计算就可以达到自适应的目的，而不必人工进行干预。

为此，本发明提出以反投影的配置参数为变量、以反投影时间为目标函数的最优化模型，然后通过遗传算法来搜索最优的配置参数，从而实现反投影的自适应计算。

本方法确定的最优化模型如下：定义策略(1)中一次反投影处理I个投影，策略(2)中block的维度为(Bx,By)、一个线程更新K个体素。本发明建立以(I,K,Bx,By)为参数、以反投影计算时间T(I,K,B_x,B_y)为代价函数的最优化模型

${(I,K,B_x,B_y)}^{opt}=\underset{(I,K,Bx,By)}{\arg }\min T(I,K,B_x,B_y),$

$s.t.\left\{\begin{array}{c}I,K,B_x,B_y\&Element;N\\ I\&Element;\&lsqb;1,min(N_p,N_{layer},N_m)\&rsqb;,K\&Element;\&lsqb;1,N_z\&rsqb;\\ B_x\&Element;\&lsqb;1,min(N_x,N_t)\&rsqb;,B_y\&Element;\&lsqb;1,min(N_y,N_t)\&rsqb;\\ B_x*B_y\&Element;\&lsqb;1,N_t\&rsqb;,\mathrm{mod}(B_x*B_y,32)=0\end{array}\right..---\left(1\right)$

式(1)中N代表自然数空间，(N_x,N_y,N_z)为重建图像沿三个坐标轴方向的体素数量，N_p为CT采集的投影数量，N_layer为GPU硬件对于多层纹理层数的上限，N_m为GPU全局存储器容量存储投影数量的上限，N_t为GPU硬件对于block中所包含线程数量的上限，约束mod(B_x*B_y,32)＝0是为了保证block中的线程数量为warp(包含32个线程)的整数倍。

最优化模型求解如下：目前，遗传算法有许多变种，但由Goldberg提出的基本遗传算法(SimpleGeneticAlgorithm，SGA)是许多算法的原型。SGA包含选择与复制算子、交叉算子和变异算子，具有遗传进化过程操作简单的优点。

本发明采用SGA来求解模型(1)的计算步骤可以描述为：

步骤1：随机生成M个个体(I,K,Bx,By)作为初始种群，并对每个个体进行编码。由于个体(I,K,Bx,By)全部是整数，因此直接采用每个参数的二进制表示作为基因编码，并将其串联在一起组成基因序列；

步骤2：编译每个个体对应的反投影核函数为动态链接库，参数最优化搜索程序调用该动态链接库并记录反投影时间，从而得到每个个体的适应度；

步骤3：当终止条件未达到时，执行下述操作：

子步骤3.1：基于轮盘赌算法执行选择与复制算子；

子步骤3.2：以概率Pc执行交叉算子；

子步骤3.3：以概率Pm执行变异算子；

子步骤3.4：验证新产生的个体是否满足模型的约束条件。如果满足，则执行下一步；否则重新执行步骤3.2和3.3，直至满足约束条件；

子步骤3.5：对每个个体的基因编码进行解码得到(I,K,Bx,By)，然后重复步骤2，计算新一代种群中每个个体的适应度；

步骤4：当终止条件达到时，输出适应度最小的个体，并对其进行解码得到最优的配置参数。

通过求解模型(1)可以得到使反投影核函数在GPU上运行时间最少的配置参数，但整个反投影运算的时间还应包括投影数据从内存传入显存和重建图像从显存传回内存的时间，因此下面介绍隐藏数据传输时间的方法：考虑到一次数据上传的时间远小于一次反投影核函数运算的时间，本发明设计了一种新的隐藏数据传输时间的流水处理，如图3所示，其在CUDA编程框架下，基于CUDA的反投影并行计算分为三个优化部分：

优化1：在第一次反投影核函数运算之前，首先完成I个投影的数据上传，然后在核函数运算的同时上传下次需要处理的投影数据；

优化2：在接下来的反投影中，GPU在对上次上传的投影数据进行反投影的同时，上传下次需要处理的投影数据；

优化3：对于最后一次反投影，将重建图像分成均等的两个部分，采用两个核函数分别完成对应的反投影，当第一个核函数完成反投影之后即开始下载数据到内存，同时第二个核函数计算另一半重建图像，待反投影结束后下载数据到内存。

本发明提出的流水处理方法优点是时序控制简单，容易实现。本方法基于RabbitCT平台，设定重建图像矩阵为5123，体素边长为0.5mm，数据类型为float。设置遗传算法的初始种群规模M＝200，交叉算子中的交叉概率Pc＝0.75，变异算子中的变异概率Pm＝0.05，最大遗传代数Ng＝500，ediff＝0.01s。最优化模型中的N_m根据各GPU的显存容量和RabbitCT中单个投影的数据量计算得到。本方法使用的GPU为NVIDIA公司的TeslaK40c，实验平台的其它参数为：2颗IntelXeonE5-26208核CPU，主频2GHz；24GB内存；Windows764位操作系统；程序开发集成环境为visualstudio2010，C语言；CUDA6.564位程序开发包。将上述参数代入本发明算法的程序运行10次，取其中最优的个体作为最优的配置参数，得到重建图像如图4所示，其中(a)-(c)图的显示窗口为[-2905HU，10300HU]，(d)图为使用VGstudio软件渲染的三维图像，在TeslaK40c上测试结果如表1所示。

表1反投影最优配置参数的计算结果

GPU	I	K	Bx	By	反投影时间/ms
						K40c	62	4	128	8	954

Claims (4)

1.一种基于遗传算法的CT图像重建反投影加速方法，其特征在于：其在CUDA编程框架下，基于CUDA的反投影并行计算策略如下：

(1)一次反投影核函数运算处理多个投影，减少重建图像体素值更新时对全局存储器的访问次数；

(2)合理设计线程网格(grid)和线程块(block)的维度和顺序来保证全局存储器的合并访问，降低重建图像体素值更新时对全局存储器的访问延迟；

(3)使用常数存储器保存反投影中的常数，比如三角函数值，减少算术指令的运算时间；

(4)使用纹理存储器实现双线性插值，减少算术指令的运算时间；

(5)约简反投影运算中的相同因子，如相同的除法，减少算术指令的运算时间。

2.根据权利要求1所述的基于遗传算法的CT图像重建反投影加速方法，其特征在于：所述定义策略(1)中一次反投影处理I个投影，策略(2)中block的维度为(Bx,By)、一个线程更新K个体素，建立以(I,K,Bx,By)为参数、以反投影计算时间T(I,K,B_x,B_y)为代价函数的最优化模型

${(I,K,B_x,B_y)}^{opt}=\underset{(I,K,Bx,By)}{\arg }\min T(I,K,B_x,B_y),$

$s.t.\left\{\begin{array}{c}I,K,B_x,B_y\&Element;N\\ I\&Element;\&lsqb;1,min(N_p,N_{layer},N_m)\&rsqb;,K\&Element;\&lsqb;1,N_z\&rsqb;\\ B_x\&Element;\&lsqb;1,min(N_x,N_t)\&rsqb;,B_y\&Element;\&lsqb;1,min(N_y,N_t)\&rsqb;\\ B_x*B_y\&Element;\&lsqb;1,N_t\&rsqb;,\mathrm{mod}(B_x*B_y,32)=0\end{array}\right..$

公式中N代表自然数空间，(N_x,N_y,N_z)为重建图像沿三个坐标轴方向的体素数量，N_p为CT采集的投影数量，N_layer为GPU硬件对于多层纹理层数的上限，N_m为GPU全局存储器容量存储投影数量的上限，N_t为GPU硬件对于block中所包含线程数量的上限，约束mod(B_x*B_y,32)＝0是为了保证block中的线程数量为warp(包含32个线程)的整数倍。

3.根据权利要求2所述的基于遗传算法的CT图像重建反投影加速方法，其特征在于：所述最优化模型采用SGA求解步骤如下：

步骤1：随机生成M个个体(I,K,Bx,By)作为初始种群，并对每个个体进行编码；由于个体(I,K,Bx,By)全部是整数，因此直接采用每个参数的二进制表示作为基因编码，并将其串联在一起组成基因序列；

步骤2：编译每个个体对应的反投影核函数为动态链接库，参数最优化搜索程序调用该动态链接库并记录反投影时间，从而得到每个个体的适应度；

步骤3：当终止条件未达到时，执行下述操作：

子步骤3.1：基于轮盘赌算法执行选择与复制算子；

子步骤3.2：以概率Pc执行交叉算子；

子步骤3.3：以概率Pm执行变异算子；

子步骤3.4：验证新产生的个体是否满足模型的约束条件；如果满足，则执行下一步；否则重新执行步骤3.2和3.3，直至满足约束条件；

子步骤3.5：对每个个体的基因编码进行解码得到(I,K,Bx,By)，然后重复步骤2，计算新一代种群中每个个体的适应度；

步骤4：当终止条件达到时，输出适应度最小的个体，并对其进行解码得到最优的配置参数。

4.根据权利要求1所述的基于遗传算法的CT图像重建反投影加速方法，其特征在于：其在CUDA编程框架下，基于CUDA的反投影并行计算分为三个优化部分：

优化1：在第一次反投影核函数运算之前，首先完成I个投影的数据上传，然后在核函数运算的同时上传下次需要处理的投影数据；

优化2：在接下来的反投影中，GPU在对上次上传的投影数据进行反投影的同时，上传下次需要处理的投影数据；

优化3：对于最后一次反投影，将重建图像分成均等的两个部分，采用两个核函数分别完成对应的反投影，当第一个核函数完成反投影之后即开始下载数据到内存，同时第二个核函数计算另一半重建图像，待反投影结束后下载数据到内存。