基于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