CN106447716A - 用于条纹投影图像的方向偏微分方程滤波方法 - Google Patents

Abstract

本发明属于光学检测和图像处理技术领域，为提出对数值求解偏微分方程提供新方法，使得偏微分方程滤波效果得到提升。本发明采用的技术方案是，用于条纹投影图像的方向偏微分方程滤波方法，步骤如下：步骤1，输入一幅带有噪声的条纹投影图u；步骤2，将条纹投影图u进行离散化：步骤3：计算条纹方向与x轴方向的夹角θ_i,j；步骤4：自适应得到时间步长Δt；步骤5：自适应得到迭代次数Nc；步骤6：对于每次迭代，求出图像u每个像素的一阶偏导数以及二阶偏导数；步骤7：求出图像u每个像素的数值解步骤8：重复步骤6和步骤7，直到达到设置的最大迭代次数Nc停止迭代，此时的数值解即为滤波图像。本发明主要应用于图像处理。

Description

用于条纹投影图像的方向偏微分方程滤波方法

技术领域

本发明属于光学检测和图像处理技术领域，涉及一种用于条纹投影轮廓术中条纹图的滤波的基于拟小波离散格式的方向偏微分方程滤波方法。

背景技术

条纹投影轮廓术(Fringe Projection Profilometry,FPP)是近年来十分受欢迎的一种三维形貌测量方法，因为它具有高灵敏度、非接触、全场观测等优点，所以被广泛应用于科研和实际工程中[1-3]。对条纹投影图像的处理分析是获取待测物体形貌信息的主要手段之一，但由于条纹投影图像中总是伴随着强烈的噪声，导致无法直接利用这些图像正确的获取相位信息，所以必须对FPP条纹图像进行滤波，提高图像的对比度，去除噪声的影响。偏微分方程(Partial differential equation,简称PDE)是一种方案灵活、处理效果良好的图像滤波技术[4,5]，尤其是方向偏微分方程，只沿着条纹方向[6,7]进行滤波，适合于条纹图像的滤波。在使用PDE时，一个重要问题是如何将PDE模型进行离散化，方程的离散化程度直接影响到滤波效果以及计算效率。

参考文献

[1]杨国光,近代光学测试技术.杭州:浙江大学出版社,1997。

[2]S.Gorthi and P.Rastogi,Fringe projection techniques:Whither weare？Opt.Lasers Eng.2010,48(2):133～140。

[3]王开福,现代光测及其图像处理.北京:科学出版社,2012:100～110。

[4]F.Zhang,W.Liu,and L.Xia,“Homomorphic partial differential equationfiltering method for electronic speckle pattern interferometry fringes basedon fringe density,”Chin.Opt.Lett.7(3),210-213(2009)。

[5]Y.Chen,C.A.Z.Barcelos,and B.A.Mair,“Smoothing and edge detectionby time-varying coupled nonlinear diffusion equations,”Comput.Vis.ImageUnd.82(2),85-100(2001)。

[6]S.Chikkerur,A.N.Cartwright,and V.Govindaraju,“Fingerprintenhancement using STFT analysis,”Pattern Recogn.198-211(2007)。

[7]L.Hong,Y.Wan,and A.Jain,“Fingerprint image enhancement:algorithmand performance evaluation,”IEEE T.Pattern Anal.20(18),777-789(1998)。

[8]W.Cai and W.Zhang,An Adaptive Spline Wavelet ADI(SW-ADI)Method forTwo-dimensional Reaction-diffusion Equations.J.Comp.Phys.1998,139:9-G～12。

发明内容

对偏微分方程离散时，差分格式的选择是非常重要的，之前已有的偏微分方程离散化用到的差分格式主要以中心差分为主，为克服现有技术的不足，本发明旨在提出一种新的基于拟小波方法的差分格式，对数值求解偏微分方程提供新方法，使得偏微分方程滤波效果得到了提升。本发明采用的技术方案是，用于条纹投影图像的方向偏微分方程滤波方法，步骤如 下：

步骤1，输入一幅带有噪声的条纹投影图u；

步骤2，将条纹投影图u进行离散化，假设条纹投影图u的大小为M×N，u_i,j为(i,j)点处的灰度值，1≤i≤M，1≤j≤N，时间步长为Δt，在方程的演化过程中，t_n＝nΔt时刻的演化图像表示为，构造方向偏微分方程的离散格式为：

