非线性图像扩散LB模型的CUDA算法设计与实现

周 明， 严壮志， 黄 彬

上海大学通信与信息工程学院，上海200072

非线性图像扩散的格子波尔兹曼(Lattice Boltzmann,LB)方法是近年发展起来的一种结合建模与数值计算的图像处理技术，并成功应用于图像去噪[1-6]、分割[7]以及修复[8]等领域.作为一种数值解法，LB方法物理意义清晰，边界条件处理简单，程序易于实现，且具有并行化和稳定性等特点，显然优于传统数值解法.尽管目前LB图像处理方法较其他非线性偏微分方程的数值实现方法计算效率有所提高，但现有的处理速度远远不能满足复杂图像(例如医学超声图像)处理的实时性要求[5,7-8].因此，发挥LB模型的边界条件处理简单、程序易于实现和并行化的特点，进一步提高计算效率是一个重要的研究课题[9].目前，由于多核处理器和GPU可以低成本使用，研究非线性扩散LB模型的并行算法及其在GPU上的实现具有重要的应用价值.

LB模型在GPU上的并行计算一般可分为通用GPU(general purpose GPU,GPGPU)和计算统一设备架构(compute unif ied device architecture,CUDA)两种，且大多用于计算流体力学领域。前者包括混合流模拟[10]和三维实时流体模拟[11]等，后者包括圆形绕流大涡模拟[12]和二维、三维方腔流模拟[9,13-14]等.这些应用都获得了较高的GPU-CPU加速比和每秒百万格子更新(million Lattice update per second,MLUPS).例如，文献[12]模拟圆形绕流大涡LB模型在nVIDIA GTX275上的计算，相对于CPU可得到18倍的加速比；文献[13-14]则报道了在nVIDIA GTX280上模拟二维和三维方腔流的LB模型计算，相对于CPU的加速比可达50倍以上，每秒百万格子更新可达300次以上.在图像处理领域，仅有少量文献报道了在GPGPU上实现LB模型的并行处理技术[6]，而通过CUDA实现图像处理LB方法的研究则尚无文献报道.

非线性图像扩散LB模型具有清晰而明确的物理解释，它以图像特征和图像先验知识作为约束条件，能兼顾自下而上和自上而下的分析方法，故在图像处理与分析领域具有很好的应用前景.而作为一种并行算法设计架构，CUDA具有易实现、可及性好和GPU-CPU加速比高等特点.为此，本文将聚焦于非线性图像扩散LB模型，研究其CUDA并行算法的设计与实现.

1 非线性图像扩散LB模型

文献[1]提出了基于LB模型的非线性扩散方法并成功地应用于图像去噪，文献[2-5]进一步改进了非线性图像扩散LB模型.本文针对文献[5]提出的非线性图像扩散LB模型，研究了CUDA算法.该模型的优点在于将图像边缘特征嵌入松弛因子中，不但具有良好的边缘保持能力，而且实现了大步长迭代，提高了计算效率.

假设图像由离散化的网格组成，每个网格的值由粒子的分布函数Ii(i=0,1,2,···,b-1)和离散速度矢量ci(i=0,1,2,···,b-1)决定.根据ci方向数b的不同，一般把二维LB模型分为D2Q5和D2Q9两类.LB模型在t时刻位于r处，沿ci方向的演化过程由迁移过程和碰撞过程组成[9]：

迁移过程

碰撞过程

式中，Ii(r,t)为粒子分布函数，(r,t)为平衡态分布函数，为发生迁移后的粒子分布函数，Δx为空间步长，Δt为时间步长，ω为松弛因子.

综合以上两个过程，可得到LB演化方程

松弛因子ω为

式中，C为迭代步长，g(▽Gσ∗I)为在大步长迭代下的边缘截止函数

式中，Gσ为方差为σ的高斯核；K为平滑阈值，用以控制扩散的范围.

在D2Q9模型中，平衡态分布函数为

松弛因子ω为

由式(3)～(8)可以看出，该LB模型的计算具有良好的局部性，非常适用于并行计算.以D2Q5为例，图像去噪的非线性扩散LB算法流程如下：

确定迭代次数N

While n＜N do

设置初始平衡态函数：

根据式(5)计算松弛因子

根据式(1)计算迁移过程

根据式(2)计算碰撞过程

更新节点分布函数

更新平衡态分布函数

n=n+1

End While

输出

2 CUDA编程模型及CUDA算法设计

2.1 CUDA编程模型

