CN101770577A - 无人机航片松材线虫病病死木的膨胀-剔除信息提取法 - Google Patents

Abstract

一种无人机航片松材线虫病病死木的膨胀-剔除信息提取法，属数字图像处理技术领域。以无人机彩色数码相片作为输入图像，经过平滑，将RGB模式转换成CIE(Luv)模式。统计直方图，以小尺度直方图用非参数聚类方法-爬峰法进行初始聚类。在图像空间选取一些同属病死林木区域的类，经合并成为种子类。在直方图内基于种子类用一6-联通三维模板持续进行数学形态学膨胀，允许有非病死林木区域像素的对应栅格进入目标类。用鼠标在图像上提取经过膨胀进入病死林木区域的非病死林木区域像素，将这些像素所在的栅格从直方图目标类中剔除，并规定这些栅格在以后的膨胀操作中不被再次膨胀。反复进行膨胀和剔除操作，直至得到满意的结果。

Description

无人机航片松材线虫病病死木的膨胀-剔除信息提取法

【技术领域】

本发明属于计算机数字图像处理技术领域，具体是一种无人机航片松材线虫病病死林木的计算机人机交互识别方法。

【背景技术】

松材线虫病是松类树种的一种毁灭性病害，是我国目前发生的最严重、最具威胁性的外来入侵森林病害。这种病在传播方式上具有跳跃式传播特性。将病死林木及时清除是目前防止此病蔓延的最有效方法。病死林木零星分布在崇山峻岭，仅仅依靠人工地面查找不仅劳动强度很大，而且很容易遗漏，给病害蔓延留下隐患。用无人机携带高清数码相机，从空中拍摄相片，然后配以地理坐标，在相片上确定病死林木位置，输入PDA用以野外导航，这是正在研究的新方法，其中的核心技术是计算机人机交互识别病死林木，这从数字图像角度看，就是图像分割，即将病死林木分割出来的技术。

图像分割的方法众多，大致可分为两大类：一类是基于边界的，一类是基于区域的。病死林木的色彩大体呈现红色，但边界不明显，同时森林图像一般都缺少明显的边界，而且基于边界的方法不能在分割的同时确定区域性质，所以不适合使用基于边界的方法。基于区域的方法是根据区域内图像的一致性进行分割的，这类方法一般又可分为基于图像空间和基于特征空间两类。基于图像空间的方法也不能在分割的同时确定区域的性质。基于特征空间的方法又可分为参数方法和非参数方法两类。参数方法需要有模型假定，就是图像区域对应的特征空间的类的分布需要有模型假定的形状，如团状等，非参数方法则不需要这样的假定。当区域特征不具备模型假定的条件，即特征分布是不规则的，就适合采用非参数方法。在采用非参数方法时，通常希望通过选择最佳的特征空间划分尺度，对特征空间进行聚类，得到和图像区域对应的类。但这通常很难得到满意的结果，因为最佳的分割尺度在很多情况下是不存在的：有的区域需要大一点，有的需要小一点，很难协调。已有方法都不能很好解决无人机航片松材线虫病病死林木的提取问题。

【发明内容】

本发明需要解决的技术问题是针对低空无人机装载的高清数码相机拍摄的数码相片，用人机交互的方式提取由于松材线虫病造成的死亡林木，提出一种无人机航片松材线虫病病死木的膨胀-剔除信息提取法。

为解决上述技术问题，本发明的技术方案是按如下步骤进行：

(1)用常规方法输入原始RGB真彩色无人机数码相片，规定病死林木区域为目标区域，其他区域为非目标区域；原始RGB模式真彩色图像中的目标区域呈现红色，灰度图中呈现亮色；

(2)用均匀模板对原始图像进行平滑运算；

(3)将图像从RGB模式转换成CIE(Luv)模式，CIE(Luv)模式下的数据拉伸到0～255并取整，CIE(Luv)模式彩色图像中目标区域呈蓝色；以下的操作均对CIE(Luv)模式数据进行；

(4)生成尺度化后的CIE(Luv)模式数据三维直方图：确定一个聚类尺度s＝2^λ，λ为整数，0≤λ≤7，原始特征空间每个分量划分等级为256，确定聚类尺度后，新的量化等级为256/s，建立三维直方图h[256/s，256/s，256/s]作为爬峰法聚类的基础，直方图中每个栅格h[i，j，k]具有一个数值，该数值就是新尺度下三个分量亮度值为[i，j，k]的像素数；

