CN111220092B - 一种光学测量中条纹滤波器的构造方法 - Google Patents

Abstract

本发明公开了一种光学测量中条纹滤波器的构造方法，该滤波器属于空间滤波器，采用一种正弦辅助的筛分计算方法，通过设置合适的滤波尺寸，直接将条纹分量从图像中分离出来，首先，在空域中对条纹图像进行二维自相关计算，提取自相关系数的横、纵中心轴，将两者进行相关性计算，得到的系数涵盖横、纵方向的共同周期信息；其次，提取相关性系数中的周期信息，并分析其变化情况和趋势，进而估算出条纹分量中对应于最大局部变形的最大截止空间周期，用该周期作为条纹滤波器的尺寸；最后，通过正弦辅助的筛分计算直接滤出含有重要信息的条纹分量，以用于后续的计算处理。本发明的有益效果在于：该方法简单、速度快、鲁棒性好，非常具有实用的潜力。

Description

一种光学测量中条纹滤波器的构造方法

技术领域

本发明涉及光学、计算机视觉和人工智能技术领域，尤其涉及一种光学测量中条纹滤波器的构造方法。

背景技术

光栅投影三维测量技术已经实现对静态或慢速运动物体的有效测量，但该技术对动态快速物体的测量还始终处于研究阶段。动态快速测量最重要的一点就是对单个条纹图像的处理，其计算方法还受计算精度、速度等的综合影响。

传统的基于傅里叶思想的方法能够实现对单帧条纹图像的有效测量，但是会受到低频的背景环境光所带来的干扰，使得频域滤波变得极为困难和不准确，而为了改善该问题，出现了各种时频分析方法，如窗口傅里叶变换、小波变换等等，可这些方法均需要较大计算量且需要预先设置参数，不适于实际应用。近年来的深度学习是一种新兴的技术手段，来实现条纹图的快速动态提取，可是目前为止，大部分基于深度学习的方法仍需要输入至少两幅不同频率的条纹图像来保证精度，如果仅用一张图像，则还是需要事先去除背景光干扰。

为了去除背景光，提取条纹分量，出现了不少基于傅里叶思想等的技术，它们仍然面临参数预置、计算量大等问题。经验模态分解近些年被大量用于条纹图的预先处理，可是这种方法通常采用先将条纹图分解成很多单模函数，然后再将分解所得的单模函数进行分类，从而提取其中的条纹分量。这种做法需要大量的迭代计算，非常浪费时间，并且实际应用情况复杂，对单模函数的组合分类很难做到非常鲁棒。本发明的条纹滤波器则是打破传统的先分解再分类的思维习惯，而是直接分析出所要提取的条纹分量的周期变动范围，找到最大的截止周期，用该周期值作为经验模态分解中筛分方法的尺寸，从而实现一次性分解，既不会迭代计算，也不再需要额外的分类组合工作，大大提高了工作效率和精度，保证了鲁棒性。

发明内容

技术问题：

条纹图像的背景光滤除很难做到自适应、快速且具有鲁棒性；基于傅里叶思想的时频分析方法容易受到频谱混叠、计算量大、需要预置很多参数等问题；基于经验模态分解的方法容易却受到先分解、再分类的思想禁锢，会产生大量的迭代计算，分类规则需要设置的极其复杂却仍难以保证鲁棒性。针对这些问题，本发明提出一种构造光学测量中条纹滤波器的方法，通过空间统计方法分析出条纹分量的统计周期范围，从而估算出滤波器的尺寸，结合基于正弦辅助的筛分过程可以一次性将条纹分量准确、快速地提取出来，方法简单、有效，且能保证鲁棒性，具有较好的应用前景。

技术方案：

一种光学测量中条纹滤波器的构造方法，基于正弦辅助的筛分计算，可自适应地设置滤波器的尺寸，从而快速、有效地直接提取出有用的条纹分量，其特征在于，具体步骤为：

步骤1：对条纹图像进行自相关计算，提取横、纵方向的自相关系数值，进一步计算两者的相关性系数，截取有用的系数部分并从中提取出代表正弦条纹周期的信息；

步骤2：分析所提取周期的变化情况和变化趋势，估算出条纹滤波器的尺寸；

步骤3：使用基于正弦辅助的筛分计算，直接滤取有用的条纹信息。

进一步的，所述步骤1对长度为r、高度为c的二维条纹图像I(x,y)进行如式(1)的二维自相关函数的计算：

其中，τ₁和τ₂为自相关系数A(τ₁,τ₂)的坐标，取值范围分别为-c+1≤τ₁≤c-1和-r+1≤τ₁≤r-1。

