CN103814286A - 修正ps-oct的测量数据的程序及安装该程序的ps-oct系统 - Google Patents

Abstract

将PS-OCT中测量的数据非线性修正，提高PS-OCT的定量解析能力，能够进行包含病变部位的病期在内的准确的定量诊断，可以作为用于计算机诊断的有效手段使用。即使在PS-OCT(1)中的光程差包含误差而成为噪音，其分布既不为正态分布、也不在真值周围对称分布的情况下，也能使用通过用蒙特卡洛仿真解析噪音特性所获得的分布变换函数变换测量数据，由此去除系统性误差，推测被噪音掩盖的真值，将PS-OCT(1)的图像修正得更清晰。

Description

修正PS-OCT的测量数据的程序及安装该程序的PS-OCT系统

技术领域

本发明涉及一种修正基于光学相干层析成像技术(OCT:Optical coherencetomography)的测量数据来进行高精度化的程序及安装该程序的OCT系统。特别地，本发明涉及一种修正偏振敏感型光学相干层析成像装置(也称为偏振敏感光学图像测量装置、偏振敏感型OCT或偏振光学干涉断层仪等，由于通常称为PS-OCT(Polarization-Sensitive OCT)，所以以下称作“PS-OCT”。)的测量数据来进行高精度化的程序及安装该程序的PS-OCT系统。

背景技术

以往，为了以非破坏方式、高分辨率获取物体的内部信息，即后方散射、反射率分布及折射率分布的微分构造，会使用OCT。

作为用于医疗领域等的非破坏式断层测量技术之一，有光学断层图像化法“光学相干层析成像技术”(OCT)(参照专利文献1)。OCT将光用作测量探头，因此具有能够测量被测量物体的反射率分布、折射率分布、分光信息、偏振信息(双折射率分布)等的优点。

基本的OCT43以迈克尔逊干涉仪为基础，用图8对其原理进行说明。从光源44射出的光通过准直透镜45被平行化后，被分束器46分离为参照光和物体光。物体光由物体臂内的物镜47会聚在被测量物体48上，在此被散射、反射后再次返回物镜47、分束器46。

而参照光通过参照臂内的物镜49后被参照镜50反射，再次通过物镜49返回分束器46。这样返回了分束器46的物体光和参照光与物体光一起向会聚透镜51入射，被会聚在光检测器52(光电二极管等)上。

OCT的光源44使用低时间相干性光(不同时刻从光源发出的光之间非常难以干涉的光)的光源。以低时间相干性光为光源的迈克尔逊型的干涉仪中，干涉信号仅在参照臂和物体臂的距离大致相等时出现。其结果是在改变参照臂和物体臂的光程差(τ)的同时由光检测器52测量干涉信号的强度时，可以得到对应于光程差的干涉信号(干涉图)。

该干涉图的形状表示被测量物体48的纵深方向的反射率分布，能够通过一维的轴方向扫描获得被测量物体48的纵深方向的构造。这样，OCT43中能够通过光程扫描测量被测量物体48的纵深方向的构造。

在这样的轴方向(A方向)的扫描以外附加横方向(B方向)的机械性扫描(B扫描)进行二维的扫描，由此获得被测量物体的二维断面图像。作为进行这种横方向的扫描的扫描装置，使用了直接移动被测量物体的结构、物体固定而移动物镜的结构、被测量物体和物镜都固定而旋转设在物镜的光瞳面附近的电流镜的角度的结构等。

以上的基本型OCT经过发展，有了扫描光源的波长来获得光谱干涉信号的波长扫描型OCT(Swept Source OCT，简称“SS-OCT”。)和使用分光器来获得光谱信号的光谱域OCT。作为后者，有傅立叶域OCT(Fourier Domain OCT，简称“FD-OCT”。参照专利文献2)以及PS-OCT(参照专利文献3)。

波长扫描型OCT通过高速波长扫描激光器改变光源的波长，使用与光谱信号同步获取的光源扫描信号将干涉信号重新排列并施加信号处理，由此获得3维光学断层图像。还有，作为改变光源波长的单元而利用了单色仪，但也可以作为波长扫描型OCT来使用。

傅立叶域OCT的特征是用光谱仪(光谱分光器)获取来自被测量物体的反射光的波长光谱，通过对该光谱强度分布进行傅立叶变换，提取实空间(OCT信号空间)上的信号，该傅立叶域OCT不需要进行纵深方向的扫描，可以通过进行x轴方向的扫描来测量被测量物体的断面构造。

PS-OCT是一种在B-扫描的同时连续调制直线偏振后的光束的偏振状态，获取试样(被检物体)具备的偏振信息，能够测量试样的更细微的构造及折射率的各向异性的光学相干断层扫描装置。

进一步详细来说，PS-OCT与傅立叶域OCT相同，用光谱分光器获取来自被测量物体的反射光的波长光谱，但是将入射光及参照光分别通过1/2波长板、1/4波长板等作为水平直线偏振光、垂直直线偏振光、45°直线偏振光、圆偏振光，将来自被测量物体的反射光和参照光叠加通过1/2波长板、1/4波长板等，例如仅使水平偏振光成分向光谱分光器入射并干涉，仅将物体光中具备特定偏振光状态的成分取出并进行傅立叶变换。该PS-OCT也不需要进行纵深方向的扫描。

现有技术文献

专利文献

专利文献1：日本特开2002-310897号公报

专利文献2：日本特开平11-325849号公报

专利文献3：日本特开2004-028970号公报

发明内容

发明所要解决的课题

一般地，在考虑将OCT应用于生物体测量时，物体的折射率分布的微分构造虽然能够以非破坏式、高分辨率获取，但在具有因纤维状构造(纤维的伸长方向等)或牙釉质不同产生的双折射所导致的偏振依存性的生物试样的测定中，分辨率降低，并且无法获取构造。