(5)爬峰法聚类：爬峰法聚类在三维特征空间——三维直方图中进行，直方图中除峰以外的每个非空栅格都有一个指针指向相邻的一个栅格，在这个方向上密度具有最大的梯度值，梯度值等于相邻栅格像素数的差值除以栅格中心间的欧氏距离，当最大梯度值多于一个时就任选一个，若一个栅格的象素数大于所有相邻栅格的象素数，它就处于局部最大值，该栅格为一个峰；峰确定以后，根据指针的指向，将所有非峰且非空的栅格归拢到各自的峰，一个峰形成一个类，从而实现特征空间的聚类；爬峰法聚类结果受聚类尺度s的影响最大，s越大，类别数越少，s越小，类别数越多；s确定的原则是：s值宜小，0≤λ≤2，以确保目标区域的像素不和非目标区域的像素聚到一起，但允许同属目标区域的像素聚成很多类；

(6)生成代码化直方图：聚类完成后将直方图中的像素数数据改写成类代码，属于同一个类的栅格赋予同一个代码，形成代码化直方图C[i/s，j/s，k/s]，同时计算各个类的平均亮度值；

(7)生成代码图像：将聚类结果映射到图像上，生成一个以代码为像素值的图像——代码图像，实现图像的初步分割，分割图像显示时用各个类的平均亮度值；代码图像生成的算法是：若一个像素的亮度值是[i，j，k]，则该像素的类代码c是

c＝C[i/s，j/s，k/s]；

(8)交互式提取种子类：用鼠标在图像上提取数个同属目标区域的小矩形，通过代码图像将这些小矩形内像素所涉及的目标类全部提取出来，将提取的这些类合并成一个类，合并后的目标类形成一个种子类，赋予他们同一个代码；

(9)设定一个膨胀模板：膨胀模板采用6-联通三维模板，即和中心栅格6个面相连的6个相邻栅格为数学形态学的膨胀对象；

(10)种子类特征空间数学形态学膨胀：在代码化直方图内对选定的种子类的直方图用6-联通三维模板反复进行数学形态学膨胀，其中一次膨胀的算法如下：膨胀模板的中心棚格遍历整个直方图，若膨胀模板的中心栅格覆盖的直方图栅格属于种子类，则将模板其他栅格覆盖的直方图非种子类栅格改写为种子类栅格；在一次膨胀中，若模板中心栅格覆盖的是本次膨胀中前面刚膨胀的种子类，则不进行膨胀操作，直接将模板中心位置移至下一栅格；一次膨胀完成后，显示膨胀后直方图对应的图像，非种子类区域不显示，这时图像种子类的区域有所扩展；这样反复数次膨胀操作，直至图像中明显出现了非目标区域像素；

(11)剔除错误膨胀部分：反复数次膨胀后，图像上明显出现的非目标区域是错误膨胀的，需予以剔除，剔除通过图像选择进行，通过人机交互的方式反复用鼠标在图像上选择错误膨胀区域的像素，将选中的图像像素映射到直方图，将对应的直方图栅格从直方图中剔除，并将落入这些剔除栅格的所有图像像素从图像目标区域中剔除；赋予被剔除部分栅格一个特殊代码，以禁止它们被再次膨胀；

(12)如果结果不满意，按步骤(10)进行，否则按步骤(13)进行；

(13)后处理：为使分割后的图像整洁，进行面积滤波，以删除一些细小的噪声区域，给定一个以像素数为单位的面积阈值；当一目标区域小于这个阈值时，就删除该区域，同时将目标区域内部小于阈值的非目标区域改为目标区域；阈值要根据图像噪声情况从小到大逐步递增，以确保被滤掉的区域是图像噪声引起的误差；

(14)输出画上边界线的目标区域图像：为了直观检验信息提取的精度，分别在原始RGB模式和CIE(Luv)模式图像上画上提取区域的边界线；

(15)以文本方式输出每个目标区域的地理坐标，本信息提取法完成。

