本發(fā)明屬于光學(xué)檢測(cè)和圖像處理技術(shù)領(lǐng)域，涉及一種用于條紋投影輪廓術(shù)中條紋圖的濾波的基于擬小波離散格式的方向偏微分方程濾波方法。

背景技術(shù)：

條紋投影輪廓術(shù)(Fringe Projection Profilometry,FPP)是近年來(lái)十分受歡迎的一種三維形貌測(cè)量方法，因?yàn)樗哂懈哽`敏度、非接觸、全場(chǎng)觀(guān)測(cè)等優(yōu)點(diǎn)，所以被廣泛應(yīng)用于科研和實(shí)際工程中[1-3]。對(duì)條紋投影圖像的處理分析是獲取待測(cè)物體形貌信息的主要手段之一，但由于條紋投影圖像中總是伴隨著強(qiáng)烈的噪聲，導(dǎo)致無(wú)法直接利用這些圖像正確的獲取相位信息，所以必須對(duì)FPP條紋圖像進(jìn)行濾波，提高圖像的對(duì)比度，去除噪聲的影響。偏微分方程(Partial differential equation,簡(jiǎn)稱(chēng)PDE)是一種方案靈活、處理效果良好的圖像濾波技術(shù)[4,5]，尤其是方向偏微分方程，只沿著條紋方向[6,7]進(jìn)行濾波，適合于條紋圖像的濾波。在使用PDE時(shí)，一個(gè)重要問(wèn)題是如何將PDE模型進(jìn)行離散化，方程的離散化程度直接影響到濾波效果以及計(jì)算效率。

參考文獻(xiàn)

[1]楊國(guó)光,近代光學(xué)測(cè)試技術(shù).杭州:浙江大學(xué)出版社,1997。

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

[3]王開(kāi)福,現(xiàn)代光測(cè)及其圖像處理.北京:科學(xué)出版社,2012:100～110。

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

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

[6]S.Chikkerur,A.N.Cartwright,and V.Govindaraju,“Fingerprint enhancement using STFT analysis,”P(pán)attern Recogn.198-211(2007)。

[7]L.Hong,Y.Wan,and A.Jain,“Fingerprint image enhancement:algorithm and 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 for Two-dimensional Reaction-diffusion Equations.J.Comp.Phys.1998,139:9-G～12。

技術(shù)實(shí)現(xiàn)要素：

對(duì)偏微分方程離散時(shí)，差分格式的選擇是非常重要的，之前已有的偏微分方程離散化用到的差分格式主要以中心差分為主，為克服現(xiàn)有技術(shù)的不足，本發(fā)明旨在提出一種新的基于擬小波方法的差分格式，對(duì)數(shù)值求解偏微分方程提供新方法，使得偏微分方程濾波效果得到了提升。本發(fā)明采用的技術(shù)方案是，用于條紋投影圖像的方向偏微分方程濾波方法，步驟如 下：

步驟1，輸入一幅帶有噪聲的條紋投影圖u；

步驟2，將條紋投影圖u進(jìn)行離散化，假設(shè)條紋投影圖u的大小為M×N，u_i,j為(i,j)點(diǎn)處的灰度值，1≤i≤M，1≤j≤N，時(shí)間步長(zhǎng)為Δt，在方程的演化過(guò)程中，t_n＝nΔt時(shí)刻的演化圖像表示為，構(gòu)造方向偏微分方程的離散格式為：

其中，θ_i,j表示條紋方向與x軸方向的夾角，u_x、u_y分別表示圖像u沿著x軸、y軸方向的一階偏導(dǎo)數(shù)，u_xx、u_yy分別表示圖像u沿著x軸、y軸方向的二階偏導(dǎo)數(shù)；表示對(duì)圖像進(jìn)行離散化，表示對(duì)圖像離散化后的差分格式；

步驟3：計(jì)算條紋方向與x軸方向的夾角θ_i,j；

步驟4：自適應(yīng)得到時(shí)間步長(zhǎng)Δt；

步驟5：自適應(yīng)得到迭代次數(shù)Nc；

步驟6：對(duì)于每次迭代，按照公式：

其中的擬小波離散格式：

其中E＝exp(-x²/(2σ²)),S＝sin(πx/Δ),C＝cos(πx/Δ)，x_i,j-x_i,k表示像素點(diǎn)之間做減法運(yùn)算，Δ表示網(wǎng)格大小，σ表示窗口大小參數(shù)，求出圖像u每個(gè)像素的一階偏導(dǎo)數(shù)u_x，u_y以及二階偏導(dǎo)數(shù)u_xx，u_xy和u_yy；

步驟7：基于步驟2中的方向偏微分方程的離散格式，求出圖像u每個(gè)像素的數(shù)值解

步驟8：重復(fù)步驟6和步驟7，直到達(dá)到設(shè)置的最大迭代次數(shù)Nc停止迭代，此時(shí)的數(shù)值 解即為濾波圖像。

步驟4自適應(yīng)得到時(shí)間步長(zhǎng)Δt包括以下步驟：

步驟4-1：輸入一幅電子散斑干涉條紋圖像u；

步驟4-2：給定初值，設(shè)置迭代初值u0為原始圖像I，迭代時(shí)間步長(zhǎng)Δt_k＝1，初始循環(huán)次數(shù)k＝0；

步驟4-3：計(jì)算兩次的濾波結(jié)果u(Δtk)和u(2Δtk)；

步驟4-4：判斷第二次濾波后的噪聲與信號(hào)的相關(guān)系數(shù)corr(u(0)-u(2Δtk),u(2Δtk))是否小于第一次濾波后的噪聲與信號(hào)的相關(guān)系數(shù)corr(u(0)-u(Δtk),u(Δtk))，若是，則執(zhí)行步驟4-5；否則設(shè)置k＝k+1，時(shí)間步長(zhǎng)在上次的基礎(chǔ)上減0.05，即賦值Δtk+1＝Δtk-0.05，重復(fù)步驟4-3和步驟4-4；

步驟4-5：輸出此時(shí)的Δtk。

步驟5自適應(yīng)得到迭代次數(shù)Nc包括以下步驟：

步驟5-1：輸入一幅電子散斑干涉條紋圖像u；

步驟5-2：固定時(shí)間步長(zhǎng)Δt為上述自適應(yīng)得到的Δt，設(shè)置迭代初值u0為初始圖像I，迭代次數(shù)n＝1；

步驟5-3：計(jì)算第一次迭代的濾波結(jié)果

步驟5-4：自n＝1開(kāi)始，計(jì)算并保存當(dāng)前濾波結(jié)果和上一次的濾波結(jié)果

步驟5-5：判斷是否同時(shí)滿(mǎn)足下述兩個(gè)條件，

條件1是當(dāng)前濾波后噪聲與信號(hào)的相關(guān)系數(shù)corr(u(0)-u(2Δt),u(2Δt))大于上一次濾波后的噪聲與信號(hào)的相關(guān)系數(shù)corr(u(0)-u(Δt),u(Δt))，

條件2是散斑系數(shù)小于0.2；

若是，則執(zhí)行步驟5-6，否則設(shè)置n＝n+1，重復(fù)步驟5-4和步驟5-5；

步驟5-6：輸出Nc＝n。

本發(fā)明的特點(diǎn)及有益效果是：

擬小波方法的離散格式相比于多點(diǎn)差分計(jì)算需要的選點(diǎn)更少，一般只需要15個(gè)點(diǎn)即可滿(mǎn)足計(jì)算精度，并且這些系數(shù)可以一次性計(jì)算出來(lái)，大大縮短計(jì)算時(shí)間。相對(duì)于中心差分格式相比，擬小波離散格式擁有更高的計(jì)算精度。

附圖說(shuō)明：

圖1噪聲方差為1.2時(shí)三種方法的濾波結(jié)果。圖中，(a)噪聲圖像，(b)SOOPDE，(c)WFF，(d)WTM。

具體實(shí)施方式

本發(fā)明基于擬小波離散的二階單方向PDE濾波方法的主要步驟：

步驟1，輸入一幅帶有噪聲的條紋投影圖u

步驟2，將圖像u進(jìn)行離散化，假設(shè)圖像u的大小為M×N，u_i,j(1≤i≤M,1≤j≤N)為(i,j)點(diǎn)處的灰度值，時(shí)間步長(zhǎng)為Δt，在方程的演化過(guò)程中，t_n＝nΔt時(shí)刻的演化圖像表示為， 構(gòu)造方向偏微分方程的離散格式為：

其中，θ_i,j表示條紋方向與x軸方向的夾角，u_x、u_y分別表示圖像u沿著x軸、y軸方向的一階偏導(dǎo)數(shù)，u_xx、u_yy分別表示圖像u沿著x軸、y軸方向的二階偏導(dǎo)數(shù)；

步驟3：計(jì)算條紋方向與x軸方向的夾角θ_i,j；

步驟4：自適應(yīng)得到時(shí)間步長(zhǎng)Δt；

步驟5：自適應(yīng)得到迭代次數(shù)N_c；

步驟6：對(duì)于每次迭代，按照公式：

其中的擬小波離散格式：

其中E＝exp(-x²/(2σ²)),S＝sin(πx/Δ),C＝cos(πx/Δ).求出圖像u每個(gè)像素的一階偏導(dǎo)數(shù)u_x，u_y以及二階偏導(dǎo)數(shù)u_xx，u_xy和u_yy；x_i,j-x_i,k表示像素點(diǎn)之間做減法運(yùn)算，Δ表示網(wǎng)格大小，σ表示窗口大小參數(shù)

步驟7：基于步驟2中的方向偏微分方程的離散格式，求出圖像u每個(gè)像素的數(shù)值解

步驟8：重復(fù)步驟6和步驟7，直到達(dá)到設(shè)置的最大迭代次數(shù)N_c停止迭代，此時(shí)的數(shù)值解即為濾波圖像。

步驟4自適應(yīng)得到時(shí)間步長(zhǎng)Δt包括以下步驟：

步驟4-1：輸入一幅電子散斑干涉條紋圖像u；

步驟4-2：給定初值，設(shè)置迭代初值u0為原始圖像I，迭代時(shí)間步長(zhǎng)Δt_k＝1，初始循環(huán)次數(shù)k＝0；

步驟4-3：計(jì)算兩次的濾波結(jié)果u(Δtk)和u(2Δtk)；

步驟4-4：判斷第二次濾波后的噪聲與信號(hào)的相關(guān)系數(shù)corr(u(0)-u(2Δtk),u(2Δtk))是否小于第一次濾波后的噪聲與信號(hào)的相關(guān)系數(shù)corr(u(0)-u(Δtk),u(Δtk))，若是，則執(zhí)行步驟4-5；否則設(shè)置k＝k+1，時(shí)間步長(zhǎng)在上次的基礎(chǔ)上減0.05，即賦值Δtk+1＝Δtk-0.05，重復(fù)步驟4-3和步驟4-4；

步驟4-5：輸出此時(shí)的Δtk。

步驟5自適應(yīng)得到迭代次數(shù)Nc包括以下步驟：

步驟5-1：輸入一幅電子散斑干涉條紋圖像u；

步驟5-2：固定時(shí)間步長(zhǎng)Δt為上述自適應(yīng)得到的Δt，設(shè)置迭代初值u0為初始圖像I，迭代次數(shù)n＝1；

步驟5-3：計(jì)算第一次迭代的濾波結(jié)果

步驟5-4：自n＝1開(kāi)始，計(jì)算并保存當(dāng)前濾波結(jié)果和上一次的濾波結(jié)果

步驟5-5：判斷是否同時(shí)滿(mǎn)足下述兩個(gè)條件，

條件1是當(dāng)前濾波后噪聲與信號(hào)的相關(guān)系數(shù)corr(u(0)-u(2Δt),u(2Δt))大于上一次濾波后的噪聲與信號(hào)的相關(guān)系數(shù)corr(u(0)-u(Δt),u(Δt))，

條件2是散斑系數(shù)小于0.2；

若是，則執(zhí)行步驟5-6，否則設(shè)置n＝n+1，重復(fù)步驟5-4和步驟5-5；

步驟5-6：輸出Nc＝n。

下面結(jié)合具體實(shí)施方式對(duì)本發(fā)明做進(jìn)一步詳細(xì)的介紹：

條紋投影輪廓術(shù)(Fringe Projection Profilometry,FPP)是近年來(lái)十分受歡迎的一種三維形貌測(cè)量方法，因?yàn)樗哂懈哽`敏度、非接觸、全場(chǎng)觀(guān)測(cè)等優(yōu)點(diǎn)，所以被廣泛應(yīng)用于科研和實(shí)際工程中。對(duì)條紋投影圖像的處理分析是獲取待測(cè)物體形貌信息的主要手段之一，但由于條紋投影圖像中總是伴隨著強(qiáng)烈的噪聲，導(dǎo)致無(wú)法直接利用這些圖像正確的獲取相位信息，所以必須對(duì)FPP條紋圖像進(jìn)行濾波，提高圖像的對(duì)比度，去除噪聲的影響。偏微分方程(Partial differential equation,簡(jiǎn)稱(chēng)PDE)是一種方案靈活、處理效果良好的圖像濾波技術(shù)，尤其是方向偏微分方程，只沿著條紋方向進(jìn)行濾波，適合于條紋圖像的濾波。在使用PDE時(shí)，一個(gè)重要問(wèn)題是如何將PDE模型進(jìn)行離散化，方程的離散化程度直接影響到濾波效果以及計(jì)算效率。本發(fā)明提出一種擬小波離散方法，對(duì)PDE濾波模型進(jìn)行離散化，同時(shí)將改進(jìn)的PDE濾波方法與窗口傅里葉方法和小波濾波方法進(jìn)行比較，驗(yàn)證本發(fā)明提出的方法的有效性。