对于这一点，PS-OCT使来自某特定部分的散射光成分和某偏振光状态的参照光干涉，在其干涉成分中突出反映偏振特性，还能获得相位延迟信息(光程差量)，其结果是能够获取在纵深方向的断面的某特定部分的偏振信息、双折射信息。所以，PS-OCT作为眼科、循环器官、呼吸器官、消化器官等广泛需要非侵袭性组织判别的临床领域中的检查装置是有效的。

这样，PS-OCT作为将组织的偏振特性以非侵袭性可视化的光学断层成像装置，适用于组织对比这种定性观察。然而，本发明者在以PS-OCT的高精度化为目的锐意研究开发过程中发现，PS-OCT观察病变部位、在病期的准确分类等定量观察中的精度未必足够，存在应该进一步改善的地方。

本发明的目的是解决在上述现有的PS-OCT中得出的应该改善的地方，课题是实现将PS-OCT中测量的数据非线性修正，提高PS-OCT的定量解析能力，能够进行病期的准确定量诊断等的程序及安装该程序的PS-OCT系统。

更具体来说，通过PS-OCT测量的光程差(相位延迟量、相位差的量)包含误差，分布在真值的周围(将其称为光程差的分布或相位差的分布)。该分布一般不对称，也包含真值不一定位于中心的非对称的情况。所以，用获取标准平均值的相位解析无法得到正确的光程差。

本发明是为了解决这种现有的PS-OCT的问题，课题是实现在光程差包含误差而成为噪音，其分布既不为正态分布，也不在真值周围对称分布的情况下，也能通过由蒙特卡洛仿真解析噪音特性所获得的分布变换函数变换测量数据，由此去除系统性误差，推测被噪音掩盖的真的值，将PS-OCT的图像修正得更清晰的程序，同时实现安装该程序的PS-OCT系统。

用于解决课题的手段

为了解决上述课题，本发明提供一种安装在处理由PS-OCT获得的测量数据的计算机上的程序，其特征在于，使计算机作为通过非线性变换函数将由PS-OCT测量的、关于光程差的相位差的分布数据变换为对称的相位差的分布数据并进行平均，由此求出真的相位量的单元来发挥作用。

为了解决上述课题，本发明提供一种安装在处理由PS-OCT获得的测量数据的计算机上的程序，其特征在于，使计算机作为以下单元来发挥作用：预先将由分别关于0～π的多个相位量的分布组成的相位差的分布数据的组针对每个不同的ESNR生成多个组并存储的单元；确定将由PS-OCT测量的相位差的分布数据变换为对称的相位差的分布数据的非线性变换函数的系数的单元；根据由PS-OCT测量的相位差的分布数据估算ESNR，并基于该估算的ESNR选择针对所述每个不同的ESNR生成的相位差的分布的多个组中的某一个组的单元；在根据所述选择的相位差的分布的组所确定的系数下，通过所述变换函数将由PS-OCT测量的相位差的分布数据变换为对称的相位差的分布数据的单元；以及通过对所述对称的相位差的分布数据进行平均来计算真的相位量的单元。

为了解决上述课题，本发明提供一种PS-OCT系统，具备PS-OCT和处理由该PS-OCT获得的测量数据的计算机，所述PS-OCT系统的特征在于，所述计算机具备输入装置、输出装置、CPU及存储装置，作为通过非线性变换函数将由PS-OCT测量的相位差的分布数据变换为对称的相位差的分布数据并进行平均，由此求出真的相位量的单元来发挥作用。

为了解决上述课题，本发明提供一种PS-OCT系统，具备PS-OCT和处理由该PS-OCT获得的测量数据的计算机，所述PS-OCT系统的特征在于，所述计算机具备输入装置、输出装置、CPU及存储装置，作为以下单元发挥作用：预先将由分别关于0～π的多个相位量的分布组成的多个相位差的分布的组针对每个不同的ESNR生成多个组并存储在存储装置中的单元；确定将由PS-OCT测量的相位差的分布数据变换为对称的相位差的分布数据的非线性变换函数的系数的单元；根据由PS-OCT测量的相位差的分布数据估算ESNR，并基于该估算的ESNR选择针对所述每个不同的ESNR生成的相位差的分布的多个组中的某一个组的单元；在根据所述选择的相位差的分布的组确定的系数下，通过所述变换函数将由PS-OCT测量的相位差的分布数据变换为对称的相位差的分布数据的单元；以及通过对所述对称的相位差的分布数据进行平均来计算真的相位量的单元。

所述非线性变换函数为后述的式(2)。

确定所述变换函数的系数的单元的结构为，针对所述非线性变换函数，从0～π为止每隔π/m改变，生成下面的m+1的联立方程式即后述的式(3)，求解该联立方程式来制作变换函数的系数b_i的表格。其中，m是划分π的次数。

求解所述联立方程式来确定所述变换函数的系数的单元的结构，也可以是数值性地求出后述式(4)的矩阵DM的伪逆矩阵DM⁺，用后述式(5)计算，唯一确定最佳的变换函数的系数b_i。

确定所述变换函数的系数的单元的结构也可以是使用变分法唯一确定最佳的变换函数的系数b_i，使得后述式(6)所示的均方误差为最小。

发明的效果

根据本发明，能够将PS-OCT中测量的数据非线性修正，提高PS-OCT的定量解析能力，因此虽然现有的PS-OCT局限于病变部位的形态观察等用途，但本发明的PS-OCT系统能够进行包含病变部位的病期在内的准确的定量诊断，可以作为用于计算机诊断的有效手段来使用。

附图说明

图1是表示本发明涉及的PS-OCT系统的全体结构的图。