本发明的有益效果是：提供一种创新性的专门针对无人机数码相片松材线虫病病死林木信息人机交互提取方法。病死林木在5到10月份之间都会出现，死亡时间的较大差异导致病死林木在图像空间的色彩差异较大，以及在特征空间中的分布状态十分复杂，即是严重非参数的。这种特殊情况使得现有方法不适合这种特殊信息的提取。本方法的整个过程都是非参数的，从前期的爬峰法聚类，到核心部分的膨胀-剔除技术，均没有任何参数模型假定。信息提取在特征空间和图像空间交叉进行，方法灵活、直观，人为控制能力强。本方法也可拓展应用于其他专题信息的提取。

【附图说明】

图1为本方法的流程图。

图2为本发明实施例中的原始图像的灰度图。

图3为本发明实施例中5＊5均匀模板平滑后的图像。

图4为本发明实施例中图像的CIE(Luv)模式图。

图5为本发明实施例中初始聚类即爬峰法聚类结果图。

图6为本发明实施例中在图像上提取种子类的示意图。

图7为本发明实施例中提取的用CIE(Luv)模式显示的种子类。

图8为本发明实施例中提取的种子类像素在直方图中的分布图。

图9为本发明特征空间数学形态学膨胀模板即6-联通三维模板。

图10为本发明实施例中用模板对种子类的直方图(图8)进行数学形态学膨胀的结果图。

图11为本发明实施例中直方图膨胀后对应的图像。

图12为本发明人机交互在图像上选择错误膨胀区域的示意图。

图13为本发明实施例中膨胀、剔除操作结束时的图像。

图14为本发明实施例中膨胀、剔除操作结束时的直方图。

图15为本发明实施例中用4联通模板、阈值分别为5和10进行两次面积滤波后的图像。

图16为本发明实施例中提取的目标区域套上边界线后的RGB模式原始图像。

图17为本发明实施例中提取的目标区域套上边界线后的CIE(Luv)模式原始图像。

图18为本发明实施例中输出的病死林木地理坐标(仅列出部分)。

【具体实施方式】

下面结合实施例并参考附图对本发明作进一步详述：

病死林木死亡时间不同，色彩有差异，刚死亡的呈现鲜红色，死亡时间久的呈现暗红色、红褐色等，本方法对这种色彩差异具有较好的适应性。本发明称图像空间中的病死林木区域为目标区域，其他区域为非目标区域，特征空间中的病死林木类为目标类，其他类为非目标类。

实施流程参见图1。

步骤1，输入无人机真彩色数码相片。用常规方法输入原始RGB真彩色无人机数码相片。规定病死林木区域为目标区域，其他区域为非目标区域。目标区域在原始RGB模式真彩色图像中呈现红色，在灰度图中呈现亮色，见图2。从纸质文件显示效果考虑附图仅取整个实施例图像的局部，大小为500＊400像素，整个实施例图像大小为3072＊2304像素；

步骤2，用均匀模板对原始图像进行平滑运算。本发明实施例所用模板为5＊5均匀模板，平滑结果见图3；

步骤3，将图像从RGB模式转换成CIE(Luv)模式。CIE(Luv)模式下的数据拉伸到0～255并取整。接下去的操作均对CIE(Luv)模式数据进行。本发明实施例CIE(Luv)模式图像见图4，CIE(Luv)模式彩色图像中目标区域呈现蓝色，在灰度图像中呈现亮色。具体算法是：令

$A=\left(\begin{array}{ccc}0.607& 0.174& 0.200\\ 0.299& 0.587& 0.114\\ 0& 0.066& 1.116\end{array}\right)$

(x₀，y₀，z₀)^T＝A(255，255，255)^T

如果x₀+15y₀+3z₀≠0，令

u₀＝4x₀/(x₀+15y₀+3z₀)

v₀＝6y₀/(x₀+15y₀+3z₀)

否则令

u₀＝0

v₀＝0

读出一个像素的红、绿、蓝(即R、G、B)亮度值，令

(x，y，z)^T＝A(R，G，B)^T

如果x+15y+3z≠0，令

u₁＝4x/(x+15y+3z)

v₁＝6y/(x+15y+3z)

否则令

u₁＝0

v₁＝0

令TMP＝y/y₀，如果TMP＜0.011，则令TMP＝0.011。令

$L=116\mathrm{TM}{P}^{\frac{1}{3}}-16$