CUDA程序的执行如图1(a)所示，即让主机(CPU)中每个核函数(kernel)按照线程网格(grid)的概念在显卡硬件(GPU)上执行；每个线程网格可以包含多个线程块(block)；每个线程块又可以包含多个线程(thread)，线程为最小的执行单元.在硬件执行的时候，如图1(b)所示，让每个线程对应一个流处理器(SP)；每个线程块对应一个多流处理器(SM)；线程网格被以线程块为单位分配到SM上执行.

图1(c)为GPU的存储器资源.在Fermi架构下，芯片上有大量线程私有的32位寄存器和线程块私有的64 kB共享内存，它们具有低延迟性.其中同一个线程块内共享内存的数据是可共享的.芯片外有全局内存即显存，虽然数据存储空间大，但对其访问有400～600个时钟周期的延迟.在全局内存中有一块区域大小为768 kB的纹理内存，纹理内存具有缓存功能.在CUDA程序中，数据是由CPU传到GPU的全局内存上.全局内存的优化思路是合并对齐访问隐藏延迟.合并对齐访问条件为：①线程访问的数据长度为4、8或16字节；②被访问地址构成一片连续内存空间；③第N个线程访问第N个全局内存地址；④起始全局内存地址对齐到所存储数据长度的16倍处.

图1 CUDA软硬件架构及GPU存储器资源[15]Figure 1 CUDA architecture and GPU memory resources[15]

2.2 CUDA算法设计

LB去噪方法的CUDA算法流程如下：

①分别在CPU和GPU上分配内存，初始化CPU和GPU上的相关数据；

②如图2所示，设置线程网格和线程块的大小，并在GPU端进行LB去噪方法的并行计算；

③当计算满足迭代次数n时结束计算，将数据拷贝回CPU做后续处理，释放GPU和CPU的内存.

图2 假设图像大小为9×9，线程网格设置为(3×9)，线程块设置为(3×1)，线程数据所对应的内存存储Figure 2 Working process of global memory when image size is 9×9,grid is 3×9,and block is 3×1

由于LB模型的计算只涉及当前节点及周围节点的信息，而每个格点的计算又相对简单，因此CUDA程序性能的瓶颈主要是GPU的数据访问.由2.1节可知，数据均存储在全局内存中，是否满足合并对齐访问条件是影响程序性能的关键.碰撞过程仅涉及当前节点的计算，数据访问容易满足合并对齐访问条件；而迁移过程需要将计算后的数据进行移动，在需要向左或向右移动的方向，数据访问的起始地址会发生偏移(不满足合并对齐访问条件的第4条)，如图3所示.

图3 迁移过程的数据访问地址发生偏移，如线程0的数据迁移到线程1上Figure 3 Diagram of shift during data access in propagation steps:the data of thread 0 shifts to thread 1

2.2.1 CUDA算法

针对上述迁移过程数据访问的特点，本文设计了以下3种CUDA算法：

算法1 借助纹理内存实现迁移过程的数据访问.纹理内存具有高速缓存，支持二维寻址，一个数据的“上下左右”数据都能读入缓存，数据读取时不需要满足合并对齐规则.因此纹理内存十分适合用于图像处理和非对齐的数据访问[15].

算法2 借助共享内存实现迁移过程的数据访问.共享内存只有1个时钟周期延迟且同一个线程块内数据可共享，因此比较适合图像处理以提高效率.

算法3 直接使用全局内存实现迁移过程的数据访问.这种算法的优点是只需一个核函数即可完成，减少了对全局内存的访问次数，且全局内存资源丰富，可有效提高多流处理器的占用率.

此外，由于在传统算法2和3中，当采用一块数组内存时，容易发生数据覆盖.如图4所示，数据由0传到1，再由1传到2，使原来1中的数据被覆盖而产生错误.为此，在本文的算法2和3中，分别采用两个相同大小的数组来储存格点的分布函数，分别设为A和B.对于奇数迭代步，数据从A中读取，经LB去噪方法处理后存储到B；对于偶数迭代步，则反过来执行.

图4 迁移过程的数据覆盖Figure 4 Data sample covering problem in propagation steps

2.2.2 CUDA算法的实现

图5具体介绍算法1的实现.首先定义纹理变量，然后使用cudaBind Texture2D()将全局内存中的图像数据绑定到纹理对象，最后使用tex2D()函数访问该纹理对象，就可利用纹理内存的特性对图像数据进行操作.

图5借助纹理内存实现迁移过程Figure 5 Propagation steps using texture memory