选取A(τ₁,τ₂)的中心横轴A(0,τ₂)和中心纵轴A(τ₁,0)，分别对它们进行差分计算以突出两者的正弦波动，即得：

计算ΔA(τ₁)和ΔA(τ₂)的相关性系数如下：

R(τ)＝∑ΔA(τ₂)ΔA(τ₁-τ),

其中τ为相关性系数R(τ)的坐标，其坐标范围为-(c+r)+1≤τ≤(c+r)-1。

对R(τ)进行差分计算以突出正弦波动，可得：

ΔR(τ)＝R(τ+1)-R(τ).

找到ΔR(τ)的最小值，记为Min[ΔR(τ)](Min[·]代表求最小值的运算符号)，以该值为起始，向右截取坐标长度为Min[(r,c)]/2的ΔR(τ)值，作为相关性最大、最能体现正弦波动的相关性系数，写作ΔR₀(τ₀)，其中τ₀代表新的坐标，取值范围为1≤τ₀≤Min[(r,c)]/2+1。

对ΔR₀(τ₀)求所有极小值点，求这些极小值点两两之间的距离作为正弦波动的统计周期值，写作P(i)，1≤i≤L，L代表统计周期值的个数。

进一步的，所述步骤2分析所提取统计周期值的变化情况和变化趋势，估算滤波器的尺寸，具体步骤为：

1)判断下式所示条件是否满足：

如果满足，则设置滤波器尺寸T为：

T＝P(1),

否则，设置滤波器尺寸为：

其中，κ为比例系数，i₁为分割第一组变化趋势一样的统计周期的坐标序数。

2)计算上式中的比例系数κ。通过下式判断统计周期值序列P(i)是否呈整体增长趋势：

如果是，则令

否则，则令κ＝1。

3)确定式(7)中的坐标序数i₁。首先，找到P(i)的第一个极值点，令该极值点的坐标序数为i₁；如果P(i)不存在极值点，则找出P(i)中最大周期值对应的序数为i₁；其次，分析P(i₁)临近的周期值情况来进一步修正确认i₁，具体做法为：

I)如果P(i₁+1)也是一个极值点，并且满足一下三种情况中的任意一个，则修改i₁＝i₁+1，三种情况为：

情况1：P(i₁)是极小值而P(i₁+1)是极大值且P(i₁+1)≥P(1)；

情况2：极大值而P(i₁+1)是极小值且P(i₁+1)≤P(1)；

情况3：从P(1)到P(i₁)都是相等值，而P(i₁+1)也是个极值或者P(i₁+1)＝P(i₁+2)。

II)如果P(i₁+1)不是极值点但P(i₁+2)≥P(1)，则修改i₁＝i₁+1。

III)如果P(i₁+1)不是极值点，但P(i₁+1)与P(i₁)相等且P(i₁)是极小值，则修改i₁＝i₁-1。

进一步的，所述步骤3使用基于正弦辅助的筛分计算，直接滤取有用的条纹信息，具体步骤为：

1.1)对I(x,y)使用快速经验模态分解法(EFEMD)进行单层尺度的分解，其中EFEMD单层分解函数的输入为I(x,y)以及步骤2所得的滤波尺寸T，输出为BIMF(x,y)，根据BIMF(x,y)的幅值和频率，设计二维正弦分布图s(x,y)＝acos(2πfx)·cos(2πfy)，其中，幅值a＝max(|BIMF(x,y)|)，max(·)为求最大值操作，频率

1.2)将二维正弦分布图s(y,x)与初始分解图BIMF(x,y)分别进行相加和相减，得到

对

进行EFEMD单层分解，输入参数为尺寸T，得到第一中间分解图

同样地，对

进行EFEMD单层分解，得到第二中间分解图

1.3)对所得到的第一中间分解图

和第二中间分解图

求平均值：

BIMF_PMS(x,y)即为最终滤得的I(x,y)所含的条纹分量，而滤除的背景光分量则为

Back(x,y)＝I(x,y)-BIMF_PMS(x,y)。

有益效果：

对比于基于傅里叶思想的条纹处理方法，本发明提供的方法可适用于复杂、非线性变化的条纹图信号，且本发明具有自适应性、计算量小；

对比于其他基于经验模态分解的条纹处理方法，本发明不需要迭代分解计算，不需要对分解结果进一步分类组合，可以直接一次性提取条纹图中的条纹分量，滤除背景光分量。此外，由于正弦辅助的设计，还会避免经验模态分解的模式混叠现象发生，从而保证筛分过程的精确性；

