基于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技术对于指纹防伪，打击指纹犯罪具有重要实用意义。在对指纹皮下组织进行三维可视化时，发现皮下乳头层形成的内指纹与表皮指纹的图案一致，如何利用相关技术提取内指纹，并进行指纹识别将是下一步工作的重点。