一种基于独立成分分析的全场光学相干显微成像去噪方法与流程

本发明属于图像处理技术领域，涉及一种基于独立成分分析的全场光学相干显微成像去噪方法。

背景技术：

2002年法国boccara团队基于时域oct(opticalcoherencetomography)原理发明的全场光学相干层析成像，不需扫描，一次性凝视提取组织内部全视场图像，可原生环境下实现亚细胞分辨率细胞检测，穿透深度大、分辨率高、速度快，是一种非侵入性、非破坏性和内源性获取生物组织图像的方法。全场光学相干层析成像技术的原理是当样品被具有低相干长度的光源照射时由样品后向散射的光与光源干涉成像，该技术利用光源的低相干性来隔离由样品深度方向上的虚拟切片后向散射的光。干涉仪的使用使得可以通过干涉现象产生表示从样品的给定切片选择性地发出的光的干涉信号，并且消除源自样品的其余部分的光。

全场光学相干层析成像可以获得具有微米量级分辨率的三维图像，使得可以观察到血管、血管壁、胶原蛋白、脂肪细胞等的大部分组织结构。同时，一次性凝视成像使其可以在几分钟内产生表面为几平方厘米，深度方向为几百微米的三维切片图像，可以潜在地用于诸如诊断癌症的医学应用。在手术期间的全场光学相干层析成像使得可以进行实时诊断，甚至在例如癌症的情况下优化手术操作，从而避免可能在切除后几天组织学分析揭示存在肿瘤细胞的情况下重复手术。

然而，由于全场光学相干层析成像系统使用了迈克尔逊干涉仪原理，装置对振动十分敏感，周围环境微米量级的振动将使样品成像时产生轴向运动伪影噪声，严重影响图像质量。目前全场光学相干层析成像系统仍主要放置在实验室气浮隔振平台上使用，阻碍了其仪器化和设备化。由于周围环境振动噪声和样品光学相干信号具有相对独立性，且环境振动噪声主要集中在低频成分，因此，可以将样品光学相干信号和噪声分别看作互相独立的两个信号，利用独立成分分析原理将环境振动产生的轴向运动伪影噪声滤除。

技术实现要素：

有鉴于此，本发明的目的在于提供一种基于独立成分分析的全场光学相干显微成像去噪方法，利用独立成分分析原理将低频振动噪声和全场光学相干信号分离，去除系统振动导致样品轴向位移产生的运动伪影噪声，从而通过图像处理算法而不需要气浮隔振实验平台滤除系统噪声。

为达到上述目的，本发明提供如下技术方案：

一种基于独立成分分析的全场光学相干显微成像去噪方法，具有包括以下步骤：

s1：获取由于系统低频振动产生运动伪影噪声的全场光学相干有噪图像xi，并对其预处理，即通过图像累加以增强有用信号并去除系统随机噪声；

s2：对预处理后的图像x1进行二维小波变换，并通过置高频通道小波系数为零后重构，获得保留低频信息的全场光学相干有噪图像x2；

s3：将图像x1和x2合并组成新的图像矩阵x3；

s4：对图像矩阵x3进行去均值和白化预处理，获得白化后的图像矩阵z；

s5：将白化后的图像矩阵作为独立成分分析算法的输入进行计算，即求解混合矩阵w，使输出矩阵y＝wz的各行向量之间相互独立；

s6：根据求解出的混合矩阵w，计算得到原始无噪干涉图像。

进一步，步骤s1中，对图像xi预处理，即图像累加，其计算公式为：

进一步，步骤s1中，通过全场光学相干显微成像装置获取全场光学相干图像，该成像装置由低相干光源(1)、成像透镜(3)、分光镜(4)、参考臂(5)、样品臂(6)、i/o驱动控制卡(7)和面阵ccd探测器(8)组成。

进一步，步骤s2中，所述二维小波变换采用daubechies8小波分解方法将原始有噪图像分解成不同频率范围的子图像集，包括ll低频系数部分、hl水平细节部分、lh垂直细节部分和hh高频系数部分。

进一步，步骤s3中，将图像x1和x2合并组成新的图像矩阵x3，具体包括：设原始光学干涉无噪图像为s，系统振动时样品轴向位移产生的伪影噪声为n，则有：

x1＝s+n

x2＝a21s+a22n

其中，{aij}为混合系数，i,j＝1,2，w为矩阵a的逆。

进一步，所述步骤s4具体包括以下步骤：

s41：去均值：以图像矩阵x3的某一行向量为对象，计算其元素的平均值，将该行向量的元素都减去此平均值，获得新向量，将矩阵其他行都按照相同方法处理，获得新矩阵x4；

s42：白化预处理：寻找一个线性白化矩阵v，使得矩阵z＝vx4的行向量相互正交，其中，x4为去均值后的矩阵，v＝d^-1/2e^t，d为矩阵x4的协方差矩阵e{x4x4^t}的特征向量的对角矩阵，矩阵d^-1/2通过简单的逐元素开方计算得到为e{x4x4^t}的特征向量的正交矩阵。

进一步，步骤s5中，求解混合矩阵w具体包括以下步骤：

s51：随机初始化向量wi(0)，其元素个数为矩阵z的行向量个数，i＝1,…,m，m为需要估计的独立成分的个数，且满足||wi(0)||＝1；

s52：迭代：对每个i＝1,…,m，令k＝1，更新wi(k)＝e{z(wi(k-1)^tz)³}-3wi(k-1)，其中，z为矩阵z的一列，e为矩阵z(wi(k-1)^tz)³的期望；

