基于FitzHugh-Nagumo反应扩散方程组的木材纹理图像处理

牛 蕾，隋振璋，张春蕊

(东北林业大学)

0 引言

因为各种木材的成长环境不同，分子结构不同，以及伤害的侵袭方式不同，其不同的切面总是呈现出不同的纹理模式.这种纹理模式与木材自身的结构属性息息相关，为了从有缺陷的图像中获得更多关于木材纹理的信息，对木材纹理的图像处理是非常必要的.木材纹理的边缘检测帮助我们了解树木成长信息，判断树木的缺陷类型和它的成长趋势.

反应扩散系统是经典的非线性动力学行为，被广泛的用来描述化学、地震、气象等中周期性变化的现象.随着人们对反应扩散系统深入的探索和学习，发现反应扩散系统可以实现对图像的处理，这就为我们研究木材纹理提供了便利和工具.该文提供的图像处理算法是建立在FitzHugh-Nagumo反应扩散方程组上的，与以往的边缘检测算法不同的是，它具有灵活的自组织时空反应［1］.Kuhnert 等用光敏 Belousov-Zhabotinsky反应扩散系统实现了图像增强和图像分割［2-3］;Kurata等用常微分的思想分析了 FitzHugh-Nagumo反应扩散方程及其稳定条件［4］，为边缘检测和图像分割奠定了理论基础;Atsushi Nomura等分析了固定值的FitzHugh-Nagumo反应扩散方程对二值图像的边缘检测［5］.该文是在文献［4-5］的基础上，改变扩散系数使得图像梯度越大扩散系数越小，来完整的保留木材纹理边缘.

1 FitzHugh-Nagumo反应扩散系统的木材纹理图像处理基础

1.1 FitzHugh-Nagumo方程的动力学性质

扩散方程描述的是在一定时间内有多少物质或者热量在空间中进行了传播，方程可以表示为

其中，u(x，t)代表物质或者热量的分布，并且u的分布是由空间向量x∈Rn和时间t定义的;D是扩散系数，它代表物质或者热量的扩散能力大小，其正负代表扩散方向;Δ2u是Laplacian算子.设u0(x)是u(x，t)的初始值，其通过与方程(1)所生成的掩膜做卷积操作可以滤除噪声增强图像［6］;经过一定的时间，空间分布U0(x)在扩散系数D的作用下，生成了输出分布:

FitzHugh-Nagumo反应扩散方程组由如下方程表示:

式中变量u和v随空间向量x∈Rn和时间t变化而变化，具有初始状态分别为u(x，t=0)=U0(x)，v(x，t=0)=V0(x).Du和Dv是扩散系数.而方程组反应项分别为:

式中a和b分别为常数，0＜ε≪1是一个极小的正常数.当扩散系数为0时，得到如下的常微分方程组(3):

利用动力系统理论，可以得到变量u和v的变化趋势(如图1所示):

图1 变量u和v的变化趋势

图1中，du/dt=0用实线表示，dv/dt=0用虚线表示，随着时间的推移，用箭头表示解(u，v)的趋近方向;其中，区域Ⅰ:du/dt＞0且dv/dt＞0;区域Ⅱ:du/dt＜0且dv/dt＞0;区域Ⅲ:du/dt＜0且dv/dt＜0;区域Ⅳ:du/dt＞0且dv/dt＜0.在微分动力学［7］中，这两条线的交点(0，0)是一个稳定平衡点，对于解(u，v)，从平衡点的某一邻域的任意初值出发，解(u，v)最终都会收敛到此平衡点.把u具有较大值的区域称为兴奋区域;而把起点及其邻域区域称为休眠区域.

从图1可以看到，反应项中的参数a是初始条件的临界值.设解(u，v)的初始值为(u0，v0)，则当u0＞a且v0≤0时，根据图像中解(u，v)的走向，解(u，v)进入兴奋区域;而当u0＜a且v0≥0时，解(u，v)快速进入休眠区域.

利用系统的这个性质，当一幅需要处理的图像作为初始值输入时，该反应扩散系统遍历这幅图像的所有像素点，在适当的临界值a和足够的时间作用下，图像像素点的灰度值划分到图1所示的两个区域，这样就扩大了不同像素点的灰度值.

1.2 FitzHugh-Nagumo反应扩散方程的数值处理

对反应扩散方程(2)，采用九点差分中的显示格式进行求解［6，8］.空间变量(x，y)和时间变量t用有限差分进行离散，空间步长和时间步长分别为 δh，δt，其中:i=［x/δh］，j=［y/δh］，k=［t/δt］，i，j和k是离散空间和时间的指数，那么变量u(x，y，t)的离散格式是uki，j，且

定义∂u/∂t，∂u/∂x，∂u/∂y的离散格式分别为Δtu，Δxu，Δyu，且

由(4)和(5)，方程组(2)的第一个方程变为:

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

方程组(2)的Neumann边界条件的离散式为:

在图像处理中，一幅图像离散化后的数字图像大小为M×N，即其有M行和N列，取i=0，1，…，M-1，j=0，1，…，N-1，那么在初值(U0，V0)和离散边界条件(7)基础上，通过用式(6)的迭代运算，最终获得随时间演化的解(u，v)，即处理后的图像.

2 对木材纹理图像处理算法的实现和结果

