易 俊，秦曉萌，岑穎珊，劉碧旺，韓定安，王茗祎，周月霞

（佛山科學技術學院 物理與光電工程學院，廣東 佛山 528200）

0 引言

圖1 FD-OCT系統(tǒng)組成Fig.1 FD-OCT system composition

頻域光學相干層析成像技術（Frequency Doman Optical Coherence Tomography, FD-OCT）是一種高軸向分辨率的新興無創(chuàng)三維光學成像技術[1]。自提出至今已經得到了迅速發(fā)展，能實現(xiàn)對生物組織的功能成像，目前在醫(yī)學領域有著廣泛的應用。近年來，隨著FD-OCT 技術的發(fā)展及超高速CMOS 線陣相機的更新，F(xiàn)D-OCT 系統(tǒng)的采集速度基本能滿足實時成像的要求。在FD-OCT 中，成像計算的關鍵步驟是插值計算。由于圖像重建過程中需要對每一線（即線陣相機的每一次曝光采集到的數(shù)據(jù)）數(shù)據(jù)都要進行插值計算并進行正向傅里葉變換。一幀完整的二維圖像至少由500 條線組成，數(shù)據(jù)吞吐量較大，從而加重了重建圖像中數(shù)據(jù)處理環(huán)節(jié)的計算負擔，在Matlab 語言下，基于CPU 算法的程序一幅由500 線組成的FD-OCT 圖像實現(xiàn)3 次樣品插值需要至少60ms 以上。其中插值計算占據(jù)了成像三之二的計算時間。因此，提高插值的計算速度，是提高FDOCT 成像速度的關鍵核心問題。

為了解決FD-OCT 系統(tǒng)在處理數(shù)據(jù)方面遇到的問題，相關研究人員提出了許多的解決方法：利用多核CPU 對FD-OCT 成像系統(tǒng)的數(shù)據(jù)進行并行處理；也有相關研究人員在FD-OCT 系統(tǒng)中增加硬件模塊來加快數(shù)據(jù)處理速度，這樣使得系統(tǒng)的成本大大地增加了[2]。隨著GPU 技術的飛速發(fā)展，它不僅在傳統(tǒng)的圖像處理方面發(fā)揮作用，也早已被用于磁共振成像，系統(tǒng)用于數(shù)據(jù)處理的加速。Watanabe 和Itagaki 將GPU 應用于FD-OCT 系統(tǒng)，利用其對FD-OCT 數(shù)據(jù)做了一維正向傅里葉變換，達到了每秒8 幀[3]。本文研究將經濟實惠的商用GPU 應用到FD-OCT 系統(tǒng)中，大幅度地加快了系統(tǒng)的數(shù)據(jù)處理速度。在每線像素點為1024，每幀1000 線的實驗條件下，成像重建速度達到52 幀每秒，實現(xiàn)了FD-OCT 的實時成像。

1 系統(tǒng)組成和工作原理

1.1 系統(tǒng)的組成

本文用到的FD-OCT 系統(tǒng)組成如圖1 所示，主要包括超輻射發(fā)光激光二極管（λ0=1310nm，Δλ=400nm），樣品臂、參考臂、光譜儀和計算機（包含數(shù)據(jù)采集卡、圖像采集卡和GPU）5 個部分。其中數(shù)據(jù)采集卡為美國國家儀器公司的PCIe-6711，圖像采集卡采用的是來自美國國家儀器公司PCIe-1429；GPU 采用了英偉達公司的Ge Force GTX1060 顯卡（NVIDA CUDA 核心數(shù)1280 個，3GB 顯存）。

1.2 系統(tǒng)的工作原理

FD-OCT 系統(tǒng)的主要構成是麥克爾遜干涉儀和光譜儀。寬帶光源的光被光纖耦合器一分為二，其中一束射向樣品臂，另一束射向參考臂；樣品臂上的后向散射光和參考臂的反射光在耦合器處相遇，產生相干光，這一部分即為麥克爾遜干涉儀。干涉光通過光譜儀分光，而后該干涉光被CCD 相機接收。根據(jù)單色光干涉理論可知，分光鏡處的光強可以表示為：