u＝13L(u₁-u₀)

v＝13L(v₁-v₀)

遍历整幅图像，找出最大、最小的L、u、v：MaxL、Maxu、Maxv、MinL、Minu、Minv。令

DeltaL₀＝MaxL-MinL，如果DeltaL₀＝0，令DeltaL₀＝1

Deltau₀＝Maxu-Minu，如果Deltau₀＝0，令Deltau₀＝1

Deltav₀＝Maxv-Minv，如果Deltav₀＝0，令Deltav₀＝1

DeltaL＝1/DeltaL₀

Deltau＝1/Deltau₀

Deltav＝1/Deltav₀

重新遍历整幅图像，读出每个像素的R、G、B亮度值，进行以下计算，得出每个像素的拉伸到0～255整数的CIE(Luv)模式的L、u、v值：

(x，y，z)^T＝A(R，G，B)^T

如果x+15y+3z≠0

u₁＝4x/(x+15y+3z)

v₁＝6y/(x+15y+3z)

否则

u₁＝0

v₁＝0

令TMP＝y/y₀，如果TMP＜0.011，则令TMP＝0.011。令

$L=255(116{\mathrm{TMP}}^{\frac{1}{3}}-16-\mathrm{MinL})*\mathrm{DeltaL}$

u＝255[13L(u₁-u₀)-Minu]*Deltau

v＝255[13L(v₁-v₀)-Minv)]*Deltav

步骤4，生成尺度化后的CIE(Luv)模式数据三维直方图。确定一个聚类尺度s＝2^λ，λ为整数，0≤λ≤7。原始特征空间每个分量划分等级为256，确定聚类尺度后，新的量化等级为256/s。建立三维直方图h[256/s，256/s，256/s]作为爬峰法聚类的基础。直方图中每个栅格h[i，j，k]具有一个数值，这个数值就是新尺度下三个分量亮度值为[i，j，k]的像素数。本发明实施例中聚类尺度为s＝2(即λ＝1)，三维直方图就是数组h[128，128，128]；

步骤5，爬峰法聚类。爬峰法聚类在三维特征空间——三维直方图中进行。直方图中每个非空栅格(除峰以外)都有一个指针指向相邻的一个栅格，在这个指针方向上密度具有最大的梯度值。梯度值等于相邻栅格像素数的差值除以相邻栅格中心间的欧氏距离。当最大梯度值多于一个时就任选一个方向作为指针方向。如果一个棚格的象素数大于所有相邻棚格的象素数，它就处于局部最大值，这个栅格就是一个峰(Peak或Mode)。峰确定以后，根据指针的指向，将所有非峰且非空的栅格归拢到各自的峰，一个峰形成一个类，从而实现特征空间的聚类。爬峰法聚类结果受聚类尺度s的影响最大。s越大，类别数越少，s越小，类别数越多。s确定的原则是：s值宜小，0≤λ≤2，以确保目标区域的像素不和非目标区域的像素聚到一起，但允许同属目标区域的像素聚成很多类。爬峰法具体算法可参考“Imagesegmentation by a parallel，non-parametric histogram basedclustering algorithm”(A.Khotanzad et al.，Pattern Recognition，23(9)：961-973)，“面向类的图像分割方法研究”(葛宏立，北京林业大学博士论文，2004)；

步骤6，生成代码化直方图。聚类完成后将直方图中的像素数数据改写成类代码，属于同一个类的栅格赋予同一个代码，形成代码化直方图C[256/s，256/s，256/s]。同时计算各个类的平均亮度值。

步骤7，生成代码图像。将聚类结果映射到图像上，生成一个以代码为像素值的图像——代码图像，实现图像的初步分割，分割图像显示时用各个类的平均亮度值。本发明实施例中聚类后对应的初始分割图像见图5，聚类尺度s为2，即原来的256个量化等级转化为128个量化等级。比较图5和图4可以看出，聚类后结果和原图没有明显区别。聚类的类别数为7118。

代码图像生成的算法是：如果一个像素的亮度值是[i，j，k]，则该像素的类代码c是

c＝C[i/s，j/s，k/s]；

