基于象元的GLCM纹理特征提取

基于象元的GLCM 纹理特征提取
1.算法原理
灰度共生矩阵(GLCM )是一个统计描述图像中的一个局部区域或整个区域相邻象元或一定间距内两象元灰度呈现某种关系的矩阵。

该矩阵中的元素值代表灰度级之间联合条件概率密度(,/,)P i j d θ，(,/,)P i j d θ表示在给定空间距离d 和方向θ时，灰度以i 为始点，出现灰度级为j 的概率（也即频数）。

算法涉及到的几个关键概念：
（1）共生矩阵的维数：在不对原图像灰度级别进行压缩的情况下，取决于原图像灰度的级数；
（2）原图像灰度级别的压缩：为了使共生矩阵的维数不至于太大（太大会降低纹理特征值的计算效率），通常先把原始图像的灰度等级进行压缩，比如从0-255的8 bit 图像压缩为0-31的 5 bit 图像，相应的共生矩阵维数从256*256降低到32*32；
（3）基准窗口（纹理窗口）：尺寸通常为奇数方阵，大小可设为3*3，5*5，7*7等； （4）移动窗口：尺寸大小由基准窗口决定
（5）移动方向和步长：移动方向通常为0度，45度，90度和135度，如图1所示：
图1 窗口移动方向
移动步长与两个因素有关：（1）移动方向；（2）基准窗口中心象元和移动窗口中心象元的距离（用像素表示）。

举例说明，如图2、3所示：
图2 基准窗口和移动窗口示例（1）
示例1：红色框表示基准窗口（大小为3*3），蓝色框表示移动窗口；移动方向为135度；移动步长dx =2，dy =2.
图3 基准窗口和移动窗口示例（2）
示例2：红色框表示基准窗口（大小为3*3），蓝色框表示移动窗口；移动方向不属于常
用的四个方向的任意一个；移动步长dx=1，dy=3.
2.算法过程详细描述
假设原图像为6 * 6大小，图像的像素如图4矩阵所示：
123456
626543
254165
413625
523651
123564
图4 图像像素的矩阵表示
下面以135度方向，移动步长dx=2，dy=2，基准窗口大小为3 * 3为例来说明共生矩
阵生成和纹理特征值计算过程。

第一步：共生矩阵的生成
图5 移动窗口和基准窗口的关系
如图5所示：
（1）先根据基准窗口大小在影像左上角选择一个3 *3大小的区域，如红色框所示；（2）然后根据移动步长dx=2，dy=2 ，方向为135度，获得生成共生矩阵的第二个区域，即为移动窗口，如蓝色框所示；
（3）以基准窗口和移动窗口（矩阵）相对相同位置的像素进行统计得到共生矩阵。

假设原始遥感影像已被压缩到6级，则生成的共生矩阵为6*6的方阵。

注意压缩采用的
具体公式为：dn dn
grayNums
dn dn
-
-
min
max min
*（该公式的原型是图像的线
性拉伸公式，详细参见文后的附录）。

dn表示原图像中某一象元值，dn
min
表原图像
所有象元值的最小值，
max
dn表原图像所有象元值的最大值，grayNums表示原图像灰度等级将要压缩到的等级，例如32。

得到的共生矩阵如图6所示：
图6 灰度共生矩阵（频率表示）
上图共生矩阵生成说明：共生矩阵的行号对应基准窗口的象元值，列号对应移动窗口的象元值。

例如红色框中的左上角像素值为1，而蓝色框中的左上角像素值为4，则共生矩阵中的（1，4）位置即为1，当然如果红色框和蓝色框的相同位置还分别出现1和4时，则共生矩阵的（1，4）位置相应的应累加1，即为2。

共生矩阵其他位置的计算以此类推。

由于纹理特征计算利用的是概率，所以还应该把共生矩阵中的值转换为概率值，方法是用各个元素值除以矩阵中所有元素的和（上图元素和为9），最后可以得到最后的方阵，即为共生矩阵。