其中，θ_i,j表示条纹方向与x轴方向的夹角，u_x、u_y分别表示图像u沿着x轴、y轴方向的一阶偏导数，u_xx、u_yy分别表示图像u沿着x轴、y轴方向的二阶偏导数；表示对图像进行离散化，表示对图像离散化后的差分格式；

步骤3：计算条纹方向与x轴方向的夹角θ_i,j；

步骤4：自适应得到时间步长Δt；

步骤5：自适应得到迭代次数Nc；

步骤6：对于每次迭代，按照公式：

其中的拟小波离散格式：

其中E＝exp(-x²/(2σ²)),S＝sin(πx/Δ),C＝cos(πx/Δ)，x_i,j-x_i,k表示像素点之间做减法运算，Δ表示网格大小，σ表示窗口大小参数，求出图像u每个像素的一阶偏导数u_x，u_y以及二阶偏导数u_xx，u_xy和u_yy；

步骤7：基于步骤2中的方向偏微分方程的离散格式，求出图像u每个像素的数值解

步骤8：重复步骤6和步骤7，直到达到设置的最大迭代次数Nc停止迭代，此时的数值 解即为滤波图像。

步骤4自适应得到时间步长Δt包括以下步骤：

步骤4-1：输入一幅电子散斑干涉条纹图像u；

步骤4-2：给定初值，设置迭代初值u0为原始图像I，迭代时间步长Δt_k＝1，初始循环次数k＝0；

步骤4-3：计算两次的滤波结果u(Δtk)和u(2Δtk)；

步骤4-4：判断第二次滤波后的噪声与信号的相关系数corr(u(0)-u(2Δtk),u(2Δtk))是否小于第一次滤波后的噪声与信号的相关系数corr(u(0)-u(Δtk),u(Δtk))，若是，则执行步骤4-5；否则设置k＝k+1，时间步长在上次的基础上减0.05，即赋值Δtk+1＝Δtk-0.05，重复步骤4-3和步骤4-4；

步骤4-5：输出此时的Δtk。

步骤5自适应得到迭代次数Nc包括以下步骤：

步骤5-1：输入一幅电子散斑干涉条纹图像u；

步骤5-2：固定时间步长Δt为上述自适应得到的Δt，设置迭代初值u0为初始图像I，迭代次数n＝1；

步骤5-3：计算第一次迭代的滤波结果

步骤5-4：自n＝1开始，计算并保存当前滤波结果和上一次的滤波结果

步骤5-5：判断是否同时满足下述两个条件，

条件1是当前滤波后噪声与信号的相关系数corr(u(0)-u(2Δt),u(2Δt))大于上一次滤波后的噪声与信号的相关系数corr(u(0)-u(Δt),u(Δt))，

条件2是散斑系数小于0.2；

若是，则执行步骤5-6，否则设置n＝n+1，重复步骤5-4和步骤5-5；

步骤5-6：输出Nc＝n。

本发明的特点及有益效果是：

拟小波方法的离散格式相比于多点差分计算需要的选点更少，一般只需要15个点即可满足计算精度，并且这些系数可以一次性计算出来，大大缩短计算时间。相对于中心差分格式相比，拟小波离散格式拥有更高的计算精度。

附图说明：

图1噪声方差为1.2时三种方法的滤波结果。图中，(a)噪声图像，(b)SOOPDE，(c)WFF，(d)WTM。

具体实施方式

本发明基于拟小波离散的二阶单方向PDE滤波方法的主要步骤：

步骤1，输入一幅带有噪声的条纹投影图u

步骤2，将图像u进行离散化，假设图像u的大小为M×N，u_i,j(1≤i≤M,1≤j≤N)为(i,j)点处的灰度值，时间步长为Δt，在方程的演化过程中，t_n＝nΔt时刻的演化图像表示为， 构造方向偏微分方程的离散格式为：

其中，θ_i,j表示条纹方向与x轴方向的夹角，u_x、u_y分别表示图像u沿着x轴、y轴方向的一阶偏导数，u_xx、u_yy分别表示图像u沿着x轴、y轴方向的二阶偏导数；

步骤3：计算条纹方向与x轴方向的夹角θ_i,j；