图2是表示本发明涉及的PS-OCT系统的进行图像处理的计算机的图。

图3是光程差(相位差、相位延迟量)的分布图(相位差的分布图)，右栏的4个图是PS-OCT测量的相位差的分布图，左栏是通过仿真生成的相位差的分布图。

图4表示本发明中的试验例1中比较仿真与测量值的比较试验的结果。

图5是表示本发明中的试验例2的试验结果的图像。

图6是图5的(c)、(d)中虚线表示的断面图像。

图7是表示本发明中的试验例3的试验结果的图像。

图8是说明以往的OCT的图。

具体实施方式

以下参照附图对本发明涉及的修正PS-OCT的测量数据的程序及安装该程序的PS-OCT系统的实施方式进行说明。

本发明涉及的PS-OCT系统具备PS-OCT和处理PS-OCT获得的图像数据的图像处理装置。图像处理装置使用通常的计算机，本发明涉及的程序作为修正PS-OCT获得的图像数据的手段发挥功能。

(PS-OCT)

PS-OCT已通过日本特许第4344829号公报等而公知，它是作为本发明的前提的技术，因此说明其概要。

本发明涉及的偏振敏感光学图像测量装置在B扫描的同时(同步)将来自光源的偏振光束(由偏振镜进行了直线偏振的光束)通过EO调制器(偏振调制器、电气光学调制器)连续调制，将这个连续调制偏振光所得的偏振光束分开，将一方作为入射光束扫描照射在试样上，获得其反射光(物体光)的同时，将另一方作为参照光，通过二者的光谱干涉进行OCT测量。

并且，该装置的特征是，其结构为通过同时由2个光检测器测定该光谱干涉成分中的垂直偏振成分(H)和水平偏振成分(V)，获得表示试样的偏振特性的琼斯矢量(H图像和V图像)。

图1是表示本发明涉及的偏振敏感光学图像测量装置的光学系统的全体结构的图。图1所示的偏振敏感光学图像测量装置1具备光源2、偏振镜3、EO调制器4、光纤耦合器(光耦合器)5、参照臂6、试样臂7、分光器8等光学构件。该偏振敏感光学图像测量装置1的光学系统中光学构件互相由光纤9结合，但也可以是不由光纤结合的类型的构造(自由空间型)。

光源2使用具有宽频带光谱的超辐射发光二极管(SLD:Super LuminessentDiode)。另外，光源2也可以是脉冲激光器。光源2上依次连接有准直透镜11、将来自光源2的光直线偏振的偏振镜3、进相轴被设置为45°方向的EO调制器(偏振调制器、电气光学调制器)4、会聚透镜13及光纤耦合器5。

EO调制器4通过将进相轴固定在45°方向，并对加在该EO调制器4上的电压进行正弦调制，连续改变进相轴和与其垂直的迟相轴之间的相位差(光程差)，由此，当从光源2发出、经偏振镜3变为(纵向)直线偏振光的光向调制器4入射时，以上述调制周期如直线偏振光→椭圆偏振光→直线偏振光……等那样被调制。EO调制器4使用市售的EO调制器即可。

光纤耦合器5上经由分支的光纤9连接了参照臂6和试样臂7。参照臂6上依次设有偏振控制器(polarization controller)10、准直透镜11、偏振镜12、会聚透镜13及参照镜(固定镜)14。参照臂6的偏振镜12用于选择即使如上述那样调制偏振状态，从参照臂6返回的光的强度也不变化的方向。该偏振镜12的方向(直线偏振光的偏振方向)的调整与偏振控制器10联合进行。

试样臂7上依次设有偏振控制器15、准直透镜11、固定镜24、电流镜16、会聚透镜13，来自光纤耦合器5的入射光束由2轴的电流镜扫描，向试样17照射。来自试样17的反射光作为物体光再次返回光纤耦合器5，与参照光叠加，作为干涉光束被送给分光器8。

分光器8具备依次连接的偏振控制器18、准直透镜11、(偏振敏感型体积相位全息)衍射光栅19、傅立叶变换透镜20、偏振分束器21及2个光检测器22、23。该实施方式中，使用线阵CCD相机（1维CCD相机）作为光检测器22、23。从光纤耦合器5送来的干涉光束通过准直透镜11被校准，通过衍射光栅19被分光为干涉光谱。

通过衍射光栅19被分光的干涉光谱光束在傅立叶变换透镜20进行傅立叶变换，被偏振分束器21分为水平及垂直成分，分别由2个线阵CCD相机(光检测器)22、23检测。这2个线阵CCD相机22、23用于检测水平及垂直偏振信号二者的相位信息，因此2个线阵CCD相机22、23必须参与同一个分光器的组成。

另外，光源2、参照臂6、试样臂7及分光器8上，分别设有偏振控制器10、15、18，它们调整从光源2送给参照臂6、试样臂7、分光器8的各个光束的初始偏振状态，进行控制使得在EO调制器4被连续调制的偏振状态对于参照光和物体光也维持振幅相互固定且相对偏振状态相互固定的关系，进一步地对于连接在光纤耦合器5上的分光器8保持振幅固定和相对偏振状态固定。

另外，当修正包含2个线阵CCD相机22、23的分光器8时，EO调制器4停止。阻挡参照光，将载玻片（slide-glass）和反射镜置于试样臂7上。该配置保证水平及垂直偏振成分的峰值位置相同。然后，来自载玻片的背面和反射镜的OCT信号由2个分光器8检测。OCT信号的峰值的相位差被监测。

该相位差在所有光轴方向的深度都应该为零。下一步，为了通过包含2个线阵CCD相机22、23的分光器8获得复谱，对信号进行加窗以及傅立叶逆变换。该相位差在所有的频率都应该为零，因此通过监视这些值可以调整2个线阵CCD相机22、23的物理位置，使得相位差最小。

