基于OCT的手指皮下組織的三維可視化

(浙江工業大學 信息工程學院，浙江 杭州 310023)

光學相干層析成像技術[1-3](Optical coherence tomograph,OCT)是一種新型的光學成像技術，它利用弱相干光干涉儀的基本原理，檢測生物組織不同深度層面對入射弱相干光的背向散射信號，能夠得到生物組織的二維或三維結構圖像。它采用紅外光掃描，穿透不深,在一般的生物組織中成像深度只有2～3 mm,但它的辨識能力比CT高了1 000 倍，分辨率可達到1～15 μm，被應用于眼、皮膚和牙齒等的高分辨率成像，可有力輔助醫生進行臨床診斷和治療。近年來，OCT被逐漸應用到手指生物特征采集[4]上，傳統的手指特征采集主要通過按壓的方式直接獲取二維的指紋圖像，然后根據指紋圖像的全局特征和細節特征來對人的身份進行鑒定，而OCT則借助弱相干光，從空氣入射到手指表皮再深入到皮下組織，除了采集表皮指紋，還采集了皮下汗腺、真皮層內指紋這些新的生物特征信息。

三維可視化[5-6]是促進對OCT采集的數據進行分析和研究的有效手段，如吳凌[7]用OCT獲得指甲邊緣的三維圖像和橘子表皮的三維圖像；萬木森[8]用OCT獲得離體牙釉質和牙本質的三維圖像；Kwon等[9]用OCT獲得了血流的三維圖像。但隨著OCT采樣速率的逐步提高，數據采集量龐大，且數據處理過程包括許多次FFT(Fast Fourier transform)變換，該算法計算復雜度高，耗時嚴重，難以達到實時成像的要求，針對這一問題，Watanabe等[10]最先提出將GPU與硬件相結合，應用于OCT系統成像，提高了數據處理速度；Jeught等[11]在不改變外界硬件的基礎上實現了利用GPU完成重采樣過程的加速。筆者利用頻域OCT系統采集得到手指樣品的原始數據，借助GPU的高速并行運算能力對OCT數據處理過程進行加速，解決了CPU串行處理方式速度慢造成無法實時成像的問題。為了能夠更加直觀準確地展示手指皮下組織的空間分布和幾何形狀，筆者利用三維可視化算法中的光線投射算法來展示手指皮下組織的三維結構，成像結果能很好地反映手指皮下組織的形狀特征和空間分布。最后分別掃描活體手指和帶假手指套的手指得到原始數據，經數據處理后進行三維可視化，比較發現：它們的成像結果都可以顯示表皮指紋，但手指皮下組織結構卻存在著明顯的差異，實驗結果證實了OCT技術在指紋特征采集和指紋防偽領域[12-13]具有廣闊的應用前景。

1 OCT系統搭建

頻域OCT系統主要由光源、樣品臂、參考臂和光譜儀等幾大部分組成。所用系統如圖1所示，寬帶光源SLD發出的低相干光(中心波長λ為848 nm，帶寬Δλ為46 nm)經4 端光纖耦合器分成2 束光，以50/50的比例分別進入參考臂和樣品臂。進入參考臂中的光經過準直透鏡準直后，平行地入射到全反射鏡上，并被平行地反射回來作為參考光。進入樣品臂中的光聚焦在樣品上，通過掃描振鏡系統對手指樣品進行掃描，被手指樣品中不同深度的散射顆粒反射成為信號光。從參考臂和樣品臂中返回的參考光和信號光在光纖耦合器中匯合，發生疊加和干涉，干涉信號從耦合器另一端出射進入光譜儀進行解析，被線陣CCD采集，通過數據采集卡送入計算機中進行一系列的數據處理，然后將處理后的數據進行2D或3D成像。

圖1 OCT系統原理結構示意圖Fig.1 OCT system principle structure diagram