步骤4：自适应得到时间步长Δt；

步骤5：自适应得到迭代次数N_c；

步骤6：对于每次迭代，按照公式：

其中的拟小波离散格式：

其中E＝exp(-x²/(2σ²)),S＝sin(πx/Δ),C＝cos(πx/Δ).求出图像u每个像素的一阶偏导数u_x，u_y以及二阶偏导数u_xx，u_xy和u_yy；x_i,j-x_i,k表示像素点之间做减法运算，Δ表示网格大小，σ表示窗口大小参数

步骤7：基于步骤2中的方向偏微分方程的离散格式，求出图像u每个像素的数值解

步骤8：重复步骤6和步骤7，直到达到设置的最大迭代次数N_c停止迭代，此时的数值解即为滤波图像。

步骤4自适应得到时间步长Δt包括以下步骤：

步骤4-1：输入一幅电子散斑干涉条纹图像u；

步骤4-2：给定初值，设置迭代初值u0为原始图像I，迭代时间步长Δt_k＝1，初始循环次数k＝0；

步骤4-3：计算两次的滤波结果u(Δtk)和u(2Δtk)；

步骤4-4：判断第二次滤波后的噪声与信号的相关系数corr(u(0)-u(2Δtk),u(2Δtk))是否小于第一次滤波后的噪声与信号的相关系数corr(u(0)-u(Δtk),u(Δtk))，若是，则执行步骤4-5；否则设置k＝k+1，时间步长在上次的基础上减0.05，即赋值Δtk+1＝Δtk-0.05，重复步骤4-3和步骤4-4；

步骤4-5：输出此时的Δtk。

步骤5自适应得到迭代次数Nc包括以下步骤：

步骤5-1：输入一幅电子散斑干涉条纹图像u；

步骤5-2：固定时间步长Δt为上述自适应得到的Δt，设置迭代初值u0为初始图像I，迭代次数n＝1；

步骤5-3：计算第一次迭代的滤波结果

步骤5-4：自n＝1开始，计算并保存当前滤波结果和上一次的滤波结果

步骤5-5：判断是否同时满足下述两个条件，

条件1是当前滤波后噪声与信号的相关系数corr(u(0)-u(2Δt),u(2Δt))大于上一次滤波后的噪声与信号的相关系数corr(u(0)-u(Δt),u(Δt))，

条件2是散斑系数小于0.2；

若是，则执行步骤5-6，否则设置n＝n+1，重复步骤5-4和步骤5-5；

步骤5-6：输出Nc＝n。

下面结合具体实施方式对本发明做进一步详细的介绍：

条纹投影轮廓术(Fringe Projection Profilometry,FPP)是近年来十分受欢迎的一种三维形貌测量方法，因为它具有高灵敏度、非接触、全场观测等优点，所以被广泛应用于科研和实际工程中。对条纹投影图像的处理分析是获取待测物体形貌信息的主要手段之一，但由于条纹投影图像中总是伴随着强烈的噪声，导致无法直接利用这些图像正确的获取相位信息，所以必须对FPP条纹图像进行滤波，提高图像的对比度，去除噪声的影响。偏微分方程(Partial differential equation,简称PDE)是一种方案灵活、处理效果良好的图像滤波技术，尤其是方向偏微分方程，只沿着条纹方向进行滤波，适合于条纹图像的滤波。在使用PDE时，一个重要问题是如何将PDE模型进行离散化，方程的离散化程度直接影响到滤波效果以及计算效率。本发明提出一种拟小波离散方法，对PDE滤波模型进行离散化，同时将改进的PDE滤波方法与窗口傅里叶方法和小波滤波方法进行比较，验证本发明提出的方法的有效性。

小波函数指的是一组有有限能量的函数,它们可以通过伸缩和平移来构造生成：

公式(1)中的参数a控制函数伸缩，参数b控制函数平移。这样通过调整尺度参数a，就可以分析函数本身的局部特性。而拟小波是对正交规范化尺度函数进行正则化处理得到的，其中正交规范化尺度函数可以定义为：

任意一个属于空间的函数f(x)都可以表示成下面的离散形式：