PS-OCT的原理：

本发明的特征如下。将来自光源2的光进行直线偏振，该直线偏振后的光束由EO调制器4进行连续的偏振状态调制。即EO调制器4通过将进相轴固定在45°方向，并对加在该EO调制器4上的电压正弦调制，连续改变进相轴和与其垂直的迟相轴之间的相位差(偏振角：光程差)，由此，当从光源2发出、经偏振镜3变为(纵向)直线偏振光的光向EO调制器4入射时，以上述调制周期被调制为直线偏振光→椭圆偏振光→直线偏振光……等那样。

然后，被直线偏振的偏振光束由EO调制器4进行连续的偏振状态调制，同时同步进行B扫描。即在1次B扫描的期间，进行多个周期EO调制器4的偏振光的连续调制。这里，1周期是指偏振角(光程差)φ在0～2π变化的期间。简而言之，该1周期的期间，来自偏振镜的光的偏振光连续调制为直线偏振光(垂直偏振光)→椭圆偏振光→直线偏振光(水平偏振光)……等那样。

这样连续调制偏振光束的偏振，同时在试样臂7上由电流镜16将入射光束在试样17上扫描来进行B扫描，在分光器8中，针对其反射光、即物体光和参照光的干涉光谱，由2个线阵CCD相机22、23检测其水平偏振成分及垂直偏振成分。由此，通过1次B扫描能够获得分别与水平偏振成分及垂直偏振成分对应的2幅A-B扫描图像。

如上所述，1次B扫描期间进行多个周期的偏振光束的偏振的连续调制，在各周期(1周期)的连续调制期间由2个线阵CCD相机22、23检测出的水平偏振成分及垂直偏振成分的偏振信息分别成为1个像素的偏振信息。1周期的连续调制中，由2个线阵CCD相机22、23与定时信号同步地检测偏振信息，1周期内可以将检测次数(测取次数)适当确定为4次、8次等。

这样，将1次B扫描内获得的2幅A-B扫描图像的数据在B扫描方向进行一维傅立叶变换。于是出现0次、1次、-1次的峰值。这里，分别提取出0次的峰值，仅使用该数据进行傅立叶逆变换，则得到H0、V0图像。同样地，分别提取出1次的峰值，仅使用该数据进行傅立叶逆变换，则得到H1、V1图像。

根据H0、H1图像，可以在表示试样17的偏振特性的琼斯矩阵的成分中求出J(1,1)及J(1,2)。并且，根据V0、V1图像，可以在表示试样17的偏振特性的琼斯矩阵的成分中求出J(2,1)及J(2,2)。

还有，这里，琼斯矩阵的成分J(1,1)、J(1,2)、J(2,1)、J(2,2)基于专利文献1的记载而记载，分别相当于本发明的后述式(1)的琼斯矩阵的成分I_0H(z)、I_1H(z)、I_0V(z)、I_0V(z)。

这样，在1次B扫描中得到包含4个偏振特性的信息。然后，若将这4个信息分别与普通的FD-OCT同样地在A扫描方向上进行傅立叶变换，则得到1次峰值具有试样17的深度方向的信息、并且分别对应于偏振特性的4幅A-B图像的测量数据。

(本发明的程序及PS-OCT系统)

通过上述结构的PS-OCT获得的测量数据被输入作为图像处理装置使用的计算机30。该计算机30是普通的计算机，如图2所示，具备输入部31、输出部32、CPU33、存储装置34及数据总线35。

本发明涉及的PS-OCT系统上安装的程序是存储在计算机30的存储装置34中的程序，其作为修正输入给计算机30的通过PS-OCT获得的图像的测量数据，将PS-OCT的图像修正得更清晰的单元，使计算机30发挥功能。通过安装该程序，本发明涉及的PS-OCT系统成为具备将图像修正得更清晰的单元的系统。

即，以往将从PS-OCT测量的原始数据(原始的偏振相位延迟数据)获得的图像直接用于诊断及其他各种用途。本发明的程序及PS-OCT系统中，对通过PS-OCT测量的原始数据以非线性函数进行修正，能够进行高精度的生物体双折射测量，该非线性函数为蒙特卡洛法(Monte-Carlo analysis)，预先使用本发明的程序在计算机30中通过仿真来设计。

以下，对本发明的程序及安装该程序的PS-OCT系统进行详细说明。

图3是光程差(相位差、相位延迟量)的分布图(相位差的分布图)。横轴为相位，纵轴为光程差产生的频度。图3的右栏是基于PS-OCT测量的光程差的相位差的分布，图3的左栏是通过仿真生成的相位差的分布。

相位差的分布，将横轴上彼此相隔π/m刻度的m+1个数据(m+1个相位差的分布)作为1组，改变ESNR生成多个组。例如，取m=60，将横轴上彼此相隔π/60刻度的61个数据作为1组，改变ESNR生成多个组。

图3中是对于2个不同的ESNR的组，分别表示关于π/2及π的相位差的分布的图。

而根据PS-OCT的光程差的噪音模型是对琼斯矩阵的各个要素上设定加法性噪音，如下面的式(1)所示。这里，S是真的光程差的值，N是噪音，I是测定值。角标0、1是测量通道(参照图1的光检测器22、23)。说明书第9页9-11行中说明的琼斯矩阵的成分对应于在真值S上加上噪音N所得的测定值I。这里z是深度方向的坐标，琼斯矩阵的各要素相当于从表面(z=0)到深度z为止的光程差等的积分量。

式(1)

