基于CUDA的快速三维医学图像分割

《基于CUDA的快速三维医学图像分割》由会员分享，可在线阅读，更多相关《基于CUDA的快速三维医学图像分割（6页珍藏版）》请在装配图网上搜索。

1、基于 CUDA 的快速三维医学图像分割孟晓林， 秦 安， 徐 建， 陈武凡， 冯前进 （南方医科大学 生物医学工程学院， 广东 广州 510515）摘要：目的： 三维分割是医学图像分析和可视化中的重要组成部分，也是医学图像分割中的一个难点。 水平集方法在三维医学图像分割中有很广阔的应用前景，但是该算法的计算量大，不能达到实时处理的要求。 针对这个问题，提出了一 种基于 CUDA 的并行加速方法。 方法： 采用 NVIDIA 公司的 GPGPU 模型 CUDA，利用图像像素的独立性和偏微分方 程求解的并发性，提高 C-V 水平集算法的分割速度。 给出了并行计算的流程图，并对 C-V 水平集算法在

2、 CUDA 上的实 现进行了详细介绍。 结果： 实现了 C-V 水平集并行加速算法，该方法在保证分割效果的前提下，具有更快的分割速度。 结论： 所提出的方法是切实可行的，实现了快速的三维医学图像分割。关键词： 水平集； 医学图像分割； CUDA； 并行图像处理中图分类号： TP391文献标识码： A文章编号： 1005-202X（2010）02-1716-05Fast 3D Medical Image Segmentation Based on CUDAMENG Xiao-lin, QIN An , XU Jian, CHEN Wu-fan, FENG Qian-jin(The Departm

3、ent of Biomedical Engineering, Southern Medical University,Guangzhou Guangdong 510515, China)Abstract:Objective: 3D segmentation is an important part of medical image analysis and visualization. It also continues to belarge challenge in the medical image segmentation. While level sets have demonstra

4、ted a great potential for 3D medical image segmentation, these algorithms have a large computational burden thus are not suitable for real time processing requirement. To solve this problem, we propose a parallel accerelated method based on CUDA. Methods: We implement C-V level set algorithm in the

5、CUDA environment which is the NVIDIAs GPGPU model.The segmentation speed can greatly improved by using independence of image pixel and concurrence of partial differential equation .The paper shows the flow chart of the parallel computing and gives the detailed introduction of the C-V level set algor

6、ithm which is implemented in the CUDA environment. Results: Realizing the C-V level set parallel accerelated algorithm. This method has faster segmentation speed while preserving the qualitative results. Conclusions: This method is viable and makes the fast 3D medical image segmentation come true.Ke

7、y words: level set; medical image segmentation; CUDA; parallel image processing医学图像三维分割也一直是图像处理领域的一个经典难题。 由于问题的重要性和困难性，从 20 世纪 70 年代起，图像分割问题就吸引了很多研究人员为之付 出了很大的努力。 在三维分割中，几何活动轮廓模型 与水平集相结合的曲线演化方法是目前广为关注的 一种图像分割方法。 该方法利用轮廓曲线的几何特 性，建立轮廓曲线运动的能量函数，通过最小化这个 能量函数， 使轮廓曲线逐渐逼近图像中目标边界，并 利用水平集函数将轮廓曲线运动方程转化成求解数 值偏

8、微分方程的问题。 在这类方法中，Chan 和 Vese 基于 Munford-Shah 分割模型提出的 C-V 模型是研究前 言图像分割是三维图像数据处理的最重要和最关键的过程之一，是体数据分析和可视化的基础，而且收稿日期：2009-11-05基金项目：广东省产学研重点项目（No. cgzhzd0714）作 者 简 介 ： 孟 晓 林 , 男, 硕 士 , 主 要 研 究 方 向 ： 医 学 图 像 处 理 。 E -mail:holymxl。通讯作者：冯前进(1974-)，男，博士，副教授，主要研究方向 ： 医学图像的热点。 与传统的基于参数形变模型和几何活动轮廓模型不同，C-V 水平集模型

9、提取物体边界时不依赖图 像的梯度， 因此适于用梯度无意义或边界模糊的图 像。水平集方法虽然在三维分割中呈现出广阔的应 用前景，但 level set 算法本身计算量太大这个缺点限 制了它的应用。 医学图像数据量很大，灰度级更丰富， 在求解 C-V 模型的偏微分方程时，运算量很大，不能 满足实时处理的要求。 针对这个问题，本文结合最近 几年发展起来的通用编程的图形处理单元（GPGPU1） 模型 CUDA 技术， 实现了 C-V 水平集算法的并行加 速。 CUDA 对海量数据的并行处理带来的加速是惊人 的， 已有许多文献对该平台与 CPU 的速度进行了详F（准，c1 ，c2 ）= v 乙 准准 准

