基于GPU的IVOCT图像管腔分割算法

任江+贾小云

摘要：为了解决传统的CPU串行算法对血管内光学相干断层扫描（IVOCT）图像进行管腔分割时耗时较多的问题，提出一种基于图形处理器（GPU）的IVOCT图像管腔分割算法。首先，分析基于CPU管腔分割算法可知算法最耗时的两部分是坐标转换算法获取极坐标图像和动态规划算法获取管腔内壁轮廓。接着，分析坐标转换算法和动态规划算法的并行性，在统一设备架构（CUDA）下实现GPU加速优化。最后，在MATLAB中进行仿真实验，定量分析GPU加速后的坐标转换算法、动态规划算法以及最终管腔分割算法的优化性能。实验表明，基于GPU的管腔分割算法较CPU串行算法达到21倍加速比，能够对IVOCT图像序列进行快速处理，基本满足了冠状动脉疾病诊断和治疗的实际需求。

关键词：CUDA；GPU；管腔分割；IVOCT图像序列；冠状动脉疾病

中图分类号：TP391.54 文献标识码：A 文章编号：1009-3044（2017）28-0198-03

1 背景

心血管病是中国居民的首位死因，心血管病患病率及死亡率仍处于上升阶段[1]。血管内光学相干成像（Intravascular Optical Coherence Tomography， IVOCT）是目前国内外较新的冠状动脉内影像技术，分辨率约为10μm，是血管内超声分辨率的十倍[2]。血管管腔区域信息在冠状动脉疾病的诊断和治疗中具有很大的价值，如评估冠状动脉的狭窄程度，获取支架植入的最佳位置，评估支架的贴壁情况等[3]，因此由IVOCT图像管腔分割算法获取血管管腔区域是冠状动脉疾病的诊断和治疗中十分重要的基础性工作。

文献[4][5]提出基于CPU的IVOCT图像管腔分割算法具有较强的鲁棒性和准确性，但是时间复杂度较高，处理时间较长。在临床实践中，通常需要对数据量很大的IVOCT图像序列进行快速处理，CPU处理速度已无法满足实际需求。近年来，通用计算能力不断增强的图形处理器（Graphics Processing Unit，GPU）已经广泛地应用到医学图像处理领域。NVIDIA提出的统一设备架构（Compute Unified Device Architecture，CUDA）提高了GPU的可编程性，极大方便了基于GPU的通用计算程序开发。本文将分析文献[4][5]提出的基于CPU的IVOCT图像的管腔分割算法，在CUDA架构下对其中耗时较多的部分进行GPU加速，提出一种基于GPU的IVOCT图像的管腔分割算法，力求大幅缩短IVOCT图像管腔分割整体耗时。

2 基于CPU的IVOCT图像管腔分割算法

文献[4][5]提出的基于CPU的IVOCT图像管腔分割算法的主要流程：利用坐标转换算法，将直角坐标系下的IVOCT图像转换到极坐标系下；在极坐标下对图像进行预处理来移除导丝和导管；对预处理后的图像进行滤波，获取能量图；在能量图中利用动态规划（Dynamic Programming，DP）算法获取管腔轮廓坐标点，从而实现管腔分割。以单张704×704的IVOCT图像为处理，在MATLAB中对基于CPU的管腔分割分割算法进行仿真实验，统计算法耗时可知单张IVOCT图像的管腔分割耗时约为2.17秒，坐标转换算法在总耗时中的占比最大，占总耗时的77%，其次是DP算法，占总耗时的20%

2.1 坐标转换算法

设IVOCT图像在直角坐标下为imR，大小为M×N，极坐标系下为imP，大小为Mp×Np。由imR到imP的转换过程中，选取imR的图像中心（Om，On）作为极点，Om=（M+1）/2，On=（N+1）/2，imP的原点对应imR的中心点，imP的x方向对应角度，y方向对应半径。对imP中的像素点（r，c），1≤r≤Mp，1≤c≤Np，通过公式（1）映射到imR中的坐标点（m，n），如果m，n都为整数，则imP（r，c）的灰度值等于imR（m，n）；如果m，n不为整数，找到与它相邻的四个像素点，双线性插值获取对应的灰度值。遍历imP中的所有的像素点，进行相同的转换操作就可以得到极坐标下的IVOCT图像。

2.2 DP算法