其中，IR和IS是參考光和樣品光的直流信號，AR為參考光的振幅，AS為樣品光的振幅，zj為等光程面的探測深度，Δφ 為相位差。當干涉光經光柵分光后，CCD 相機的每個線陣單元接收到的強度信號可以表示為：

由于光柵方程：

m 為衍射級次，λ 為衍射波長，φ 為入射角，θ 為衍射角以及d 為光柵常數(shù)。根據(jù)該光柵方程對入射光譜進行分光[4]。干涉光通過光柵分為不同波長的光，之后衍射到線陣CCD 上。因此，相機記錄到的信號是一個以波長λ 為參變量的光強信號，根據(jù)干涉公式（2），實際需要獲得的信號是以波矢k 為參變量的信號。由于將波長信號插值到波矢信號是FD-OCT 成像的關鍵步驟，需要對相機采集到的光譜信號進行插值，獲得真正的干涉信號。占據(jù)了成像的三分之二的計算時間，因而提高插值的計算速度是成像的關鍵核心問題。本文應用了基于GPU 的CUDA 語言編寫了一套三次樣條插值算法，加速了整個插值計算過程。

1.3 三次樣條插值

三次樣條插值（Cubic Spline Interpolation），是一種應用廣泛的插值函數(shù) 。三次樣條插值函數(shù)是由多段的低次多項式組成，并且它的曲線非常光滑；得到插值點處的函數(shù)值及提供（n-1）個邊界節(jié)點處的導數(shù)信息就的三次樣條插值函數(shù)[5]。

1.3.1 三次樣條插值函數(shù)的定義[6-8]

圖2 三次樣條插值求解流程圖Fig.2 Three-spline interpolation solution flowchart

假設有以下兩個節(jié)點x:a=x0＜x1＜...＜xn=b，y:y0y1...yn。樣條曲線是一個分段定義的函數(shù)。在n+1 個數(shù)據(jù)點中，共有n 個區(qū)間，三次樣條插值函數(shù)滿足以下條件：1）在每個區(qū)間[xi,xi+1]上，S(x)都是三次多項式；2）滿足S(xi)=yi；3）S(x)的導數(shù)S'(x)及二階導數(shù)S"(x)在區(qū)間[a,b]上都是連續(xù)的，即S(x)是光滑的曲線[9]。所以n 個三次多項式分段可以寫作：

其中ai，bi，ci，di代表曲線擬合的4 個系數(shù)。

1.3.2 三次樣條函數(shù)的求解

求解三次樣條插值函數(shù)是利用n+1 個數(shù)據(jù)節(jié)(x0,y0),(x1,y1),(x2,y2),…,(xn,yn)先計算步hi=xi+1-xi，令mi=Si"(xi)，由插值條件Si(xi+1)=yi+1以及Si(x)在樣點xi處二階導數(shù)的條件，增加自然邊界條件：

得到如下矩陣方程：

解矩陣方程（7），求出二次微分值mi[10]，從而求得系數(shù)ai，bi，ci，di并將系數(shù)代入獲得擬合曲線，獲得待求插值點的數(shù)值，最終完成三次樣條插值函數(shù)求解。

2 基于GPU并行計算的數(shù)據(jù)處理模塊的實現(xiàn)

2.1 CUDA：并行計算架構

CUDA（Compute Unified Device Architecture）是一種并行計算架構[11]。CUDA 將C 語言作為編程語言，并提供大量的計算指令，使得效率更高的密集數(shù)據(jù)計算解決方案能夠在GPU 突出計算能力的基礎上建立起來[12]。