总体而言，本发明计算量小、鲁棒性好，可处理绝大部分现实生活中出现的条纹图情况。

附图说明

图1为本发明的整体算法流程示意图。

图2为一个仿真的复杂条纹图像。

图3为用本发明方法所过滤分解得到的条纹分量。

图4为用本发明方法滤出的背景光分量。

图5为本发明对实际采集的光学测量条纹图像处理所滤出的背景光分量。

具体实施方式

为详细说明本发明的技术内容、构造特征、所实现目的及其效果，下面将结合本发明实施例中附图对本发明详予说明。

如图1所示，本发明的一种光学测量中条纹滤波器的构造方法，基于正弦辅助的筛分计算，可自适应地设置滤波器的尺寸，从而快速、有效地直接提取出有用的条纹分量，方法包括：

步骤1：对条纹图像进行自相关计算，提取横、纵方向的自相关系数值，进一步计算两者的相关性系数，截取有用的系数部分并从中提取出代表正弦条纹周期的信息；

步骤2：分析所提取周期的变化情况和变化趋势，估算出条纹滤波器的尺寸；

步骤3：使用基于正弦辅助的筛分计算，直接滤取有用的条纹信息。

以上是发明的总体流程步骤，具体实施时，可按下述步骤进行：

1、针对步骤1的具体实施方式为：

对长度为r、高度为c的二维条纹图像首先进行简单的噪声滤波处理(可任意选择去噪方法)，去噪之后的条纹图像I(x,y)如图2所示，对其进行如下式的二维自相关函数的计算：

其中，τ₁和τ₂为自相关系数A(τ₁,τ₂)的坐标，取值范围分别为-c+1≤τ₁≤c-1和-r+1≤τ₁≤r-1。

选取A(τ₁,τ₂)的中心横轴A(0,τ₂)和中心纵轴A(τ₁,0)，分别对它们进行差分计算以突出两者的正弦波动，即得：

计算ΔA(τ₁)和ΔA(τ₂)的相关性系数如下：

R(τ)＝∑ΔA(τ₂)ΔA(τ₁-τ),

其中τ为相关性系数R(τ)的坐标，其坐标范围为-(c+r)+1≤τ≤(c+r)-1。

对R(τ)进行差分计算以突出正弦波动，可得：

ΔR(τ)＝R(τ+1)-R(τ).

找到ΔR(τ)的最小值，记为Min[ΔR(τ)](Min[·]代表求最小值的运算符号)，以该值为起始，向右截取坐标长度为Min[(r,c)]/2的ΔR(τ)值，作为相关性最大、最能体现正弦波动的相关性系数，写作ΔR₀(τ₀)，其中τ₀代表新的坐标，取值范围为1≤τ₀≤Min[(r,c)]/2+1。

对ΔR₀(τ₀)求所有极小值点，求这些极小值点两两之间的距离作为正弦波动的统计周期值，写作P(i)，1≤i≤L，L代表统计周期值的个数。

2、针对步骤2的具体实施方式为：

1)判断下式所示条件是否满足：

如果满足，则设置滤波器尺寸T为：

T＝P(1),

否则，设置滤波器尺寸为：

其中，κ为比例系数，i₁为分割第一组变化趋势一样的统计周期的坐标序数。

2)计算式(7)的比例系数κ。通过下式判断统计周期值序列P(i)是否呈整体增长趋势：

如果是，则令

否则，则令κ＝1。

3)确定式(7)中的坐标序数i₁。首先，找到P(i)的第一个极值点，令该极值点的坐标序数为i₁；如果P(i)不存在极值点，则找出P(i)中最大周期值对应的序数为i₁；其次，分析P(i₁)临近的周期值情况来进一步修正确认i₁，具体做法为：

I)如果P(i₁+1)也是一个极值点，并且满足一下三种情况中的任意一个，则修改i₁＝i₁+1，三种情况为：

情况1：P(i₁)是极小值而P(i₁+1)是极大值且P(i₁+1)≥P(1)；

情况2：极大值而P(i₁+1)是极小值且P(i₁+1)≤P(1)；

情况3：从P(1)到P(i₁)都是相等值，而P(i₁+1)也是个极值或者P(i₁+1)＝P(i₁+2)。

II)如果P(i₁+1)不是极值点但P(i₁+2)≥P(1)，则修改i₁＝i₁+1。

III)如果P(i₁+1)不是极值点，但P(i₁+1)与P(i₁)相等且P(i₁)是极小值，则修改i₁＝i₁-1。

