任江+賈小云

摘要：為了解決傳統(tǒng)的CPU串行算法對(duì)血管內(nèi)光學(xué)相干斷層掃描（IVOCT）圖像進(jìn)行管腔分割時(shí)耗時(shí)較多的問(wèn)題，提出一種基于圖形處理器（GPU）的IVOCT圖像管腔分割算法。首先，分析基于CPU管腔分割算法可知算法最耗時(shí)的兩部分是坐標(biāo)轉(zhuǎn)換算法獲取極坐標(biāo)圖像和動(dòng)態(tài)規(guī)劃算法獲取管腔內(nèi)壁輪廓。接著，分析坐標(biāo)轉(zhuǎn)換算法和動(dòng)態(tài)規(guī)劃算法的并行性，在統(tǒng)一設(shè)備架構(gòu)（CUDA）下實(shí)現(xiàn)GPU加速優(yōu)化。最后，在MATLAB中進(jìn)行仿真實(shí)驗(yàn)，定量分析GPU加速后的坐標(biāo)轉(zhuǎn)換算法、動(dòng)態(tài)規(guī)劃算法以及最終管腔分割算法的優(yōu)化性能。實(shí)驗(yàn)表明，基于GPU的管腔分割算法較CPU串行算法達(dá)到21倍加速比，能夠?qū)VOCT圖像序列進(jìn)行快速處理，基本滿足了冠狀動(dòng)脈疾病診斷和治療的實(shí)際需求。

關(guān)鍵詞：CUDA；GPU；管腔分割；IVOCT圖像序列；冠狀動(dòng)脈疾病

中圖分類號(hào)：TP391.54 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1009-3044（2017）28-0198-03

1 背景

心血管病是中國(guó)居民的首位死因，心血管病患病率及死亡率仍處于上升階段[1]。血管內(nèi)光學(xué)相干成像（Intravascular Optical Coherence Tomography， IVOCT）是目前國(guó)內(nèi)外較新的冠狀動(dòng)脈內(nèi)影像技術(shù)，分辨率約為10μm，是血管內(nèi)超聲分辨率的十倍[2]。血管管腔區(qū)域信息在冠狀動(dòng)脈疾病的診斷和治療中具有很大的價(jià)值，如評(píng)估冠狀動(dòng)脈的狹窄程度，獲取支架植入的最佳位置，評(píng)估支架的貼壁情況等[3]，因此由IVOCT圖像管腔分割算法獲取血管管腔區(qū)域是冠狀動(dòng)脈疾病的診斷和治療中十分重要的基礎(chǔ)性工作。

文獻(xiàn)[4][5]提出基于CPU的IVOCT圖像管腔分割算法具有較強(qiáng)的魯棒性和準(zhǔn)確性，但是時(shí)間復(fù)雜度較高，處理時(shí)間較長(zhǎng)。在臨床實(shí)踐中，通常需要對(duì)數(shù)據(jù)量很大的IVOCT圖像序列進(jìn)行快速處理，CPU處理速度已無(wú)法滿足實(shí)際需求。近年來(lái)，通用計(jì)算能力不斷增強(qiáng)的圖形處理器（Graphics Processing Unit，GPU）已經(jīng)廣泛地應(yīng)用到醫(yī)學(xué)圖像處理領(lǐng)域。NVIDIA提出的統(tǒng)一設(shè)備架構(gòu)（Compute Unified Device Architecture，CUDA）提高了GPU的可編程性，極大方便了基于GPU的通用計(jì)算程序開(kāi)發(fā)。本文將分析文獻(xiàn)[4][5]提出的基于CPU的IVOCT圖像的管腔分割算法，在CUDA架構(gòu)下對(duì)其中耗時(shí)較多的部分進(jìn)行GPU加速，提出一種基于GPU的IVOCT圖像的管腔分割算法，力求大幅縮短IVOCT圖像管腔分割整體耗時(shí)。

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

文獻(xiàn)[4][5]提出的基于CPU的IVOCT圖像管腔分割算法的主要流程：利用坐標(biāo)轉(zhuǎn)換算法，將直角坐標(biāo)系下的IVOCT圖像轉(zhuǎn)換到極坐標(biāo)系下；在極坐標(biāo)下對(duì)圖像進(jìn)行預(yù)處理來(lái)移除導(dǎo)絲和導(dǎo)管；對(duì)預(yù)處理后的圖像進(jìn)行濾波，獲取能量圖；在能量圖中利用動(dòng)態(tài)規(guī)劃（Dynamic Programming，DP）算法獲取管腔輪廓坐標(biāo)點(diǎn)，從而實(shí)現(xiàn)管腔分割。以單張704×704的IVOCT圖像為處理，在MATLAB中對(duì)基于CPU的管腔分割分割算法進(jìn)行仿真實(shí)驗(yàn)，統(tǒng)計(jì)算法耗時(shí)可知單張IVOCT圖像的管腔分割耗時(shí)約為2.17秒，坐標(biāo)轉(zhuǎn)換算法在總耗時(shí)中的占比最大，占總耗時(shí)的77%，其次是DP算法，占總耗時(shí)的20%

