基于FitzHugh-Nagumo反應(yīng)擴(kuò)散方程組的木材紋理圖像處理

牛 蕾，隋振璋，張春蕊

(東北林業(yè)大學(xué))

0 引言

因?yàn)楦鞣N木材的成長(zhǎng)環(huán)境不同，分子結(jié)構(gòu)不同，以及傷害的侵襲方式不同，其不同的切面總是呈現(xiàn)出不同的紋理模式.這種紋理模式與木材自身的結(jié)構(gòu)屬性息息相關(guān)，為了從有缺陷的圖像中獲得更多關(guān)于木材紋理的信息，對(duì)木材紋理的圖像處理是非常必要的.木材紋理的邊緣檢測(cè)幫助我們了解樹(shù)木成長(zhǎng)信息，判斷樹(shù)木的缺陷類(lèi)型和它的成長(zhǎng)趨勢(shì).

反應(yīng)擴(kuò)散系統(tǒng)是經(jīng)典的非線(xiàn)性動(dòng)力學(xué)行為，被廣泛的用來(lái)描述化學(xué)、地震、氣象等中周期性變化的現(xiàn)象.隨著人們對(duì)反應(yīng)擴(kuò)散系統(tǒng)深入的探索和學(xué)習(xí)，發(fā)現(xiàn)反應(yīng)擴(kuò)散系統(tǒng)可以實(shí)現(xiàn)對(duì)圖像的處理，這就為我們研究木材紋理提供了便利和工具.該文提供的圖像處理算法是建立在FitzHugh-Nagumo反應(yīng)擴(kuò)散方程組上的，與以往的邊緣檢測(cè)算法不同的是，它具有靈活的自組織時(shí)空反應(yīng)［1］.Kuhnert 等用光敏 Belousov-Zhabotinsky反應(yīng)擴(kuò)散系統(tǒng)實(shí)現(xiàn)了圖像增強(qiáng)和圖像分割［2-3］;Kurata等用常微分的思想分析了 FitzHugh-Nagumo反應(yīng)擴(kuò)散方程及其穩(wěn)定條件［4］，為邊緣檢測(cè)和圖像分割奠定了理論基礎(chǔ);Atsushi Nomura等分析了固定值的FitzHugh-Nagumo反應(yīng)擴(kuò)散方程對(duì)二值圖像的邊緣檢測(cè)［5］.該文是在文獻(xiàn)［4-5］的基礎(chǔ)上，改變擴(kuò)散系數(shù)使得圖像梯度越大擴(kuò)散系數(shù)越小，來(lái)完整的保留木材紋理邊緣.

1 FitzHugh-Nagumo反應(yīng)擴(kuò)散系統(tǒng)的木材紋理圖像處理基礎(chǔ)

1.1 FitzHugh-Nagumo方程的動(dòng)力學(xué)性質(zhì)

擴(kuò)散方程描述的是在一定時(shí)間內(nèi)有多少物質(zhì)或者熱量在空間中進(jìn)行了傳播，方程可以表示為

其中，u(x，t)代表物質(zhì)或者熱量的分布，并且u的分布是由空間向量x∈Rn和時(shí)間t定義的;D是擴(kuò)散系數(shù)，它代表物質(zhì)或者熱量的擴(kuò)散能力大小，其正負(fù)代表擴(kuò)散方向;Δ2u是Laplacian算子.設(shè)u0(x)是u(x，t)的初始值，其通過(guò)與方程(1)所生成的掩膜做卷積操作可以濾除噪聲增強(qiáng)圖像［6］;經(jīng)過(guò)一定的時(shí)間，空間分布U0(x)在擴(kuò)散系數(shù)D的作用下，生成了輸出分布:

FitzHugh-Nagumo反應(yīng)擴(kuò)散方程組由如下方程表示:

式中變量u和v隨空間向量x∈Rn和時(shí)間t變化而變化，具有初始狀態(tài)分別為u(x，t=0)=U0(x)，v(x，t=0)=V0(x).Du和Dv是擴(kuò)散系數(shù).而方程組反應(yīng)項(xiàng)分別為:

式中a和b分別為常數(shù)，0＜ε?1是一個(gè)極小的正常數(shù).當(dāng)擴(kuò)散系數(shù)為0時(shí)，得到如下的常微分方程組(3):