3、对步骤3的具体实施方式为：

1.1)对I(x,y)使用快速经验模态分解法(EFEMD)进行单层尺度的分解，其中EFEMD单层分解函数的输入为I(x,y)以及步骤2所得的滤波尺寸T，输出为BIMF(x,y)。

根据BIMF(x,y)的幅值和频率，设计二维正弦分布图s(x,y)＝acos(2πfx)·cos(2πfy)，其中，幅值a＝max(|BIMF(x,y)|)，max(·)为求最大值操作，频率

具体地，EFEMD的单层分解步骤如下：

步骤2.1：对I(x,y)进行镜面延拓，产生I_e(x,y)，具体地，基于Matlab的镜面延拓步骤如下：

步骤2.1.1：确定延拓尺寸为es＝round[min(c,l)/2]；

步骤2.1.2：提取图像Im(x,y)宽度为es的外框，并求其镜面作为外围延拓值，即：

水平顶部延拓值为ext1＝flipud[Im(1:es,1:c)]，其中flipud(·)为Matlab中自带的矩阵上下翻转函数；水平底部延拓值为ext2＝flipud[Im(l-es+1:l,1:c)]；

左边延拓值为ext3＝fliplr[Im(1:l,1:es)]，其中fliplr(·)为Matlab中自带的矩阵左右翻转函数；右边延拓值为ext4＝fliplr[Im(1:l,c-es+1:c)]；

左上角延拓值为ext5＝fliplr{rot90[Im(1:es,1:es)]}，其中rot90(·)为Matlab自带函数以实现矩阵逆时针旋转90度；右上角延拓值为ext6＝Im(1:es,c-es+1:c)^T，(·)^T为矩阵求转置；左下角延拓值为ext7＝Im(l-es+1:l,1:es)^T；右下角延拓值为ext8＝fliplr{rot90[Im(l-es+1:l,c-es+1:c)]}；

步骤2.1.3：将上述步骤所有的延拓值与Im(x,y)组合起来，最终得到镜面延拓后的结果为：

步骤2.2：对I_e(x,y)分别进行两次形态学膨胀操作，即

MaxD(x,y)＝imdilate(I_e(x,y),ms)

和

MinD(x,y)＝-imdilate(-I_e(x,y),ms)

其中ms＝strel(‘disk’,AvDistPre)，strel(·)为Matlab工具箱中自带函数，用以创建指定形状的结构元素，imdilate(·)为Matlab中自带的形态学膨胀函数，这些函数亦可用其他工具箱，其输入参数也可自行改变；求I_e(x,y)的极大值分布为：

MaxM(x,y)＝～(I_e(x,y)-MaxD(x,y))

I_e(x,y)的极小值分布为：

MinM(x,y)＝～(I_e(x,y)-MinD(x,y))

其中“～”为“逻辑非”操作，即将矩阵中元素为0的值置为1、非0的值置为0。

步骤2.3：对步骤2.2中的两个形态学膨胀结果求平均，并对得到的平均值进行二维卷积平滑操作：

其中mC＝fspecial(“disk”,round(T/2))，conv2(·)为Matlab工具箱中的二维卷积操作函数，fspecial(·)为Matlab工具箱中建立预定义滤波算子的函数，这些函数亦可采用其他工具箱具有相同功能的函数；

步骤2.4：将步骤2.3中得到的平滑结果从I_e(x,y)中减去：

B(x,y)＝I_e(x,y)-SM(x,y)

依照步骤2.1将B(x,y)中延拓区域的值去除，得到的结果BIMF(x,y)为EFEMD的最终输出。

1.2)将二维正弦分布图s(y,x)与初始分解图BIMF(x,y)分别进行相加和相减，

得到

对

进行EFEMD单层分解，输入参数为尺寸T，得到第一中间分解图

同样地，对

进行EFEMD单层分解，输入参数为尺寸T，得到第二中间分解图

1.3)对所得到的第一中间分解图

和第二中间分解图

求平均值：

BIMF_PMS(x,y)即为最终滤得的I(x,y)所含的条纹分量，如图3所示；而滤除的背景光分量则为Back(x,y)＝I(x,y)-BIMF_PMS(x,y)，如图4所示。

Claims (1)

1.一种光学测量中条纹滤波器的构造方法，其特征在于，具体步骤为：

步骤1：对条纹图像进行自相关计算，提取横、纵方向的自相关系数值，进一步计算两者的相关性系数，截取有用的系数部分并从中提取出代表正弦条纹周期的信息；