图6 具体介绍算法2的实现，本文选取了粒子向右方向的迁移.首先将最右端的粒子分布值移到最左端，其次将其余的粒子依次向右移，这样在迁移过程就满足合并对齐访问条件.经过上述处理，虽然满足了合并对齐访问，但线程块两端的数据没有正确的迁移.因此，需要另一个核函数来处理线程块两端的数据，将每个线程块最左端的数据移到右边线程块的最左端，从而正确地完成了粒子的向右迁移.同样，算法2也可以通过粒子向左方向的迁移来实现.

图6 借助共享内存实现粒子向右方向的迁移Figure 6 Propagation steps from left to right using shared memory

算法3的实现按照算法流程即可.数据传输到GPU后，根据非线性扩散LB算法流程实施：即先计算松弛因子ω，接着根据式(1)完成迁移过程以及式(2)完成碰撞过程，再更新Ii(r,t)和，满足迭代次数后停止计算.

3 实验及讨论

3.1 实验平台

为了验证CUDA算法在图像处理中的可行性，本文搭建一个实验平台.平台包括Xeon W3550的CPU和8 GB内存.为了分析不同GPU的性能，本文分别使用Quadro4000和GT550M两款GPU在此平台上做实验.程序编译运行环境为CUDA4.0及VS2010.两款GPU均为Fermi架构.

在正式实验之前，本文首先做一个前置实验来选取合适的线程块尺寸(线程块将对GPU计算效率有影响，见图2). 在前置实验中选取D2Q9(一种常见的二维LB模型)为对象.其中针对算法1，为了满足blockDim.y≥2(线程块在纵坐标的尺寸，因为这个算法使用二维纹理)，选用(32×4)，(64×2)，(64×4)，(32×8)，(128×2)等5种线程块尺寸对512×512大小的图像做实验.针对算法2和算法3，一般将blockDim.y设置为1[9,14]，选取(32×1)，(64×1)，(128×1)，(256×1)等4种线程块尺寸对512×512大小的图像做实验.

图7给出了本文CUDA算法每次迭代时间与不同线程块尺寸的关系.从图7中可以看出，不同的线程块尺寸，程序性能有差异；算法1～3的最优线程块尺寸分别为(64×2)、(128×1)、(128×1).因此，这3个线程块将被运用到正式的实验中去.

图7 CUDA算法每次迭代时间与不同线程块尺寸的关系Figure 7 Relationships between caculation time of CUDA algorithms and different blocks

3.2 实验方法

为了验证基于LB模型的CUDA算法在图像处理中的有效性和效率，本文分别针对图像去噪质量、计算效率和真实医学图像处理效果进行3个实验.

3.2.1 评价图像去噪质量的实验方法

本实验的目的是定量评价分别利用CPU和GPU实现非线性图像扩散去噪的质量.实验采用Lena图像作为标准图像来测试.在该标准图像中加入均值为0和方差为0.01、0.03、0.05、0.07、0.09的高斯噪声.采用图像的峰值信噪比(PSNR)作为评价标准.PSNR定义为

式中，I(r)为未加入噪声的原始图像，ρ(r,N)为去噪后的图像，m和n分别为图像的行数和列数.实验中分别用CPU和GPU对上述加噪图像进行处理，以处理后图像的PSNR为指标对图像去噪质量进行评价，实验结果见图8.

3.2.2 评价计算效率的实验方法

本实验用来比较分析上述CUDA算法的性能及两款GPU实现的程序性能.实验采用256×256、512×512、1 024×1 024大小的图片做处理.实验中针对上述3种大小的图片，分别用CPU和两款GPU对它们进行处理，以每次迭代运行时间和MLUPS值为指标，对CUDA算法的计算效率进行评价，实验结果见表1.

3.2.3 真实图像去噪验证的实验方法

本实验进一步利用上海肿瘤医院超声科提供的肿瘤图像(256×256大小)，来验证基于LB模型的CUDA算法在真实医学图像中的处理效果.为便于比较，将CUDA算法与非线性扩散去噪差分方法中应用广泛的加性算子分裂算法(AOS)进行实验对比.实验中取步长C=2，迭代次数为8，边缘截止函数阈值K=4.实验以两种算法处理图像后是否存在伪纹理和对图像边缘的保持能力为评级指标，处理后图像的细节见图9.

3.3 实验结果及讨论

3.3.1 图像去噪质量的结果

图8为CUDA算法处理图像后获得的PSNR随噪声方差增加的变化曲线.从图8中可以看出，3种CUDA算法的PSNR随噪声方差增加有相同的下降趋势.由于CPU在硬件执行、浮点计算精度方面与GPU存在差异，CPU和GPU上实现的PSNR略有偏差，但总体保持一致，具有相同的去噪效果[15].