$\left[\begin{array}{cc}{I}_{0,H}\left(z\right)& I_{1,H}\left(z\right)\\ I_{0,V}\left(z\right)& I_{1,V}\left(z\right)\end{array}\right]=\left[\begin{array}{cc}{S}_{0,H}\left(z\right)& S_{1,H}\left(z\right)\\ S_{0,V}\left(z\right)& S_{1,V}\left(z\right)\end{array}\right]+\left[\begin{array}{cc}{N}_{0,H}\left(z\right)& N_{1,H}\left(z\right)\\ N_{0,V}\left(z\right)& N_{1,V}\left(z\right)\end{array}\right]$

通常，使用S/N比(也称作SN比、SNR等)作为表示测量数据的噪音的评价值，而本发明中使用ESNR作为更有效地表示噪音的评价值。本说明书中，ESNR表示为γ，以下式定义。

1/γ=1/4(1/SNR_S,0+1/SNR_S,1+1/SNR_r,0+1/SNR_r,1)

而该式中，s是PS-OCT的样本臂，r是参照臂。0和1是次数。

为了将多个独立的SNR以单一值表示，导入了ESNR。角标0、1是测量通道(参照图1的光检测器22、23)。SNR_S,0、SNR_S,1、SNR_r,0、SNR_r,1分别通过角标的组合，表示从样本臂向0通道贡献的SNR、从样本臂向1通道贡献的SNR、从参照臂向0通道贡献的SNR、从参照臂向1通道贡献的SNR。

如图3右栏从上起第2个图所示，PS-OCT测量的光程差(相位延迟量、相位差)包含误差，分布在真值(π/2)的周围。而当ESNR变小时，则如图3右栏从上起第1个图所示，误差的幅度变大。

另外，该相位差的分布如图3右栏从上起第3、第4个图所示，并不一定对称，也包含真值不一定位于中心的非对称情况。

简而言之，也包含光程差包含误差而成为噪音、相位差的分布不为正态分布、不在真的相位量的周围对称分布的情况。这样的情况下，在获取标准平均值的相位解析中无法获得正确的光程差。

本发明涉及的程序是作为通过使用由蒙特卡洛仿真解析噪音特性所获得的分布变换函数变换PS-OCT测量的图像的测量数据，去除系统性误差、推测被噪音掩盖的真的相位量、将PS-OCT的图像修正得更清晰的单元，使计算机30发挥功能的程序，本发明的系统是安装该程序并进行上述修正的PS-OCT系统。

原理：

本发明涉及的程序作为修正根据PS-OCT测量的图像的测量数据的实测光程差的分布(相位差的分布)，求出真的相位量δ_T(真的光程差的值)的单元，使计算机30发挥功能，首先对其原理(算法)说明如下。

本发明中，使用非线性函数f，从根据PS-OCT测量的图像的测量数据的实测的相位量δ_M变换为对称的相位差的分布δ_E。简而言之，对称的相位差的分布δ_E用实测的相位量δ_M的幂函数δ_E=f(δ_M)表示。本说明书中将该函数称为变换函数，用下面的式(2)表示。这意味着为了将测量值的非对称分布δ_M变换为以正确值作为叠加中心的分布δ_E所进行的非线性映射。

式(2)

${\mathrm{\δ}}_E=\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{b}_i\left(\mathrm{\γ}\right){\mathrm{\δ}}_M^i=b_0\left(\mathrm{\γ}\right)+b_1\left(\mathrm{\γ}\right){\mathrm{\δ}}_M+b_2\left(\mathrm{\γ}\right){\mathrm{\δ}}_M^2+...+b_n\left(\mathrm{\γ}\right){\mathrm{\δ}}_M^n$

该式(2)中，b_i为变换函数的系数。最佳的变换函数的系数b_i(i为0、1、2、……n)是ESNRγ的函数，如下这样求解联立方程式来确定，或使用最小二乘法确定，预先制作b_i(γ)的表格。

为此，将相位值从0到π每隔π/m改变，根据蒙特卡洛法求出δ_M。系数b_i通过求解式(3)所示的联立方程式来确定。式(3)中，<δ_M，π/m>表示每隔π/m，例如表示每隔π/60仿真所得的、从0到π的光程差的仿真值δ_M。这里，m为例如1、2、……60。还有，式(3)中的<>的意思是表示由多次蒙特卡洛仿真得到的值的平均值。具体来说，将式(3)如式(4)所示以矩阵改写，求出系数b_i(γ)。

式(3)

$\left\{\begin{array}{c}0=b_0\left(\mathrm{\&gamma;}\right)+b_1\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,0}>+b_2\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,0}^2>+...+b_n\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,0}^n>\\ \frac{\mathrm{\&pi;}}{m}=b_0\left(\mathrm{\&gamma;}\right)+b_1\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,\mathrm{\&pi;}/m}>+b_2\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,\mathrm{\&pi;}/m}^2>+...+b_n\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,\mathrm{\&pi;}/m}^n>\\ .\\ .\\ .\\ \mathrm{\&pi;}=b_0\left(\mathrm{\&gamma;}\right)+b_1\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,\mathrm{\&pi;}}>+b_2\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,\mathrm{\&pi;}}^2>+...+b_n\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,\mathrm{\&pi;}}^n>\end{array}\right.$

式(4)

D_T=D_M·B

这里，

D_T=[O,n/1,...,n/m,...,n]^T,B=[b_o(Y),...,b_n(γ)]^T,

而为了确定系数，数值性求出通过式(4)的仿真得到的δ_M,π/m ⁿ的n×m的矩阵DM的伪逆矩阵DM⁺，用下面的式(5)计算求出系数。还有，式(4)中的<>的意思是表示由多次蒙特卡洛仿真得到的值的平均值。

式(5)

B≡D_M ⁺·D_T

使用最小二乘法确定系数的手段可以使用变分法，为使下面的式(6)所示