其中表示是在Hilbert空间的Paley-Wiener重构核。这里{x_k}是定义在无穷区间(-∞,∞)的一组离散空间点的坐标。为了实现计算，这里对a赋值，一般a＝π/ΔΔ为单元网格大小。

但是正交规范化尺度函数δ_Δ,σ(x)在坐标空间上没有明显的局部特性，在文献[8]中，Wei[8]等人建议对它们进行正则化处理(regularization procedure)。即引入一个Gauss函数作为正则因子:

Gauss函数定义为：

R_σ(x)＝exp[-x²/(2σ²)]，σ＞0 (4)

其中σ为窗口大小参数。

Gauss正则化之后的尺度函数可以定义为：

δ_α,σ(x)＝δ_α(x)R_σ(x). (5)

Gauss正则化样本尺度函数不再是精确的正交规范化小波尺度函数，因为它不满足正交规范化小波尺度函数的条件,因此，我们把Gauss正则化样本尺度函数δ_Δ,σ(x)称为拟尺度函数，由拟尺度函数生成的小波就相应地称为拟小波。

因此对于一个函数f(x)，通过δ_Δ,σ(x)可以被唯一重构为：

由于上式求和运算定义在无限区间(-∞,∞)，这在实际计算中是做不到的。因此，为了计算需要，有必要把计算域定义在有限区间。由于我们已经引进Gauss正则因子，它具有很好的快速衰减特性，即局域特性。在实际计算中，只需要在格点附近取2W(15左右)个计算点x_k，就可以达到机器精度。因此，上式可进一步写成：

上式中的上标(n)表示对空间坐标x的n阶导数，2W+1为计算带宽，它通常远远小于整个计算域的网格点数。公式(8)中，如果△和σ给定，在整个计算域中则只有一个核出现，这对数值计算来说是十分经济和有利的、我们把公式(8)称为拟小波数值离散格式。

应用拟小波离散方法对二阶偏微分方程进行离散化：

其中偏微分方程滤波模型为：

将图像大小为M×N，对于1<i<M,1<j<N的任意一点时刻n的二阶导数为：把空间x坐标均匀等分,单元网格大小记为Δx＝1/w(w是化分区间所采用的单元网格总数)，网格点坐标记为x_i＝i·Δx，(i＝0,1,2,…,M),于是有x_i-x_i+k＝-kΔx。

方程(5-25)的拟小波离散格式可以写成如下格式：

其中的拟小波离散格式：

其中E＝exp(-x²/(2σ²)),S＝sin(πx/Δ),C＝cos(πx/Δ).

这里[-W,+W]是计算带宽,在计算中这意味着计算只在很小的带宽上进行,极大减小了数值计算量，并且这些系数计算只与单元网格尺度大小有关,因此当网格分布确定后,这些系数可一次性计算出来,然后把它们保存起来供计算机调用,不需要重复计算。通过上述拟小波离散和最速下降法，可以实现对二阶单方向偏微分方程的数值求解，从而得到图像的滤波结果。

下面给出实验结果图。

我们采用图1-(a)模拟FPP条纹图作为实验图像，为了验证算法的滤波性能，采用服从均匀分布的随机噪声，幅值大小分别为0.2、0.6、0.8、1.0、1.2、1.5、2.0。这里我们仅展示方差大小为1.2的模拟噪声图，如图1-(a)所示，图1-(a)为一幅像素大小为512×512的模拟条纹图，模拟公式为：

I(x,y)＝a(x,y)+b(x,y)cos(φ(x,y)+2πf₀(x+y))+NOISE； (17)

其中相位为：

Peaks为Matlab中自带函数，其中，Re{}代表实部，条纹的载频频率为f₀＝1/16，背景强度a(x,y)为0.5φ(x,y)，调制强度b(x,y)为1，高斯随机噪声的方差为1.2。

从图中可以清晰的看到这幅图像含有大量噪声，并且条纹密集。采用我们改进的二阶单方向偏微分方程(SOOPDE)方法，窗傅里叶滤波方法(WFT)，小波变换滤波(WTM)方法，对图1-(a)做滤波处理。其中图1-(b),(c),(d),表示为SOOPDE,WFT和WTM三种方法分别对图1-(a)的滤波结果，在进行处理时各种算法的最优参数是通过反复试验得出来的并且参数设置展示在表1-1中。