获取管腔内壁轮廓可视为在能量图中搜索一条从第1行到最后一行累积能量最小的路径，由于轮廓的连续性，第1行到第i行的最小累积能量路径与第1行到i-1行的最小累积能量路径相关，想要得到第1行到最后一行的路径就必须依赖于第一行到中间行的最小累积能量路径，是一种典型动态规划问题，可以从局部最优的最小累积能量路径一步步递推到全局最优的路径。用公式（2）定义累积累积能量函数E：

E（i，j）表示第一行到i行j列像素點的累积能量，Min为最小值函数，j-m≤ j*≦j+m，m为相关系数，用来确定搜索邻域的宽度，本文中取m=2。e（i，j）为该像素点的能量，由滤波获取。通过上述累积能量函数获取第一行到每个像素点的累积能量后，从最后一行开始反向搜索使得累积能量全局最小的路径，最终获取管腔内壁轮廓的坐标点。

3 基于GPU的管腔分割算法

在实际应用中，需要对100-270帧的IVOCT图像序列进行快速处理，基于CPU的分割算法显然无法满足冠状动脉疾病诊断和治疗的实际需求。因此，我们提出一种GPU加速优化方案，利用GPU强大的并行计算能力优化坐标转换算法和DP算法，从而缩短IVOCT图像管腔分割的耗时。下面介绍针对坐标转换算法和DP算法的GPU优化方案。

3.1 坐标转换算法的GPU加速优化

通过2.1的介绍可知，IVOCT图像由直角坐标系转换到极坐标系的过程中，imP的每个像素都需要与imR进行坐标映射，然后在imR上双线性插值获取灰度值。每个像素执行的操作具有高度的并行性而且不存在数据依赖，适合在CUDA架构下进行并行计算。为imP中的每个像素都分配一个CUDA线程，将像素执行的坐标映射和双线性插值操作写在核函数函数中，每个线程并行执行Kernel函数完成计算就可以得到imP.。CUDA架构下实现的主要步骤：endprint

1） 初始化GPU设备，为imR，imP分配GPU存储空间，将存放CPU端内存中的imR数据拷贝到GPU中，dev_imR，dev_imP分别指向imR，imP显存空间的首地址。CPU端计算imR的中心点（Om，On），作为参数传入核函数。

2） 确定Block和Thread数目。对于单张IVOCT图像，选用二维线程网格和二维线程块，线程块中的线程数设置为16×16，线程网格中线程块个数为（ceil（Np/16））×（ceil（Mp/16）），Mp，Np分别是imP的高度和宽度，ceil函数是向上取整函数，确保有足够多的线程块来存放与imP像素个数相同的线程。对于IVOCT图像序列，选用三维线程网格和二维线程块，线程块的大小为16×16，线程网格中线程块的个数（ceil（Np/16））×（ceil（Mp/16））×F，F表示图像序列的帧数。

3） 启动核函数，GPU端的每个CUDA线程执行核函数完成计算。单张IVOCT图像执行核函数ImToPolar_2D，imP中每个像素点与CUDA线程一一对应，每个线程负责对应像素点的计算。ImToPolar_2D伪代码如下：

4） 将显存中的计算结果传回CPU端内存中，释放开辟的GPU显存空间

3.2 DP算法的GPU加速优化

由2.2可知，计算每个像素点的累积能量时，除第一行外，每行像素点的累积能量的计算都依赖前一行的计算结果，不同行中的像素点计算累积能量时存在数据依赖关系，但是对于在同一行不同列的像素，累积能量的计算不存在数据依赖关系，具有并行性，适合在CUDA架构下进行并行计算。对于M×N能量图e，基于GPU并行计算能量累积矩阵E时，将计算像素点累积能量的操作写入核函数中，执行M-1个核函数，依次计算第二行到第M行各行像素的累积能量，每个核函数开辟N个线程，这些线程并行计算对应行中各列像素点的累积能量。CUDA架构下实现的主要步骤：1）初始化GPU设备，在CPU端创建大小为M×N的累积能量矩阵E，将能量图e中第一行数据复制给E的第一行，接着将CPU端内存中存放的E和e复制到GPU端的显存中，dev_e，dev_E分别指向e和E显存空间的首地址。2）确定Block和Thread数目。选用一维线程网格和一维线程块，线程块中的线程数设置为16×16，每个线程网格中包含ceil（N/16）个线程块。3）依次启动M-1个内核函数DPTrack_GPU，每个CUDA线程并行执行核函数完成计算。DPTrack_GPU的伪代码如下：