頻域OCT系統主要通過掃描振鏡系統完成對手指樣品的橫向掃描，首先對手指樣品上的一個點進行一次軸向掃描(A-SCAN)，CCD將采集到一條包含深度信息的線性光譜，由于手指皮下組織對光的散射和吸收，本系統的成像深度大約為3 mm。然后對手指樣品進行橫向掃描(B-SCAN)，獲得多條A-SCAN組成的深度信息光譜，可以通過數據處理重構出一幅二維圖像，如果依次往后重復地進行橫向掃描可以得到一個體掃描數據，這個體數據包含了被測手指樣品的三維結構信息。筆者所用的頻域OCT系統的軸向分辨率是7 μm，橫向分辨率是11 μm，CCD的像素點為2 048 個，則一個A-SCAN獲得的像素個數為2 048，一個B-SCAN包含500 個A-SCAN信號，即對手指樣品某一方向進行了500 次深度數據的獲取，掃描速率是36 klines/s，每個B-SCAN大小為500 lines×2 048 pixels/line×2 bytes/pixel，共計1.9 MB數據，300 個B-SCAN可以組成一個大小為570 MB的體數據，數據量很大，這將造成數據處理過程耗時嚴重。

2 基于GPU并行的數據處理

2.1 OCT數據處理

本系統的數據處理過程如圖2所示,除了對500 條光譜求平均，數據處理是一條線一條線進行的，每條線包括數據截取、插值、減直流、FFT變換、取模和取對數這些操作。數據截取是指一個A-SCAN有2 048 個像素，只截取400～1 850的1 451 個像素；由于光譜儀測量到的均勻干涉信號光譜I(λ)經空間轉換后變成了非均勻的波數空間函數I(k),而FFT變換要求輸入數據是均勻的，需要通過插值變換將1 451 個數據變成4 096 個均勻采樣的數據；減直流是為了去除參考臂和樣品臂返回光的直流信號以及樣品不同反射層返回光之間的干涉信號，它們屬于需要去除的背景光；FFT變換可以解析出數據的深度信息，是數據處理過程中最重要的一個步驟；FFT變換之后，對前500 個數據取模得到光功率譜的幅值，由于CCD采集的數據會存在圖像灰度值偏低的問題，取對數可以對圖像進行銳化。

圖2 OCT數據處理流程圖Fig.2 OCT data processing flow chart

OCT系統的數據采集量大，數據處理過程中的計算量也很大，一個體掃描數據(500×2 048×300)需要執行1 500 次數據截取，1 501 次插值(500 條光譜求平均后也需要插值)，1 500 次減直流，1 500 次FFT變換，1 500 次取模和1 500 次取對數。其中，一次數據截取要對1 451 個數據進行計算，一次插值、減直流和FFT要分別對4 096 個數據進行計算，一次取模和取對數要分別對500 個數據進行計算。所有的計算中，FFT算法復雜度最高，耗時最嚴重，如果采樣傳統的CPU串行處理方式，將難以滿足實時成像的需求。鑒于一個B-SCAN中的500 個A-SCAN相互獨立，且在進行數據截取、插值、減直流、FFT變換、取模和取對數這些計算時，大量單個的數據元可以獨立完成計算，因此適合在GPU中并行處理。

2.2 GPU數據處理

CUDA(Computer unified device architecture)是NVIDIA在2006 年為GPU上的并行編程提供的編程API和架構。在CUDA模型中，CPU作為主機(Host),負責內存和顯存的分配、數據傳輸以及內核的啟動，GPU作為協處理器或者設備(Device)，負責并行部分的大量密集型數據的計算。