CUDA 構架分為兩部分：Host 和Device。通常而言，Host 指的是主機，Device 指的是GPU[13]。在CUDA 構架中，CPU 作為主機，負責執(zhí)行程序的串行部分；GPU 作為協(xié)助處理器主要負責對密集且大量的數(shù)據(jù)進行并行計算[14]。用CUDA 編寫好的程序叫做核（kernel）函數(shù)。CUDA 允許程序員定義稱為核的C 語言函數(shù)或使用CUDA 指令，在調用此類函數(shù)時，它將由N 個不同的CUDA 線程并行執(zhí)行N 次，執(zhí)行核的每個線程都會被分配一個獨特的線程ID，可通過內置的threadIdx 變量在內核中訪問此ID[15]。在 CUDA 程序中，主程序在調用GPU 內核之前，務必對核執(zhí)行配置操作，以此確定線程塊數(shù)和每個線程塊中的線程數(shù)以及共享內存大小。

2.2 數(shù)據(jù)處理流程

在FD-OCT 系統(tǒng)處理過程中，首先，CCD 相機采集到的原始數(shù)據(jù)經過圖像采集卡傳輸?shù)紺PU 內存上，然后再傳輸?shù)紾PU 顯存上面進行數(shù)據(jù)處理；其中數(shù)據(jù)處理的主要過程包括對采集數(shù)據(jù)進行去背景、色散補償、波長空間（λ空間）到波數(shù)空間（k 空間）的轉換、插值、FFT 變換、取模和取對數(shù)。最后，把處理好的數(shù)據(jù)再傳輸?shù)紺PU 內存上并由計算機進行顯示。由于插值過程計算復雜且系統(tǒng)采集到的圖像數(shù)據(jù)每列是相互獨立的，處理的時候也是分開一列一列數(shù)據(jù)進行處理的，所以可以基于GPU 強大的并行計算能力對數(shù)據(jù)進行插值計算，并把CPU 計算插值部分代碼改寫成在GPU 上執(zhí)行的kemel 函數(shù)，從而加快數(shù)據(jù)處理的速度，使成像時間減少。

2.3 程序設計

首先，定義整個數(shù)據(jù)處理過程中需要的變量，并給他們分配對應大小的內存空間；其中Host 端用malloc 函數(shù)開空間，Device 端用CUDA 內部指令cudaMalloc 開空間；利用cudaMemcpy 將原始數(shù)據(jù)以及運算所需中間變量從主機內存復制到設備內存。根據(jù)分配的線程空間借用對應ID 進行數(shù)據(jù)的調用和處理，其中包括調整數(shù)據(jù)類型以及去除背景噪聲。其次，由于CUDA 編程沒有自帶插值的庫函數(shù)，本文編寫了一套三次樣條插值算法并用CUDA 語言編寫成在GPU 上運行的kernel 函數(shù)。

圖2 為三次樣條插值函數(shù)求解的流程圖，首先計算出步長，然后通過計算得到三對角矩陣，利用CUDA 自帶的庫函數(shù)求解三對角矩陣[16]；最后用求解矩陣得到的結果去計算多項式擬合系數(shù)，并將系數(shù)代入得到擬合曲線，獲得待求插值點的數(shù)值，完成三次樣條插值求解。

按照三次樣條插值的原理以及流程圖，以下為本文編

寫的基于GPU 三次樣條插值具體程序代碼：

//計算步長

__global__ void Kernel001(float *K, float *ks, int *index,float *h, float *dx, float *x31, float *xn)

{

h[idx] = K[idx] - ks[index[idx]]；

}

//得到三對角矩陣

__global__ void Kernel002(float *dx, float *x31, float *xn,float *A, float *B, float *C)

{

int idx = threadIdx.x + blockDim.x * blockIdx.x；

C[idx] = x31[0]； B[idx] = dx[1]； A[idx] = 0；

C[idx] = dx[idx - 1]； B[idx] = 2 * (dx[idx] + dx[idx - 1])；

A[idx] = dx[idx]；

C[idx] = 0； B[idx] = dx[idx - 2]；

A[idx] = xn[0]；

}

//三對角矩陣求解

cusparseCreate(&cusparseH)；