利用動(dòng)力系統(tǒng)理論，可以得到變量u和v的變化趨勢(shì)(如圖1所示):

圖1 變量u和v的變化趨勢(shì)

圖1中，du/dt=0用實(shí)線(xiàn)表示，dv/dt=0用虛線(xiàn)表示，隨著時(shí)間的推移，用箭頭表示解(u，v)的趨近方向;其中，區(qū)域Ⅰ:du/dt＞0且dv/dt＞0;區(qū)域Ⅱ:du/dt＜0且dv/dt＞0;區(qū)域Ⅲ:du/dt＜0且dv/dt＜0;區(qū)域Ⅳ:du/dt＞0且dv/dt＜0.在微分動(dòng)力學(xué)［7］中，這兩條線(xiàn)的交點(diǎn)(0，0)是一個(gè)穩(wěn)定平衡點(diǎn)，對(duì)于解(u，v)，從平衡點(diǎn)的某一鄰域的任意初值出發(fā)，解(u，v)最終都會(huì)收斂到此平衡點(diǎn).把u具有較大值的區(qū)域稱(chēng)為興奮區(qū)域;而把起點(diǎn)及其鄰域區(qū)域稱(chēng)為休眠區(qū)域.

從圖1可以看到，反應(yīng)項(xiàng)中的參數(shù)a是初始條件的臨界值.設(shè)解(u，v)的初始值為(u0，v0)，則當(dāng)u0＞a且v0≤0時(shí)，根據(jù)圖像中解(u，v)的走向，解(u，v)進(jìn)入興奮區(qū)域;而當(dāng)u0＜a且v0≥0時(shí)，解(u，v)快速進(jìn)入休眠區(qū)域.

利用系統(tǒng)的這個(gè)性質(zhì)，當(dāng)一幅需要處理的圖像作為初始值輸入時(shí)，該反應(yīng)擴(kuò)散系統(tǒng)遍歷這幅圖像的所有像素點(diǎn)，在適當(dāng)?shù)呐R界值a和足夠的時(shí)間作用下，圖像像素點(diǎn)的灰度值劃分到圖1所示的兩個(gè)區(qū)域，這樣就擴(kuò)大了不同像素點(diǎn)的灰度值.

1.2 FitzHugh-Nagumo反應(yīng)擴(kuò)散方程的數(shù)值處理

對(duì)反應(yīng)擴(kuò)散方程(2)，采用九點(diǎn)差分中的顯示格式進(jìn)行求解［6，8］.空間變量(x，y)和時(shí)間變量t用有限差分進(jìn)行離散，空間步長(zhǎng)和時(shí)間步長(zhǎng)分別為 δh，δt，其中:i=［x/δh］，j=［y/δh］，k=［t/δt］，i，j和k是離散空間和時(shí)間的指數(shù)，那么變量u(x，y，t)的離散格式是uki，j，且

定義?u/?t，?u/?x，?u/?y的離散格式分別為Δtu，Δxu，Δyu，且

由(4)和(5)，方程組(2)的第一個(gè)方程變?yōu)?

式中，r=Duδt/δh2.

方程組(2)的Neumann邊界條件的離散式為:

在圖像處理中，一幅圖像離散化后的數(shù)字圖像大小為M×N，即其有M行和N列，取i=0，1，…，M-1，j=0，1，…，N-1，那么在初值(U0，V0)和離散邊界條件(7)基礎(chǔ)上，通過(guò)用式(6)的迭代運(yùn)算，最終獲得隨時(shí)間演化的解(u，v)，即處理后的圖像.

2 對(duì)木材紋理圖像處理算法的實(shí)現(xiàn)和結(jié)果

2.1 FitzHugh-Nagumo離散模型圖像增強(qiáng)算法的實(shí)現(xiàn)

按照以上的理論分析，對(duì)于一個(gè)不帶擴(kuò)散項(xiàng)的常微系統(tǒng)(3)，或者r=0的離散系統(tǒng)(6)模型來(lái)說(shuō)，該模型可以用來(lái)增強(qiáng)圖像，拉伸模糊圖像的對(duì)比度.當(dāng)a為一個(gè)固定值時(shí)，該模型可以用來(lái)處理二值圖像;那么當(dāng)一幅圖像為灰度圖像時(shí)，可以借鑒文獻(xiàn)［9］的思想，利用擴(kuò)散項(xiàng)在時(shí)間上的迭代，引入變量