的均方误差R²为最小，而通过蒙特卡洛法唯一确定最佳的变换函数的系数b_i(γ)。这里也将相位值从0到π为止每隔π/m改变，例如取m=60，每隔π/60改变，通过蒙特卡洛法求出仿真值δ_M。k表示将相位渐次测量的测量次数，为0、1、2、……m，例如m=60等。

式(6)

$R^2=\underset{k=0}{\overset{m}{\mathrm{\&Sigma;}}}{[{\mathrm{\&delta;}}_{T,k}-(b_0\left(\mathrm{\&gamma;}\right)+b_1\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,k}>+b_2\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,k}^2>+...+b_n\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,k}^n>)]}^2$

最佳的变换函数的系数b_i确定后，表示对称的相位差的分布δ_E的变换函数δ_E=f(δ_M)即确定，代入PS-OCT实测的相位量δ_M，则求出δ_E。

然后，通过将δ_E如下面的式(7)所示进行平均(式(7)中的<>是平均的意思)，可以求出真的相位量δ_T。按为了求出一个相位量而进行的实测次数进行平均。

式(7)

6_T=<δ_M>

结构：

本发明的原理如上，根据该原理，在本发明的PS-OCT的系统中，本发明涉及的程序使作为图像处理装置使用的计算机30作为以下的单元发挥功能。

(1)相位差的分布的生成单元

相位差的分布的生成单元，是预先将如图3的左栏所示的光程差的相位差的分布进行仿真并生成多个的单元。

具体来说，生成的相位差的分布为将横轴按每π/m、例如m=60即按π/60划分，将61个数据(61个相位差的分布)作为1组，改变ESNR来将相位差的分布的组生成多组。这些多组数据生成后被存储在存储装置34中。

(2)变换函数的系数的确定单元

变换函数的系数的确定单元是针对多个相位差的分布的组，分别确定变换函数的系数b_i的单元。

该变换函数的系数b_i如上所述，求解式(3)的联立方程式来确定，或者在式(6)中使用最小二乘法来确定。b_i被确定后生成其表格。该变换函数的系数b_i的表格被存储在存储装置34中。

(3)相位差的分布的选择单元

相位差的分布的选择单元是在由PS-OCT测量的图像的实测相位差的分布、相位量的修正之前，根据其测量数据及系统设计等估算ESNR的值，针对上述(1)所示的多个相位差的分布(如图3的左栏的相位差的分布)的组，确定使用哪个组。

(4)基于变换函数的相位差的分布的变换单元

变换函数的计算单元是在所述存储装置34中存储的变换函数的系数b_i或其表格中，使用与由选择单元选择的相位差的分布的组对应的系数b_i或其表格，代入式(2)的b_i中，根据PS-OCT测量的实测的光程差的相位量δ_M求出表示对称的相位差的分布的变换函数δ_E=f(δ_M)的单元。

(5)真的相位量的计算单元

真的相位量的计算单元是根据表示如上所述求出的对称的相位差的分布的δ_E=f(δ_M)，基于所述式(7)进行平均，由此计算出真的相位量δ_T的单元。

(试验例1)

本发明者为了评价本发明的程序及安装该程序的PS-OCT系统，进行了将仿真值与实测值比较的比较试验。在该比较试验中，试样(测量对象物)为玻璃板、1/8波长板及1/4波长板，分别产生的光程差的真值为0、π/2及π。

还有，仿真值是如上所述，预先用蒙特卡洛仿真确定b_i的系数，使用由此获得的系数变换该仿真数据所得的值。以上述系数变换实测数据后为实测值。蒙特卡洛仿真中由于混入了随机噪音、且b_i在有限的项数截止，因此虽然是仿真却并不完全一致。

该试验中，为了进行基于PS-OCT的测量，使用了1.3μm的探头。该PS-OCT在空气中的深度解析为8.3μm，具有30,000行/秒的测定速度。通过在试样的前方放置可以改变光强度的非偏振滤光片，改变ESNR来进行测量。

图4表示比较仿真与试验的比较试验的结果。虚线表示真的光程差的值(0、π/2及π)。实线(b)、(c)、(e)是不进行本发明的修正，仅单纯将光程差的实测值取平均而得到的。实线(a)、(d)、(f)是本发明的仿真的结果，用于制作系数的仿真数为65536个，这是用于求系数的数据，也是用于得到结果而使用的相位数据。

位于实线(b)、(c)、(e)及(a)、(d)、(f)附近的●■＋标记是分别根据试验得到的光程差的实测值及将该实测值通过本发明修正所得的试验修正值，这里使用64个测定数据。仿真与试验修正值在绝对值及ESNR依存性等方面基本一致。

(试验例2)

本发明者为了评价本发明，进行了将通常用于评价PS-OCT的鸡胸肉(ex-vivo，离体)作为试样进行测量的比较试验。

该试验中，对于试样的10mm的部分，进行了由512条A线组成的B扫描。简而言之，是将在A方向(深度方向)的测量中获得的A线通过B扫描(与A方向垂直的B方向的扫描)获得512条。而且该测量对于试样在同一地点进行128次。小于8db的ESNR不稳定，因此去除。

图5是表示试验例2的试验结果的图像图。图5(a)表示OCT强度图像，图5(b)表示将说明书第10页21-23行的ESNR的值用图像浓度表示的图像，图5(c)表示单纯平均相位图像，图5(d)表示通过本发明修正的最优化的图像。

在试样的较深部分，图5(c)的单纯平均相位图像中对比度降低，但图5(d)则获得了对比度不太降低的结果。

图6是图5的(c)、(d)中以虚线表示的断面图像，细线表示ESNR的值，粗虚线表示光程差。图5(c)、(d)中，ESNR的值(以细线表示)在试样越往内部(较深部分)变得越低，简而言之表示信号强度变小、噪音变大。