GPU中需要完成的主要有數據截取、插值、減直流、FFT變換、取模和取對數這些計算，將它們改寫成Kernel函數，一個Kernel函數對應一個線程網格(Grid)，線程網格由大量線程塊(Block)組成,線程塊由大量線程(Thread)構建。具體的線程布局如表1所示，Grid(x,y,1)用于定義Grid的維度和尺寸，即一個Grid有多少個Block，Block(x,y,z)用于定義Thread的維度和尺寸，即一個Block有多少個Thread。在線程布局時，定義的Grid和Block都是一維的，線程中的線程數量都是一個束(Warp)的整數倍，即32的整數倍，因為GPU的執行是以Warp為單位進行調度，如果Block中的線程數不設為32的整數倍，則最后一個Warp中的部分線程是沒有用的。線程布局完成后，每個數據元都有唯一對應線程和線程塊索引。其中，FFT運算部分采用CUDA自帶的CUFFT庫實現，首先用CufftHandle plan定義一個句柄plan，然后利用Cufftplan1d函數創建一個一維快速傅立葉變換句柄，設置好plan后，運行CufftExexcC2C函數，即可在GPU上完成FFT運算，由于plan可以根據輸入數據的大小預先配置好內存和計算資源，在運算時處理器可達到最佳性能。具體程序設計如下：1) CPU上讀取OCT采集到的原始bin數據到內存；2) 利用CUDA中的cudaMalloc函數在GPU上分配顯存，由于筆者將300 個B-SCAN依次在GPU中并行處理，故分配的顯存大小由一個B-SCAN的數據量決定；3) 利用CUDA中的cudaMemcpy函數將CPU上準備好的數據拷貝至顯存，每次拷貝的數據量大小為一個B-SCAN的數據量；4) CPU上調用Kernel函數，GPU上執行并行計算，完成數據截取、插值、減直流、FFT變換、取模和取對數計算，其中，FFT運算部分采用CUDA自帶的CUFFT庫來實現；5) 利用CUDA中的cudaMemcpy函數將GPU顯存中的結果拷貝至CPU，300 個B-SCAN的計算結果存儲在三維紋理數組中，得到一個體數據，為后文的三維可視化做準備。

表1 線程布局Table 1 Thread layout

2.3 數據處理時間測試結果

以Microsoft Visual Studio 2010集成CUDA 9.1為開發環境，GPU為NVIDIA GeForce GTX 1050 Ti,計算能力6.1。現將筆者提出的GPU并行處理方式與Matlab處理方式和CPU串行處理方式進行對比。處理一個體掃描數據(由300 個B-SCAN數據構成，共570 M)的時間測試結果如表2所示。其中Matlab處理方式耗時最長，總共需要45 000 ms。CPU串行處理方式中用到了FFTW庫，FFTW是由麻省理工學院超級計算技術組開發的一套離散傅里葉變換(DFT)的計算庫，被稱為是最快的FFT，由于數據處理過程中要執行很多次FFT變換，耗時嚴重，運用FFTW庫以后可大大縮減計算時間，整個處理過程總共需要7 700 ms。GPU并行處理方式運用了GPU強大的并行計算能力，FFT運算部分采用CUDA自帶的CUFFT庫來實現，它由FFTW庫發展而來，同樣具有很快的計算速率，除此之外，數據截取、插值、減直流、取模和取對數等操作都改寫成Kernel函數，在GPU上并行實現，進一步提升了處理速度，最后總共耗時2 476 ms。比較發現：筆者提出的方法耗時最短，速度是Matlab處理方式的18 倍，是CPU串行處理方式(運用FFTW庫)的3 倍。OCT數據的實時三維成像要求OCT系統剛采集完一個體掃描數據，計算機屏幕上就能顯示樣品的三維圖像，由于所用OCT系統的掃描速率是36 klines/s,采集完一個體掃描數據(300 個B-SCAN，共150 000 個A-SCAN)大約需要4 167 ms，如果一邊采集數據一邊對每個B-SCAN進行數據處理，數據處理的時間將會被隱藏，這種情況下，Matlab處理和CPU串行處理時間都大于采集時間，而且所用系統后期還會改進，采集速率會有進一步的提升，因此不滿足要求，而GPU并行處理時間小于采集時間，能滿足實時三維成像的要求。

表2 不同方式下數據處理時間對比Table 2 Comparison of data processing time in different ways

3 OCT數據的三維可視化

處理后的OCT數據被保存到三維紋理數組中，一個三維紋理數組由300 個二維數組構成，它們是手指的二維斷層圖像序列，圖3是其中的一幅二維斷層圖像，從淺到深依次是表皮層、真皮層。從表皮層中可以看到角質層和汗腺，角質層的圖像灰度更高，汗腺與手指表面的汗孔相連，真皮層最厚，在它與表皮層的分界處呈波浪狀的是乳頭層。二維圖像只能反映手指的截面信息，而三維圖像可以反映手指皮下組織的三維空間結構，較之二維圖像更加直觀準確，現利用移動立方體算法、最大密度投影算法和光線投射算法對手指皮下組織進行三維可視化，并比較渲染效果的差異。

