-
基于SWT变换的NSCT分解的低频图像能更好保留图像完整信息,将适应于高信噪比的偏微分方程的非线性扩散滤波应用在低频图像,使低频图像在去噪的同时也保持图像结构完整,最高尺度高频图像更好表达图像细节,利用快速非局部均值滤波的高效和冗余性,在保持高频图像细节的同时也减少滤波时间, 最后为了表达完整高低频图像信息,进行小波低频最大值融合。
-
离散平稳小波具有冗余性,更能捕捉图像的特征与细节,它在每两个滤波器系数间插入零值来实现滤波器的延展,具有平移不变性,可避免Gibbs振荡的问题。分解公式为:
$ \left\{ \begin{array}{l} {c_{j,k}} = \sum\limits_q {{h_{0,j}}\left( {q - 2k} \right){c_{j - 1,q}}} \\ {d_{j,k}} = \sum\limits_q {{h_{1,j}}\left( {q - 2k} \right){d_{j - 1,q}}} \end{array} \right. $
(1) 式中, cj, k为近似部分系数;dj, k为细节系数;j为离散化尺度系数,k为离散化平移系数。h0, j=〈φ1, 0, φ0, k〉,h1, j=〈Ψ1, 0 , Ψ0, k〉;h0, j和h1, j为分别在h0, h1中间插入2j-1个0;Ψ1, 0和Ψ0, k是小波函数;φ1, 0和φ0, k是尺度函数;q=0, 1, …Q-1,Q为信号长度。
-
轮廓波(contourlet)变换由于不具有平移不变性,容易造成频谱混叠,所以将非下采样思想引入contourlet变换,其中包括非下采样塔式滤波器组(nonsubsampled pyramid, NSP)和非下采样方向滤波器组(nonsubsampled directional filter bank, NSDFB)。NSCT变换步骤为:(1)将原图像进行NSP分解,得到低频和高频两部分; (2)用NSDFB分解将高频图像分成若干方向; (3)组合两个滤波器,实现NSCT分解。
-
PERONA和MALIK提出PM方程理论[20],后人经过改进,Alvarez扩散模型代替传统扩散模型,该模型沿着边缘扩散,可以更好地保护边缘[21]。Alvarez模型为:
$ \left\{ \begin{array}{l} \frac{{\partial u}}{{\partial t}} = \left| {\nabla u} \right|{\rm{div}}\left( {\frac{{\nabla u}}{{\left| {\nabla u} \right|}}} \right)\\ u\left( {x,y,0} \right) = I\left( {x,y} \right) \end{array} \right. $
(2) 式中, u为图像灰度,▽u为图像灰度的梯度,t为时间参量,I(x, y)是初始图像像素点(x, y)处的灰度值,u(x, y, 0)为初始状态下图像灰度,K为梯度阈值。
利用Tukey’s biweight函数代替传统扩散函数,使扩散具有自适应性,数学模型为:
$ g\left( x \right) = \left\{ \begin{array}{l} x \times \left[ {1 - {{\left( {\frac{x}{K}} \right)}^2}} \right],\left( {\left| x \right| \le K} \right)\\ 0,\left( {{\rm{other}}} \right) \end{array} \right. $
(3) 将Alvarez扩散模型与g函数结合,得到新的数学模型为:
$ \frac{{\partial u}}{{\partial t}} = g\left( {\left| {\nabla {G_\varepsilon } * u} \right|} \right)\left| {\nabla u} \right|{\rm{div}}\left( {\frac{{\nabla u}}{{\left| {\nabla u} \right|}}} \right) $
(4) 式中, Gε为高斯核,ε为尺度参量,*为卷积符号。
利用中心差法,将新的数学模型转化图像扩散项为:
$ \begin{array}{*{20}{c}} {\left| {\nabla u} \right|{\rm{div}}\left( {\frac{{\nabla u}}{{\left| {\nabla u} \right|}}} \right) = g\left( {\left| {\nabla {G_\varepsilon } * u} \right|} \right) \cdot }\\ {\frac{{{I_{xx}}I_y^2 - 2{I_{xy}}{I_x}{I_y} + {I_{yy}}I_x^2}}{{I_x^2 + I_y^2}}} \end{array} $
(5) 离散形式为:
$ \begin{array}{*{20}{c}} {{I_{i,o}}^{p + 1} = {I_{i,o}}^p + \Delta t \times }\\ {\left[ {g\left( {\left| {\nabla {G_\varepsilon } * u} \right|} \right) \cdot \frac{{{I_{xx}}I_y^2 - 2{I_{xy}}{I_x}{I_y} + {I_{yy}}I_x^2}}{{I_x^2 + I_y^2}}} \right]} \end{array} $
(6) 式中, Δt为迭代时间步长,p为迭代次数,Ii, o为像素点(i, o)的灰度值,Ix和Ixx为初始图像I对x的1阶偏导数和2阶偏导数,Iy和Iyy为初始图像I对y的1阶偏导数和2阶偏导数,Ixy为初始图像I首先对x的1阶偏导数,再对y的1阶偏导数。
-
传统非局部均值滤波利用图像冗余信息和自相似性质,在滤波同时,保留图像边缘细节及纹理信息,获得明显效果[22]。
它的主要思想是当前像素的估计值由与它具有相似邻域结构的像素加权平均得到,突破局部信息的限制,在整个图像进行寻找对应像素值。具体实现办法是每个像素的权值用以它为中心的图像区域与以当前像素为中心的图像区域之间的相似程度来确定,相似程度由高斯加权欧氏距离得到,最后由权值获得加权平均值,进而估计当前像素值。令含有噪声图像为:v(c)={u(c), c∈A},其中A是噪声图像,条纹图去噪图像像素加权平均值为:
$ v\left( c \right) = \sum\limits_{f \in c} {w\left( {c,f} \right)u\left( f \right)} $
(7) 式中, w(c, f)为权值,它是以分别计算c和f为中心的图像区域相似性来得到的。数学模型为:
$ \begin{array}{*{20}{c}} {w\left( {c,f} \right) = }\\ {\frac{1}{{z\left( c \right)}}\exp \left[ { - \frac{{\left\| {v\left( {{N_c}} \right) - v\left( {{N_f}} \right)} \right\|_{2,\sigma }^2}}{h}} \right]} \end{array} $
(8) 式中, σ为噪声标准偏差; 0≤ w(c, f) ≤1且∑w(c, f)=1;h为控制指数函数衰减的去噪参量; N为图像区域; z(c)为标准化常数; v(Nc)是以c为中心的区域,v(Nf)是以f为中心的区域。$\left\| {v\left( {{N_c}} \right) - v\left( {{N_f}} \right)} \right\|_{2, \sigma }^2 $为高斯加权欧氏距离,用来计算v(Nc)和v(Nf)的相似程度。
FROMENT等人提出构造一个关于像素差值的积分图像,从而降低两邻域相似性计算时间,进而降低算法复杂度,缩短算法时间。设图像共有M个像素点,搜索窗口大小为D×D,两邻域窗口大小相等为d×d,两邻域窗口在搜索窗口中滑动,算法复杂度为O(d2),对于搜索窗口需计算D2个像素点相似性,所以整个图像算法复杂度为O(MD2d2)。构造像素差值积分图像:
$ {S_l}\left( x \right) = \sum\limits_{{z_1} \le {x_1},{z_2} \le {x_2}} {{s_l}\left( z \right),\left( {x = \left( {{x_1},{x_2}} \right)} \right)} $
(9) 式中, $ {S_l}\left( x \right) = {\left\| {V\left( x \right) - V\left( {x + l} \right)} \right\|^2}$,所以计算两邻域V(x)和V(y)(y=x+l)欧氏距离,在常量时间内完成时,算法复杂度降到O(MD2):
$ \begin{array}{*{20}{c}} {{{\left\| {V\left( x \right) - V\left( y \right)} \right\|}^2} = \frac{1}{{{d^2}}}\left[ {{S_l}\left( {{x_1} + {d_{\rm{r}}},{x_2} + {d_{\rm{r}}}} \right) + } \right.}\\ {{S_l}\left( {{x_1} - {d_{\rm{r}}} - 1,{x_2} - {d_{\rm{r}}} - 1} \right) - {S_l}\left( {{x_1} + {d_{\rm{r}}},} \right.}\\ {\left. {\left. {{x_2} - {d_{\rm{r}}} - 1} \right) - {S_l}\left( {{x_1} - {d_{\rm{r}}} - 1,{x_2} + {d_{\rm{r}}}} \right)} \right]} \end{array} $
(10) 式中, l为两邻域平移量,Sl(x)为区域积分图像,sl(z)为区域点的积分图像,dr表示邻域窗口半径,x1为行像素值,x2为列像素值,d2为邻域窗口像素点个数。
-
(1) 采用SWT变换分解ESPI条纹图,分解成近似图像; (2)采用NSCT变换将近似图像分解成低频图像和高频细节图像; (3)选取低频图像和最高尺度高频细节图像; (4)对低频图像进行非线性扩散滤波,和高频细节图像进行快速非局部均值滤波,最后进行小波图像低频最大值融合。
图 1为选取ESPI条纹图基于本文算法步骤处理图像。本例中采用SWT分解中sym小波基单层分解,NSCT采用maxflat非下采样塔式分解和dmaxflat7非下采样方向滤波器组,尺度为2,方向为1,向量形式为[0, 0]。对低频图像采用基于偏微分方程非线性扩散滤波,采用梯度阈值k=40;迭代次数p=40的3次非线性扩散滤波。对高频细节图像采用基于快速非局部均值滤波,邻域窗口半径dr=2;搜索窗口半径Dr=5;控制指数函数衰减的去噪参量h=10,最后采用小波低频最大值融合图像。
-
滤波图像在视觉评价效果的基础上,还需要参量定量评价,本文中采用散斑指数Is和峰值信噪比(peak signal-to-noise-ratio, PSNR)PSNR定量评价指标。
Is表示为:
$ {I_{\rm{s}}} = \frac{1}{{ab}}\sum\limits_i^a {\sum\limits_o^b {\frac{{\sigma \left( {i,o} \right)}}{{\mu \left( {i,o} \right)}}} } $
(11) 式中, a和b表示图像大小,μ(i, o)和σ(i, o)分别为像素点(i, o)灰度值平均值和标准差。Is表示噪声抑制能力,参量越小,抑制能力越强。
PSNR表示为:
$ {P_{{\rm{SNR}}}} = 10\lg \frac{{{{255}^2}}}{{\frac{1}{{ab}}\sum\limits_{i = 1}^a {\sum\limits_{o = 1}^b {\left[ {I\left( {i,o} \right) - {I_0}\left( {i,o} \right)} \right]} } }} $
(12) 式中, I0(i, o)和I(i, o)分别为像素点(i, o)前后图像灰度值。PSNR表征图像信息保持能力,PSNR越大,表示图像信息保持越完整。
-
为验证本文中算法有效性,分别对计算机模拟和实验ESPI条纹图,进行小波软阈值、curvelet硬阈值和快速傅里叶变换滤波3种滤波算法对比。在计算机模拟图添加随机散斑噪声。对计算机模拟图和实验图采用小波基为sym,进行单层分解,NSCT采用maxflat非下采样塔式分解和dmaxflat7非下采样方向滤波器组。模拟图像大小a和b分别为500,如图 2所示,ESPI条纹图相位公式为:
$ \begin{array}{*{20}{c}} {\phi \left( {i,o} \right) = 30\exp \left[ { - \frac{{{{\left( {i - 0.1a} \right)}^2} + {{\left( {o - 0.2b} \right)}^2}}}{{40000}}} \right] + }\\ {20{{\left( {\frac{{o - 0.5a}}{{200}}} \right)}^2} + }\\ {50\exp \left[ { - \frac{{{{\left( {i - 0.2a} \right)}^2} + {{\left( {o - 0.2b} \right)}^2}}}{{30000}}} \right]} \end{array} $
(13) 对低频图像采用3次非线性扩散滤波,对高频图像采用2次快速非局部均值滤波,最后完成小波低频最大值图像融合。
如图 3所示,ESPI条纹图相位公式为:
$ \begin{array}{*{20}{c}} {\phi \left( {i,o} \right) = 90\exp \left[ { - \frac{{{{\left( {i - 0.2a} \right)}^2} + {{\left( {o - 0.5b} \right)}^2}}}{{30000}}} \right] - }\\ {70\exp \left[ { - \frac{{{{\left( {i - 0.45a} \right)}^2} + {{\left( {o - 0.75b} \right)}^2}}}{{18000}}} \right] - }\\ {60\exp \left[ { - \frac{{{{\left( {i - 0.8a} \right)}^2} + {{\left( {o - 0.25b} \right)}^2}}}{{10000}}} \right]} \end{array} $
(14) 对低频图像采用3次非线性扩散滤波,对高频图像采用2次快速非局部均值滤波,最后完成小波低频最大值图像融合。
图 4为实验得到ESPI减条纹图像,它是周边固支的圆盘受集中力作用。对低频图像采用1次非线性扩散滤波,对高频图像采用1次快速非局部均值滤波,最后完成小波低频最大值图像融合。
从表 1可以看出,图 2和图 3中,本文中算法PSNR值分别为11.5070dB和8.2018dB, 相较其它3种算法,PSNR值最大,对细节保护能力最强。从表 2可以看出,对于模拟和实验图,本文中算法的散斑指数为0.41121,0.38043,0.35362,散斑指数最小,噪声点存留较少。从视觉上来看,curvelet方法有明显“交错划痕”现象,小波方法有明显噪声,快速傅里叶变换滤波会造成较密条纹处的模糊,本文中算法处理图像效果较清晰。综上所述,在计算机模拟和实验条纹图里,本文中算法从综合评价上来看,滤波效果最佳,有明显优势。
Table 1. Comparison of parameters PSNR of simulated ESPI stripe patterns
Table 2. Comparison of Is parameter of simulation and experiment ESPI stripe patterns
改进非下采样轮廓波在散斑条纹中的滤波处理
Improvement of filtering processing of nonsubsampled contourlet transform in speckle stripes
-
摘要: 为了改进传统多尺度变换滤波在电子散斑干涉(ESPI)条纹图中去噪效果和边缘细节保护不理想问题,提出改进非下采样轮廓波(NSCT)滤波算法。采用离散平稳小波变换和NSCT变换模型,联合非线性扩散和改进的快速非局部均值滤波算法,进行了理论分析和实验验证,取得了将本文中算法应用于模拟和实验ESPI条纹图滤波效果定量分析的数据。结果表明,本文中的算法在模拟ESPI条纹图和实验图相比其它算法散斑指数最小分别为0.41121,0.38043,0.35362,对应峰值信噪比最大;该算法在提升去噪能力的同时,能够更好地恢复条纹细节信息。研究结果为以后应用多尺度变换滤波在ESPI条纹图打下了基础。Abstract: In order to improve the denoising effect and the edge detail protection of traditional multi-scale transform filtering in the fringe pattern of electronic speckle pattern interferometry (ESPI), an improved nonsubsampling contourlet transform (NSCT) filtering algorithm was proposed.Discrete stationary wavelet transform and NSCT model were adopted with nonlinear diffusion and improved fast non local mean filtering algorithm.The algorithm was applied to the quantitative analysis of simulated and experimental ESPI fringe pattern filtering results.After theoretical analysis and experimental verification, the results show that, the minimum speckle index of the algorithm in this paper is 0.41121, 0.38043 and 0.35362 respectively in simulated and experimental ESPI fringe pattern.The peak signal-to-noise ratio is the largest.This algorithm improves the ability of denoising and recovers the details of stripes.The results lay the foundation for the application of multi-scale transform filtering in ESPI fringe pattern.
-
Table 1. Comparison of parameters PSNR of simulated ESPI stripe patterns
Table 2. Comparison of Is parameter of simulation and experiment ESPI stripe patterns
-
[1] XU W J, TANG Ch.Two parabolic-hyperbolic oriented partial diffe-rential equations for denoising in electronic speckle pattern interferometry fringes[J].Applied Optics, 2015, 54(15):4720-4726. doi: 10.1364/AO.54.004720 [2] LI B Y, TANG Ch. General filtering method for electronic speckle pa-ttern interferometry fringe images with various densities based on variational image decomposition[J].Applied Optics, 2017, 56(16):4843-4853. doi: 10.1364/AO.56.004843 [3] ZHOU Q L, TANG Ch, WANG L L. A daptive oriented PDEs filter-ing methods based on new controlling speedfunction for discontinuous optical fringe patterns[J].Optics and Lasers in Engineering, 2018, 100:111-117. doi: 10.1016/j.optlaseng.2017.07.018 [4] RONG W X, QIAN X F, LIU B, et al. Algorithm of in-plane displacement measured by speckle photography based on phase of Fourier transform[J].Laser Technology, 2017, 41(4):473-478(in Chin-ese). [5] WANG M, WEN Y W, CHEN Zh B.Semi-smooth Newton method for total variation noise removal[J].Laser Technology, 2017, 41(2):289-295(in Chinese). [6] GU G Q, WANG K F, YAN X J.Electron speckle interference image processing based on homomorphic filtering[J].Laser Technology, 2010, 34(6):750-752(in Chinese). [7] LANG Ch P, YANG R H. High quality infrared video image restoration algorithm based on the improved ridgelet transform[J].Laser Technology, 2015, 39(2):247-251(in Chinese). [8] LI X, ZHANG W M.Study on speckle noise filter in the ESPI image based on curvelet transform[J]. Laser & Infrared, 2011, 41(9):1045-1048(in Chinese). [9] REN H W. Key techniques of phase retrieval for ESPI and experimental research on dynamic thermal deformation of circuit system[D].Tianjin: Tianjin University, 2013: 47-51(in Chinese). [10] XU W J, TANG Ch. Combination of oriented partial differential equation and shearlet transform for denoising in electronic speckle pattern interferometry fringe patterns[J].Applied Optics, 2017, 56(10):2843-2850. doi: 10.1364/AO.56.002843 [11] WU Y Q, YIN J, CAO Zh Q. Fusion of infrared thermal wave image based on nonsubsampledshearlet transform and weighted nonnegative matrix factorization[J].Acta Photonica Sinica, 2014, 43(10):1010001(in Chinese). doi: 10.3788/gzxb [12] LIU L, JIA Z, YANG J. A medical image enhancement method using adaptive thresholding in NSCT domain combined unsharp masking[J]. Medical Device & Sensors, 2015, 25(3):199-205. [13] LI C J, YANG H X, CAI Y Y, et al. Image denoising algorithm based on nonsubsampled contourlet transform and bilateral filtering[C]//2015 International Conference on Computer Information Systems and Industrial Applications.Bangkok: Atlantis Press, 2015: 666-669. [14] TIAN X, JIAO L, GUO K. An affinity-based algorithm in nonsubsampled contourlet transform domain:application to synthetic aperture radar image denoising[J].Journal of Signal Processing Systems, 2016, 83(3):373-388. doi: 10.1007/s11265-015-1024-2 [15] CHEN P, ZHAN Y, JIA Z, et al. Remote sensing image change detection based on NSCT-HMT model and its application[J]. Sensors, 2017, 17(6):1-15. doi: 10.1109/JSEN.2017.2656005 [16] GUO H, LI X J.Non-subsampled contourlet transform algorithm to promote detail information capturing ability[J].Journal of Image and Graphics, 2012, 17(8):920-922(in Chinese). [17] STARCK J L, CANDES E J, DONOHO D L, et al. The curvelet transform for image denoising[J]. IEEE Transactions on Image Processing, 2002, 11(6):670-684. doi: 10.1109/TIP.2002.1014998 [18] LIU Y, LI Y B, LI A. NSCT underwater image threshold denoising algorithm based on double filtering[J].Journal of Harbin Engineering University, 2013, 34(2):252-255(in Chinese). [19] FROMENT J. Parameter-free fast pixelwise non-local means denoising[J]. Image Processing on Line, 2014(4):300-326. [20] PERONA P, MALIK J. Scale-space and edge detection using anisotropic diffusion[J].IEEE Transactions on Pattern Analysis and Machine Intelligence, 1990, 12(7):629-639. doi: 10.1109/34.56205 [21] ALVAREZ L, LIONS P L, MOREL J M.Image selective smoothing and edge detection by nonlinear diffusion[J].Siam Journal on Numerical Analysis, 1992, 29(3):845-866. doi: 10.1137/0729052 [22] BUADES A, COLL B, MOREL J M. A non-local algorithm for image denoising[C]//2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition.San Diego, USA: IEEE Xplore, 2005, 2(7): 60-65.