对于光程差，使用本发明的图6(d)所示的图比基于单纯平均的图6(c)所示的图显示得更为明确。特别地，图6(d)表示的图可以看出距离0及π附近较近的值被重现。由此评价为本发明在PS-OCT的图像的解析中是有效的。

另外，不进行变换的处理(参照图5(c))中，0及π附近的相位数据缺失，相位变化失去连续性，而在通过本发明处理后的结果(参照图5(d))中，可以看出0及π附近的相位也被重现，深度方向的光程差的绝对值得到了正确的测量。

(试验例3)

本发明者为了评价本发明，通过PS-OCT对人的眼底(in-vivo，体内)进行了测量。使用1.0μm的探头光束，对于11.0μm的眼底深度，通过具有30,000行/秒的测定速度的PS-OCT进行该测量。

图7表示由该试验获得的图像。图7(a)是之前的单纯平均处理图像，图7(b)是进行了本发明的修正处理的图。如图7(b)所示，可以认为根据本发明，在相位量为0的表面附近被准确地重现，得到更接近真实的相位量、相位差的分布。

以上，对本发明涉及的程序及安装该程序的PS-OCT系统的实施方式进行了说明，显然本发明并不局限于这样的实施方式，当然在权利要求书中记载的技术事项的范围内具有各种各样的实施方式。

产业上的可利用性

以往，PS-OCT局限于病变部位的形态观察用途，而本发明涉及的程序及安装该程序的PS-OCT系统由于具有上述的结构，在眼科、循环器官、呼吸器官、消化器官等广泛需要非侵袭性组织判别的临床领域，能够进行包含病变部位的病期在内的准确的定量诊断，可以作为用于计算机诊断的有效手段使用。

例如，通过置入前眼部PS-OCT装置来使用，能够对青光眼手术后的手术瘢痕进行非接触、非侵袭性观察，其结果是眼科治疗的方针能够大大提高效率。

符号说明

1 PS-OCT(偏振敏感光学图像测量装置)

2 光源

3、12   偏振镜

4 EO调制器(偏振调制器、电气光学调制器)

5 光纤耦合器(光耦合器)

6 参照臂

7 试样臂

8 分光器

9 光纤

10、15、18 偏振控制器(polarization controller)

11   准直透镜

13   会聚透镜

14   参照镜(固定镜)

16   电流镜

17   试样

19   衍射光栅

20   傅立叶变换透镜

21   偏振分束器

22、23  光检测器(线阵CCD相机)

24   固定镜

30   计算机(图像处理装置)

31   输入部

32   输出部

33   CPU

34   存储装置

35   数据总线

Claims (12)

1.一种安装在处理由PS-OCT获得的测量数据的计算机上的程序，其特征在于，

使计算机作为通过非线性变换函数将由PS-OCT测量的、关于光程差的相位差的分布数据变换为对称的相位差的分布数据并进行平均，由此求出真的相位量的单元来发挥作用。

2.一种安装在处理由PS-OCT获得的测量数据的计算机上的程序，其特征在于，

使计算机作为以下单元来发挥作用：

预先将由分别关于0～π的多个相位量的分布组成的相位差的分布数据的组针对每个不同的ESNR生成多个组并存储的单元；

确定将由PS-OCT测量的相位差的分布数据变换为对称的相位差的分布数据的非线性变换函数的系数的单元；

根据由PS-OCT测量的相位差的分布数据估算ESNR，并基于该估算的ESNR选择针对所述每个不同的ESNR生成的相位差的分布的多个组中的某一个组的单元；

在根据所述选择的相位差的分布的组所确定的系数下，通过所述变换函数将由PS-OCT测量的相位差的分布数据变换为对称的相位差的分布数据的单元；以及

通过对所述对称的相位差的分布数据进行平均来计算真的相位量的单元。

3.根据权利要求2所述的程序，其特征在于，

所述非线性变换函数是下面的式(2)，

式(2)

${\mathrm{\&delta;}}_E=\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}{b}_i\left(\mathrm{\&gamma;}\right){\mathrm{\&delta;}}_M^i=b_0\left(\mathrm{\&gamma;}\right)+b_1\left(\mathrm{\&gamma;}\right){\mathrm{\&delta;}}_M+b_2\left(\mathrm{\&gamma;}\right){\mathrm{\&delta;}}_M^2+...+b_n\left(\mathrm{\&gamma;}\right){\mathrm{\&delta;}}_M^n$

其中，式(2)中，δE为对称的相位差的分布，δM为由PS-OCT测量的相位量，b_i(i=0、1、2、3、……n)为变换函数的系数，γ为ESNR。

4.根据权利要求3所述的程序，其特征在于，

确定所述变换函数的系数的单元的结构为，针对所述非线性变换函数，从0到π为止每隔π/m改变，生成下面的m+1的联立方程式即式(3)，求解该联立方程式来制作变换函数的系数b_i的表格，

式(3)

$\left\{\begin{array}{c}0=b_0\left(\mathrm{\&gamma;}\right)+b_1\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,0}>+b_2\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,0}^2>+...+b_n\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,0}^n>\\ \frac{\mathrm{\&pi;}}{m}=b_0\left(\mathrm{\&gamma;}\right)+b_1\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,\mathrm{\&pi;}/m}>+b_2\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,\mathrm{\&pi;}/m}^2>+...+b_n\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,\mathrm{\&pi;}/m}^n>\\ .\\ .\\ .\\ \mathrm{\&pi;}=b_0\left(\mathrm{\&gamma;}\right)+b_1\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,\mathrm{\&pi;}}>+b_2\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,\mathrm{\&pi;}}^2>+...+b_n\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,\mathrm{\&pi;}}^n>\end{array}\right.$

其中，式(3)中的m是划分π的次数。

5.根据权利要求4所述的程序，其特征在于，