cusparseSgtsvInterleavedBatch(cusparseH,0, pixel, d_dl0,d_d0, d_du0, d_b, line*pic, d_work)；

以上兩條函數(shù)其專門用于求解三對角矩陣源于CUDA自帶的庫函數(shù)CUSPARSE LIBRARY。第一條用于創(chuàng)建句柄，第二條負責運算求解矩陣。調用GPU 自帶函數(shù)求解矩陣，極大地提高了數(shù)據(jù)處理的速度以及精確度。

//用求解矩陣結果計算多項式擬合系數(shù)

_global__ void Kernel0061(float *divdif, float *s, float *dx,float *Y, float *dd, float *cc, float *bb, float *aa)

{

int idx = threadIdx.x + blockDim.x * blockIdx.x；

int idy = threadIdx.y + blockDim.y * blockIdx.y；

int idz = threadIdx.z + blockDim.z * blockIdx.z；

if (idy ＜ line && idx ＜ pixel - 1 && idz ＜ pic)

{ aa[idx + (pixel - 1)*idy + idz*s_frame] = Y[idx +pixel*idy + idz*frame]；

bb[idx + (pixel - 1)*idy + idz*s_frame] = s[idx + pixel*idy+ idz*frame]；

cc[idx + (pixel - 1)*idy + idz*s_frame] = 2 * (divdif[idx +(pixel - 1)*idy + idz*s_frame] - s[idx + pixel*idy + idz*frame]) /dx[idx] - (s[idx + 1 + pix-el*idy + idz*frame] - divdif[idx + (pixel- 1)*idy + idz*s_frame]) / dx[idx]；

dd[idx + (pixel - 1)*idy + idz*s_frame] = ((s[idx + 1 +pixel*idy + idz*frame] - divdif[idx + (pixel - 1)*idy + idz*s_frame]) / dx[idx] - (divdif[idx + (pixel - 1)*idy + idz*s_frame] -s[idx + pixel*idy + idz*frame]) / dx[idx]) / dx[idx]；

}

}

//系數(shù)代入得到擬合曲線，獲得待求插值點數(shù)值

__global__ void Kernel007(float *h, int *index, float *Y,float *dd, float *cc, float *bb, float *aa)

{ int idx = threadIdx.x + blockDim.x * block-Idx.x；

int idy = threadIdx.y + blockDim.y * blockIdx.y；

int idz = threadIdx.z + blockDim.z * blockIdx.z；

if (idx ＜ pixel && idy ＜ line && idz ＜ pic) {

Y[idx + pixel*idy + idz*frame] = h[idx] * (h[idx] * (h[idx] *dd[index[idx] + idy*(pixel - 1) + idz*s_frame] + cc[index[idx] +idy*(pixel - 1) + idz*s_frame]) + bb[index[idx] + idy*(pixel - 1)+ idz*s_frame]) + aa[index[idx] + idy*(pixel - 1) + idz*s_frame]；

}

}

然后，應用CUDA 自帶的cuFFT 庫函數(shù)對插值后的數(shù)據(jù)做快速傅里葉變化，并改寫成kernel 函數(shù)；同時添加對快速傅里葉變換后數(shù)據(jù)進行取模和取對數(shù)處理的核函數(shù)，實現(xiàn)FD-OCT 系統(tǒng)數(shù)據(jù)在GPU 完成所有數(shù)據(jù)處理環(huán)節(jié)。最后，將最終處理完的數(shù)據(jù)再次用cudaMemcpy 從顯存復制到內存并在計算機上顯示和釋放所有顯存以及內存空間。

3 程序驗證以及實驗結果分析

本系統(tǒng)采用Microsoft Visual Studio2015 中集成CUDA Toolkit 64 bit 9.2 和Nvidia Driver for Win-dows7 64bit 為開發(fā)環(huán)境，分別用基于CPU 模式的Matlab 語言和基于GPU 模式的CUDA 語言對采集到的人體表皮圖像數(shù)據(jù)進行數(shù)據(jù)處理得到圖像。Matlab 語言是基于CPU 數(shù)據(jù)處理單元的串行數(shù)據(jù)處理方法，而CUDA 語言則用到了GPU 突出的并行計算能力來實現(xiàn)成像。