由于模拟图自身具有已知的理想背景和条纹，此处我们采用峰值信噪比(PSNR)对每一种模型的滤波结果进行定量研究，其中峰值信噪比(PSNR)的定义形式为：

表1-2展示对于图1-(a)三种方法处理获得的滤波图像的峰值信噪比以及滤波时间。三种方法对应程序的计算环境为：MATLAB R2014b，Intel Core i3-2100@2.9GHz CPU，运行内存2G。

接下来对于不同噪声等级下的模拟条纹图进行三种方法的去噪性能比较分析，结果如表1-3所示。

表1-1三种方法参数设置

从图1，表1-1、表1-2和表1-3可以看出，窗傅里叶变换方法能获得很好的滤波效果但是由于需要大量的参数设置，因此需要花费大量的时间，不利于提高运算效率。小波变换方 法花费时间最少，但是滤波质量与本文提出的方法有一定的差距。本文提出的基于拟小波离散格式的方向PDE滤波方法在得到高质量滤波结果前提下，需要更少的迭代次数，这大大节省了计算时间。

表1-2三种方法的PSNR及其计算时间

表1-3不同噪声等级下三种滤波方法的PSNR及其计算时间

尽管上面结合图对本发明进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨的情况下，还可以作出很多变形，这些均属于本发明的保护之内。

Claims (3)

1.一种用于条纹投影图像的方向偏微分方程滤波方法，其特征是，步骤如下：

步骤1，输入一幅带有噪声的条纹投影图u；

步骤2，将条纹投影图u进行离散化，假设条纹投影图u的大小为M×N，u_i,j为(i,j)点处的灰度值，1≤i≤M，1≤j≤N，时间步长为Δt，在方程的演化过程中，t_n＝nΔt时刻的演化图像表示为，构造方向偏微分方程的离散格式为：

$u_{i,j}^{n+1}=u_{i,j}^n+\mathrm{\&Delta;}t\&CenterDot;({\left(u_{xx}\right)}_{i,j}^n{\cos }^2{\mathrm{\&theta;}}_{i,j}+{\left(u_{yy}\right)}_{i,j}^n{\sin }^2{\mathrm{\&theta;}}_{i,j}+2{\left(u_{xy}\right)}_{i,j}^n{\mathrm{sin\&theta;}}_{i,j}{\mathrm{cos\&theta;}}_{i,j})$

其中，θ_i,j表示条纹方向与x轴方向的夹角，u_x、u_y分别表示图像u沿着x轴、y轴方向的一阶偏导数，u_xx、u_yy分别表示图像u沿着x轴、y轴方向的二阶偏导数；表示对图像进行离散化，表示对图像离散化后的差分格式；

步骤3：计算条纹方向与x轴方向的夹角θ_i,j；

步骤4：自适应得到时间步长Δt；

步骤5：自适应得到迭代次数Nc；

步骤6：对于每次迭代，按照公式：

$\begin{array}{cc}{\left(u_t\right)}_{i,j}^n=\frac{u_{i,j}^{n+1}-u_{i,j}^n}{\mathrm{\&Delta;}t}& {\left(u_{xx}\right)}_{i,j}^n=\underset{j=1}{\overset{N}{\&Sigma;}}\underset{k=-w}{\overset{+w}{\&Sigma;}}{\mathrm{\&delta;}}_{\mathrm{\&Delta;},\mathrm{\&sigma;}}^{\left(2\right)}(x_{i,j}-x_{k,j})u\left(x_{k,j}\right)\end{array}$

${\left(u_{yy}\right)}_{i,j}^n=\underset{i=1}{\overset{M}{\&Sigma;}}\underset{k=-w}{\overset{+w}{\&Sigma;}}{\mathrm{\&delta;}}_{\mathrm{\&Delta;},\mathrm{\&sigma;}}^{\left(2\right)}(x_{i,j}-x_{i,k})u\left(x_{i,k}\right)$

${\left(u_{xy}\right)}_{i,j}^n=\frac{\&part;x}{\&part;}{\left(u_y\right)}_{i,j}^n$

${\left(u_y\right)}_{i,j}^n=\underset{i=1}{\overset{M}{\&Sigma;}}\underset{k=-w}{\overset{+w}{\&Sigma;}}{\mathrm{\&delta;}}_{\mathrm{\&Delta;},\mathrm{\&sigma;}}^{\left(1\right)}(x_{i,j}-x_{i,k})u\left(x_{i,k}\right)$

