(浙江工業(yè)大學(xué) 信息工程學(xué)院,浙江 杭州 310023)
光學(xué)相干層析成像技術(shù)[1-3](Optical coherence tomograph,OCT)是一種新型的光學(xué)成像技術(shù),它利用弱相干光干涉儀的基本原理,檢測生物組織不同深度層面對入射弱相干光的背向散射信號,能夠得到生物組織的二維或三維結(jié)構(gòu)圖像。它采用紅外光掃描,穿透不深,在一般的生物組織中成像深度只有2~3 mm,但它的辨識能力比CT高了1 000 倍,分辨率可達(dá)到1~15 μm,被應(yīng)用于眼、皮膚和牙齒等的高分辨率成像,可有力輔助醫(yī)生進(jìn)行臨床診斷和治療。近年來,OCT被逐漸應(yīng)用到手指生物特征采集[4]上,傳統(tǒng)的手指特征采集主要通過按壓的方式直接獲取二維的指紋圖像,然后根據(jù)指紋圖像的全局特征和細(xì)節(jié)特征來對人的身份進(jìn)行鑒定,而OCT則借助弱相干光,從空氣入射到手指表皮再深入到皮下組織,除了采集表皮指紋,還采集了皮下汗腺、真皮層內(nèi)指紋這些新的生物特征信息。
三維可視化[5-6]是促進(jìn)對OCT采集的數(shù)據(jù)進(jìn)行分析和研究的有效手段,如吳凌[7]用OCT獲得指甲邊緣的三維圖像和橘子表皮的三維圖像;萬木森[8]用OCT獲得離體牙釉質(zhì)和牙本質(zhì)的三維圖像;Kwon等[9]用OCT獲得了血流的三維圖像。但隨著OCT采樣速率的逐步提高,數(shù)據(jù)采集量龐大,且數(shù)據(jù)處理過程包括許多次FFT(Fast Fourier transform)變換,該算法計算復(fù)雜度高,耗時嚴(yán)重,難以達(dá)到實時成像的要求,針對這一問題,Watanabe等[10]最先提出將GPU與硬件相結(jié)合,應(yīng)用于OCT系統(tǒng)成像,提高了數(shù)據(jù)處理速度;Jeught等[11]在不改變外界硬件的基礎(chǔ)上實現(xiàn)了利用GPU完成重采樣過程的加速。筆者利用頻域OCT系統(tǒng)采集得到手指樣品的原始數(shù)據(jù),借助GPU的高速并行運算能力對OCT數(shù)據(jù)處理過程進(jìn)行加速,解決了CPU串行處理方式速度慢造成無法實時成像的問題。為了能夠更加直觀準(zhǔn)確地展示手指皮下組織的空間分布和幾何形狀,筆者利用三維可視化算法中的光線投射算法來展示手指皮下組織的三維結(jié)構(gòu),成像結(jié)果能很好地反映手指皮下組織的形狀特征和空間分布。最后分別掃描活體手指和帶假手指套的手指得到原始數(shù)據(jù),經(jīng)數(shù)據(jù)處理后進(jìn)行三維可視化,比較發(fā)現(xiàn):它們的成像結(jié)果都可以顯示表皮指紋,但手指皮下組織結(jié)構(gòu)卻存在著明顯的差異,實驗結(jié)果證實了OCT技術(shù)在指紋特征采集和指紋防偽領(lǐng)域[12-13]具有廣闊的應(yīng)用前景。
頻域OCT系統(tǒng)主要由光源、樣品臂、參考臂和光譜儀等幾大部分組成。所用系統(tǒng)如圖1所示,寬帶光源SLD發(fā)出的低相干光(中心波長λ為848 nm,帶寬Δλ為46 nm)經(jīng)4 端光纖耦合器分成2 束光,以50/50的比例分別進(jìn)入?yún)⒖急酆蜆悠繁邸_M(jìn)入?yún)⒖急壑械墓饨?jīng)過準(zhǔn)直透鏡準(zhǔn)直后,平行地入射到全反射鏡上,并被平行地反射回來作為參考光。進(jìn)入樣品臂中的光聚焦在樣品上,通過掃描振鏡系統(tǒng)對手指樣品進(jìn)行掃描,被手指樣品中不同深度的散射顆粒反射成為信號光。從參考臂和樣品臂中返回的參考光和信號光在光纖耦合器中匯合,發(fā)生疊加和干涉,干涉信號從耦合器另一端出射進(jìn)入光譜儀進(jìn)行解析,被線陣CCD采集,通過數(shù)據(jù)采集卡送入計算機(jī)中進(jìn)行一系列的數(shù)據(jù)處理,然后將處理后的數(shù)據(jù)進(jìn)行2D或3D成像。
圖1 OCT系統(tǒng)原理結(jié)構(gòu)示意圖Fig.1 OCT system principle structure diagram
頻域OCT系統(tǒng)主要通過掃描振鏡系統(tǒng)完成對手指樣品的橫向掃描,首先對手指樣品上的一個點進(jìn)行一次軸向掃描(A-SCAN),CCD將采集到一條包含深度信息的線性光譜,由于手指皮下組織對光的散射和吸收,本系統(tǒng)的成像深度大約為3 mm。然后對手指樣品進(jìn)行橫向掃描(B-SCAN),獲得多條A-SCAN組成的深度信息光譜,可以通過數(shù)據(jù)處理重構(gòu)出一幅二維圖像,如果依次往后重復(fù)地進(jìn)行橫向掃描可以得到一個體掃描數(shù)據(jù),這個體數(shù)據(jù)包含了被測手指樣品的三維結(jié)構(gòu)信息。筆者所用的頻域OCT系統(tǒng)的軸向分辨率是7 μm,橫向分辨率是11 μm,CCD的像素點為2 048 個,則一個A-SCAN獲得的像素個數(shù)為2 048,一個B-SCAN包含500 個A-SCAN信號,即對手指樣品某一方向進(jìn)行了500 次深度數(shù)據(jù)的獲取,掃描速率是36 klines/s,每個B-SCAN大小為500 lines×2 048 pixels/line×2 bytes/pixel,共計1.9 MB數(shù)據(jù),300 個B-SCAN可以組成一個大小為570 MB的體數(shù)據(jù),數(shù)據(jù)量很大,這將造成數(shù)據(jù)處理過程耗時嚴(yán)重。
本系統(tǒng)的數(shù)據(jù)處理過程如圖2所示,除了對500 條光譜求平均,數(shù)據(jù)處理是一條線一條線進(jìn)行的,每條線包括數(shù)據(jù)截取、插值、減直流、FFT變換、取模和取對數(shù)這些操作。數(shù)據(jù)截取是指一個A-SCAN有2 048 個像素,只截取400~1 850的1 451 個像素;由于光譜儀測量到的均勻干涉信號光譜I(λ)經(jīng)空間轉(zhuǎn)換后變成了非均勻的波數(shù)空間函數(shù)I(k),而FFT變換要求輸入數(shù)據(jù)是均勻的,需要通過插值變換將1 451 個數(shù)據(jù)變成4 096 個均勻采樣的數(shù)據(jù);減直流是為了去除參考臂和樣品臂返回光的直流信號以及樣品不同反射層返回光之間的干涉信號,它們屬于需要去除的背景光;FFT變換可以解析出數(shù)據(jù)的深度信息,是數(shù)據(jù)處理過程中最重要的一個步驟;FFT變換之后,對前500 個數(shù)據(jù)取模得到光功率譜的幅值,由于CCD采集的數(shù)據(jù)會存在圖像灰度值偏低的問題,取對數(shù)可以對圖像進(jìn)行銳化。
圖2 OCT數(shù)據(jù)處理流程圖Fig.2 OCT data processing flow chart
OCT系統(tǒng)的數(shù)據(jù)采集量大,數(shù)據(jù)處理過程中的計算量也很大,一個體掃描數(shù)據(jù)(500×2 048×300)需要執(zhí)行1 500 次數(shù)據(jù)截取,1 501 次插值(500 條光譜求平均后也需要插值),1 500 次減直流,1 500 次FFT變換,1 500 次取模和1 500 次取對數(shù)。其中,一次數(shù)據(jù)截取要對1 451 個數(shù)據(jù)進(jìn)行計算,一次插值、減直流和FFT要分別對4 096 個數(shù)據(jù)進(jìn)行計算,一次取模和取對數(shù)要分別對500 個數(shù)據(jù)進(jìn)行計算。所有的計算中,F(xiàn)FT算法復(fù)雜度最高,耗時最嚴(yán)重,如果采樣傳統(tǒng)的CPU串行處理方式,將難以滿足實時成像的需求。鑒于一個B-SCAN中的500 個A-SCAN相互獨立,且在進(jìn)行數(shù)據(jù)截取、插值、減直流、FFT變換、取模和取對數(shù)這些計算時,大量單個的數(shù)據(jù)元可以獨立完成計算,因此適合在GPU中并行處理。
CUDA(Computer unified device architecture)是NVIDIA在2006 年為GPU上的并行編程提供的編程API和架構(gòu)。在CUDA模型中,CPU作為主機(jī)(Host),負(fù)責(zé)內(nèi)存和顯存的分配、數(shù)據(jù)傳輸以及內(nèi)核的啟動,GPU作為協(xié)處理器或者設(shè)備(Device),負(fù)責(zé)并行部分的大量密集型數(shù)據(jù)的計算。
GPU中需要完成的主要有數(shù)據(jù)截取、插值、減直流、FFT變換、取模和取對數(shù)這些計算,將它們改寫成Kernel函數(shù),一個Kernel函數(shù)對應(yīng)一個線程網(wǎng)格(Grid),線程網(wǎng)格由大量線程塊(Block)組成,線程塊由大量線程(Thread)構(gòu)建。具體的線程布局如表1所示,Grid(x,y,1)用于定義Grid的維度和尺寸,即一個Grid有多少個Block,Block(x,y,z)用于定義Thread的維度和尺寸,即一個Block有多少個Thread。在線程布局時,定義的Grid和Block都是一維的,線程中的線程數(shù)量都是一個束(Warp)的整數(shù)倍,即32的整數(shù)倍,因為GPU的執(zhí)行是以Warp為單位進(jìn)行調(diào)度,如果Block中的線程數(shù)不設(shè)為32的整數(shù)倍,則最后一個Warp中的部分線程是沒有用的。線程布局完成后,每個數(shù)據(jù)元都有唯一對應(yīng)線程和線程塊索引。其中,F(xiàn)FT運算部分采用CUDA自帶的CUFFT庫實現(xiàn),首先用CufftHandle plan定義一個句柄plan,然后利用Cufftplan1d函數(shù)創(chuàng)建一個一維快速傅立葉變換句柄,設(shè)置好plan后,運行CufftExexcC2C函數(shù),即可在GPU上完成FFT運算,由于plan可以根據(jù)輸入數(shù)據(jù)的大小預(yù)先配置好內(nèi)存和計算資源,在運算時處理器可達(dá)到最佳性能。具體程序設(shè)計如下:1) CPU上讀取OCT采集到的原始bin數(shù)據(jù)到內(nèi)存;2) 利用CUDA中的cudaMalloc函數(shù)在GPU上分配顯存,由于筆者將300 個B-SCAN依次在GPU中并行處理,故分配的顯存大小由一個B-SCAN的數(shù)據(jù)量決定;3) 利用CUDA中的cudaMemcpy函數(shù)將CPU上準(zhǔn)備好的數(shù)據(jù)拷貝至顯存,每次拷貝的數(shù)據(jù)量大小為一個B-SCAN的數(shù)據(jù)量;4) CPU上調(diào)用Kernel函數(shù),GPU上執(zhí)行并行計算,完成數(shù)據(jù)截取、插值、減直流、FFT變換、取模和取對數(shù)計算,其中,F(xiàn)FT運算部分采用CUDA自帶的CUFFT庫來實現(xiàn);5) 利用CUDA中的cudaMemcpy函數(shù)將GPU顯存中的結(jié)果拷貝至CPU,300 個B-SCAN的計算結(jié)果存儲在三維紋理數(shù)組中,得到一個體數(shù)據(jù),為后文的三維可視化做準(zhǔn)備。
表1 線程布局Table 1 Thread layout
以Microsoft Visual Studio 2010集成CUDA 9.1為開發(fā)環(huán)境,GPU為NVIDIA GeForce GTX 1050 Ti,計算能力6.1?,F(xiàn)將筆者提出的GPU并行處理方式與Matlab處理方式和CPU串行處理方式進(jìn)行對比。處理一個體掃描數(shù)據(jù)(由300 個B-SCAN數(shù)據(jù)構(gòu)成,共570 M)的時間測試結(jié)果如表2所示。其中Matlab處理方式耗時最長,總共需要45 000 ms。CPU串行處理方式中用到了FFTW庫,F(xiàn)FTW是由麻省理工學(xué)院超級計算技術(shù)組開發(fā)的一套離散傅里葉變換(DFT)的計算庫,被稱為是最快的FFT,由于數(shù)據(jù)處理過程中要執(zhí)行很多次FFT變換,耗時嚴(yán)重,運用FFTW庫以后可大大縮減計算時間,整個處理過程總共需要7 700 ms。GPU并行處理方式運用了GPU強(qiáng)大的并行計算能力,F(xiàn)FT運算部分采用CUDA自帶的CUFFT庫來實現(xiàn),它由FFTW庫發(fā)展而來,同樣具有很快的計算速率,除此之外,數(shù)據(jù)截取、插值、減直流、取模和取對數(shù)等操作都改寫成Kernel函數(shù),在GPU上并行實現(xiàn),進(jìn)一步提升了處理速度,最后總共耗時2 476 ms。比較發(fā)現(xiàn):筆者提出的方法耗時最短,速度是Matlab處理方式的18 倍,是CPU串行處理方式(運用FFTW庫)的3 倍。OCT數(shù)據(jù)的實時三維成像要求OCT系統(tǒng)剛采集完一個體掃描數(shù)據(jù),計算機(jī)屏幕上就能顯示樣品的三維圖像,由于所用OCT系統(tǒng)的掃描速率是36 klines/s,采集完一個體掃描數(shù)據(jù)(300 個B-SCAN,共150 000 個A-SCAN)大約需要4 167 ms,如果一邊采集數(shù)據(jù)一邊對每個B-SCAN進(jìn)行數(shù)據(jù)處理,數(shù)據(jù)處理的時間將會被隱藏,這種情況下,Matlab處理和CPU串行處理時間都大于采集時間,而且所用系統(tǒng)后期還會改進(jìn),采集速率會有進(jìn)一步的提升,因此不滿足要求,而GPU并行處理時間小于采集時間,能滿足實時三維成像的要求。
表2 不同方式下數(shù)據(jù)處理時間對比Table 2 Comparison of data processing time in different ways
處理后的OCT數(shù)據(jù)被保存到三維紋理數(shù)組中,一個三維紋理數(shù)組由300 個二維數(shù)組構(gòu)成,它們是手指的二維斷層圖像序列,圖3是其中的一幅二維斷層圖像,從淺到深依次是表皮層、真皮層。從表皮層中可以看到角質(zhì)層和汗腺,角質(zhì)層的圖像灰度更高,汗腺與手指表面的汗孔相連,真皮層最厚,在它與表皮層的分界處呈波浪狀的是乳頭層。二維圖像只能反映手指的截面信息,而三維圖像可以反映手指皮下組織的三維空間結(jié)構(gòu),較之二維圖像更加直觀準(zhǔn)確,現(xiàn)利用移動立方體算法、最大密度投影算法和光線投射算法對手指皮下組織進(jìn)行三維可視化,并比較渲染效果的差異。
圖3 手指二維斷層圖像Fig.3 Finger 2D tomographic image
現(xiàn)采用MC算法、MIP算法和Ray-casting算法分別對手指皮下組織進(jìn)行三維可視化,圖4展示了不同算法下任意視角的手指皮下組織三維結(jié)構(gòu),可以發(fā)現(xiàn):手指皮下組織主要由表皮指紋、汗腺和真皮層3 部分構(gòu)成。
圖4 不同算法下手指皮下組織可視化效果的比較Fig.4 Visualization of subcutaneous tissues by different algorithms
汗腺是手指皮下組織中非常典型的組織結(jié)構(gòu),筆者主要依據(jù)汗腺來比較可視化效果。汗腺的繪制結(jié)果如圖5(a~c)所示,與MC算法和MIP算法相比,Ray-casting算法繪制的圖像具有明顯的層次性,更清晰地反映了皮下汗腺的空間分布與結(jié)構(gòu)特征。
圖5 不同算法下皮下汗腺可視化效果的比較Fig.5 Visualization of subcutaneous sweat glands under different algorithms
MC算法在求等值面的過程中,將等值面與立方體棱邊的交點簡單相連得到三角面片,最后將三角面片連接得到逼近的等值面,但這個等值面只是一個近似面,有時會不準(zhǔn)確,如果數(shù)據(jù)場的數(shù)據(jù)量很少且間隙很大,最終繪制的圖像誤差將會非常大。但筆者的數(shù)據(jù)來自于OCT系統(tǒng),分辨率很高,在這種情況下,立方體中的三角面片非常微小,最終提取的等值面不會有太大誤差,所以圖5(a)所示的皮下汗腺的形狀及分布基本與實際符合,但由于MC算法通常只適用于表面特征分明的器官和組織,而汗腺的形態(tài)特征并不分明,所以最后繪制的圖像不能清晰地反映汗腺的細(xì)節(jié)特征。
MIP算法實現(xiàn)較簡單,繪制速度優(yōu)于另外兩種算法,但這也導(dǎo)致繪制效果較差。圖5(b)中的汗腺疊加在一起,幾乎觀察不到汗腺的細(xì)節(jié)特征,看起來更像是一張X光片。由于血管造影對細(xì)節(jié)要求不高,所以該算法經(jīng)常應(yīng)用于三維血管結(jié)構(gòu)的顯示。
從圖5(c)可以看出:Ray-casting算法繪制的圖像清晰地展示了皮下汗腺的形狀特征和層次分布,具有非常好的形狀感知和深度感知,這是因為Ray-casting算法考慮了所有體素對結(jié)果圖像的貢獻(xiàn),且OCT數(shù)據(jù)具有很高的分辨率,特征信息豐富,所以繪制的圖像充分表達(dá)了體數(shù)據(jù)中原有的信息,圖像質(zhì)量比MC算法和MIP算法高很多。
現(xiàn)將一些細(xì)節(jié)汗腺放大作進(jìn)一步的比較分析。圖6(a,b)為MC算法下繪制的細(xì)節(jié)汗腺,圖6(c,d)分別為MIP算法、Ray-casting算法下繪制的細(xì)節(jié)汗腺。從圖6(a)中可以發(fā)現(xiàn)有的汗腺被斬斷了,即不具有連通性,這是由于MC算法存在面二義性問題,當(dāng)值為1的角點和值為0的角點分別位于對角線的兩端時,會有兩種可能的連接方式,這個面就叫做二義性面,如果把它作為兩個立方體的公共面,且兩個立方體在公共面上采用了不同的連接方式,那么最后構(gòu)造的等值面就會出現(xiàn)“空洞”現(xiàn)象,就像圖6(a)中一樣,汗腺發(fā)生了斷裂;圖6(b)中有一根汗腺呈螺旋狀,由于OCT數(shù)據(jù)的高分辨率,它跟實際的形狀已十分接近,但是與圖6(d)中的螺旋汗腺相比,它顯得毫無立體感,這是由MC算法本身的缺陷造成的,在獲取等值面的過程中,將等值面與立方體棱邊的交點簡單連接成三角面片,這樣必然會丟失掉一些細(xì)節(jié)信息,比如立方體內(nèi)的環(huán)形結(jié)構(gòu)和臨界點(等值面改變的點),就像圖6(b)中的螺旋狀汗腺一樣,螺旋處信息單薄,導(dǎo)致形狀不夠逼真;圖6(c)中的汗腺缺乏形狀感知和深度感知,可以發(fā)現(xiàn)能看到3 根汗腺的基本輪廓但不清晰,而且感受不到汗腺的前后順序,反而覺得它們處于同一平面,這是因為MIP算法只投射光線上的最大密度值點,這樣會導(dǎo)致圖像中缺少深度信息,甚至一些位于后面的組織由于密度大于前面的組織,顯示時出現(xiàn)后面的組織遮擋前面組織的情況,給用戶造成混淆。
圖6 汗腺細(xì)節(jié)放大圖Fig.6 Sweat gland detail larger image
盡管Ray-casting算法繪圖質(zhì)量很高,但由于體素數(shù)量和投射光線數(shù)量過多,導(dǎo)致Ray-casting算法繪制速度緩慢,但隨著國內(nèi)外研究者提出了許多關(guān)于Ray-casting算法的加速技術(shù)[14-16],這個問題漸漸得到解決,基本可以實現(xiàn)實時交互操作。綜上所述,最終選擇Ray-casting算法來實現(xiàn)手指皮下組織的三維可視化。
利用頻域OCT系統(tǒng)對活體手指和帶假手指套(2,3 mm厚)的手指分別進(jìn)行掃描,掃描面積為5 mm×3 mm,得到3 組實驗數(shù)據(jù)。為了清晰地觀察到活體指紋和偽造指紋皮下組織的三維結(jié)構(gòu)并進(jìn)行對比,現(xiàn)用光線投射算法對活體指紋和偽造指紋3 組實驗數(shù)據(jù)進(jìn)行三維可視化,其中,1 組活體指紋數(shù)據(jù)來自于活體手指,2 組偽造指紋數(shù)據(jù)來自于帶假手指套的手指。實驗開發(fā)環(huán)境為Visual Studio 2010,算法實現(xiàn)使用GLSL。
構(gòu)建傳輸函數(shù)時,為表皮指紋和皮下組織設(shè)置不同的不透明度,從而區(qū)分出手指的不同部位。圖7(a)是活體指紋體數(shù)據(jù)的三維可視化側(cè)視圖,除了可以看到手指表皮指紋的紋路外,還可以看到手指皮下汗腺的分布;圖7(b)是活體指紋體數(shù)據(jù)的三維可視化仰視圖,可以發(fā)現(xiàn)汗腺被一條條的紋路連接起來,這些紋路實際上是由乳頭層的上邊沿構(gòu)成,形成了手指的內(nèi)指紋,具有跟表皮指紋相同的紋理結(jié)構(gòu),當(dāng)表皮指紋被偽造或損壞時,內(nèi)指紋不會受到影響[17-19]。由此可見:OCT系統(tǒng)可以用于指紋防偽,而且相對傳統(tǒng)的表皮指紋識別具有更高的準(zhǔn)確性和抗干擾性。圖7(c,d)分別為偽造指紋體數(shù)據(jù)的三維可視化側(cè)視圖、仰視圖,假手指套的厚度為2 mm,觀察發(fā)現(xiàn)假手指套表面分布有脊線和谷線,與活體指紋十分相似,而傳統(tǒng)的指紋識別方法就是基于脊線和谷線的細(xì)節(jié)特征來進(jìn)行識別的,如果不法分子秘密獲取他人的指紋,用硅膠等廉價的材料制造偽造指紋,進(jìn)行違法犯罪活動,就會給社會帶來極大的安全隱患,進(jìn)一步觀察發(fā)現(xiàn),偽造指紋皮下沒有汗腺,再往下又隱約可以看到活體指紋,說明OCT系統(tǒng)掃描到了假手指套和一部分活體手指。圖7(e,f)分別為偽造指紋體數(shù)據(jù)的三維可視化側(cè)視圖、仰視圖,假手指套的厚度為3 mm,由于筆者所用OCT系統(tǒng)掃描深度大約為3 mm,此時掃描不到活體手指,所以偽造指紋下面只能看到一些硅膠氣泡,觀察不到活體指紋。綜上所述,OCT技術(shù)在指紋特征采集和指紋防偽領(lǐng)域具有廣闊的應(yīng)用前景。
圖7 活體指紋和偽造指紋皮下組織的三維可視化對比圖Fig.7 3D visualization of live and forged fingerprints subcutaneous tissue
對頻域OCT系統(tǒng)進(jìn)行了介紹,闡釋了該系統(tǒng)的數(shù)據(jù)處理過程,并對數(shù)據(jù)處理過程進(jìn)行了GPU加速。最后結(jié)合光線投射算法,對活體指紋和偽造指紋皮下組織進(jìn)行三維可視化,比較發(fā)現(xiàn):活體指紋皮下分布有豐富的汗腺,而偽造指紋皮下沒有組織結(jié)構(gòu),該實驗結(jié)果證明了OCT技術(shù)對于指紋防偽,打擊指紋犯罪具有重要實用意義。在對指紋皮下組織進(jìn)行三維可視化時,發(fā)現(xiàn)皮下乳頭層形成的內(nèi)指紋與表皮指紋的圖案一致,如何利用相關(guān)技術(shù)提取內(nèi)指紋,并進(jìn)行指紋識別將是下一步工作的重點。