10、dxdydz2乙乙H 准准 准dxdydzI - c1+12准准准 准准dxdydzI - c21-H+2（2）由欧拉-拉格朗日方法求解上式，得到如下式（3）（6）的偏微分方程式6：乙I（x，y，z）H 准准 准dxdydzc （准） =（3）1乙 H 准准 准dxdydz乙I（x，y）（1- H 准准 准）dxdydz细的比较，该方法一方面保证了 C-V 模型本身的分c （准） =（4）2乙割效果；另一方面，通过并行计算提高图像数据处理的速度，缩短图像处理的时间，以达到实时处理的需 要2,3。1 Level set 方法和 C-V 模型三维分割时，level set 方法的基本思想是将随时

11、间演化的闭合曲面 C（p，t）隐含地表达为同样随时间 演化的四维连续函数 准（x，y，z，t）的一个具有相同函 数值的同值曲面，通常是准（x，y，z，t） = 0，称为零水 平集，而 准（x，y，z，t）称为水平集函数。 初始化时，level set 函数定义为 准（x，y，z，t） = d。 其中 d 是图像中任 意一点（x，y，z）到初始化水平集的最短距离 ， 函数的 符号取决于该点在曲面的内部还是外部，一般定义在曲面的内部点的距离为负值， 外部为正值，t 表示时 间。 Level set 图像分割方法中演化曲面的停止是依靠 图像的梯度信息，理论上当演化曲面到达物体的边界 时，梯度变化大，

12、曲面停止演化；然而当图像边界模糊 时，因为梯度变化不明显，演化曲面会在边缘处继续 运动，出现越界现象4,5。C-V 模型是基于水平集思想和 Munford-Shah 模 型，不利用梯度信息，而是通过最小化能量函数式（1） 来演化曲面。 当图像中只有目标和背景两类分片光滑 区域时，设图像 I（x，y，z）的全部图像域 被闭合边界 C 分割为目标 1（C 的内部）和背景 2（C 的外部）两个同质区域， 两个区域的平均灰度值分别为 c1 和 c2，建立如式（1）的能量函数：（1- H 准准 准）dxdydz= 准准 准准vdiv 准准-准I-c2 准准（5）准坠准22准I-c1 准+12坠tarct

13、g 准准准准准1+准12准H 准准 准=准准2准准准（6）准准 准准 准= 1 准准2 + 准2准准准准2 统一计算设备框架 CUDACUDA （compute unified device architecture） 是 NVIDIA 公司为通用并行计算开发的在 GPU 上计算 的统一设备构架， 这个架构可以将 GPU 视为一个并 行数据计算的设备， 对所进行的计算进行分配和管 理。 它是一个完整的 GPGPU 解决方案，提供了硬件 的直接访问接口，而不必像传统方式一样必须依赖图形 API 接口来实现 GPU 的访问。在架构上采用了一种全新的计算体系结构来使用 GPU 提 供 的 硬 件 资

14、源，从而给大规模的数据计算应用提供了一种比 CPU更加强大的计算能力。 CUDA API 是一种 C 语言的扩 展，并且提供一般 DRAM 内存寻址方式，从而给程序 员提供了最大的编程灵活性。在 CUDA 框架下，一个支持 CUDA 的 GPU 作为CPU 的协处理器， 适用于可以分解为单任务多数据（SIMD）并行模式的算法，他有自己的静态存储器、片 上共享存储器以及由一系列流处理单元 （streaming processor）整合而成的流式多处理器。 就编程模式而 言，GPU 执行的所有任务都有用户为 GPU 分配的一乙2I - c1F（C，c1 ，c2 ）= 1dxdydz1乙2I - c

15、2+2dxdydz +vC（1）系列线程模块（Block）中的线程（Threads）完成。各个2其中，1，20 分别为各项加权系数。 式中第一项和第二项为保真项，当闭合轮廓曲线 C 位于目标的边界 时，此两项和为最小值。 第三项为正则约束项，是闭合轮廓曲线 C 的长度约束，使轮廓曲线平滑。 将上式用 水平集的方法表示为式（2）：线程以阵列的形式组成了线程模块，而所有的线程模块又以阵列的形式组成一个最终并行计算的 网 格（Grid）。 一旦通过一系列 CUDA API 函数指定网格和 线程模块的维数大小，CUDA API 中的内嵌变量就会 据此来为这张网格中的每一个线程指定一个索引号，准用户通过