小波函數(shù)指的是一組有有限能量的函數(shù),它們可以通過(guò)伸縮和平移來(lái)構(gòu)造生成：

公式(1)中的參數(shù)a控制函數(shù)伸縮，參數(shù)b控制函數(shù)平移。這樣通過(guò)調(diào)整尺度參數(shù)a，就可以分析函數(shù)本身的局部特性。而擬小波是對(duì)正交規(guī)范化尺度函數(shù)進(jìn)行正則化處理得到的，其中正交規(guī)范化尺度函數(shù)可以定義為：

任意一個(gè)屬于空間的函數(shù)f(x)都可以表示成下面的離散形式：

其中表示是在Hilbert空間的Paley-Wiener重構(gòu)核。這里{x_k}是定義在無(wú)窮區(qū)間(-∞,∞)的一組離散空間點(diǎn)的坐標(biāo)。為了實(shí)現(xiàn)計(jì)算，這里對(duì)a賦值，一般a＝π/ΔΔ為單元網(wǎng)格大小。

但是正交規(guī)范化尺度函數(shù)δ_Δ,σ(x)在坐標(biāo)空間上沒(méi)有明顯的局部特性，在文獻(xiàn)[8]中，Wei[8]等人建議對(duì)它們進(jìn)行正則化處理(regularization procedure)。即引入一個(gè)Gauss函數(shù)作為正則因子:

Gauss函數(shù)定義為：

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

其中σ為窗口大小參數(shù)。

Gauss正則化之后的尺度函數(shù)可以定義為：

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

Gauss正則化樣本尺度函數(shù)不再是精確的正交規(guī)范化小波尺度函數(shù)，因?yàn)樗粷M(mǎn)足正交規(guī)范化小波尺度函數(shù)的條件,因此，我們把Gauss正則化樣本尺度函數(shù)δ_Δ,σ(x)稱(chēng)為擬尺度函數(shù)，由擬尺度函數(shù)生成的小波就相應(yīng)地稱(chēng)為擬小波。

因此對(duì)于一個(gè)函數(shù)f(x)，通過(guò)δ_Δ,σ(x)可以被唯一重構(gòu)為：

由于上式求和運(yùn)算定義在無(wú)限區(qū)間(-∞,∞)，這在實(shí)際計(jì)算中是做不到的。因此，為了計(jì)算需要，有必要把計(jì)算域定義在有限區(qū)間。由于我們已經(jīng)引進(jìn)Gauss正則因子，它具有很好的快速衰減特性，即局域特性。在實(shí)際計(jì)算中，只需要在格點(diǎn)附近取2W(15左右)個(gè)計(jì)算點(diǎn)x_k，就可以達(dá)到機(jī)器精度。因此，上式可進(jìn)一步寫(xiě)成：

