遥感影像纹理特征是描述影像中像素间空间关系的统计特征，常用于地物分类、目标识别和变化检测等遥感应用中。常见的纹理特征计算方式包括灰度共生矩阵（GLCM）、灰度差异矩阵（GLDM）、灰度不均匀性矩阵（GLRLM）等。这些方法通过对像素灰度值的统计分析，揭示了影像中像素间的空间分布规律，从而提取出纹理特征。
常用的纹理特征包括
1. 对比度（Contrast）： 描述图像中灰度级之间的对比程度。
2. 同质性（Homogeneity）： 描述图像中灰度级分布的均匀程度。
3. 熵（Entropy）： 描述图像中灰度级分布的不确定性。
4. 方差（Variance）： 描述图像中灰度级的离散程度。
5. 相关性（Correlation）： 描述图像中像素间的灰度值相关程度。
6. 能量(Energy)： 是灰度共生矩阵各元素值的平方和。
7. 均值（Mean）： 描述图像中每个灰度级出现的平均频率。
8. 不相似度（Dissimilarity）： 描述图像中不同灰度级之间的差异程度。
这些纹理特征可以通过GLCM等方法计算得到，用于描述遥感影像的纹理信息，对于遥感影像的分类和分析具有重要意义。
在ENVI中可以直接使用工具箱的工具直接计算得到输出影像
此外也可以使用python的skimage.feature库，这个库的greycoprops函数只有contrast、dissimilarity、homogeneity、ASM、energy和correlation这几个特征没有均值、方差和熵。因此下面对这个库的原函数进行修正，添加了方差等计算的greycoprops函数
def greycoprops(P, prop='contrast'):
"""Calculate texture properties of a GLCM.
Compute a feature of a grey level co-occurrence matrix to serve as
a compact summary of the matrix. The properties are computed as
follows:
- 'contrast': :math:`\\sum_{i,j=0}^{levels-1} P_{i,j}(i-j)^2`
- 'dissimilarity': :math:`\\sum_{i,j=0}^{levels-1}P_{i,j}|i-j|`
- 'homogeneity': :math:`\\sum_{i,j=0}^{levels-1}\\frac{P_{i,j}}{1+(i-j)^2}`
- 'ASM': :math:`\\sum_{i,j=0}^{levels-1} P_{i,j}^2`
- 'energy': :math:`\\sqrt{ASM}`
- 'correlation':
.. math:: \\sum_{i,j=0}^{levels-1} P_{i,j}\\left[\\frac{(i-\\mu_i) \\
(j-\\mu_j)}{\\sqrt{(\\sigma_i^2)(\\sigma_j^2)}}\\right]
Each GLCM is normalized to have a sum of 1 before the computation of texture
properties.
Parameters
----------
P : ndarray
Input array. `P` is the grey-level co-occurrence histogram
for which to compute the specified property. The value
`P[i,j,d,theta]` is the number of times that grey-level j
occurs at a distance d and at an angle theta from
grey-level i.
prop : {'contrast', 'dissimilarity', 'homogeneity', 'energy', \
'correlation', 'ASM'}, optional
The property of the GLCM to compute. The default is 'contrast'.
Returns
-------
results : 2-D ndarray
2-dimensional array. `results[d, a]` is the property 'prop' for
the d'th distance and the a'th angle.
References
----------
.. [1] The GLCM Tutorial Home Page,
http://www.fp.ucalgary.ca/mhallbey/tutorial.htm
Examples
--------
Compute the contrast for GLCMs with distances [1, 2] and angles
[0 degrees, 90 degrees]
>>> image = np.array([[0, 0, 1, 1],
... [0, 0, 1, 1],
... [0, 2, 2, 2],
... [2, 2, 3, 3]], dtype=np.uint8)
>>> g = greycomatrix(image, [1, 2], [0, np.pi/2], levels=4,
... normed=True, symmetric=True)
>>> contrast = greycoprops(g, 'contrast')
>>> contrast
array([[0.58333333, 1. ],
[1.25 , 2.75 ]])
"""
check_nD(P, 4, 'P')
(num_level, num_level2, num_dist, num_angle) = P.shape
if num_level != num_level2:
raise ValueError('num_level and num_level2 must be equal.')
if num_dist <= 0:
raise ValueError('num_dist must be positive.')
if num_angle <= 0:
raise ValueError('num_angle must be positive.')
# normalize each GLCM
P = P.astype(np.float32)
glcm_sums = np.apply_over_axes(np.sum, P, axes=(0, 1))
glcm_sums[glcm_sums == 0] = 1
P /= glcm_sums
# create weights for specified property
I, J = np.ogrid[0:num_level, 0:num_level]
if prop == 'contrast':
weights = (I - J) ** 2
elif prop == 'dissimilarity':
weights = np.abs(I - J)
elif prop == 'homogeneity':
weights = 1. / (1. + (I - J) ** 2)
elif prop in ['ASM', 'energy', 'correlation', 'entropy', 'mean', 'variance']:
pass
else:
raise ValueError('%s is an invalid property ?' % (prop))
# compute property for each GLCM
if prop == 'energy':
asm = np.apply_over_axes(np.sum, (P ** 2), axes=(0, 1))[0, 0]
results = np.sqrt(asm)
elif prop == 'ASM':
results = np.apply_over_axes(np.sum, (P ** 2), axes=(0, 1))[0, 0]
elif prop == 'entropy':
results = np.apply_over_axes(np.sum, -P * np.log10(P+0.00001), axes=(0, 1))[0, 0]
elif prop == 'correlation':
results = np.zeros((num_dist, num_angle), dtype=np.float64)
I = np.array(range(num_level)).reshape((num_level, 1, 1, 1))
J = np.array(range(num_level)).reshape((1, num_level, 1, 1))
diff_i = I - np.apply_over_axes(np.sum, (I * P), axes=(0, 1))[0, 0]
diff_j = J - np.apply_over_axes(np.sum, (J * P), axes=(0, 1))[0, 0]
std_i = np.sqrt(np.apply_over_axes(np.sum, (P * (diff_i) ** 2),
axes=(0, 1))[0, 0])
std_j = np.sqrt(np.apply_over_axes(np.sum, (P * (diff_j) ** 2),
axes=(0, 1))[0, 0])
cov = np.apply_over_axes(np.sum, (P * (diff_i * diff_j)),
axes=(0, 1))[0, 0]
# handle the special case of standard deviations near zero
mask_0 = std_i < 1e-15
mask_0[std_j < 1e-15] = True
results[mask_0] = 1
# handle the standard case
mask_1 = mask_0 == False
results[mask_1] = cov[mask_1] / (std_i[mask_1] * std_j[mask_1])
elif prop in ['contrast', 'dissimilarity', 'homogeneity']:
weights = weights.reshape((num_level, num_level, 1, 1))
results = np.apply_over_axes(np.sum, (P * weights), axes=(0, 1))[0, 0]
elif prop == 'mean':
I = np.array(range(num_level)).reshape((num_level, 1, 1, 1))
results = np.apply_over_axes(np.sum, (P * I), axes=(0, 1))[0, 0]
elif prop == 'variance':
I = np.array(range(num_level)).reshape((num_level, 1, 1, 1))
mean = np.apply_over_axes(np.sum, (P * I), axes=(0, 1))[0, 0]
results = np.apply_over_axes(np.sum, (P * (I - mean) ** 2), axes=(0, 1))[0, 0]
return results
具体影像分析时还需要考虑灰色关联矩阵计算的角度、步长、灰度级数和窗口大小。以一张多光谱影像为例，相关性使用了greycoprops，其他的特征使用公式计算，实际上导入上面的新greycoprops函数后其他特征都可以用函数计算，具体代码如下，输出结果为多光谱各个波段的纹理特征的多通道影像：
from osgeo import gdal, osr
import os
import numpy as np
import cv2
from skimage import data
from skimage.feature import greycoprops
def export_tif(out_tif_name, var_lat, var_lon, NDVI, bandname=None):
# 判断数组维度
if len(NDVI.shape) > 2:
im_bands, im_height, im_width = NDVI.shape
else:
im_bands, (im_height, im_width) = 1, NDVI.shape
LonMin, LatMax, LonMax, LatMin = [min(var_lon), max(var_lat), max(var_lon), min(var_lat)]
# 分辨率计算
Lon_Res = (LonMax - LonMin) / (float(im_width) - 1)
Lat_Res = (LatMax - LatMin) / (float(im_height) - 1)
driver = gdal.GetDriverByName('GTiff')
out_tif = driver.Create(out_tif_name, im_width, im_height, im_bands, gdal.GDT_Float32) # 创建框架
# 设置影像的显示范围
# Lat_Res一定要是-的
geotransform = (LonMin, Lon_Res, 0, LatMax, 0, -Lat_Res)
out_tif.SetGeoTransform(geotransform)
# 获取地理坐标系统信息，用于选取需要的地理坐标系统
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326) # 定义输出的坐标系为"WGS 84"，AUTHORITY["EPSG","4326"]
out_tif.SetProjection(srs.ExportToWkt()) # 给新建图层赋予投影信息
# 数据写出
if im_bands == 1:
out_tif.GetRasterBand(1).WriteArray(NDVI)
else:
for bands in range(im_bands):
# 为每个波段设置名称
if bandname is not None:
out_tif.GetRasterBand(bands + 1).SetDescription(bandname[bands])
out_tif.GetRasterBand(bands + 1).WriteArray(NDVI[bands])
# 生成 ENVI HDR 文件
hdr_file = out_tif_name.replace('.tif', '.hdr')
with open(hdr_file, 'w') as f:
f.write('ENVI\n')
f.write('description = {Generated by export_tif}\n')
f.write('samples = {}\n'.format(im_width))
f.write('lines = {}\n'.format(im_height))
f.write('bands = {}\n'.format(im_bands))
f.write('header offset = 0\n')
f.write('file type = ENVI Standard\n')
f.write('data type = {}\n'.format(gdal.GetDataTypeName(out_tif.GetRasterBand(1).DataType)))
f.write('interleave = bsq\n')
f.write('sensor type = Unknown\n')
f.write('byte order = 0\n')
band_names_str = ', '.join(str(band) for band in bandname)
f.write('band names = {%s}\n'%(band_names_str))
del out_tif
def fast_glcm(img, vmin=0, vmax=255, levels=8, kernel_size=5, distance=1.0, angle=0.0):
'''
Parameters
----------
img: array_like, shape=(h,w), dtype=np.uint8
input image
vmin: int
minimum value of input image
vmax: int
maximum value of input image
levels: int
number of grey-levels of GLCM
kernel_size: int
Patch size to calculate GLCM around the target pixel
distance: float
pixel pair distance offsets [pixel] (1.0, 2.0, and etc.)
angle: float
pixel pair angles [degree] (0.0, 30.0, 45.0, 90.0, and etc.)
Returns
-------
Grey-level co-occurrence matrix for each pixels
shape = (levels, levels, h, w)
'''
mi, ma = vmin, vmax
ks = kernel_size
h,w = img.shape
# digitize
bins = np.linspace(mi, ma+1, levels+1)
gl1 = np.digitize(img, bins) - 1
# make shifted image
dx = distance*np.cos(np.deg2rad(angle))
dy = distance*np.sin(np.deg2rad(-angle))
mat = np.array([[1.0,0.0,-dx], [0.0,1.0,-dy]], dtype=np.float32)
gl2 = cv2.warpAffine(gl1, mat, (w,h), flags=cv2.INTER_NEAREST,
borderMode=cv2.BORDER_REPLICATE)
# make glcm
glcm = np.zeros((levels, levels, h, w), dtype=np.uint8)
for i in range(levels):
for j in range(levels):
mask = ((gl1==i) & (gl2==j))
glcm[i,j, mask] = 1
kernel = np.ones((ks, ks), dtype=np.uint8)
for i in range(levels):
for j in range(levels):
glcm[i,j] = cv2.filter2D(glcm[i,j], -1, kernel)
# 灰色关联矩阵归一化，之后结果与ENVI相同
glcm = glcm.astype(np.float32)
glcm_sums = np.apply_over_axes(np.sum, glcm, axes=(0, 1))
glcm_sums[glcm_sums == 0] = 1
glcm /= glcm_sums
return glcm
def fast_glcm_texture(img, vmin=0, vmax=255, levels=8, ks=5, distance=1.0, angle=0.0):
'''
calc glcm texture
'''
h,w = img.shape
glcm = fast_glcm(img, vmin, vmax, levels, ks, distance, angle)
mean = np.zeros((h,w), dtype=np.float32)
var = np.zeros((h,w), dtype=np.float32)
cont = np.zeros((h,w), dtype=np.float32)
diss = np.zeros((h,w), dtype=np.float32)
homo = np.zeros((h,w), dtype=np.float32)
asm = np.zeros((h,w), dtype=np.float32)
ent = np.zeros((h,w), dtype=np.float32)
cor = np.zeros((h, w), dtype=np.float32)
for i in range(levels):
for j in range(levels):
mean += glcm[i,j] * i
homo += glcm[i,j] / (1.+(i-j)**2)
asm += glcm[i,j]**2
cont += glcm[i,j] * (i-j)**2
diss += glcm[i,j] * np.abs(i-j)
ent -= glcm[i, j] * np.log10(glcm[i, j] + 0.00001)
for i in range(levels):
for j in range(levels):
var += glcm[i, j] * (i - mean)**2
greycoprops(glcm, 'correlation') # 上面计算的几个特征均可这样写
cor[cor == 1.0] = 0.
homo[homo == 1.] = 0
asm[asm == 1.] = 0
return [mean, var, cont, diss, homo, asm, ent, cor]
# 遥感影像
image = r'..\path\to\your\file\image.tif'
# 文件输出路径
out_path = r'..\output\path'
dataset = gdal.Open(image) # 读取数据
adfGeoTransform = dataset.GetGeoTransform()
nXSize = dataset.RasterXSize # 列数
nYSize = dataset.RasterYSize # 行数
im_data = dataset.ReadAsArray(0, 0, nXSize, nYSize) # 读取数据
im_data[im_data==65536] = 0
var_lat = [] # 纬度
var_lon = [] # 经度
for j in range(nYSize):
lat = adfGeoTransform[3] + j * adfGeoTransform[5]
var_lat.append(lat)
for i in range(nXSize):
lon = adfGeoTransform[0] + i * adfGeoTransform[1]
var_lon.append(lon)
result = []
band_name=[]
for i, img in enumerate(im_data):
img = np.uint8(255.0 * (img - np.min(img))/(np.max(img) - np.min(img))) # normalization
fast_result = fast_glcm_texture(img, vmin=0, vmax=255, levels=32, ks=3, distance=1.0, angle=0.0)
result += fast_result
variable_names = ['mean', 'variance', 'contrast', 'dissimilarity', 'homogeneity', 'ASM', 'entropy', 'correlation']
for names in variable_names:
band_name.append('band_'+str(i+1)+'_'+names)
name = os.path.splitext(os.path.basename(image))[0]
export_tif(os.path.join(out_path, '%s_TF.tif' % name), var_lat, var_lon, np.array(result), bandname=band_name)
print('over')
总结
目前熵值特征计算结果与ENVI有差异，不过只是数值差异，使用影像打开的结果显示一致，说明只是值的范围差异，不影响分析，其他特征均与ENVI的计算结果一致