步骤8，交互式提取种子类。如图6所示，用鼠标在图像上提取数个同属目标区域的小矩形。通过代码图像将这些小矩形内像素所涉及的目标类全部提取出来。本实施例提取结果如图7所示，共提取79类。将提取的这些类合并成一个类，合并后的目标类形成一个种子类，赋予他们同一个代码。种子类的三维直方图见图8；

步骤9，设定一个膨胀模板。膨胀模板采用6-联通三维模板，即和中心栅格6个面相连的6个相邻栅格为数学形态学的膨胀对象，如图9所示；

步骤10，种子类特征空间数学形态学膨胀。在代码化直方图内对选定的种子类的直方图用6-联通三维模板(见图9)反复进行数学形态学膨胀。其中一次膨胀的算法如下：膨胀模板的中心栅格遍历整个直方图，如果膨胀模板的中心栅格覆盖的直方图栅格属于种子类，则将模板其他栅格覆盖的直方图非种子类栅格改写为种子类栅格。在一次膨胀中，如果模板中心栅格覆盖的是本次膨胀中前面刚膨胀的种子类则不进行膨胀操作，直接将模板中心位置移至下一栅格。一次膨胀完成后，显示膨胀后直方图对应的图像，非种子类区域不显示，这时图像种子类的区域(即目标区域)有所扩展。这样反复数次膨胀操作，直至图像中明显出现了非目标区域像素。本实施例数次膨胀后的直方图见图10，对应的图像见图11；

步骤11，剔除错误膨胀部分。反复数次膨胀后，图像上明显出现的非目标区域是错误膨胀的，如图11中的暗色部分，需予以剔除。剔除通过图像选择进行。用人机交互的方式在图像上反复用鼠标选择错误膨胀区域的像素，如图12所示，将选中的图像像素映射到直方图，将对应的直方图栅格从直方图中剔除，并将落入这些剔除栅格的所有图像像素从图像目标区域中剔除。赋予被剔除部分栅格一个特殊代码，以禁止它们被再次膨胀；

步骤12，如果结果不满意，转步骤10，否则停止膨胀-剔除的重复操作，转步骤13。本发明实施例膨胀-剔除操作后的图像见图13，对应的直方图见图14；

步骤13，后处理。为使分割后的图像看上去整洁，进行面积滤波，以删除一些细小的噪声区域。给定一个以像素数为单位的面积阈值，当一目标区域小于这个阈值时，就删除该区域。同时目标区域内部小于阈值的非目标区域改为目标区域。阈值要根据图像噪声情况从小到大逐步递增，以确保被滤掉的区域是图像噪声引起的误差。图15为本发明实施例中用4-联通模板、阈值分别为5和10进行两次面积滤波后的图像；

步骤14，输出画上边界线的目标区域图像。为直观检验信息提取的精度，分别在原始RGB模式和CIE(Luv)模式图像上画上提取的目标区域的边界线，分别参见图16、17；

步骤15，以文本方式输出每个目标区域的地理坐标，见图18。本信息提取法完成。

上述步骤10描述的膨胀和步骤11描述的剔除，是本发明首次提出的，它们一起构成本发明的核心，姑且称之为膨胀-剔除信息提取法。

本发明可拓展应用于下列情况：无人机之外的其他飞行器作为搭载工具，如卫星、有人驾驶飞机、气球等；数码相机之外的其他传感器获取的彩色图像；松材线虫病病死木之外的其他专题信息提取。只要彩色图像的空间分辨率高于50cm，即图像一个像素代表的地面面积小于50cm＊50cm。

Claims (1)

1.一种无人机航片松材线虫病病死木的膨胀-剔除信息提取法，其特征在于按如下步骤进行：

(1)用常规方法输入原始RGB真彩色无人机数码相片，规定病死林木区域为目标区域，其他区域为非目标区域；原始RGB模式真彩色图像中的目标区域呈现红色，灰度图中呈现亮色；

(2)用均匀模板对原始图像进行平滑运算；

(3)将图像从RGB模式转换成CIE(Luv)模式，CIE(Luv)模式下的数据拉伸到0～255并取整，CIE(Luv)模式彩色图像中目标区域呈蓝色；以下的操作均对CIE(Luv)模式数据进行；