2.1 坐標(biāo)轉(zhuǎn)換算法

設(shè)IVOCT圖像在直角坐標(biāo)下為imR，大小為M×N，極坐標(biāo)系下為imP，大小為Mp×Np。由imR到imP的轉(zhuǎn)換過(guò)程中，選取imR的圖像中心（Om，On）作為極點(diǎn)，Om=（M+1）/2，On=（N+1）/2，imP的原點(diǎn)對(duì)應(yīng)imR的中心點(diǎn)，imP的x方向?qū)?yīng)角度，y方向?qū)?yīng)半徑。對(duì)imP中的像素點(diǎn)（r，c），1≤r≤Mp，1≤c≤Np，通過(guò)公式（1）映射到imR中的坐標(biāo)點(diǎn)（m，n），如果m，n都為整數(shù)，則imP（r，c）的灰度值等于imR（m，n）；如果m，n不為整數(shù)，找到與它相鄰的四個(gè)像素點(diǎn)，雙線性插值獲取對(duì)應(yīng)的灰度值。遍歷imP中的所有的像素點(diǎn)，進(jìn)行相同的轉(zhuǎn)換操作就可以得到極坐標(biāo)下的IVOCT圖像。

2.2 DP算法

獲取管腔內(nèi)壁輪廓可視為在能量圖中搜索一條從第1行到最后一行累積能量最小的路徑，由于輪廓的連續(xù)性，第1行到第i行的最小累積能量路徑與第1行到i-1行的最小累積能量路徑相關(guān)，想要得到第1行到最后一行的路徑就必須依賴于第一行到中間行的最小累積能量路徑，是一種典型動(dòng)態(tài)規(guī)劃問(wèn)題，可以從局部最優(yōu)的最小累積能量路徑一步步遞推到全局最優(yōu)的路徑。用公式（2）定義累積累積能量函數(shù)E：

E（i，j）表示第一行到i行j列像素點(diǎn)的累積能量，Min為最小值函數(shù)，j-m≤ j*≦j+m，m為相關(guān)系數(shù)，用來(lái)確定搜索鄰域的寬度，本文中取m=2。e（i，j）為該像素點(diǎn)的能量，由濾波獲取。通過(guò)上述累積能量函數(shù)獲取第一行到每個(gè)像素點(diǎn)的累積能量后，從最后一行開(kāi)始反向搜索使得累積能量全局最小的路徑，最終獲取管腔內(nèi)壁輪廓的坐標(biāo)點(diǎn)。

3 基于GPU的管腔分割算法

在實(shí)際應(yīng)用中，需要對(duì)100-270幀的IVOCT圖像序列進(jìn)行快速處理，基于CPU的分割算法顯然無(wú)法滿足冠狀動(dòng)脈疾病診斷和治療的實(shí)際需求。因此，我們提出一種GPU加速優(yōu)化方案，利用GPU強(qiáng)大的并行計(jì)算能力優(yōu)化坐標(biāo)轉(zhuǎn)換算法和DP算法，從而縮短IVOCT圖像管腔分割的耗時(shí)。下面介紹針對(duì)坐標(biāo)轉(zhuǎn)換算法和DP算法的GPU優(yōu)化方案。

3.1 坐標(biāo)轉(zhuǎn)換算法的GPU加速優(yōu)化

通過(guò)2.1的介紹可知，IVOCT圖像由直角坐標(biāo)系轉(zhuǎn)換到極坐標(biāo)系的過(guò)程中，imP的每個(gè)像素都需要與imR進(jìn)行坐標(biāo)映射，然后在imR上雙線性插值獲取灰度值。每個(gè)像素執(zhí)行的操作具有高度的并行性而且不存在數(shù)據(jù)依賴，適合在CUDA架構(gòu)下進(jìn)行并行計(jì)算。為imP中的每個(gè)像素都分配一個(gè)CUDA線程，將像素執(zhí)行的坐標(biāo)映射和雙線性插值操作寫在核函數(shù)函數(shù)中，每個(gè)線程并行執(zhí)行Kernel函數(shù)完成計(jì)算就可以得到imP.。CUDA架構(gòu)下實(shí)現(xiàn)的主要步驟：endprint