4） 将显存中累积能量的计算结果传回CPU端的内存中，接着在CPU端对累积能量矩阵从最后一行开始反向搜索使得累积能量全局最小的路径，最终获取管腔内壁轮廓。最后释放开辟的GPU显存空间。

3 实验结果与分析

实验硬件平台为Intel Core i5-4460 CPU，内核主频3.20GHz，内存为8GB。GPU型号为NVIDIA Geforce GTX 770，1536个CUDA处理核心，显存位宽256位，显存带宽224.3GB/s，显存为2GB。软件平台为Microsoft Windows 7操作系统，Matlab2015b，CUDA Toolkit7.0，CUDA Driver 378.92。

实验采用的IVOCT图像来自FD-OCT（C7-XR）系统。为了进行多组数据的对比实验，对原始图像进行了预处理，通过裁剪获取6组图像数据，图像大小分别为128×128，256×256，512×512，1024×1024，2048×2048，4096×4096。表1为坐标转换算法和DP算法的CPU与GPU实现耗时对比。

由表1可知基于GPU实现的两种并行优化算法较CPU都有良好的加速效果。随着图像大小的增加，两种算法的加速比也随之增大。图2为两种算法的加速比增长对比图。

如图2所示，当图像大小从128×128增加到256×256时，坐标转换算法加速比曲线斜率相对较小，加速比增长趋势较为平缓。从256×256增加到1024×1024时加速比曲线斜率增大，加速比增长较快，加速效果十分明显。但从1024×1024增加到4096×4096时，加速比曲线仍然处于上升趋势，但是上升趋势明显变平缓。出现这种增长趋势是因为图像大小较小时，处理的数据量较小，CPU可以在相对较短的時间内完成计算。同时，GPU在初始化和主机端和设备端数据通信需要消耗一定的时间，这部分时间在GPU总运行时间中占比较大。然而随着图像大小的逐渐增大，CPU的计算能力趋于饱和，GPU额外开销时间在GPU总运行时间的占比越来越小，GPU相对于CPU的并行计算能力优势凸显，故该阶段加速比增长快。但是，随着图像大小的进一步增大，GPU的计算能力也逐渐趋于饱和状态，加速比增长放缓。DP算法的加速比增长趋势类似于直角坐标系到极坐标系转换算法，但是由于GPU的计算能力在图像大小为4096×4096时尚未趋于饱和，所以加速比上升趋势未出现放缓的迹象。对271帧707×704的IVOCT图像序列，实验测得CPU算法总耗时约为588秒，本文提出的GPU算法总耗时约为28秒，达到21倍的加速比。

3 结束语

本文通过对基于CPU的IVOCT图像管腔分割算法中耗时最多的两部分进行GPU加速优化，提出了一种基于GPU的IVOCT图像管腔分割算法。实验表明，随着图像的增大，坐标转换算法和 DP算法的加速比也随之增大。与CPU串行算法相比，基于GPU的IVOCT管腔分割算法能够对IVOCT图像序列进行快速处理，加速比达到21倍，基本满足了冠状动脉疾病诊断和治疗的实际需求。

参考文献：

[1] 陈伟伟， 高润霖， 刘力生， 等. 《中国心血管病报告2015》概要[J]. 中国循环杂志， 2016， 31（6）：617-622.

[2] 郭军， 陈韵岱， 田峰， 等. 光学相干断层成像与血管内超声在冠状动脉介入诊疗中的应用[J]. 中国医学影像学杂志， 2012， 20（11）：866-870.

[3] Antonios K， Jurgen L， Karen W， et al. Optical coherence tomography： potential clinical applications[J]. Current Cardiovascular Imaging Reports， 2012， 5（4）：206-220.

[4] Wang Z， Kyono H， Bezerra H G， et al. Automatic segmentation of intravascular optical coherence tomography images for facilitating quantitative diagnosis of atherosclerosis[C]//Proc. SPIE. 2011， 7889： 78890N.

[5] Cao Y H，Jin Q H， Chen Y D， et al.Automatic identification of side branch and main vascular measurements in intracascular optical coherence tomography images[C]//Biomedical Imaging （ISBI 2017）， 2017 IEEE 14th International Symposium on. IEEE， 2017： 608-611.endprint