并行各向异性扩散算法与实时医学图像增强技术

秦 安，孟晓林，冯前进，陈武凡

（南方医科大学 广东省医学图像处理重点实验室， 广东 广州510515）

并行各向异性扩散算法与实时医学图像增强技术

秦 安，孟晓林，冯前进，陈武凡

（南方医科大学 广东省医学图像处理重点实验室， 广东 广州510515）

各向异性扩散算法在去除图像噪声同时能保持重要的边缘和局部细节，因此在医学图像增强中得到广泛应用。但医学图像数据分辨率和灰度级都很高，求解各向异性扩散的偏微分方程时运算量很大，在许多实时系统中应用该类算法存在显著迟滞。本文采用CUDA架构，利用像素的独立性和偏微分方程求解的并发性，实现了并行各向异性扩散算法，在保证效果的同时显著降低了处理时间，满足了实时性的需求。

各向异性扩散方程；图像增强；CUDA；医学图像

Abstract:Anisotropic diffusion can remove noises while preserving image boundaries, and has been widely used in medical image enhancement.But typical medical images have high resolution and grayscale, and the computation of the partial differential equation is huge and the time cost has blocked its applications in real-time medical imaging system. This paper implements anisotropic diffusion image enhancement under CUDA architecture. By taking use of the independence of image pixel and concurrence of partial differential equation, parallel computing in GPU can greatly reduce the computing time while keeping the good performance of image smoothing.

Key words:parallel anisotropic diffusion equation; image enhancement; CUDA; medical image

随着成像设备在空间和时间分辨率上突飞猛进的发展，高维海量的医学图像数据自动增强和分析技术成为研究热点。图像增强是医学图像分析和重建的重要基础,它是后续处理（如分割、重建）的关键预处理。图像的平滑和细节保持是一对矛盾关系，常用的低通滤波在降低噪声的同时，产生图像边界的模糊，而人对图像的高频成分是很敏感的。基于偏微分方程的图像平滑在这个领域获得了广泛的重视，因为相对于传统的图像平滑方法，偏微分方程图像处理方法在平滑的同时，能很好地保留图像原有的边缘特征，有效地解决了图像平滑和细节保持的矛盾。在基于偏微分方程的图像处理方法中，各向异性扩散模型是其中重要的分支。由于各向异性扩散模型所具有的优良特性，已广泛应用于图像处理的许多领域，包括图像增强、图像分割中[1,2]。

医学图像尤其是三维影像数据量大，灰度级比自然图像丰富，在求解各向异性扩散的偏微分方程时运算量很大，不能满足实时系统的要求。本文采用P-M各向异性扩散方程，利用离散像素的相对独立性和求解偏微分方程的并发性，结合最近几年发展起来的通用编程的图形处理单元（GPGPU）模型CUDA技术，采用并行算法来提高图像数据处理的速度。该方法在保证滤除噪声的同时，保持了图像原有的边缘特性，并通过并行计算极大地缩短图像处理时间,以达到实时系统的需求[3]。

1 各向异性扩散算法

Perona和Malik[4]提出了基于偏微分方法的各向异性扩散方程(P-M方程)，即在非线性尺度空间，根据图像不同方向上的梯度确定扩散系数，其滤除噪声和保留细节的性能有了较大地提升，其方程如下:

其中，k是扩散系数中的一个参数，由它控制其扩散强度[6,8]。

P-M方程的离散化形式如下:

其中q为离散化图像中像素的位置， uq为图像中像素q处的值。t为迭代的次数，λ为决定扩散速率的常数，Nq为像素q的邻域，通常，二维图像取上、下、左、右四邻域，要取得更好的效果，可以取周围的八邻域[5]。

由于各向异性扩散采用了图像在不同方向上的梯度的单调递减函数作为扩散系数，在同质区域内部，它的扩散系数较大，可以有效地滤除同质区域内的噪声，而在边缘区域，灰度值变化剧烈，梯度较大，扩散系数就较小，能够保留图像的边缘信息[9]。

2 GPU和CUDA简介

通用编程的图形处理单元（GPGPU）将通用计算问题重构为图形渲染的方法来解决，对非图形应用的有效性支持依然有限，限制了解决方案在GPU上的可移植性。CUDA(Compute Unified Device Architecture)是NVIDIA公司为通用并行计算开发的在GPU上计算的统一设备构架，这个架构可以将GPU视为一个并行数据计算的设备，对所进行的计算进行分配和管理。它是一个完整的GPGPU解决方案，而提供了最大的编程灵活性，允许程序员更好地利用平行的权力来进行GPU的通用计算。

图 1 CUDA编程模型

从硬件层来看，GPU是组多处理器，每个多处理器使用单一指令多数据架构（SIMD）。从编程角度看，CUDA架构如图1所示，GPU作为主CPU的一个协处理器，可以被视为执行海量并行线程的计算设备。主机CPU进行必要的预处理，之后将要进行高度并行处理的数据拷贝到GPU内存，对其进行必要的线程分配，然后调用设备函数对数据进行并行处理。

3 并行各向异性平滑的实现

P-M方程的具体并行实现步骤如下：