1） 初始化GPU設(shè)備，為imR，imP分配GPU存儲(chǔ)空間，將存放CPU端內(nèi)存中的imR數(shù)據(jù)拷貝到GPU中，dev_imR，dev_imP分別指向imR，imP顯存空間的首地址。CPU端計(jì)算imR的中心點(diǎn)（Om，On），作為參數(shù)傳入核函數(shù)。

2） 確定Block和Thread數(shù)目。對(duì)于單張IVOCT圖像，選用二維線程網(wǎng)格和二維線程塊，線程塊中的線程數(shù)設(shè)置為16×16，線程網(wǎng)格中線程塊個(gè)數(shù)為（ceil（Np/16））×（ceil（Mp/16）），Mp，Np分別是imP的高度和寬度，ceil函數(shù)是向上取整函數(shù)，確保有足夠多的線程塊來(lái)存放與imP像素個(gè)數(shù)相同的線程。對(duì)于IVOCT圖像序列，選用三維線程網(wǎng)格和二維線程塊，線程塊的大小為16×16，線程網(wǎng)格中線程塊的個(gè)數(shù)（ceil（Np/16））×（ceil（Mp/16））×F，F(xiàn)表示圖像序列的幀數(shù)。

3） 啟動(dòng)核函數(shù)，GPU端的每個(gè)CUDA線程執(zhí)行核函數(shù)完成計(jì)算。單張IVOCT圖像執(zhí)行核函數(shù)ImToPolar_2D，imP中每個(gè)像素點(diǎn)與CUDA線程一一對(duì)應(yīng)，每個(gè)線程負(fù)責(zé)對(duì)應(yīng)像素點(diǎn)的計(jì)算。ImToPolar_2D偽代碼如下：

4） 將顯存中的計(jì)算結(jié)果傳回CPU端內(nèi)存中，釋放開(kāi)辟的GPU顯存空間

3.2 DP算法的GPU加速優(yōu)化

由2.2可知，計(jì)算每個(gè)像素點(diǎn)的累積能量時(shí)，除第一行外，每行像素點(diǎn)的累積能量的計(jì)算都依賴前一行的計(jì)算結(jié)果，不同行中的像素點(diǎn)計(jì)算累積能量時(shí)存在數(shù)據(jù)依賴關(guān)系，但是對(duì)于在同一行不同列的像素，累積能量的計(jì)算不存在數(shù)據(jù)依賴關(guān)系，具有并行性，適合在CUDA架構(gòu)下進(jìn)行并行計(jì)算。對(duì)于M×N能量圖e，基于GPU并行計(jì)算能量累積矩陣E時(shí)，將計(jì)算像素點(diǎn)累積能量的操作寫入核函數(shù)中，執(zhí)行M-1個(gè)核函數(shù)，依次計(jì)算第二行到第M行各行像素的累積能量，每個(gè)核函數(shù)開(kāi)辟N個(gè)線程，這些線程并行計(jì)算對(duì)應(yīng)行中各列像素點(diǎn)的累積能量。CUDA架構(gòu)下實(shí)現(xiàn)的主要步驟：1）初始化GPU設(shè)備，在CPU端創(chuàng)建大小為M×N的累積能量矩陣E，將能量圖e中第一行數(shù)據(jù)復(fù)制給E的第一行，接著將CPU端內(nèi)存中存放的E和e復(fù)制到GPU端的顯存中，dev_e，dev_E分別指向e和E顯存空間的首地址。2）確定Block和Thread數(shù)目。選用一維線程網(wǎng)格和一維線程塊，線程塊中的線程數(shù)設(shè)置為16×16，每個(gè)線程網(wǎng)格中包含ceil（N/16）個(gè)線程塊。3）依次啟動(dòng)M-1個(gè)內(nèi)核函數(shù)DPTrack_GPU，每個(gè)CUDA線程并行執(zhí)行核函數(shù)完成計(jì)算。DPTrack_GPU的偽代碼如下：

4） 將顯存中累積能量的計(jì)算結(jié)果傳回CPU端的內(nèi)存中，接著在CPU端對(duì)累積能量矩陣從最后一行開(kāi)始反向搜索使得累積能量全局最小的路徑，最終獲取管腔內(nèi)壁輪廓。最后釋放開(kāi)辟的GPU顯存空間。

3 實(shí)驗(yàn)結(jié)果與分析