16、索引号来控制线程对存储于 GPU 上的数据读写或运算。 而同时，GPU 内部根据线程模块的分配 情况以及存储器的使用情况，使得尽可能多的流式多处理器同时工作， 顺序完成程序上分配好的并行线 程。在实际操作时，用户必须通过 CUDA API 将需要 并行处理的数据从 CPU 载入 GPU 的静态存储器中， 并进行线程分配，编写 GPU 的核函数（kernel）， 最后型 也 是 dim3， 指 定 每 个 Block 的 维 度 和 大 小 ，Db.x*Db.y*Db.z 等于各块的线程数量。 dim3 类型是一种 整形向量类型，基于用于指定维度的 uint37。 如图 1 所示，GPU/Dev

17、ice 中 Grid 里每一个 Block 里的每一 个 Thread 处理一个像素点， 分配的 Thread 的总数等于图像体素的数目，每一个 Thread 并行地对相应的每一个体素处理。 Kernel 2 用于式（6）的并行计算，计算结果存在 gh 和 gs 中，通过 cudaMemcpy（）把 gh 中 的数据拷入到 CPU 中用于对式（3）、（4）的串行处理。 Kernel 3 用于对式（5）进行并行处理。完成一定的迭代次数后把存在 g f 中的结果数据拷回 CPU 中。 其零水平集就是图像中待分割目标 的边界。并行计算流程如图 2。4 实验结果本文实验采用了 NVIDIA 公司为通用

18、并行计算 开发的在 GPU 上计算的统一设备构架 CUDA 来对 图像数据进行并行处理。 所用的硬件配置如下：CPU 为 Pentium （R）4 3.00 GHz； 内存为 1.0 GB；GPU 为 NVIDIA GeForce 9800 GTX（显存 512 MB）。 实验数 据为剑桥大学手术规划实验室网站中提供的脑瘤数 据，大小为 256256124，灰度级为 16 位8。 参数选择 如下：时间步长 t = 0.000001，空间离散步长 h = 1，再将处理后的结果放回 CPU 中。如图 1 所示。3 CUDA 上的实现实现的具体步骤如下：将图像数据读入 CPU 中。CUDA 的编程模

19、型在 GPU 上分配空间 gI，g f，gh，gs，空间的大小和图像的大小相同，gI 为存放图像的空间，g f 为存放水平集 准 的空间，gh 为存放 H（准）的空间，gs 为存放 （准）的空间，并调用函数 cudaMemcpy() 将主机内存 中 的 数 据 拷 入 gI 中 。 如 图 1 所 示 ， 通 过 函 数 cud- aMemcpy() 将 CPU 中的数据拷入到 GPU/Device 中。 cudaMemcpy 可以在 CPU 和 GPU 之间进行数据的拷 贝， 也可以在主机与主机之间，GPU 和 GPU 之间进行，通过设置该函数的最后一个参数决定进行何种操作。定义三个 Ker

20、nel 函数。 Kernel 1 用于求解符号 距 离 函 数 ， 即 t = 0 时 刻 的 水 平 集 函 数 。 闭 合 曲 面 C （p，t）初始化为一个球，球的半径根据不同的应用确 定。 kernel1 并行计算图像中各个体素点到球面的距离，分配的线程总数为图像所包含的体素的数目，一个线程处理一个体素数据。 每个线程之 间进行并行的处理。 Dg 的类型是 dim3，指定 Grid 的维度和大小，Dg.x*Dg.y 等于启动的块数量。Db 的类图 1 CUDA 编程模型Fig.1 CUDA programming model图 2 并行计算流程图Fig.2 The flow chart

21、 of parallel computing图 3 实验效果图Fig.3 Experimental results表 1 计算时间对比表Tab.1 Comparison of the computing time图像数据 CPU 计算时间（s）GPU 计算时间(s)GPU 相对 CPU 的时间加速比256256124523171= 1，2 =1， =1，v =0.000012552。 图 3 为迭代 12 次的分割结果，初始轮廓为半径为 6 的一个球。 表 1 为 相同的条件下 CPU 和 GPU 上的计算时间对比。从图 3 可以看出， 本文的基于 C-V 水平集的三 维医学图像分割方法在脑瘤