(1)将图像数据读入主机内存,在设备内存分配空间并调用函数 cudaMemcpy()，将主机内存中的数据拷入设备内存分配的空间中。cudaMemcpy控制主机和GPU之间进行数据交换。

(2)程序进入设备函数Anisotropic_ kernel<<< Dg, Db >>>进行各个像素点的并行计算，分配的线程总数为图像所包含的像素的数目，每个线程按公式(3)处理一个像素数据。每个线程之间进行并行处理，直到完成指定的迭代次数。如图1所示，GPU/Device中Grid里每一个Block里的每一个Thread处理一个像素点，分配的Thread的总数等于图像像素的数目，每一个Thread并行地对对应的每一个像素进行迭代计算，直到完成指定的迭代次数[8]。

(3)把结果通过函数cudaMemcpy()将设备内存上的结果数据拷回主机内存中进行显示或写入磁盘文件中。

4 实验结果及分析

实验设置如下：硬件：CPU为Pentium 4（3.00GHz）；内存为1.0GB；GPU为NVIDIA GeForce 9800 GTX（显存512MB）。参数：噪声均值为0，方差为0.01的高斯噪声；迭代步长λ都取为0.08，除MRI图像的扩散参数k2取200外，其它图像都取100，迭代次数均取25。峰值信噪比的计算公式为：

表 1 实验图像数据列表

表 2 迭代过程中的时间及信噪比对比表

表 3 不同大小三维医学影像运行时间及加速比

图 2 不同图像实验效果图

图 3 DR图像增强效果

图 4 不同图像尺寸运行时间及信噪比曲线示意图

图 5 三维医学图像并行各向异性滤波结果体绘制

表2显示了不同的迭代次数的运行时间及信噪比变化。可以看出并行方法的处理时间相对串行大大缩减，且信噪比随着迭代次数增大。

图2和图3显示了本文算法增强效果。图3是高分辨率DR图像增强效果图，对局部细节进行了放大。可以看出，本文算法具有良好的去噪效果。

图4是不同图像大小下（4倍递增），串行和并行的运行时间及信噪比曲线图，该图显示，并行算法时间在图像数据增加是处理的时间基本保持不变，信噪比也是保持稳定，这也充分证明了CUDA在高分辨率数据并行处理中的快速优势。

图5和表3显示了在三维医学影像上采用本文算法进行图像增强的加速和体绘制效果比较。从表2看出在海量三维影像增强处理中，并行算法的仍然接近实时（小于1s）；而且从体绘制效果看，增强后三维影像的细节和边界更为清晰。

5 结论

由于医学图像数据量很大，灰度级丰富，求解偏微分方程时差分运算量大，各向异性扩散增强算法很难满足医学实时系统实时交互的要求。本文采用了CUDA架构，采用并行计算的方式对算法进行了加速。实验证明，基于CUDA加速的各项异性扩散算法在保证图像增强效果的同时，显著减少了运算时间，达到实时系统的要求。作者下一步将结合体绘制技术，对三维体数据滤波算法进行加速研究。

[1] 杨新.图像偏微分方程的原理与应用[M].上海:上海交通大学出版社,2003:68-73.

[2] Ceccarelli M,De Simone V,Murli A.Well-posed anisotropic diffusion for image denoising[J].IEEE Proceedings VISP,2002, 149(4):244-252.

[3] 程群,田有先.基于各向异性扩散方程的并行图像去噪研究[J].计算机工程与应用,2008,22(44):186-188.

[4] Perona P,Malik J. Scale space and edge detection using anisotropicdiffusion[J].IEEE Transon PAMI,1990,12(7):629-639.

[5] 蒋先刚.基于各向异性扩散的图像平滑及在三维重构预处理中的应用[J].计算机应用,2007,27(1):249-251.

[6] Yongjian Yu, Acton S.T. Speckle reducing anisotropic diffusion[J]. IEEE Transactions on Image Processing, 2002,11(11):1260-1270.

[7] NVIDIA CUDA Compute Unified Device Architecture Programming Guide[Z/OL].(2008-07-06)[2009-07-01].http:// cuda.csdn.net/upload/bczn2.0.doc.

[8] Catte F,lions P L,Morel J.Image selective smoothing and edge detection by nonlinear diffusion[J]. SIAM Journal on Numerical Analysis,1992,29(1):182-193.

[9] S.Tabik,E.M.Garzon,I.Garcla.Implementation of Anisotropic Nonlinear Diffusion for Filtering 3D Images in Structural Biology on SMP Clusters[J].Parallel Computing,2006,33(11):727-734.

Parallel Anisotropic Diffusion Algorithm and Real-time Medical Image Enhancement Technology

QIN An,MENG Xiao-lin, FENG Qian-jin, CHEN Wu-fan
(Key Labs of Medical Image Processing in Guangdong, Southern Medical University, Guangzhou Guangdong 510515,China)

TP391

A

10.3969/j.issn.1674-1633.2010.03.008

1674-1633(2010)03-0029-03

2009-07-02

国家自然科学基金重点项目(30730036)；广东产学研项目（cgzhzd0714）。

本文作者：秦安，博士，主要从事医学图像分析与可视化研究。

冯前进，博士，副教授。

作者邮箱：qa.free@gmail.com