3.1 判斷依據(jù)

為了比較兩種模式成像速度的快慢，可以通過比較兩種模式下處理數(shù)據(jù)和顯示每幀圖總用時。因此，對于兩種模式，在設計程序時添加了幀率計算及顯示和記錄，以及通過CUDA 的計時函數(shù)和Matlab 的計時函數(shù)可得出數(shù)據(jù)處理每個環(huán)節(jié)所用時間以及總用時，這方便比較兩者的成像速度。

3.2 實驗結果及分析

本文用B 掃描模式對人體表皮進行了成像實驗（采集到的每幀圖像數(shù)據(jù)大小為1000 線×1024 像素/線×12 字節(jié)/像素），實驗成像效果如圖3 所示，圖(a)為基于CPU模式的Matlab 語言B 掃描成像截面圖，圖(b)為基于GPU模式的CUDA 語言B 掃描成像截面圖，圖(c)對應于(a)圖中白色虛線的強度圖，圖(d)對應于(b)圖中白色虛線的強度圖。由圖3 可知，基于CPU 模式的Matlab 語言和基于GPU 模式的CUDA 語言處理數(shù)據(jù)得到的圖像質量無明顯差異，且經過強度圖對比，數(shù)據(jù)相同。

圖3 人體表皮B掃描實驗圖(1000Lines)Fig.3 Human epidermis B scan experiment (1000Lines)

圖4 不同大小圖片在基于CPU模式的Matlab語言和基于GPU模式的CUDA語言的成像時間對比Fig.4 Comparison of image times for pictures of different sizes in the CPU-based matlab language and the GPU-based CUDA language

除此之外，本文還改變了每幀圖像的線數(shù)，得到兩種模式的成像時間如圖4 所示。從圖中可以得出，在不同線數(shù)每幀圖的情況下，基于GPU 模式的CUDA 語言的成像時間都要比基于CPU 模式的Matlab 要短15 倍左右。特別是隨著數(shù)據(jù)量的變大，應用GPU 進行插值計算以及成像的速度更快，用時更短；而基于CPU 的Matlab 語言則隨著數(shù)據(jù)量的增加變得更慢，進一步說明了把GPU 應用到FD-OCT系統(tǒng)中的插值計算以及并行數(shù)據(jù)處理極大地提高了FDOCT 成像的速度。

4 結論

本文在沒有改變任何硬件的FD-OCT 系統(tǒng)的基礎上，將基于GPU 的CUDA 語言應用到系統(tǒng)的數(shù)據(jù)處理過程中，特別是針對三次樣條插值巨大且重復的數(shù)據(jù)計算，很好地利用了GPU 突出的并行計算能力，將插值計算以及需要并行計算的數(shù)據(jù)處理過程用CUDA 進行改寫，使得整個的成像時間縮短了一個數(shù)量級。在國內已有的相關研究中，劉巧艷，劉銳，朱珊珊等人也把GPU 應用到OCT 成像系統(tǒng)的數(shù)據(jù)處理中，取得了不錯的成績。本文著重是應用GPU 對整個數(shù)據(jù)處理過程的插值計算部分進行加速和優(yōu)化以達到整個數(shù)據(jù)處理過程加速，速度比單一的基于CPU 的Matlab語言處理速度加快了將近15 倍。所以說明GPU 加速模式很好地解決了傳統(tǒng)基于CPU 的Matlab 語言插值計算慢以及不能實時成像的問題，實現(xiàn)了FD-OCT 系統(tǒng)實時的2D 成像。但是，整個程序的圖像顯示環(huán)節(jié)沒有實現(xiàn)并行化處理，整個成像速度還有提升的空間，需要進一步探索，更大程度地發(fā)揮GPU 的突出并行計算能力。本文研究為FD-OCT系統(tǒng)三維實時成像打下了堅實的基礎。