實(shí)驗(yàn)硬件平臺(tái)為Intel Core i5-4460 CPU，內(nèi)核主頻3.20GHz，內(nèi)存為8GB。GPU型號(hào)為NVIDIA Geforce GTX 770，1536個(gè)CUDA處理核心，顯存位寬256位，顯存帶寬224.3GB/s，顯存為2GB。軟件平臺(tái)為Microsoft Windows 7操作系統(tǒng)，Matlab2015b，CUDA Toolkit7.0，CUDA Driver 378.92。

實(shí)驗(yàn)采用的IVOCT圖像來(lái)自FD-OCT（C7-XR）系統(tǒng)。為了進(jìn)行多組數(shù)據(jù)的對(duì)比實(shí)驗(yàn)，對(duì)原始圖像進(jìn)行了預(yù)處理，通過(guò)裁剪獲取6組圖像數(shù)據(jù)，圖像大小分別為128×128，256×256，512×512，1024×1024，2048×2048，4096×4096。表1為坐標(biāo)轉(zhuǎn)換算法和DP算法的CPU與GPU實(shí)現(xiàn)耗時(shí)對(duì)比。

由表1可知基于GPU實(shí)現(xiàn)的兩種并行優(yōu)化算法較CPU都有良好的加速效果。隨著圖像大小的增加，兩種算法的加速比也隨之增大。圖2為兩種算法的加速比增長(zhǎng)對(duì)比圖。

如圖2所示，當(dāng)圖像大小從128×128增加到256×256時(shí)，坐標(biāo)轉(zhuǎn)換算法加速比曲線斜率相對(duì)較小，加速比增長(zhǎng)趨勢(shì)較為平緩。從256×256增加到1024×1024時(shí)加速比曲線斜率增大，加速比增長(zhǎng)較快，加速效果十分明顯。但從1024×1024增加到4096×4096時(shí)，加速比曲線仍然處于上升趨勢(shì)，但是上升趨勢(shì)明顯變平緩。出現(xiàn)這種增長(zhǎng)趨勢(shì)是因?yàn)閳D像大小較小時(shí)，處理的數(shù)據(jù)量較小，CPU可以在相對(duì)較短的時(shí)間內(nèi)完成計(jì)算。同時(shí)，GPU在初始化和主機(jī)端和設(shè)備端數(shù)據(jù)通信需要消耗一定的時(shí)間，這部分時(shí)間在GPU總運(yùn)行時(shí)間中占比較大。然而隨著圖像大小的逐漸增大，CPU的計(jì)算能力趨于飽和，GPU額外開(kāi)銷時(shí)間在GPU總運(yùn)行時(shí)間的占比越來(lái)越小，GPU相對(duì)于CPU的并行計(jì)算能力優(yōu)勢(shì)凸顯，故該階段加速比增長(zhǎng)快。但是，隨著圖像大小的進(jìn)一步增大，GPU的計(jì)算能力也逐漸趨于飽和狀態(tài)，加速比增長(zhǎng)放緩。DP算法的加速比增長(zhǎng)趨勢(shì)類似于直角坐標(biāo)系到極坐標(biāo)系轉(zhuǎn)換算法，但是由于GPU的計(jì)算能力在圖像大小為4096×4096時(shí)尚未趨于飽和，所以加速比上升趨勢(shì)未出現(xiàn)放緩的跡象。對(duì)271幀707×704的IVOCT圖像序列，實(shí)驗(yàn)測(cè)得CPU算法總耗時(shí)約為588秒，本文提出的GPU算法總耗時(shí)約為28秒，達(dá)到21倍的加速比。

3 結(jié)束語(yǔ)

本文通過(guò)對(duì)基于CPU的IVOCT圖像管腔分割算法中耗時(shí)最多的兩部分進(jìn)行GPU加速優(yōu)化，提出了一種基于GPU的IVOCT圖像管腔分割算法。實(shí)驗(yàn)表明，隨著圖像的增大，坐標(biāo)轉(zhuǎn)換算法和 DP算法的加速比也隨之增大。與CPU串行算法相比，基于GPU的IVOCT管腔分割算法能夠?qū)VOCT圖像序列進(jìn)行快速處理，加速比達(dá)到21倍，基本滿足了冠狀動(dòng)脈疾病診斷和治療的實(shí)際需求。

參考文獻(xiàn)：

[1] 陳偉偉， 高潤(rùn)霖， 劉力生， 等. 《中國(guó)心血管病報(bào)告2015》概要[J]. 中國(guó)循環(huán)雜志， 2016， 31（6）：617-622.

[2] 郭軍， 陳韻岱， 田峰， 等. 光學(xué)相干斷層成像與血管內(nèi)超聲在冠狀動(dòng)脈介入診療中的應(yīng)用[J]. 中國(guó)醫(yī)學(xué)影像學(xué)雜志， 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