上式中的上標(biāo)(n)表示對(duì)空間坐標(biāo)x的n階導(dǎo)數(shù)，2W+1為計(jì)算帶寬，它通常遠(yuǎn)遠(yuǎn)小于整個(gè)計(jì)算域的網(wǎng)格點(diǎn)數(shù)。公式(8)中，如果△和σ給定，在整個(gè)計(jì)算域中則只有一個(gè)核出現(xiàn)，這對(duì)數(shù)值計(jì)算來(lái)說(shuō)是十分經(jīng)濟(jì)和有利的、我們把公式(8)稱(chēng)為擬小波數(shù)值離散格式。

應(yīng)用擬小波離散方法對(duì)二階偏微分方程進(jìn)行離散化：

其中偏微分方程濾波模型為：

將圖像大小為M×N，對(duì)于1<i<M,1<j<N的任意一點(diǎn)時(shí)刻n的二階導(dǎo)數(shù)為：把空間x坐標(biāo)均勻等分,單元網(wǎng)格大小記為Δx＝1/w(w是化分區(qū)間所采用的單元網(wǎng)格總數(shù))，網(wǎng)格點(diǎn)坐標(biāo)記為x_i＝i·Δx，(i＝0,1,2,…,M),于是有x_i-x_i+k＝-kΔx。