圖3 手指二維斷層圖像Fig.3 Finger 2D tomographic image

現采用MC算法、MIP算法和Ray-casting算法分別對手指皮下組織進行三維可視化，圖4展示了不同算法下任意視角的手指皮下組織三維結構，可以發現：手指皮下組織主要由表皮指紋、汗腺和真皮層3 部分構成。

圖4 不同算法下手指皮下組織可視化效果的比較Fig.4 Visualization of subcutaneous tissues by different algorithms

汗腺是手指皮下組織中非常典型的組織結構，筆者主要依據汗腺來比較可視化效果。汗腺的繪制結果如圖5(a～c)所示，與MC算法和MIP算法相比，Ray-casting算法繪制的圖像具有明顯的層次性，更清晰地反映了皮下汗腺的空間分布與結構特征。

圖5 不同算法下皮下汗腺可視化效果的比較Fig.5 Visualization of subcutaneous sweat glands under different algorithms

MC算法在求等值面的過程中，將等值面與立方體棱邊的交點簡單相連得到三角面片，最后將三角面片連接得到逼近的等值面，但這個等值面只是一個近似面，有時會不準確，如果數據場的數據量很少且間隙很大，最終繪制的圖像誤差將會非常大。但筆者的數據來自于OCT系統，分辨率很高，在這種情況下，立方體中的三角面片非常微小，最終提取的等值面不會有太大誤差，所以圖5(a)所示的皮下汗腺的形狀及分布基本與實際符合，但由于MC算法通常只適用于表面特征分明的器官和組織，而汗腺的形態特征并不分明，所以最后繪制的圖像不能清晰地反映汗腺的細節特征。

MIP算法實現較簡單，繪制速度優于另外兩種算法，但這也導致繪制效果較差。圖5(b)中的汗腺疊加在一起，幾乎觀察不到汗腺的細節特征，看起來更像是一張X光片。由于血管造影對細節要求不高，所以該算法經常應用于三維血管結構的顯示。

從圖5(c)可以看出：Ray-casting算法繪制的圖像清晰地展示了皮下汗腺的形狀特征和層次分布，具有非常好的形狀感知和深度感知,這是因為Ray-casting算法考慮了所有體素對結果圖像的貢獻，且OCT數據具有很高的分辨率，特征信息豐富，所以繪制的圖像充分表達了體數據中原有的信息，圖像質量比MC算法和MIP算法高很多。

現將一些細節汗腺放大作進一步的比較分析。圖6(a，b)為MC算法下繪制的細節汗腺，圖6(c，d)分別為MIP算法、Ray-casting算法下繪制的細節汗腺。從圖6(a)中可以發現有的汗腺被斬斷了，即不具有連通性，這是由于MC算法存在面二義性問題，當值為1的角點和值為0的角點分別位于對角線的兩端時，會有兩種可能的連接方式，這個面就叫做二義性面，如果把它作為兩個立方體的公共面，且兩個立方體在公共面上采用了不同的連接方式，那么最后構造的等值面就會出現“空洞”現象，就像圖6(a)中一樣，汗腺發生了斷裂；圖6(b)中有一根汗腺呈螺旋狀，由于OCT數據的高分辨率，它跟實際的形狀已十分接近，但是與圖6(d)中的螺旋汗腺相比，它顯得毫無立體感，這是由MC算法本身的缺陷造成的，在獲取等值面的過程中，將等值面與立方體棱邊的交點簡單連接成三角面片，這樣必然會丟失掉一些細節信息，比如立方體內的環形結構和臨界點(等值面改變的點)，就像圖6(b)中的螺旋狀汗腺一樣，螺旋處信息單薄，導致形狀不夠逼真；圖6(c)中的汗腺缺乏形狀感知和深度感知，可以發現能看到3 根汗腺的基本輪廓但不清晰，而且感受不到汗腺的前后順序，反而覺得它們處于同一平面，這是因為MIP算法只投射光線上的最大密度值點，這樣會導致圖像中缺少深度信息，甚至一些位于后面的組織由于密度大于前面的組織，顯示時出現后面的組織遮擋前面組織的情況，給用戶造成混淆。