其中的拟小波离散格式：

${\mathrm{\&delta;}}_{\mathrm{\&Delta;},\mathrm{\&sigma;}}^{\left(1\right)}\left(x\right)=E(\frac{C}{x}-\frac{S\mathrm{\&Delta;}}{{\mathrm{\&pi;x}}^2}-\frac{S\mathrm{\&Delta;}}{{\mathrm{\&pi;\&sigma;}}^2}),x\&NotEqual;0$

${\mathrm{\&delta;}}_{\mathrm{\&Delta;},\mathrm{\&sigma;}}^{\left(1\right)}\left(x\right)=0,x=0$

${\mathrm{\&delta;}}_{\mathrm{\&Delta;},\mathrm{\&sigma;}}^{\left(2\right)}\left(x\right)=E\&lsqb;C(-\frac{2}{{\mathrm{\&sigma;}}^2}-\frac{2}{x^2})\&rsqb;+E\&lsqb;S(\frac{\mathrm{\&Delta;}}{{\mathrm{\&pi;x\&sigma;}}^2}+\frac{\mathrm{\&Delta;}x}{{\mathrm{\&pi;\&sigma;}}^4}-\frac{\mathrm{\&pi;}}{\mathrm{\&Delta;}x}+\frac{2\mathrm{\&Delta;}}{{\mathrm{\&pi;x}}^3})\&rsqb;,x\&NotEqual;0$

${\mathrm{\&delta;}}_{\mathrm{\&Delta;},\mathrm{\&sigma;}}^{\left(2\right)}\left(x\right)=0.x=0$

其中E＝exp(-x²/(2σ²)),S＝sin(πx/Δ),C＝cos(πx/Δ)，x_i,j-x_i,k表示像素点之间做减法运算，Δ表示网格大小，σ表示窗口大小参数，求出图像u每个像素的一阶偏导数u_x，u_y以及二阶偏导数u_xx，u_xy和u_yy；

步骤7：基于步骤2中的方向偏微分方程的离散格式，求出图像u每个像素的数值解

步骤8：重复步骤6和步骤7，直到达到设置的最大迭代次数Nc停止迭代，此时的数值解即为滤波图像。

2.如权利要求1所述的用于条纹投影图像的方向偏微分方程滤波方法，其特征是，步骤4自适应得到时间步长Δt包括以下步骤：

步骤4-1：输入一幅电子散斑干涉条纹图像u；

步骤4-2：给定初值，设置迭代初值u0为原始图像I，迭代时间步长Δt_k＝1，初始循环次数k＝0；

步骤4-3：计算两次的滤波结果u(Δtk)和u(2Δtk)；

步骤4-4：判断第二次滤波后的噪声与信号的相关系数corr(u(0)-u(2Δtk),u(2Δtk))是否小于第一次滤波后的噪声与信号的相关系数corr(u(0)-u(Δtk),u(Δtk))，若是，则执行步骤4-5；否则设置k＝k+1，时间步长在上次的基础上减0.05，即赋值Δtk+1＝Δtk-0.05，重复步骤4-3和步骤4-4；

步骤4-5：输出此时的Δtk。

3.如权利要求1所述的用于条纹投影图像的方向偏微分方程滤波方法，其特征是，步骤5自适应得到迭代次数Nc包括以下步骤：

步骤5-1：输入一幅电子散斑干涉条纹图像u；

步骤5-2：固定时间步长Δt为上述自适应得到的Δt，设置迭代初值u0为初始图像I，迭代次数n＝1；

步骤5-3：计算第一次迭代的滤波结果

步骤5-4：自n＝1开始，计算并保存当前滤波结果和上一次的滤波结果

步骤5-5：判断是否同时满足下述两个条件，

条件1是当前滤波后噪声与信号的相关系数corr(u(0)-u(2Δt),u(2Δt))大于上一次滤波后的噪声与信号的相关系数corr(u(0)-u(Δt),u(Δt))，

条件2是散斑系数小于0.2；

若是，则执行步骤5-6，否则设置n＝n+1，重复步骤5-4和步骤5-5；

步骤5-6：输出Nc＝n。