步骤2：分析所提取周期的变化情况和变化趋势，估算出条纹滤波器的尺寸；

步骤3：使用基于正弦辅助的筛分计算，直接滤取有用的条纹信息；

所述步骤1对长度为r、高度为c的二维条纹图像I(x,y)进行如式(1)的二维自相关函数的计算：

其中，τ₁和τ₂为自相关系数A(τ₁,τ₂)的坐标，取值范围分别为-c+1≤τ₁≤c-1和-r+1≤τ₂≤r-1；

选取A(τ₁,τ₂)的中心横轴A(0,τ₂)和中心纵轴A(τ₁,0)，分别对它们进行差分计算以突出两者的正弦波动，即得：

计算ΔA(τ₁)和ΔA(τ₂)的相关性系数如下：

R(τ)＝∑ΔA(τ₂)ΔA(τ₁-τ), (3)

其中τ为相关性系数R(τ)的坐标，其坐标范围为-(c+r)+1≤τ≤(c+r)-1；

对R(τ)进行差分计算以突出正弦波动，可得：

ΔR(τ)＝R(τ+1)-R(τ). (4)

找到ΔR(τ)的最小值，记为Min[ΔR(τ)]，Min[·]代表求最小值的运算符号，以该值为起始，向右截取坐标长度为Min[(r,c)]/2的ΔR(τ)值，作为相关性最大、最能体现正弦波动的相关性系数，写作ΔR₀(τ₀)，其中τ₀代表新的坐标，取值范围为1≤τ₀≤Min[(r,c)]/2+1；

对ΔR₀(τ₀)求所有极小值点，求这些极小值点两两之间的距离作为正弦波动的统计周期值，写作P(i)，1≤i≤L，L代表统计周期值的个数；

所述步骤2分析所提取统计周期值的变化情况和变化趋势，估算滤波器的尺寸，具体步骤为：

1)判断式(5)所示条件是否满足：

如果满足，则设置滤波器尺寸T为：

T＝P(1), (6)

否则，设置滤波器尺寸为：

其中，κ为比例系数，i₁为分割第一组变化趋势一样的统计周期的坐标序数；

2)计算式(7)中的比例系数κ，通过式(8)判断统计周期值序列P(i)是否呈整体增长趋势：

如果满足，则令

否则，则令κ＝1；

3)确定式(7)中的坐标序数i₁，首先，找到P(i)的第一个极值点，令该极值点的坐标序数为i₁；如果P(i)不存在极值点，则找出P(i)中最大周期值对应的序数为i₁；其次，分析P(i₁)临近的周期值情况来进一步修正确认i₁，具体做法为：

I)如果P(i₁+1)也是一个极值点，并且满足一下三种情况中的任意一个，则修改i₁＝i₁+1，三种情况为：

情况1：P(i₁)是极小值而P(i₁+1)是极大值且P(i₁+1)≥P(1)；

情况2：P(i₁)是极大值而P(i₁+1)是极小值且P(i₁+1)≤P(1)；

情况3：从P(1)到P(i₁)都是相等值，而P(i₁+1)也是个极值或者P(i₁+1)＝P(i₁+2)；

II)如果P(i₁+1)不是极值点但P(i₁+2)≥P(1)，则修改i₁＝i₁+1；

III)如果P(i₁+1)不是极值点，但P(i₁+1)与P(i₁)相等且P(i₁)是极小值，则修改i₁＝i₁-1；

所述步骤3使用基于正弦辅助的筛分计算，直接滤取有用的条纹信息，具体步骤为：

1.1)对条纹图像I(x,y)使用快速经验模态分解法进行单层尺度的分解，其中快速经验模态分解法单层分解函数的输入为条纹图像I(x,y)以及步骤2所得的滤波尺寸T，输出为BIMF(x,y)，根据BIMF(x,y)的幅值和频率，设计二维正弦分布图s(x,y)＝acos(2πfx)·cos(2πfy)，其中，幅值a＝max(|BIMF(x,y)|)，max(·)为求最大值操作，频率

1.2)将二维正弦分布图s(y,x)与初始分解图BIMF(x,y)分别进行相加和相减，得到

对

进行快速经验模态分解法单层分解，输入参数为尺寸T，得到第一中间分解图

同样地，对

进行快速经验模态分解法单层分解，得到第二中间分解图

1.3)对所得到的第一中间分解图

和第二中间分解图

求平均值：

BIMF_PMS(x,y)即为最终滤得的I(x,y)所含的条纹分量，而滤除的背景光分量则为Back(x,y)＝I(x,y)-BIMF_PMS(x,y)。

Images