圖6 汗腺細節放大圖Fig.6 Sweat gland detail larger image

盡管Ray-casting算法繪圖質量很高，但由于體素數量和投射光線數量過多，導致Ray-casting算法繪制速度緩慢，但隨著國內外研究者提出了許多關于Ray-casting算法的加速技術[14-16]，這個問題漸漸得到解決，基本可以實現實時交互操作。綜上所述，最終選擇Ray-casting算法來實現手指皮下組織的三維可視化。

4 三維可視化的應用與討論

利用頻域OCT系統對活體手指和帶假手指套(2，3 mm厚)的手指分別進行掃描，掃描面積為5 mm×3 mm，得到3 組實驗數據。為了清晰地觀察到活體指紋和偽造指紋皮下組織的三維結構并進行對比，現用光線投射算法對活體指紋和偽造指紋3 組實驗數據進行三維可視化，其中，1 組活體指紋數據來自于活體手指，2 組偽造指紋數據來自于帶假手指套的手指。實驗開發環境為Visual Studio 2010，算法實現使用GLSL。

構建傳輸函數時，為表皮指紋和皮下組織設置不同的不透明度，從而區分出手指的不同部位。圖7(a)是活體指紋體數據的三維可視化側視圖，除了可以看到手指表皮指紋的紋路外，還可以看到手指皮下汗腺的分布；圖7(b)是活體指紋體數據的三維可視化仰視圖，可以發現汗腺被一條條的紋路連接起來，這些紋路實際上是由乳頭層的上邊沿構成，形成了手指的內指紋，具有跟表皮指紋相同的紋理結構，當表皮指紋被偽造或損壞時，內指紋不會受到影響[17-19]。由此可見：OCT系統可以用于指紋防偽，而且相對傳統的表皮指紋識別具有更高的準確性和抗干擾性。圖7(c，d)分別為偽造指紋體數據的三維可視化側視圖、仰視圖，假手指套的厚度為2 mm，觀察發現假手指套表面分布有脊線和谷線，與活體指紋十分相似，而傳統的指紋識別方法就是基于脊線和谷線的細節特征來進行識別的，如果不法分子秘密獲取他人的指紋，用硅膠等廉價的材料制造偽造指紋，進行違法犯罪活動，就會給社會帶來極大的安全隱患，進一步觀察發現，偽造指紋皮下沒有汗腺，再往下又隱約可以看到活體指紋，說明OCT系統掃描到了假手指套和一部分活體手指。圖7(e，f)分別為偽造指紋體數據的三維可視化側視圖、仰視圖，假手指套的厚度為3 mm，由于筆者所用OCT系統掃描深度大約為3 mm，此時掃描不到活體手指，所以偽造指紋下面只能看到一些硅膠氣泡，觀察不到活體指紋。綜上所述，OCT技術在指紋特征采集和指紋防偽領域具有廣闊的應用前景。

圖7 活體指紋和偽造指紋皮下組織的三維可視化對比圖Fig.7 3D visualization of live and forged fingerprints subcutaneous tissue

5 結 論

對頻域OCT系統進行了介紹，闡釋了該系統的數據處理過程，并對數據處理過程進行了GPU加速。最后結合光線投射算法，對活體指紋和偽造指紋皮下組織進行三維可視化，比較發現：活體指紋皮下分布有豐富的汗腺，而偽造指紋皮下沒有組織結構，該實驗結果證明了OCT技術對于指紋防偽，打擊指紋犯罪具有重要實用意義。在對指紋皮下組織進行三維可視化時，發現皮下乳頭層形成的內指紋與表皮指紋的圖案一致，如何利用相關技術提取內指紋，并進行指紋識別將是下一步工作的重點。