高速光学相干层析三维扫描系统成像软件的设计实现

单位代码 10006 学 号 10101004 分 类 号 R318 毕业设计(论文)高速光学相干层析三维扫描系统成像软件的设计实现 学院名称生物与医学工程学院 专业名称光学成像系统 学生姓名吴玮莹 指导教师李 昂 2014年 5月 论文封面书脊高速光学相干层析三维扫描系统成像软件的设计实现 吴玮莹 北京航空航天大学 四号黑体字 四号黑体字 小四号黑体字北京航空航天大学本科毕业设计（论文）任务书、毕业设计（论文）题目：高速光学相干层析三维扫描系统成像软件的设计实现 、毕业设计（论文）使用的原始资料（数据）及设计技术要求： 使用光学相干层析技术（OCT）对葡萄三维结构进行研究，相关设备平台包括由近红外光谱仪，宽频光源，系统控制模块，自动对焦系统，系统控制界面搭建而成的OCT系统，OCT技术探测的精度能够较好地达到其内部结构的探测（OCT系统三维扫描速度预计为25frame/sec，横向分辨率可达1520um，深度分辨率可达10um）。 对于葡萄三维结构的探测所需的成像系统需要极高的软硬件要求，其中，本课题主要针对软件系统进行设计。对于软件设计的要求主要包括：1.基于多核CPU的并行OCT解调引擎实现：对葡萄三维结构进行并行光学相干断层扫描，解调并进行成像，包括基于Visual Studio开发环境下Intel MKL函数库的使用进行波数域线性化及FFT等信号处理内容；2.葡萄三维结构图像的三维重建：综合使用栅格化、光线投射、辐射着色、光线跟踪等算法，完成由模型到图像的转化,调用VTK、ImageJ 3D Viewer作为接口进行交互，实现图像三维重建；3.数据留盘实现：对葡萄解调图像所得进行存储，使用多线程高级编程对所获数据进行处理；4.用户友好化软件交互界面设计。 、毕业设计（论文）工作内容：20143下半月：选取合适的OCT系统波数域线性优化算法，对OCT光学成像系统软件进行结构分析，完成软件工程需求分析 20144上半月：基于Intel MKL函数库对波数域数据进行处理，完成A扫后由波长域信息向深度信息的算法转换，与VTK、ImageJ 3D Viewer进行交互，使用多线程接口同时对数据的输入输出流进行处理 20144 下半月：对软件系统进行测试和优化，设计较为友好的用户界面，完善软件的相关功能； 20145 上半月：对软件系统进行测试和完善，协调配合OCT硬件系统对软件设计进行调试和改善；对三维结构样品过程进行样品实验成像； 20145下半月：对样品实验结果进行验收，撰写毕业设计论文，整理参考文献和实验结果； 20146上半月：修改论文，答辩。 、主要参考资料：1. Huang D, Swanson E A, Lin C P, et al. Optical coherence tomographyJ. Science, 1991, 254(5035): 1178-1181.2. Yaqoob Z, Wu J, Yang C. Spectral domain optical coherence tomography: a better OCT imaging strategyJ. Biotechniques, 2005. 北京航空航天大学毕业设计(论文)第 VIII 页 3.Yazdanfar S, Kulkarni M, Izatt J. High resolution imaging of in vivo cardiac dynamics using color Doppler optical coherence tomographyJ. Optics Express, 1997, 1(13): 424-431.4. Yazdanfar S, Rollins A M, Izatt J A. Imaging and velocimetry of the human retinal circulation with color Doppler optical coherence tomographyJ. Optics Letters, 2000, 25(19): 1448-1450.5. Puliafito C A, Hee M R, Lin C P, et al. Imaging of macular diseases with optical coherence tomographyJ. Ophthalmology, 1995, 102(2): 217-229.6. Rollins A M, Ung-Arunyawee R, Chak A, et al. Real-time in vivo imaging of human gastrointestinal ultrastructure by use of endoscopic optical coherence tomography with a novel efficient interferometer designJ. Optics letters, 1999, 24(19): 1358-1360.7. Sivak Jr M V, Kobayashi K, Izatt J A, et al. High-resolution endoscopic imaging of the GI tract using optical coherence tomographyJ. Gastrointestinal endoscopy, 2000, 51(4): 474-479.8. Tearney G J, Bouma B E. Atherosclerotic plaque characterization by spatial and temporal speckle pattern analysisJ. Optics letters, 2002, 27(7): 533-535.9. Fujimoto J G. Optical coherence tomography for ultrahigh resolution in vivo imagingJ. Nature biotechnology, 2003, 21(11): 1361-1367.10. Bouma B E, Tearney G J, Yabushita H, et al. Evaluation of intracoronary stenting by intravascular optical coherence tomographyJ. Heart, 2003, 89(3): 317-320.11. Wang X J, Milner T E, de Boer J F, et al. Characterization of dentin and enamel by use of optical coherence tomographyJ. Applied Optics, 1999, 38(10): 2092-2096.12. Brandenburg R, Haller B, Hauger C. Real-time in vivo imaging of dental tissue by means of optical coherence tomography (OCT)J. Optics communications, 2003, 227(4): 203-211.13. Yelbuz T M, Choma M A, Thrane L, et al. Optical coherence tomography a new high-resolution imaging technology to study cardiac development in chick embryosJ. Circulation, 2002, 106(22): 2771-2774.14. Perea-Gomez A, Camus A, Moreau A, et al. Initiation of gastrulation in the mouse embryo is preceded by an apparent shift in the orientation of the anterior-posterior axisJ. Current biology, 2004, 14(3): 197-207.15. Boppart S A, Bouma B E, Brezinski M E, et al. Imaging developing neural morphology using optical coherence tomographyJ. Journal of Neuroscience Methods, 1996, 70(1): 65-72. 生物与医学工程学 院（系）生物医学工程 专业类 101012 班学生 吴玮莹 毕业设计（论文）时间： 自 年 月 日至 年 月 日答辩时间： 年 月 日 成绩 指导教师： 兼职教师或答疑教师（并指出所负责部分）： 教研室主任 注：任务书应该附在已完成的毕业设计（论文）的首页。本人声明我声明，本论文及其研究工作是由本人在导师指导下独立完成的，在完成论文时所利用的一切资料均已在参考文献中列出。 作者： 吴玮莹签字：时间：2014年 5 月高速光学相干层析扫描系统成像软件的设计实现 学 生：吴玮莹 指导教师：李 昂摘 要 本论文描述了一种高速光学相干层析三维成像软件系统的实现方法，该成像软件针对谱域光学相干层析术（SD-OCT）的数据类型。其数据处理过程主要包括两个部分，分别为数据解调过程以及三维重建过程。其中，数据解调过程包括数据重排、快速傅里叶变换以及LOG灰度级压缩与数据存储，其主要开发工具为Intel MKL数学核心函数库以及Visual Studio 2012运行环境。解调图像获得了较高的信噪比以及较快的成像速率，并对数据重排过程中不同的插值算法进行了对比，对LOG灰度级压缩的相关参数进行了优化；三维重建过程主要包括基于VTK视觉化工具函式库的体绘制光线投射法的使用以及基于Imagej 3D Viewer相关三维重建接口的调用，并对二者的三维重建效果进行了对比与分析。 关键词：光学相干层析成像，解调，三维重建 High-speed optical coherence tomography imaging software system designAuthor : Wu Weiying Tutor : LI AngAbstract This paper describes a data processing method for high-speed optical coherence tomography imaging software system, which aims at Spectral Domain Optical Coherence Tomography (SD-OCT). The data processing procedure can be divided into two parts, the data demodulation process and the 3D reconstruction process. The data demodulation process involves k-domain remapping of the original spectralscopic data, Fast Fourier Transform (FFT) of the remapped data, gray-scale manipulation of the FFT result and restoration of the demodulation result, using the Intel Math Kernel Library and Visual Studio 2012 as the development tools. The demodulation results are with relatively high signal noise ratio and imaging speed, with different interpolation algorithms use in the remapping process and optimization for the relative parameters in the LOG gray-scale manipulation. The 3D reconstruction process involves the use of Volume Rendering Ray Casting Algorithms in the Visualization Toolkit (VTK) and the calling for the interface card of the 3D reconstruction based on the Imagej 3D Viewer in JAVA. Also, the comparison between the two 3D reconstruction methods are discussed.Key words：Optical coherence tomography, demodulation,3D reconstruction目录1 绪论1 1.1课题研究背景及现实意义1 1.2研究内容42光学相干成像基本原理5 2.1OCT系统基本理论5 2.2OCT相干成像72.2.1 时域光学相干成像72.2.2 频域光学相干成像73 OCT成像算法原理10 3.1 OCT成像数据解调103.1.1 直流分量滤除103.1.2 数据重排113.1.3 傅里叶变换133.1.4 LOG灰度压缩133.1.5 光谱仪校正与色散补偿14 3.2 OCT图像三维重建153.2.1 OCT三维数据场153.2.3 OCT三维重建面绘制算法173.2.3 OCT三维重建体绘制算法17 3.3 OCT图像存储184 OCT成像软件开发19 4.1软件系统开发工具194.1.1 Intel MKL数学核心函数库194.1.2 三维重建软件平台204.1.3 VTK（Visualization Toolkit）视觉化工具函式库204.1.4 ImageJ 3D Viewer22 4.2 软件系统构成235 OCT成像算法验证25 5.1数据来源25 5.2 二维解调算法结果265.2.1不同插值算法数据重排结果比较265.2.2 二维解调OCT伪彩色成像结果275.2.3 色散补偿前后结果比较285.2.4 灰度级压缩参数优化295.2.4 解调效率评估30 5.3 三维重建结果315.3.1 VTK三维重建结果315.3.2 Imagej 3D Viewer三维重建结果32结论35致谢36参考文献38附录41 核心代码构成41 附加代码构成89北京航空航天大学毕业设计(论文)第 123 页1 绪论1.1 课题研究背景及现实意义光学技术在生物医学工程领域的应用具有很长的历史，自从18世纪显微镜问世后，光学技术就成为生物学领域不可或缺的工具。1960年激光问世之后，一种新型的手术工具成为医疗领域的新宠，而光纤技术的发展更是使得允许直接对人体内部器官进行深入观察的内窥镜操作成为可能。在现代化的医疗实验室，光学技术的使用能够较好地进行组织样本的化学分析以及血细胞数量和大小的测量。OCT光学相干层析术10年前由MIT的Fujimoto的团队提出以来1-2，不断地发展成熟，目前已成为重要的临床成像方法。OCT光学相干层析能够在微米级分辨率对图像进行补偿，而其成像所要求的光学流动水平较低，灵活性较高，能够使用于敏感组织区域，如人眼等部位的成像；而其扫描探测器体积较小，能够在医疗领域人体胃肠道以及主要的血管里面进行操作。对于一个光学相干层析系统，其基本功能能够保证对目标一定深度的结构图像进行成像，而更多复杂的光学相干断层扫描则能够提供附加的功能信息，如流动（通过多普勒光学相干层析扫描），组织结构布局（通过双折射光学相干层析扫描），特定对照的空间分布（通过分子对照光学相干层析扫描）等的成像3-4。OCT光学相干层析术是一种能够提供高清活体生物内部微结构的横截面的典型成像技术，该技术深深根植于早期对白光相干性的研究所催生的光相干域反射仪的发展，这是一种一维的光学测量技术。尽管OCDR光相干域反射仪原先主要应用于光学电缆和网络元件查错，但它在探索眼部和其他生物组织中的潜质渐渐地浮出水面。例如，光学相干层析术在黄斑裂孔疾病发展过程中能够提供较为清晰的深度成像轮廓，对视网膜稀释与分离具有良好的成像效果，而这种成像效果对于那种表层的图像成像方法则是难以企及的。光学相干层析术在消化道应用中亦具有极高的实用性，其结构信息提供了一种对疾病诊断较为有效的生理相关信息的间接测量方法。OCT光学层析术在人体其他部位的应用包括对肿瘤诊断的内窥镜GI神经纤维束的成像6-7，对血管系统中脂质成形的血小板进行监督和评估8-9，监督经皮冠状动脉介入治疗10，对与牙齿退化有关的结构伤害进行追踪探测11-12等。OCT光学相干层析术在生命科学中的应用则朝鸡胚心脏发育13，鼠胚胎原肠胚形成14，以及有爪蟾蜍的神经形态学成像方向发展15。OCT光学相干层析术的光学分辨能力是通过实现光源宽带短暂的相干性，使得OCT扫描仪能够获得不同深度的组织结构微观图像，远远地超过了传统光学领域以及共焦望远镜所能达到的精度。对透明组织要求探查深度需要大于2cm，包括眼睛和青蛙胚胎。对于皮肤和其他高度分散的组织，OCT光学相干断层扫描技术可以对微小的血管组织和其他表皮层以下1-2mm深度的结构进行成像。OCT具有可进行高频超声成像的优点，该技术具有极强的竞争力，能够实现较大的深度探测，而其硬件系统的要求则相对简单。OCT的深度成像范围可达6mm，而其一个A扫的轴向分辨率可达5um，与超声波、共焦显微等技术相比，它在较高的分辨率以及较大的成像深度范围之间达到了较好的权衡，与其他成像技术相比，OCT的性能如图1.1所示。 图1.1 成像性能比较随着OCT光学相干层析术的发展，其发展趋势逐渐由原来的时域光学相干层析术向谱域光学相干层析术转变。谱域相干层析术(SD-OCT)相对于时域相干层析术(TD-OCT)能够提供更高的感光度，而高感光度则能够提高OCT光学相干层析术的获取样本信息的速率，获得更高的深度洞穿率，从而提高不同的功能光学相干层析扫描方法的敏感度。谱域相干层析术(SD-OCT)所使用的线性感光阵列具有极高的扫描速度，其采集数据的速度可达几十兆赫兹。对于如此高的数据采集速度，需要较高的数据处理速度来与之匹配，这就要求OCT光学相干层析术的软件系统具有较高的性能。对于一个OCT光学相干层析三维成像软件系统，其主要任务为快速高效地完成对硬件接口裸数据的解调以及对二维图像的三维显示。其中，对于OCT裸数据的解调，Intel MKL(英特尔数学核心函数库)封装了大量高度优化的函数例程，能够快速高效地完成OCT解调所需的样条插值、傅里叶变换等计算过程，能够较好地满足SD-OCT的成像要求，其处理数据的性能（以FFT为例）与传统算法相比如图1.2所示。 图1.2 Intel MKL与传统算法FFT性能比较 对于OCT光学相干层析术二维解调图像的三维显示，VTK（Visualization Toolkit）提供了一个开源、多语言、并行处理的图形处理开发函数库，可对输入的二维图像进行快速高效的三维显示及重建,ImageJ 3D Viewer提供了基于JAVA的三维重建接口，便于高效精确地对其进行调用。基于Intel MKL以及VTK的三维图像显示、ImageJ 3D Viewer接口的调用与基于MFC的API相互嵌合，可满足光学相干层析实时精确的成像要求。1.2 研究内容本论文第二章对OCT光学相干成像基本原理进行了简要的阐述和介绍，包括OCT系统的基本理论以及时域光学相干层析术（TD-OCT）与频域光学相干层析术(SD-OCT)的相关原理；第三章则对OCT成像算法原理进行了阐述，包括直流分量滤除的相关原理、数据重排所使用的三次样条插值算法原理、快速傅里叶变换、LOG灰度级压缩以及光谱仪校正和色散补偿的相关原理、OCT三维数据场基本原理以及三维重建面绘制、体绘制相关算法基本原理等；第四章则对本软件所使用的相关函数库以及本软件系统的构成进行了介绍，包括Intel MKL数学核心函数库、VTK视觉化工具函式库、ImageJ 3D Viewer插件相关函数的使用以及三维成像软件的构成与流程等；第五章则对所使用的算法进行了实验验证，通过对葡萄OCT数据使用上述方法进行二维解调、三维重建，论证了本成像软件的可行性以及实用性。 本论文的研究内容包括：OCT数据解调算法的编写、不同插值算法的比较、基于CPU的OCT解调算法的优化、灰度级压缩参数的优化、VTK三维重建体绘制算法的优化、ImageJ 3D Viewer 三维显示JAVA接口的调用、VTK三维重建体绘制算法的优化与ImageJ 3D Viewer三维显示效果的对比。2 光学相干成像基本原理2.1 OCT系统基本理论图2.1展示了基本OCT光学相干成像系统的基本组成。该系统的中心元件为一个光学干涉仪，其照明由宽带光源提供。 图2.1：OCT系统模块 OCT光学相干断层扫描仪将一个宽带场分成参考场Er和样本场Es，其中样本场通过光学扫描和透镜对组织表面以下一定深度的点进行聚焦。从组织回扫之后，标准化的样本场Es 与Er在图像探测器的表面混合。假设图像探测器捕捉了从参考臂与样本臂的光源发射出来的所有的光，则光学探测器的光强为 Id=0.5Ir+Is+Re (2.1)其中Ir和Is分别为从干涉仪的参考臂和样本臂发射来的平均光强，该公式的第二个等式决定于光学的时延，它由参考镜的位置所决定，表示了携带组织结构信息的干涉边缘的振幅。干涉边缘的特征（或者无任何边缘形成）则决定于Es 与Er所匹配的时空特征。因此，干涉仪起到了一个交叉相关器的作用：经过在探测器表面叠加后的相干信号的振幅的聚集提供了对交叉相关信号的测量方法。为了使交叉相关信号从光强的成分中分离出来，须使用不同的技术用来对进行调整。 基于生物组织可表现为一个理想镜对光束不作改变的特征，相关振幅决定于时间相干特征，根据Re=|G|cos20+ (2.2)其中c为光速，0=c/0是光源的中心频率，而G()是它的关于变量()的复合时间相关方程。根据Wiener-Khinchin定理，G()与光源的功率谱密度S()有关,如等式2.3所示。G=0Sexp-j2d (2.3)光源的发射谱的形状和宽度在OCT光学相干断层扫描技术中是较为重要的变量，因为它们对不同光路干涉仪敏感度的影响较大。对于一个OCT成像系统，宽频的光源是可以实现的，从而可产生短暂的时间或者空间的光学干涉。而S()与G()之间的关系都可以清晰地用高斯函数来表示： S()G() (2.4)而对于S()，有S=2ln2/exp-4ln2-02 (2.5)对于G()，有G=exp-2ln22exp(-j20) (2.6)在这些等式中，半功率带宽表示光学频率定义域内光源的频率宽度,相应的相关宽度的测量方法可从等式2.6中衍生而来，其相关长度可由如下等式得到：lc=2cln21 (2.7) 0.4402 (2.8)其中是相干方程中波长最大值测量中的宽度,在OCT光学相干层析术相关文献中，等式2.8是其最为普遍的定义。2.2 OCT相干成像2.2.1 时域光学相干成像一个经典的时间光学相干断层扫描系统如图2.2所示。对于这样一个系统，一个A扫是按照时间的序列进行的，称之为时域光学相干层析成像(TD-OCT)，TDOCT的深度扫描原来是通过在参考臂安装压电陶瓷，压电陶瓷的变化驱动反射镜使参考臂光程发生变化，实现对深度扫描的探测。由于在TD-OCT中，A扫的速度具有局限性，因而早期的OCT速度较为缓慢。图2.2：基于时间以及基于光谱的光学相干断层扫描系统示意图。M为镜子，ODL为光学延迟 在一个时域光学相干层析成像系统中，其参考臂通常为不固定的，可在一定的样本深度范围内对一定距离范围进行扫描，该距离与深度范围相等，表示为xdepth。对于一个给定的扫描持续时间T，我们可以得到一个时域光学相干成像系统将消耗lcxdepthT时间来获取从任何样本表面到达的相干信号。2.2.2 频域光学相干成像 随着OCT成像技术的发展，谱域光学相干层析成像（SD-OCT）逐渐进入人们的视线，它大大地提高了一个A扫的速度，在成像的信噪比和灵敏度上也都有了较大的提升。TD-OCT中，所探测信号是编码在时间域中宽带光源以波数为基准其各个成分经所测样品深度调制后所得的和；而在SD-OCT中，所获信号是按频率域分布的，需通过傅里叶变换被解调出来。对于这样一个系统，其干涉信号在参考臂和样品臂上可分别表示成：Ik=S(k)ARexpik2r+0ASzexp-ik2(r+nz)dz2 （2.9） 谱域光学相干层析扫描系统组成如图2.2(A)所示。该系统使用一个宽带光源，通常为一个超辐射发光二极管或者一个锁模激光器，它与一个2*2的光纤组成了一个Micelson干涉仪。在探测臂中，一个低损耗的光谱仪从参考臂或者样本臂返回光的干涉形式的光谱，而该基于栅格的光谱仪通过一个带有波长的公式对干涉信号进行测量。通常情况下，在k空间，光谱数据将被重新调节和取样，在此之前，则通过傅里叶变换或A扫获得样本深度范围（如图2.2，B-D）。 图2.2：基于光谱仪的光谱光学相干断层扫描系统，宽带激光器能够作为一个超辐射发光二极管或者锁模激光器。SMF为单模纤维，G为衍射栅栏，PC为计算机，CCD为电荷耦合设备。 其中，k为波数，S(k)为光源的功率谱，2r为参考臂往返光程，n为所测样本折射率。其中，为了简化论述，假设AR(参考臂反射对光的反射率)为1；AS(z)则表示样本在不同z值的散射率信息，解调出AS(z)即可实现对样品进行深度断层成像。假设样品散射率如下：ASz=ASz+AS(-z) (2.10)其对样品表面，散射率呈对称分布，带入式(2.14)中，可化为：Ik=Sk1+20ASzcos2knzdz+00ASzASzcos2nkz-zdzdz=Sk1+20ASzcos2knzdz+14-ASzASzcos2nkz-zdzdz (2.11)上式中第一项为信号的直流偏置，第二项为参考臂和样品臂干涉信号，第三项为样本臂信号的自相干；左右两边进行傅里叶逆变换，可得：F-1Ik=F-1S(k)z+12ASz+18ACASz (2.12)其中S(k)对ASz调制，从而表现为所使用光源光谱对系统深度分辨率造成的影响。为了提高信噪比，由(2.12)式可知，须抑制第一项DC噪声和第三项AC噪声。在OCT系统中，AC噪声主要分布样品表面位置，强度较弱，因此较为容易滤除。对于直流噪声，可直接在所测信号中减除直流分量进行去除。由上，我们可得在样品不同深度z的散射率信息ASz，从而用于光学成像。对于SD-OCT技术，其最大优势为不须进行深度扫描，从而大大提高了数据采样的速度，其成像速度较TD-OCT提高了一个数量级。缺点则为探测器阵列价格较高、动态范围较小；另一个缺点则为物体结构自相关宽度为其实际宽度的两倍，从而浪费了一半的信号带宽，对于1024阵列的探测器，只有512个阵列得到有效的利用。3 OCT成像算法原理3.1 OCT成像数据解调3.1.1 直流分量滤除 OCT读出的信号为按波长域线性排列的光谱强度信号，而SD-OCT的深度信息则是由基于波数域线性排列的光谱强度信号的傅氏变换，对于一个SD-OCT系统，对于由样本臂和参考臂得到的裸数据，首先需要进行直流分量滤除、数据重排、快速傅里叶变换、LOG灰度级压缩等步骤，从而获得一个B扫（截面成像）的二维截图图像，其基本原理基于公式3.1： FszFTASK （3.1） 其中，Fsz为所测样本在不同深度z反射光的光强值，ASK为光谱仪所测的对应于波数K的反射光的光强值，两者之间为傅里叶变换对的关系。 对于直流分量滤除，由于光谱仪的输出信号可表示为: Ik=S(k)F-1z+12ASz+18ACASz （3.2） 由上式可得，反射光的直流分量是由于参考臂的反射光作用造成的，所以，在数据解调的过程中，只需在所测光谱仪信号的基础上减去参考臂测得的信号即可。在测量信号之前，遮挡住样本臂，对参考臂的反射光光强信号进行测量，然后进行样本臂光谱强度信号测量即可，无需进行重复测量，OCT输出信号如图4.1所示。 (a)参考臂反射光谱信号 (b)光谱仪输出信号 (c)滤除直流分量之后的交流信号图3.1 OCT输出信号3.1.2 数据重排 在进行直流分量滤除后，须对所得波长域数据进行数据重排，对于数据重排算法，常见的包括近邻插值、线性插值以及三次样条插值等，其中，使用近邻插值以及线性插值无疑会引入较大误差，降低所得图像的信噪比，故而本算法中，使用的是三次样条插值(图3.2)。 三次样条插值使用的插值方法为使用分段多项式对特定点进行拟合，使用多段低阶多项式降低插值误差，故而可以避免使用高阶多项式产生的龙格现象。 三次样条插值主要包括自然插值、Akima插值、Bessel插值等多种方法16-17，其中，这些插值方法均有以下插值形式(公式3.3)：s(x) = w0xfi-1+w1xfi+w2xfi+1+w3xfi+2 （3.3） 由插值条件，我们可得， 图3.2 三次样条插值示意图 sxi=f(xi) （3.4） 而所取样条是相互连接的，故而有 si-1xi=si(xi)(i=1,2n-1) (3.5) 对于若干三次多项式，其连接处两次连续可导，故有： si-1xi=si(xi)(i=1,2n-1) (3.6) 对于给定n个多项式，由于每一个三次多项式都需要4个条件来确定其唯一性，故而共需要4n个条件，公式3.4给出了n+1个条件，而内部数据点则共给出了n-1个条件，共给出4n-2个条件，而另外两个边界条件选择的不同则成为不同三次样条插值方法的诱因。 其中，对于不同的插值方法，其边界条件如表3.1所示： 表3.1 不同插值方法边界条件插值算法 边界条件自然插值si-12xi=si2(xi)Akima插值si=wi+1xi-1,xif+wi-1xi,xi+1fwi+1+wi-1,where wi=|xi,xi+1f-xi-1,xif|Bessel插值si=(xixi-1,xif+xi-1xi,xi+1f)/(xi+xi+1) 3.1.3 傅里叶变换 将数据从波长域重采样至波数域后，对所得结果进行快速傅里叶变换即可得到相应的复数据，对其实部虚部计算结果的幅度值，可得到相应的灰度信息。3.1.4 LOG灰度压缩 经傅里叶变换后所得的二维数组数量级较大，须对其进行灰度级压缩，由于其动态范围较大，直接线性映射到0-255灰度空间将会造成较大的失真，故而在本实验中，采用先对所得傅里叶变换结果进行LOG变换，其映射公式如式3.7。gxi=a1*log(gxi-a2) (3.7)其中，a1与a2与图像本身光强有较大关系，须通过不同取值的灰度直方图来确定其最优取值。另外，为了增强对比度，本实验还对其进行了伪彩色显示，将其映射到RGB空间进行显示。3.1.5 光谱仪校正与色散补偿本实验室所采用的OCT光谱仪为Bayspec OCT engine，其波长域范围为132565nm，在OCT系统中，由于样本臂与参考臂色散的不匹配以及光谱仪与线性感光阵列的对应问题，因此，须对光谱仪进行标定，对读入光谱的波长域进行校正。校正方法为：遮挡系统的参考臂，使其反射的信号为零；在样本臂中置放厚度约为0.1mm的盖玻片，使其上下表面反射的光信号发生干涉。假设两束光所经过的介质基本相同，那么干涉信号的相位谱就应当满足线性： IACk=Acos(2kl) （3.8） 根据以上标准，我们可设定光谱仪分光之后不同像素点所对应的波长关系，当所获得的相位谱最接近线性分布时，认为所设定的对应关系最接近实际情况。假设实际波长与不同像素点满足如下多项式关系： n=Nn-n0+0+c1n-n02+c2n-n03 n=1,2,1024 (3.9)图3.3 光谱仪校正信号示意图上式中，表示所覆盖的波长范围，为了简化论述，设为1300nm；N表示感光阵列像素点数，所使用系统其值为1024；0为中心波长，n0则表示使二次和三次分式误差为0的像素点。式3.9中最后两项各自用以拟合分光系统二次和三次的非线性误差。由于实际实验中，系统三次的非线性误差相对于其他分式非常小，因此采用二次模型就可以达到拟合的精度。因此，经简化，我们得到： n=1301024n-n0+0+c1n-n02 n=1,2,1024 (3.10)因此，公式中共有n0，0，c1三个参数需进行最优拟合估计。以相位曲线和直线误差的平方和作为误差评判标准，求得当n0=680， 0=1346，c1=2.9367e-5时所得误差最小。此时，盖玻片的上下表面反射光干涉信号相位曲线与直线最为接近。(c)不同像素点波长曲线(b)最优化相位曲线(a)误差-c1的变化曲线线图3.4 光谱仪最优化结果经校正后，我们可得到不同像素点对应的光谱波长关系为：n=1301024n-680+1346+2.9367e-5n-6802 n=1,2,1024 (3.11)3.2 OCT图像三维重建3.2.1 OCT三维数据场 三维数据场也称体数据场，在一个三维空间中，根据构成体数据的离散数据间的数据关系，可以分成结构化数据、非结构化数据、结构化与非结构化混合数据三类。 其中，结构化数据为可以在逻辑上组织成一个三维数组的离散空间数据，每个数据都有自己的行号、列号、层号，该空间离散数据每个元素之间具有一个三维数组每个元素间的逻辑关系。按照三维数组的均匀属性，又可将结构化数据分为规则网格结构化数据与非规则网格结构化数据，对于规则网格结构化数据，其数据组成是由均匀网格构成，其中一个网格即为一个元素，称之为体素（图3.5）。对于OCT系统所采集到的数据，其组成是均一的，故而属于这个范畴。 对于所采得的数据，其函数值是按点分布的，对于一个体素而言，假设所采得的值分别位于其八个顶点，而在三个方向上，不同像素点的间距分别设为X、Y、Z，对于八个顶点，我们设其索引值分别为（i,j,k）、（i+1,j,k），(i+1,j+1,k),（i，j+1,k），（i,j+1,k+1）,(i+1,j,k+1),(i,j+1,k+1),(i+1,j+1,k+1)。对于索引值为(i，j，k)的点，它位于空间(xi,yj,zk)处，所对应的函数值为f(xi,yj,zk)，而对于空间中其他点，其函数值则通过八个顶点的插值来获得。对于OCT数据场，它实际上是一个标量场，是一个与方向无关的三次线性插值模型，故可以通过对八个顶点进行三次线性插值获得。 图3.5 体素示意图 OCT系统所采集的数据解调后可看作一系列二维图像的叠加，对于一系列的二维图像，可构成一个三维数据场，如图3.6所示。 图3.6 二维图像构成的三维数据场对于医学上所说的三维重建，即为对所给的三维数据场进行可视化处理，其主要算法可分为两类：面绘制算法和体绘制算法。其中，面绘制算法为在三维空间的数据场构造出所需的曲面或者多边形平面，对所构造出的几何图元使用计算机图形学的方法进行绘制。体绘制算法则并不构造中间的几何图元，而是对三维空间体素按照一定的规律赋予不透明度值和颜色值，随即产生二维图像，它实际上是对图像的重新采样与合成。3.2.3 OCT三维重建面绘制算法 三维重建面绘制算法可分为两类：基于轮廓线的面绘制算法与基于体素的面绘制算法。 其中，基于轮廓线的面绘制算法又包括全局优化法、局部优化法以及曲面拟合法等。全局优化法最初由上世纪七十年代Kepple等人首次采用18，优化参数为由生成的轮廓表面之间相互连接所得的体积，而Fuches等人则采用生成轮廓表面积为优化参数19，全局优化法计算量大且精度有限，不能解决其多轮廓线拓扑重构的问题。局部优化法则包括最大体积法、相邻轮廓线同步前进法等，相对于全局优化法，计算量大大地减小，极大地提高了三维重构的效率。曲面拟合法则包括双三次B样条曲面法、非均匀有理B样条曲面法等，其曲面重构问题得以较好地解决，但是计算量偏大，也不能较好地解决多轮廓线拓扑重构问题。而基于体素的面绘制算法则包括立方体算法、移动立方体算法（Marching Cube）20、剖分立方体算法等，移动立方体算法（MC）应用较为广泛，它能够较好地解决多轮廓线拓扑重构问题，但是由于移动立方体算法需要设置小立方体逐个对体素进行遍历，因而执行效率较低，而且生成大量零散的小面片，对存储空间要求较高。而且在体