方程(5-25)的擬小波離散格式可以寫(xiě)成如下格式：

其中的擬小波離散格式：

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

這里[-W,+W]是計(jì)算帶寬,在計(jì)算中這意味著計(jì)算只在很小的帶寬上進(jìn)行,極大減小了數(shù)值計(jì)算量，并且這些系數(shù)計(jì)算只與單元網(wǎng)格尺度大小有關(guān),因此當(dāng)網(wǎng)格分布確定后,這些系數(shù)可一次性計(jì)算出來(lái),然后把它們保存起來(lái)供計(jì)算機(jī)調(diào)用,不需要重復(fù)計(jì)算。通過(guò)上述擬小波離散和最速下降法，可以實(shí)現(xiàn)對(duì)二階單方向偏微分方程的數(shù)值求解，從而得到圖像的濾波結(jié)果。

下面給出實(shí)驗(yàn)結(jié)果圖。

我們采用圖1-(a)模擬FPP條紋圖作為實(shí)驗(yàn)圖像，為了驗(yàn)證算法的濾波性能，采用服從均勻分布的隨機(jī)噪聲，幅值大小分別為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中自帶函數(shù)，其中，Re{}代表實(shí)部，條紋的載頻頻率為f₀＝1/16，背景強(qiáng)度a(x,y)為0.5φ(x,y)，調(diào)制強(qiáng)度b(x,y)為1，高斯隨機(jī)噪聲的方差為1.2。

從圖中可以清晰的看到這幅圖像含有大量噪聲，并且條紋密集。采用我們改進(jìn)的二階單方向偏微分方程(SOOPDE)方法，窗傅里葉濾波方法(WFT)，小波變換濾波(WTM)方法，對(duì)圖1-(a)做濾波處理。其中圖1-(b),(c),(d),表示為SOOPDE,WFT和WTM三種方法分別對(duì)圖1-(a)的濾波結(jié)果，在進(jìn)行處理時(shí)各種算法的最優(yōu)參數(shù)是通過(guò)反復(fù)試驗(yàn)得出來(lái)的并且參數(shù)設(shè)置展示在表1-1中。

由于模擬圖自身具有已知的理想背景和條紋，此處我們采用峰值信噪比(PSNR)對(duì)每一種模型的濾波結(jié)果進(jìn)行定量研究，其中峰值信噪比(PSNR)的定義形式為：

表1-2展示對(duì)于圖1-(a)三種方法處理獲得的濾波圖像的峰值信噪比以及濾波時(shí)間。三種方法對(duì)應(yīng)程序的計(jì)算環(huán)境為：MATLAB R2014b，Intel Core i3-2100@2.9GHz CPU，運(yùn)行內(nèi)存2G。

接下來(lái)對(duì)于不同噪聲等級(jí)下的模擬條紋圖進(jìn)行三種方法的去噪性能比較分析，結(jié)果如表1-3所示。

表1-1三種方法參數(shù)設(shè)置

從圖1，表1-1、表1-2和表1-3可以看出，窗傅里葉變換方法能獲得很好的濾波效果但是由于需要大量的參數(shù)設(shè)置，因此需要花費(fèi)大量的時(shí)間，不利于提高運(yùn)算效率。小波變換方 法花費(fèi)時(shí)間最少，但是濾波質(zhì)量與本文提出的方法有一定的差距。本文提出的基于擬小波離散格式的方向PDE濾波方法在得到高質(zhì)量濾波結(jié)果前提下，需要更少的迭代次數(shù)，這大大節(jié)省了計(jì)算時(shí)間。

表1-2三種方法的PSNR及其計(jì)算時(shí)間

表1-3不同噪聲等級(jí)下三種濾波方法的PSNR及其計(jì)算時(shí)間

盡管上面結(jié)合圖對(duì)本發(fā)明進(jìn)行了描述，但是本發(fā)明并不局限于上述的具體實(shí)施方式，上述的具體實(shí)施方式僅僅是示意性的，而不是限制性的，本領(lǐng)域的普通技術(shù)人員在本發(fā)明的啟示下，在不脫離本發(fā)明宗旨的情況下，還可以作出很多變形，這些均屬于本發(fā)明的保護(hù)之內(nèi)。