2.1 FitzHugh-Nagumo离散模型图像增强算法的实现

按照以上的理论分析，对于一个不带扩散项的常微系统(3)，或者r=0的离散系统(6)模型来说，该模型可以用来增强图像，拉伸模糊图像的对比度.当a为一个固定值时，该模型可以用来处理二值图像;那么当一幅图像为灰度图像时，可以借鉴文献［9］的思想，利用扩散项在时间上的迭代，引入变量

利用a(x)取代固定的a值.得到改进的图像增强算法.该算法步骤如下:

步骤1:取亮度图像u(x，t=0)并归一化到［0，1］;

步骤2:准备好初始条件U0(x)=u(x，t=0)和V0(x)=0;

步骤3:计算a(x);

步骤4:把a(x)代入r=0的离散系统(6)中，通过迭代运算，最后得到边缘检测图.

2.2 FitzHugh-Nagumo离散模型图像增强算法的处理结果及分析

按照2.1算法步骤，对不清晰的有活节的木材亮度图像(图2(a)，181×188)用Matlab7.0软件进行增强处理，这里选择b=1，ε=10-3，Du=0.01，Dv=0.25，δh=1/6，δt=10-3，a(x)=A(U0;1，2，0.3).截取不同演化次数即在时域上的演化图像及其对应的灰度直方图，结果如图2所示.

图2 不同时刻得到的木材纹理图像以及对应的灰度直方图

由图2可知该算法在演化过程中对原活节木材纹理图像进行了有效的增强，其对应的灰度直方图也从原来呈现的尖锐状变化到平滑、均匀状.

该文利用Michelson对比度公式来计算图2中对应图像的对比度

式中，Ⅰmax和Ⅰmin分别表示图像中最亮的亮度和最暗的亮度.相应的图像对比度为(a)CM=1，(b)CM=1.0043，(c)CM=1.0147，(d)CM=1.0250.这个结果表明，图像对比度得到拉伸，纹理图像得到增强.演化结果与我们在1.2中的理论分析相一致:纹理图像在很短的时间内得到增强以后，图像灰度值向稳定值0趋近，在图像亮度上呈暗亮度.

2.3 FitzHugh-Nagumo反应扩散系统对木材纹理图像的边缘检测算法

对于带有扩散项的FitzHugh-Nagumo反应扩散系统模型，参数a的取法，对检测一幅图像的边缘检测同样重要;当a为固定常数时，该模型可以实现二值图像的边缘检测［5］;当参数a取a(x)=A(U0;D，T)［9］来代替固定的a值，该模型可以实现亮度图像的边缘检测;但是这种做法也不尽完美，在该模型作用下图像产生真假双边缘［10］，且随着反应扩散系数的增加，双边缘之间的距离拉长，甚至会产生覆盖边缘的现象，使图像边缘模糊.

针对这一缺陷，该文通过添加一个调整项［11］来改进该算法.这个调整项是一个关于梯度的单调递减的非负函数，取值范围为(0，1］，它能灵活的调整扩散系数，进而调整双边缘的距离，该调整项为:

式中，p为常数，它主要用来控制h的下降速率，

其中，借用Sobel算子的两组3×3的矩阵，分别将之与图像U0作平面卷积，即可分别得出横向及纵向的亮度差分近似值Gx及Gy:

同时，为了得到边缘的稳定的稳态解，该反应扩散系统必须是在强抑制Du≪Dv作用下.该文算法如下:Δ Δ

具体步骤为:

步骤1:取亮度图像u(x，t=0)并归一化到［0，1］;

步骤2:准备好初始图像U0(x)=u(x，t=0)和V0(x)=0;

步骤3:用式(8)和(9)计算得到a(x);

步骤4:把a(x)代入到反应扩散离散系统(10)中，通过式(10)和(11)的离散形式的迭代运算，最后得到边缘检测图.

2.4 FitzHugh-Nagumo反应扩散系统对木材纹理边缘检测结果及分析

按照2.3分析的算法步骤，对有死节的木材亮度图像(图3(a)，183×190)用Matlab 7.0软件进行纹理边缘检测得到图3.

图3

对比图3(a)(b)(c)三幅纹理图像，图3(b)是对图(a)用经典的Sobel算子得到的边缘检测图，很明显，其纹理边缘的线条连续，但是漏检线条太多;图3(c)是用先前的算法得到的边缘检测图，很明显突出了死节内部纹理，周围的部分纹理走势模糊，死节外部纹理线条呈间断、点状;图3(d)是用在先前算法基础上得到的本文算法得到的边缘检测图，与图3(c)比较，死节内部纹理更完整，外部纹理走势更清晰，线条连续性有提高;而与图3(b)比较，死节内部纹理清晰，得到的边缘覆盖较大，缺点就是线条不连续.

总体而言，该文边缘检测算法对木材纹理的边缘检测在一定程度上提高了完整性和清晰性，方便对木材纹理缺陷进行检测;但是在连续性上还有待提高.

［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］ 冈萨雷斯，伍兹著，阮秋琦，等译，数字图像处理［M］.北京:电子工业出版社.2011.

［7］ 陆启韶，彭临平，杨卓琴.常微分方程与动力系统［M］.北京:北京航空航天大学出版社，2009.

［8］ 李荣华，刘播.微分方程数值解法［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的医学图像增强技术研究［J］.2011:40-41.