22、分割中能取得较好的分割效果。 从表 1 可知，本文在 CUDA 上实现的加速算 法的计算时间仅为 CPU 上串行处理时间的 1/17。5 结束语基于 C-V 水平集的分割方法在三维分割中有很 广阔的应用前景，但由于其计算量大，很难达到实时 处理的要求。 本文采用 CUDA 模型，采用并行计算的方式对算法进行了加速，实验证明，基于 CUDA 加速的 C-V 分割模型能很有效地应用于医学图像的三维 分割中，并且在保证分割效果的同时，极大的降低了图像处理的时间，满足了医学图像实时处理的要求。参考文献：1 Schenke S, Wuensche B, Denzler J. GPU-based volu

23、me segmentationJ. Proc of IVCNZ,2005.2 Joshua E.C, Aaron E.L, Ross T.W. An interactive, GPU -based level set segmentation tool for 3D medical images J. Medical Image Anal - ysis, 2004,8:217-231.3 Aaron E.L, Joshua E.C, Ross T.W. Interactive ，GPU -based level setsfor 3D segmentation J. Medical Image

24、Computing and Computer-Assisted Intervention, 2003:564-572.4 Oliver Klar. Interactive GPU -based segmentation of large medical volume data with level sets M. Masters Thesis,VRVIS, 2006.5 A.E.Lefohn, R.T.Whitake. GPU -based, three dimensional level set solver with curvature flow R. Technical report,

25、school of Computing,University of Utah, 2002.6 王大凯主 编. 图像处理的偏微分方程方法 M. 北 京 ： 科 学 出 版 社 ，2008：103-108.7 NVIDIA .CUDA Compute Unified Device Architecture ProgrammingGuide M. 8 http:/www.spl.harvard.edu/pages/Software.（上接第 1707 页）中华放射肿瘤学杂志, 2004,13(2):113-116.7 Jang SY, Liu HH, Mohan R. Underestimation

26、 of low-dose radiation in treatment planning of intensity -modulated radiotherapy. Int J Radiat Oncol Biol Phys, 2008,71(5):1537-46.8 Nioutsikou E, Bedford JL, Christian JA, Brada M, Webb S. Segmenta -tion of IMRT plans for radical lung radiotherapy delivery with the step-and-shoot technique. Med Ph

27、ys, 2004,31(4):892-901.9 傅卫华,胡逸民,戴建荣. 侧向电子失衡对肺部肿瘤放射治疗计划设计的影响J. 中国生物医学工程学报, 2004，23(1):49-54.10 张富利,郑明民,高军茂,等. 6MV 与 15MV X 线在肺癌调强放疗中的剂量学比较J. 2008,25(5):801-803.11 Pirzkall A, Carol M, Lohr F, et al. Comparison of intensity modulated radiotherapy with conventional conformal radiotherapy for complex -

28、 shaped tumors. Int J Radiat Oncol Biol Phys, 2000,48(5):1371-1380.12 Galvin JM, Ezzell G, Eisbrauch A, et al. Implementing IMRT in clin -ical practice: A joint document of the American Society for Thera- peutic radiology and oncology and the American Association of Physicists in Medicine. Int J Rad

29、iat Oncol Biol Phys, 2004,58 (5):1616-1634.（上接第 1711 页）技术的剂量学比较J 中国医学物理学杂志， 2008，25(5)：781-7846 郭小毛, 梅欣, 朱国培, 等 食管癌常规放疗与三维适形放疗比较及4 种不同设野技术的剂量学研究 J 中国癌症杂志 ， 2005，15(5)：462-4657 武新虎，李丹明 食管癌术后三种放疗方案的对比研究J 医学研 究生学报， 2007，20(8)8 Cominos M，Mosleh Shirazi MA，Tait D， et al Quantification and re- duction of

30、cardiac dose in radical radiotherapy for oesophageal can- cer J Br J Radiol，2005，78 (936)：1069-10749 Fu WH，Wang LH，Zhou ZM，et al Comparison of conformal and in-tensity - modulated techniques for simultaneous integrated boost ra -diotherapy of esophageal carcinomal J World J Gastroenterol，2004，10(8)：

31、1098-110210 蒋治勤 食管癌三维适形与束流调强放射治疗的研究 D 复旦 大学，2006.11 Graham MV，Purdy JA，Emami B，et alClinical dose -volume his- togram analysis for pneumonitis after 3D treatment for non - small2cell lung cancer (NSCLC) J Int J RadiatOncol B iolPhys，1999，45(2)：323-329file:/D|/我的资料/Desktop/新建文本文档.txtAppliance Error (configuration_error)Your request could not be processed because of a configuration error: Could not connect to LDAP server.For assistance, contact your network support team.file:/D|/我的资料/Desktop/新建文本文档.txt2012-07-12 20:42:52