利用a(x)取代固定的a值.得到改進(jìn)的圖像增強(qiáng)算法.該算法步驟如下:

步驟1:取亮度圖像u(x，t=0)并歸一化到［0，1］;

步驟2:準(zhǔn)備好初始條件U0(x)=u(x，t=0)和V0(x)=0;

步驟3:計(jì)算a(x);

步驟4:把a(bǔ)(x)代入r=0的離散系統(tǒng)(6)中，通過(guò)迭代運(yùn)算，最后得到邊緣檢測(cè)圖.

2.2 FitzHugh-Nagumo離散模型圖像增強(qiáng)算法的處理結(jié)果及分析

按照2.1算法步驟，對(duì)不清晰的有活節(jié)的木材亮度圖像(圖2(a)，181×188)用Matlab7.0軟件進(jìn)行增強(qiáng)處理，這里選擇b=1，ε=10-3，Du=0.01，Dv=0.25，δh=1/6，δt=10-3，a(x)=A(U0;1，2，0.3).截取不同演化次數(shù)即在時(shí)域上的演化圖像及其對(duì)應(yīng)的灰度直方圖，結(jié)果如圖2所示.

圖2 不同時(shí)刻得到的木材紋理圖像以及對(duì)應(yīng)的灰度直方圖

由圖2可知該算法在演化過(guò)程中對(duì)原活節(jié)木材紋理圖像進(jìn)行了有效的增強(qiáng)，其對(duì)應(yīng)的灰度直方圖也從原來(lái)呈現(xiàn)的尖銳狀變化到平滑、均勻狀.

該文利用Michelson對(duì)比度公式來(lái)計(jì)算圖2中對(duì)應(yīng)圖像的對(duì)比度

式中，Ⅰmax和Ⅰmin分別表示圖像中最亮的亮度和最暗的亮度.相應(yīng)的圖像對(duì)比度為(a)CM=1，(b)CM=1.0043，(c)CM=1.0147，(d)CM=1.0250.這個(gè)結(jié)果表明，圖像對(duì)比度得到拉伸，紋理圖像得到增強(qiáng).演化結(jié)果與我們?cè)?.2中的理論分析相一致:紋理圖像在很短的時(shí)間內(nèi)得到增強(qiáng)以后，圖像灰度值向穩(wěn)定值0趨近，在圖像亮度上呈暗亮度.

2.3 FitzHugh-Nagumo反應(yīng)擴(kuò)散系統(tǒng)對(duì)木材紋理圖像的邊緣檢測(cè)算法

對(duì)于帶有擴(kuò)散項(xiàng)的FitzHugh-Nagumo反應(yīng)擴(kuò)散系統(tǒng)模型，參數(shù)a的取法，對(duì)檢測(cè)一幅圖像的邊緣檢測(cè)同樣重要;當(dāng)a為固定常數(shù)時(shí)，該模型可以實(shí)現(xiàn)二值圖像的邊緣檢測(cè)［5］;當(dāng)參數(shù)a取a(x)=A(U0;D，T)［9］來(lái)代替固定的a值，該模型可以實(shí)現(xiàn)亮度圖像的邊緣檢測(cè);但是這種做法也不盡完美，在該模型作用下圖像產(chǎn)生真假雙邊緣［10］，且隨著反應(yīng)擴(kuò)散系數(shù)的增加，雙邊緣之間的距離拉長(zhǎng)，甚至?xí)a(chǎn)生覆蓋邊緣的現(xiàn)象，使圖像邊緣模糊.