(4)生成尺度化后的CIE(Luv)模式数据三维直方图：确定一个聚类尺度s＝2^λ，λ为整数，0≤λ≤7，原始特征空间每个分量划分等级为256，确定聚类尺度后，新的量化等级为256/s，建立三维直方图h[256/s，256/s，256/s]作为爬峰法聚类的基础，直方图中每个栅格h[i，j，k]具有一个数值，该数值就是新尺度下三个分量亮度值为[i，j，k]的像素数；

(5)爬峰法聚类：爬峰法聚类在三维特征空间——三维直方图中进行，直方图中除峰以外的每个非空栅格都有一个指针指向相邻的一个栅格，在这个方向上密度具有最大的梯度值，梯度值等于相邻栅格像素数的差值除以栅格中心间的欧氏距离，当最大梯度值多于一个时就任选一个，若一个栅格的象素数大于所有相邻栅格的象素数，它就处于局部最大值，该栅格为一个峰；峰确定以后，根据指针的指向，将所有非峰且非空的栅格归拢到各自的峰，一个峰形成一个类，从而实现特征空间的聚类；爬峰法聚类结果受聚类尺度s的影响最大，s越大，类别数越少，s越小，类别数越多；s确定的原则是：s值宜小，0≤λ≤2，以确保目标区域的像素不和非目标区域的像素聚到一起，但允许同属目标区域的像素聚成很多类；

(6)生成代码化直方图：聚类完成后将直方图中的像素数数据改写成类代码，属于同一个类的栅格赋予同一个代码，形成代码化直方图C[i/s，j/s，k/s]，同时计算各个类的平均亮度值；

(7)生成代码图像：将聚类结果映射到图像上，生成一个以代码为像素值的图像——代码图像，实现图像的初步分割，分割图像显示时用各个类的平均亮度值；代码图像生成的算法是：若一个像素的亮度值是[i，j，k]，则该像素的类代码c是

c＝C[i/s，j/s，k/s]；

(8)交互式提取种子类：用鼠标在图像上提取数个同属目标区域的小矩形，通过代码图像将这些小矩形内像素所涉及的目标类全部提取出来，将提取的这些类合并成一个类，合并后的目标类形成一个种子类，赋予他们同一个代码；

(9)设定一个膨胀模板：膨胀模板采用6-联通三维模板，即和中心栅格6个面相连的6个相邻栅格为数学形态学的膨胀对象；

(10)种子类特征空间数学形态学膨胀：在代码化直方图内对选定的种子类的直方图用6-联通三维模板反复进行数学形态学膨胀，其中一次膨胀的算法如下：膨胀模板的中心栅格遍历整个直方图，若膨胀模板的中心栅格覆盖的直方图栅格属于种子类，则将模板其他栅格覆盖的直方图非种子类栅格改写为种子类栅格；在一次膨胀中，若模板中心栅格覆盖的是本次膨胀中前面刚膨胀的种子类，则不进行膨胀操作，直接将模板中心位置移至下一栅格；一次膨胀完成后，显示膨胀后直方图对应的图像，非种子类区域不显示，这时图像种子类的区域有所扩展；这样反复数次膨胀操作，直至图像中明显出现了非目标区域像素；

(11)剔除错误膨胀部分：反复数次膨胀后，图像上明显出现的非目标区域是错误膨胀的，需予以剔除，剔除通过图像选择进行，通过人机交互的方式反复用鼠标在图像上选择错误膨胀区域的像素，将选中的图像像素映射到直方图，将对应的直方图栅格从直方图中剔除，并将落入这些剔除栅格的所有图像像素从图像目标区域中剔除；赋予被剔除部分栅格一个特殊代码，以禁止它们被再次膨胀；

(12)如果结果不满意，按步骤(10)进行，否则按步骤(13)进行；

(13)后处理：为使分割后的图像整洁，进行面积滤波，以删除一些细小的噪声区域，给定一个以像素数为单位的面积阈值；当一目标区域小于这个阈值时，就删除该区域，同时将目标区域内部小于阈值的非目标区域改为目标区域；阈值要根据图像噪声情况从小到大逐步递增，以确保被滤掉的区域是图像噪声引起的误差；

(14)输出画上边界线的目标区域图像：为了直观检验信息提取的精度，分别在原始RGB模式和CIE(Luv)模式图像上画上提取区域的边界线；

(15)以文本方式输出每个目标区域的地理坐标，本信息提取法完成。

Images