图8 CUDA算法处理图像后的PSNR与噪声方差的关系Figure 8 Relationships between PSNR of CUDA algorithms and noise variance

3.3.2 计算效率的结果

表1为D2Q9在CPU和两款GPU上实现的每次迭代运行时间和MLUPS值，每行为图像大小，每列为CPU和GPU上实现的每次迭代运行时间和MLUPS值，其中GPU分为Quadro4000和GT550M，两款GPU下又分为3种CUDA算法(算法1～3).本文中所提到的加速比为CPU和GPU上的运行时间比，如GPU-CPU加速比为CPU和GPU的运行时间比.

表1 D2Q9在CPU和两款GPU上实现的每次迭代运行时间和MLUPS值Table 1 Calculation time per iteration and MLUPS value of D2Q9 amongst two GPUs and CPU

由表1可以看出，在Quadro4000上实现，算法1的GPU-CPU加速比为70+，MLUPS值为300+；算法2的GPU-CPU加速比最低也有40+，MLUPS值最低也有200+；算法3的GPU-CPU加速比可达90+，MLUPS值可达400+.在GT550M上实现，算法3的GPU-CPU加速比为30+，MLUPS值为150+；算法1和算法2的MLUPS值也有110+.本文的MLUPS值低于流体LB方法[9,14]，原因是LB去噪方法在碰撞过程需进行梯度和滤波运算，数据访问相对复杂.即便如此，D2Q9在两款GPU上MLUPS值最低为110+，最高可达400+；而文献[6]在GPGPU上实现D2Q9的LB去噪方法，其MLUPS值仅为7.相比之下，本文CUDA上实现的计算效率得到了大幅度的提高.

从表1中还可以看出，3种算法的效率由高到低依次是算法3、算法1、算法2.原因如下：算法1和2虽然解决了全局内存的非合并对齐访问，但算法1中需要纹理绑定，且多个核函数增加对全局内存的访问次数，影响程序性能.算法2中需要使用if语句来判断线程块两端数据的迁移和边界处理，这样就使线程束内的线程走向不同的分支，所需要的时间将是不同分支之和，分支严重影响效率.随着计算规模的增大，线程块对共享存储器使用率增大使得加速比逐渐提高.算法3更高效，原因是使用全局内存可使本去噪方法在同一个核函数内执行，避免了多个核函数对全局内存数据的多次访问；同时线程块不受资源限制，可使SM中SP的占用率更高.

3.3.3 真实图像去噪验证结果

图9为CUDA算法与AOS处理超声图像后的细节对比，图9(a)和9(d)为原超声图像，图9(b)和9(e)为AOS处理后的图像，图9(c)和9(f)为CUDA算法(LB去噪方法)处理后的图像.由图9可以发现，AOS处理后的图像会产生伪纹理(图9(b)和9(e)箭头所示处)，而CUDA算法不会出现伪纹理(图9(c)和9(f)箭头所示处)，具有良好的边缘保持能力.

3.3.4 关于两款GPU性能对比的结果和讨论

从理论上来说，GPU计算性能主要由单精度计算能力和带宽决定[16-17].在不考虑带宽等因素的影响下，GPU计算性能为

式中，µ为换算系数，是一个常数；P为计算性能；f为峰值单精度浮点值.此外，当显存带宽提升一倍，约可提升30%的计算性.表2为实验中选取的两款GPU的主要技术指标.由表2可以得出Quadro4000和GT550M的峰值单精度浮点值和带宽之比分别为1.71和3，于是可推出Quadro4000计算性能约为GT550M的2.74倍，与两款GPU的流处理器数之比相当.

表2 Quadro4000和GT550M主要技术指标Table 2 Main specif ications of Quadro4000 and GT550M

图9 CUDA算法与AOS处理超声图像后的细节对比Figur e 9 Comparisons of CUDA algorithm with AOS in detail of ultrasound images

由表1可以看出，算法1和3的GPU-GPU加速比在2.7～3.1之间，和上述推论吻合；算法2偏低的原因是if语句使束内线程走向不同分支而串行执行，对程序的并行性产生了不确定性因素.因此，在保证程序并行性的前提下，综合考虑GPU的单精度浮点计算能力和带宽，不同GPU的计算性能和流处理器数量成正比.

4 结语