s53：wi(k)归一化处理：wi(k)＝wi(k)/wi(k)，确保分离结果具有单位能量；

s54：如果|wi(k)|趋近于1，则输出wi(k)，否则，令k＝k+1，并返回到步骤s52；

s55：上述迭代过程结束后，得到的wi(k)再取转置为混合矩阵w的一个行向量，从而求得混合矩阵w。

本发明的有益效果在于：本发明基于独立成分分析原理，利用小波变换产生一个保留低频伪影噪声信息的虚拟图像，和原图像一起作为独立成分分析算法的两个输入进行图像和噪声的分离，从而实现伪影噪声滤除。本发明仅依靠图像处理算法而不需要气浮隔振实验平台滤除系统噪声。

本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书来实现和获得。

附图说明

为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作优选的详细描述，其中：

图1为本发明基于独立成分分析的全场光学相干显微成像去噪方法的流程原理图；

图2为全场光学相干层析成像装置原理示意图；

图3为含运动伪影噪声的全场光学相干图像；

图4为去除的运动伪影噪声图像；

图5为去除噪声后的全场光学相干图像；

附图标记：1-低相干光源，2-会聚透镜，3-成像透镜，4-分光镜，5-参考臂，6-样品臂，7-i/o驱动控制卡，8-面阵ccd探测器，9-补偿玻片，10-高数值孔径物镜。

具体实施方式

以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需要说明的是，以下实施例中所提供的图示仅以示意方式说明本发明的基本构想，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。

其中，附图仅用于示例性说明，表示的仅是示意图，而非实物图，不能理解为对本发明的限制；为了更好地说明本发明的实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；对本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。

本发明实施例的附图中相同或相似的标号对应相同或相似的部件；在本发明的描述中，需要理解的是，若有术语“上”、“下”、“左”、“右”、“前”、“后”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此附图中描述位置关系的用语仅用于示例性说明，不能理解为对本发明的限制，对于本领域的普通技术人员而言，可以根据具体情况理解上述术语的具体含义。

请参阅图1～图5，图1为本发明基于独立成分分析的全场光学相干显微成像去噪方法，该去噪方法具有包括以下步骤：

s1：获取由于系统低频振动产生运动伪影噪声的全场光学相干有噪图像xi，并对面阵ccd探测器获取的全场光学相干有噪图像进行预处理，即通过图像累加以增强有用信号并去除系统随机噪声，其计算公式为：其中，n值越大，图像x1的信噪比越高，但所需图像采集时间越长。在本实施例中，折中取n＝10，计算后得到的图像x1，如图3所示(白色箭头为伪影噪声)。

s2：对预处理后的图像x1进行二维小波变换，并通过置高频通道hh小波系数为零，获得保留低频信息的全场光学干涉有噪图像x2。

其中，二维小波变换采用daubechies8小波分解方法将原始有噪图像分解成不同频率范围的子图像集，包括ll低频系数部分、hl水平细节部分、lh垂直细节部分和hh高频系数部分。

s3：将图像x1和x2合并组成新的图像矩阵x3，具体包括：设原始光学干涉无噪图像为s，系统振动时样品轴向位移产生的伪影噪声为n，则有：

x1＝s+n

x2＝a21s+a22n

其中，{aij}为混合系数，i,j＝1,2，w为矩阵a的逆。

s4：对图像矩阵x3进行去均值和白化预处理，获得白化后的图像矩阵z；具体包括以下步骤：

s41：去均值：以图像矩阵x3的某一行向量为对象，计算其元素的平均值，将该行向量的元素都减去此平均值，获得新向量，将矩阵其他行都按照相同方法处理，获得新矩阵x4；

s42：白化预处理：寻找一个线性白化矩阵v，使得矩阵z＝vx4的行向量相互正交，其中，x4为去均值后的矩阵，v＝d^-1/2e^t，d为矩阵x4的协方差矩阵e{x4x4^t}的特征向量的对角矩阵，矩阵d^-1/2通过简单的逐元素开方计算得到为e{x4x4^t}的特征向量的正交矩阵。

s5：将白化后的图像矩阵作为独立成分分析算法的输入进行计算，即求解混合矩阵w，使输出矩阵y＝wz的各行向量之间相互独立。求解混合矩阵w具体包括以下步骤：

s51：随机初始化向量wi(0)，其元素个数为矩阵z的行向量个数，i＝1,…,m，m为需要估计的独立成分的个数，且满足||wi(0)||＝1；

s52：迭代：对每个i＝1,…,m，令k＝1，更新wi(k)＝e{z(wi(k-1)^tz)³}-3wi(k-1)，其中，z为矩阵z的一列，e为矩阵z(wi(k-1)^tz)³的期望；

s53：wi(k)归一化处理：wi(k)＝wi(k)/wi(k)，确保分离结果具有单位能量；

s54：如果|wi(k)|趋近于1，则输出wi(k)，否则，令k＝k+1，并返回到步骤s52；

s55：上述迭代过程结束后，得到的wi(k)再取转置为混合矩阵w的一个行向量，从而求得混合矩阵w。

s6：根据求解出的混合矩阵w，计算得到原始无噪图像，如图5所示。

如图2所示，通过全场光学相干显微成像装置获取全场干涉图像，该成像装置由低相干光源1、成像透镜3、分光镜4、参考臂5、样品臂6、i/o驱动控制卡7和面阵ccd探测器8组成。参考臂5和样品臂6中均采用和补偿玻片9和高数值孔径物镜10。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。