求解所述联立方程式来确定所述变换函数的系数的单元的结构为，数值性地求出式(4)的矩阵DM的伪逆矩阵DM⁺，用式(5)计算，唯一确定最佳的变换函数的系数b_i，

式(4)

D_T=D_M·B

这里，

D_T=[O,n/1,...,n/m,...,n]^T,B=[b_o(Y),...,b_n(γ)]^T,

式(5)

B≡D_M ⁺·D_T

其中，式(4)中的m是划分π的次数。

6.根据权利要求3所述的程序，其特征在于，

确定所述变换函数的系数的单元的结构为，使用变分法唯一确定最佳的变换函数的系数b_i，使得下面的式(6)所示的均方误差为最小，

式(6)

$R^2=\underset{k=0}{\overset{m}{\mathrm{\&Sigma;}}}{[{\mathrm{\&delta;}}_{T,k}-(b_0\left(\mathrm{\&gamma;}\right)+b_1\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,k}>+b_2\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,k}^2>+...+b_n\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,k}^n>)]}^2$

其中，式(6)中的m是划分π的次数。

7.一种PS-OCT系统，具备PS-OCT和处理由该PS-OCT获得的测量数据的计算机，所述PS-OCT系统的特征在于，

所述计算机具备输入装置、输出装置、CPU及存储装置，作为通过非线性变换函数将由PS-OCT测量的相位差的分布数据变换为对称的相位差的分布数据并进行平均，由此求出真的相位量的单元来发挥作用。

8.一种PS-OCT系统，具备PS-OCT和处理由该PS-OCT获得的测量数据的计算机，所述PS-OCT系统的特征在于，

所述计算机具备输入装置、输出装置、CPU及存储装置，

作为以下单元发挥作用：

预先将由分别关于0～π的多个相位量的分布组成的多个相位差的分布的组针对每个不同的ESNR生成多个组并存储在存储装置中的单元；

确定将由PS-OCT测量的相位差的分布数据变换为对称的相位差的分布数据的非线性变换函数的系数的单元；

根据由PS-OCT测量的相位差的分布数据估算ESNR，并基于该估算的ESNR选择针对所述每个不同的ESNR生成的相位差的分布的多个组中的某一个组的单元；

在根据所述选择的相位差的分布的组确定的系数下，通过所述变换函数将由PS-OCT测量的相位差的分布数据变换为对称的相位差的分布数据的单元；以及

通过对所述对称的相位差的分布数据进行平均来计算真的相位量的单元。

9.根据权利要求8所述的PS-OCT系统，其特征在于，

所述非线性变换函数为下面的式(2)，

式(2)

${\mathrm{\&delta;}}_E=\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}{b}_i\left(\mathrm{\&gamma;}\right){\mathrm{\&delta;}}_M^i=b_0\left(\mathrm{\&gamma;}\right)+b_1\left(\mathrm{\&gamma;}\right){\mathrm{\&delta;}}_M+b_2\left(\mathrm{\&gamma;}\right){\mathrm{\&delta;}}_M^2+...+b_n\left(\mathrm{\&gamma;}\right){\mathrm{\&delta;}}_M^n$

其中，式(2)中，δE为对称的相位差的分布，δM为由PS-OCT测量的相位量，b_i(i=0、2、3、……n)为变换函数的系数，γ为ESNR。

10.根据权利要求9所述的PS-OCT系统，其特征在于，

确定所述变换函数的系数的单元的结构为，针对所述非线性变换函数，从0～π为止每隔π/m改变，生成下面的m+1的联立方程式即式(3)，求解该联立方程式来制作变换函数的系数b_i的表格，

式(3)

$\left\{\begin{array}{c}0=b_0\left(\mathrm{\&gamma;}\right)+b_1\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,0}>+b_2\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,0}^2>+...+b_n\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,0}^n>\\ \frac{\mathrm{\&pi;}}{m}=b_0\left(\mathrm{\&gamma;}\right)+b_1\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,\mathrm{\&pi;}/m}>+b_2\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,\mathrm{\&pi;}/m}^2>+...+b_n\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,\mathrm{\&pi;}/m}^n>\\ .\\ .\\ .\\ \mathrm{\&pi;}=b_0\left(\mathrm{\&gamma;}\right)+b_1\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,\mathrm{\&pi;}}>+b_2\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,\mathrm{\&pi;}}^2>+...+b_n\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,\mathrm{\&pi;}}^n>\end{array}\right.$

其中，式(3)中的m是划分π的次数。

11.根据权利要求10所述的PS-OCT系统，其特征在于，

求解所述联立方程式来确定变换函数的系数的单元的结构为，数值性地求出式(4)的矩阵DM的伪逆矩阵DM⁺，用式(5)计算，唯一确定最佳的变换函数的系数b_i，

式(4)

D_T=D_M·B

这里，

D_T=[O,n/1,...,n/m,...,n]^T,B=[b_o(Y),...,b_n(γ)]^T,

式(5)

B≡D_M ⁺·D_T

其中，式(4)中的m是划分π的次数。

12.根据权利要求9所述的PS-OCT系统，其特征在于，

确定所述变换函数的系数的单元的结构是使用变分法唯一确定最佳的变换函数的系数b_i，使得下面的式(6)所示的均方误差为最小，

式(6)

$R^2=\underset{k=0}{\overset{m}{\mathrm{\&Sigma;}}}{[{\mathrm{\&delta;}}_{T,k}-(b_0\left(\mathrm{\&gamma;}\right)+b_1\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,k}>+b_2\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,k}^2>+...+b_n\left(\mathrm{\&gamma;}\right)<{\mathrm{\&delta;}}_{M,k}^n>)]}^2$

其中，式(6)中的m是划分π的次数。

Images