本文针对非线性图像扩散LB模型研究其并行实现的3种CUDA算法，并通过实验验证了三种算法在图像处理中的可行性.验证指标包括处理后图像的PSNR，加速比和MLUPS值.实验结果表明，在保证去噪质量的前提下，CUDA并行算法比已有的CPU算法计算效率大幅提高.其中直接使用全局内存实现迁移过程的CUDA算法(算法3)在效率上明显优于借助纹理内存(算法1)和共享内存(算法2)的CUDA算法.本文还通过医学超声图像去噪，验证CUDA算法(LB去噪方法)在大步长迭代下具有良好的保持边缘能力.此外，本文进一步使用GT550M和Quadro4000做实验对比，通过以单精度浮点峰值性能和带宽为指标，验证了在保证程序并行性的前提下，GPU的计算性能与其流处理器的数目成正比.

[1]JAw ERTHB,LINP,SINZINGERE.Lattice Boltzmann models for anisotropic diffusion of images[J].Journal of Mathematical Imaging and Vision,1999,11(3):231-237.

[2]陈玉，严壮志，钱跃竑.基于格子波尔兹曼模型的图像去噪[J].电子学报，2009,37(3):574-580.

CHENYu,YANZhuangzhi,QIANYuehong.The Lattice Boltzmann method based image denoising[J].Acta Electronica Sinica,2009,37(3):574-580.(in Chinese)

[3]CHANGQ S,YANG T.A Lattice Boltzmann method for image denoising[J].IEEE Transactions on Image Processing,2009,12(18):2797-2802.

[4]ZHANGW H,SHI B C.Application of Lattice Boltzmann method to image f iltering[J].Journal of Mathematical Imaging and Vision,2012,43:135-142.

[5]王志强，严壮志，钱跃竑.图像非线性扩散去噪的格子波尔兹曼方法[J].应用科学学报，2010,28(4):367-373.

WANG Zhiqiang,YAN Zhuangzhi,QIAN Yuehong.Nonlinear diffusion for image denoising using Lattice Boltzmann method[J].Journal of Applied Sciences,2010,28(4):367-373.(in Chinese)

[6]ZHAOY.Lattice Boltzmann based PDE solver on the GPU[J].The Visual Computer,2008,24:323-333.

[7]WANG Z Q,YAN Z Z,CHEN G.Lattice Boltzmann method of active contour for image segmentation[C]//Sixth International Conference in Image and Graphics(ICIG),2011:338-343.

[8]张蕊，严壮志，刘玮.图像修复的格子波尔兹曼方法[J].电子测量技术，2011,34(3):46-65.

ZHANGRui,YANZhuangzhi,LIUWei.Lattice Boltzmann method based image inpainting[J].Electronic Measurement Technology,2011,34(3):46-65.(in Chinese)

[9]KUZNIKF,OBRECHTC,RUSAOUENG.LBMbased flowsimulationusingGPUcomputingprocessor[J].Computers&MathematicswithApplications,2010,59(7):2380-2392.

[10]朱红斌，刘学慧，柳有权.基于LatticeBoltzmann模型的液混合流模拟[J].计算机学报，2006,29:2071-2079.

ZHUHongbin,LIUXuehui,LIUYouquan.Binary mixturessimulationbasedonLatticeBoltzmann method[J].ChineseJournalofComputers,2006,29:2071-2079.(inChinese)[

11]柳有权，刘学慧，吴恩华.基于GPU带有复杂边界的三维实时流体模拟[J].软件学报，2006,17:568-576.

LIUYouquan,LIUXuehui,WUEnhua.Real-time3D fluidsimulationonGPUwithcomplexobstacles[J].JournalofSoftware,2006,17:568-576.(inChinese)

[12]ZHOUH,MOGY,WUF.GPUimplementation ofLatticeBoltzmannmethodforflowswithcurved boundaries[J].ComputerMethodsinAppliedMechanicsandEngineering,2012,225-228:65-73.

[13]OBRECHTC,KUZNIKF,TOURANCHEAUB.Anew approachtotheLatticeBoltzmannmethodfor graphicsprocessingunits[J].Computers&MathematicswithApplications,2011,61:3628-363.

[14]TOLKEJ.ImplementationofaLatticeBoltzmann kernelusingthecomputeunifieddevicearchitecture developedbynVidia[J].ComputingandVisualizationinScience,2010,13:29-39.

[15]陈曙晖，熊淑华.大规模并行处理编程实战[M].北京：清华大学出版社，2010:106-116.

[16]http://www.realworldtech.com/amd-nvidia-gpu-per formance/.

[17]http://www.realworldtech.com/gpu-memory-bandw idth/.