如图7所示：
图7 灰度共生矩阵（概率表示）
第二步：基于共生矩阵的纹理特征值计算
以熵（Entropy）纹理特征值为例来说明，Entropy主要用来衡量图像的无序性。

当图像的纹理极不一致时，灰度共生矩阵中各元素的值将偏小，这将意味该纹理具有较大的熵值。

例如，如果一个窗口内像素的灰度值呈现完全的随机分布(白噪声)，这样窗口内图像的直方图分布将呈现一个常数函数，也就是说当所有灰度值对出现的频率相等时，Entropy达到最
大。

熵的表达式如下所示：
*ln (,)k k
quant quant i=0
j=0
Entropy =p(i,j)p i j -∑
∑
其中k quant 表示原图像灰度压缩的等级（共生矩阵的维数），(,)p i j 表示共生矩阵中元素的值，也就是概率值。

通过上式可以计算出一个熵值，把这个熵值赋给红色框区域的中心象元即可完成该基准窗口（纹理窗口）的纹理特征值计算。

至此，算法已经计算得出了一个像素的纹理特征值，接下来我们要计算其他像素的纹理特征值。

图8显示的是下一步要计算的纹理特征值所依据的基准窗口和移动窗口，而计算过程与上面描述的过程类似。

图8 基准窗口和移动窗口的下一位置
总的来说，基准窗口的移动是逐个像素移动的，而移动窗口的移动是根据基准窗口和移动方向、移动步长大小得到的。

Haralick(1973)根据灰度共生矩阵定义了14种纹理特征，常用的用于提取遥感图像中纹理信息的特征统计量主要有：均值（Mean ）、标准差（Standard deviation ）、方差（Variance ）、同质性(Homogeneity)、对比度(Contrast)、不相似性（Dissimilarity ）、熵（Entropy ）、角度二阶矩（或称为能量）(Angular Second Moment)、相关性（Correlation ）。

上文已经详细介绍了熵的计算过程，其他纹理特征统计量的计算与熵纹理特征值计算类似，不同之处仅仅在于公式的表达。

下面列出了除熵统计量以外的其他纹理统计量计算公式：
*k k
quant quant i=0
j=0
Mean =
p(i,j)i ∑∑
2*()k k
quant quant i=0
j=0
Variance =
p(i,j)i Mean -∑∑
Std =
22*()k k
quant quant i=0
j=0
Contrast =
p(i,j)i j -∑∑
*k k
quant quant i=0
j=0
Dissimilarity =
p(i,j)i j -∑∑
2()k k
quant quant i=0
j=0
ASM Energy =
p(i,j)∑∑
2
()*()*(,)k k
quant quant i=0
j=0
i Mean j Mean p i j Correlation =
Variance
--∑∑
3.算法特殊情况处理
由于算法中的每次计算都是针对基准窗口的中心象元进行纹理特征赋值，所以算法实现过程涉及到边界元素的处理，下面提供几种边界元素的处理的方法：
方法1：把所有边界元素，也就是不能作为基准窗口中心象元的元素，将这些元素的最终的纹理特征值全部赋值为，ENVI 软件采用该种方法处理，计算效率较好；
方法2：根据处理的窗口大小和步长大小，用0值或者边界值扩展源图像的边界，使源图像扩大维数，并且将添加的边界象元值也参加计算，最终使原来的边界元素也能成为中心元素进行计算，但最终处理结果仍保持原始图像尺寸，处理效率较方法1差。

4.算法实现和试验结果
采用MapGIS70作为二次开发平台实现了基于GLCM 的纹理特征提取算法。

界面如图9下所示：
图9 纹理特征提取界面
图10为一副三波段QuickBird 数据纹理特征熵的计算结果示例：
图10 纹理特征熵提取结果
图11 原始遥感影像图
附：图像的线性拉伸
图像的线性拉伸是将范围为[a，b]数字图像的灰度f（x，y）变换为范围为[c，d]的灰度g（x，y），如图12所示。

图12灰度输入输出变换简图
线性拉伸的变换数学公式为：。