針對(duì)這一缺陷，該文通過(guò)添加一個(gè)調(diào)整項(xiàng)［11］來(lái)改進(jìn)該算法.這個(gè)調(diào)整項(xiàng)是一個(gè)關(guān)于梯度的單調(diào)遞減的非負(fù)函數(shù)，取值范圍為(0，1］，它能靈活的調(diào)整擴(kuò)散系數(shù)，進(jìn)而調(diào)整雙邊緣的距離，該調(diào)整項(xiàng)為:

式中，p為常數(shù)，它主要用來(lái)控制h的下降速率，

其中，借用Sobel算子的兩組3×3的矩陣，分別將之與圖像U0作平面卷積，即可分別得出橫向及縱向的亮度差分近似值Gx及Gy:

同時(shí)，為了得到邊緣的穩(wěn)定的穩(wěn)態(tài)解，該反應(yīng)擴(kuò)散系統(tǒng)必須是在強(qiáng)抑制Du?Dv作用下.該文算法如下:Δ Δ

具體步驟為:

步驟1:取亮度圖像u(x，t=0)并歸一化到［0，1］;

步驟2:準(zhǔn)備好初始圖像U0(x)=u(x，t=0)和V0(x)=0;

步驟3:用式(8)和(9)計(jì)算得到a(x);

步驟4:把a(bǔ)(x)代入到反應(yīng)擴(kuò)散離散系統(tǒng)(10)中，通過(guò)式(10)和(11)的離散形式的迭代運(yùn)算，最后得到邊緣檢測(cè)圖.

2.4 FitzHugh-Nagumo反應(yīng)擴(kuò)散系統(tǒng)對(duì)木材紋理邊緣檢測(cè)結(jié)果及分析

按照2.3分析的算法步驟，對(duì)有死節(jié)的木材亮度圖像(圖3(a)，183×190)用Matlab 7.0軟件進(jìn)行紋理邊緣檢測(cè)得到圖3.

圖3

對(duì)比圖3(a)(b)(c)三幅紋理圖像，圖3(b)是對(duì)圖(a)用經(jīng)典的Sobel算子得到的邊緣檢測(cè)圖，很明顯，其紋理邊緣的線(xiàn)條連續(xù)，但是漏檢線(xiàn)條太多;圖3(c)是用先前的算法得到的邊緣檢測(cè)圖，很明顯突出了死節(jié)內(nèi)部紋理，周?chē)牟糠旨y理走勢(shì)模糊，死節(jié)外部紋理線(xiàn)條呈間斷、點(diǎn)狀;圖3(d)是用在先前算法基礎(chǔ)上得到的本文算法得到的邊緣檢測(cè)圖，與圖3(c)比較，死節(jié)內(nèi)部紋理更完整，外部紋理走勢(shì)更清晰，線(xiàn)條連續(xù)性有提高;而與圖3(b)比較，死節(jié)內(nèi)部紋理清晰，得到的邊緣覆蓋較大，缺點(diǎn)就是線(xiàn)條不連續(xù).

總體而言，該文邊緣檢測(cè)算法對(duì)木材紋理的邊緣檢測(cè)在一定程度上提高了完整性和清晰性，方便對(duì)木材紋理缺陷進(jìn)行檢測(cè);但是在連續(xù)性上還有待提高.

［1］ Zaikin A N，Zhabotinsky A M.Concentration wave propagation in two-dimensional liquid-phase self-oscillating system［J］.Nature，1970(225):535–537.

［2］ Kuhnert L.A new optical photochemical memory device in a light sensitive chemical active medium ［J］.Nature，1986(319):393-394.

［3］ Kuhnert L，Agladze K I，Krinsky V I.Image processing using light-sensitive chemical waves［J］.Nature，1989(337):244-247.

［4］ Kurata N，Kitahata H，Mahara H，et al.Stationary pattern formation in a discrete excitable system with strong inhibitory coupling［J］.Physical Review E，2009(79):056-203.

［5］ Nomura A，Ichikawa M，Miike H，et al.Realizing visual functions with the reaction-diffusion mechanism［J］.Journal of the Physical Society of Japan，2003(72):2385-2395.

［6］ 岡薩雷斯，伍茲著，阮秋琦，等譯，數(shù)字圖像處理［M］.北京:電子工業(yè)出版社.2011.

［7］ 陸啟韶，彭臨平，楊卓琴.常微分方程與動(dòng)力系統(tǒng)［M］.北京:北京航空航天大學(xué)出版社，2009.

［8］ 李榮華，劉播.微分方程數(shù)值解法［M］.北京:高等教育出版社，2009.

［9］ Nomura A，Ichikawa M，Sianipar R H，et al.Edge detection with reaction-diffusion equations having a local average threshold［J］.Pattern Recognition and Image Analysis，2008(18):289-299.

［10］Nomura A，Ichikawa M，Okada K，et al.Edge detection algorithm inspired by pattern formation processes of reaction-diffusion systems［J］.International Journal Of Circuits，Systems And Signal Processing，2011(5):105-114.

［11］荊萍.基于PDE的醫(yī)學(xué)圖像增強(qiáng)技術(shù)研究［J］.2011:40-41.