216 lines
7.0 KiB
Python
216 lines
7.0 KiB
Python
|
# coding: utf-8
|
|||
|
|
# The code is written by Linghui
|
|||
|
|
|
|||
|
|
import numpy as np
|
|||
|
|
from skimage import data
|
|||
|
|
from matplotlib import pyplot as plt
|
|||
|
|
import get_glcm
|
|||
|
|
import time
|
|||
|
|
from PIL import Image
|
|||
|
|
import spectral
|
|||
|
|
import os
|
|||
|
|
import warnings
|
|||
|
|
warnings.filterwarnings('ignore')
|
|||
|
|
|
|||
|
|
# 使用get_glcm模块进行GLCM计算
|
|||
|
|
|
|||
|
|
def main():
|
|||
|
|
pass
|
|||
|
|
|
|||
|
|
|
|||
|
|
def load_hyperspectral_data(image_path, band_index=None):
|
|||
|
|
"""
|
|||
|
|
加载高光谱数据并可选地提取指定波段
|
|||
|
|
|
|||
|
|
参数:
|
|||
|
|
- image_path: 高光谱图像文件路径
|
|||
|
|
- band_index: 要提取的波段索引，如果为None则返回整个图像
|
|||
|
|
|
|||
|
|
返回:
|
|||
|
|
- img: 处理后的图像数据
|
|||
|
|
- h, w: 图像尺寸
|
|||
|
|
"""
|
|||
|
|
print(f'正在读取高光谱图像: {image_path}')
|
|||
|
|
|
|||
|
|
try:
|
|||
|
|
# 使用spectral库读取ENVI格式
|
|||
|
|
img_obj = spectral.open_image(image_path)
|
|||
|
|
hyperspectral_data = img_obj.load()
|
|||
|
|
print(f'使用spectral库读取成功，数据形状: {hyperspectral_data.shape}')
|
|||
|
|
except Exception as e:
|
|||
|
|
print(f'spectral库读取失败: {e}')
|
|||
|
|
raise ValueError(f"无法读取图像文件: {image_path}")
|
|||
|
|
|
|||
|
|
# 检查数据维度
|
|||
|
|
if len(hyperspectral_data.shape) != 3:
|
|||
|
|
raise ValueError(f"高光谱数据应该是3维的，当前形状: {hyperspectral_data.shape}")
|
|||
|
|
|
|||
|
|
h, w, bands = hyperspectral_data.shape
|
|||
|
|
print(f'图像尺寸: {h}x{w}, 波段数: {bands}')
|
|||
|
|
|
|||
|
|
# 如果指定了波段索引，提取该波段
|
|||
|
|
if band_index is not None:
|
|||
|
|
if band_index < 0 or band_index >= bands:
|
|||
|
|
raise ValueError(f"波段索引 {band_index} 超出范围 [0, {bands-1}]")
|
|||
|
|
|
|||
|
|
print(f'提取波段 {band_index}...')
|
|||
|
|
img = hyperspectral_data[:, :, band_index].astype(np.float32)
|
|||
|
|
|
|||
|
|
# 归一化到0-255范围
|
|||
|
|
img_min, img_max = img.min(), img.max()
|
|||
|
|
if img_max > img_min:
|
|||
|
|
img = np.uint8(255.0 * (img - img_min) / (img_max - img_min))
|
|||
|
|
else:
|
|||
|
|
img = np.zeros_like(img, dtype=np.uint8)
|
|||
|
|
|
|||
|
|
print(f'提取的波段数据范围: [{img_min:.2f}, {img_max:.2f}] -> [0, 255]')
|
|||
|
|
else:
|
|||
|
|
# 如果没有指定波段，返回第一个波段作为示例
|
|||
|
|
print('未指定波段，使用第一个波段...')
|
|||
|
|
img = hyperspectral_data[:, :, 0].astype(np.float32)
|
|||
|
|
img_min, img_max = img.min(), img.max()
|
|||
|
|
if img_max > img_min:
|
|||
|
|
img = np.uint8(255.0 * (img - img_min) / (img_max - img_min))
|
|||
|
|
else:
|
|||
|
|
img = np.zeros_like(img, dtype=np.uint8)
|
|||
|
|
|
|||
|
|
return img, h, w
|
|||
|
|
|
|||
|
|
|
|||
|
|
if __name__ == '__main__':
|
|||
|
|
|
|||
|
|
main()
|
|||
|
|
|
|||
|
|
start = time.time()
|
|||
|
|
|
|||
|
|
print('---------------0. Parameter Setting-----------------')
|
|||
|
|
nbit = 64 # gray levels
|
|||
|
|
mi, ma = 0, 255 # max gray and min gray
|
|||
|
|
slide_window = 7 # sliding window
|
|||
|
|
# step = [2, 4, 8, 16] # step
|
|||
|
|
# angle = [0, np.pi/4, np.pi/2, np.pi*3/4] # angle or direction
|
|||
|
|
step = [2]
|
|||
|
|
angle = [0]
|
|||
|
|
|
|||
|
|
# 指定波段索引（可以修改）
|
|||
|
|
band_index = 25 # 要处理的波段索引，从0开始
|
|||
|
|
|
|||
|
|
print('-------------------1. Load Data---------------------')
|
|||
|
|
# 修改为高光谱图像路径
|
|||
|
|
image_path = r"C:\Program Files\Spectronon3\_internal\examples\leaf_small.bip.hdr"
|
|||
|
|
# image_path = r"./test.tif" # 备用的TIFF图像路径
|
|||
|
|
|
|||
|
|
# 检查文件是否存在
|
|||
|
|
if not os.path.exists(image_path):
|
|||
|
|
print(f"警告: 指定的图像文件不存在: {image_path}")
|
|||
|
|
print("请确保文件路径正确，或使用其他可用的图像文件")
|
|||
|
|
exit(1)
|
|||
|
|
|
|||
|
|
# 加载高光谱数据
|
|||
|
|
img, h, w = load_hyperspectral_data(image_path, band_index)
|
|||
|
|
print('------------------2. 计算GLCM特征---------------------')
|
|||
|
|
img = img.squeeze()
|
|||
|
|
# 使用get_glcm计算GLCM特征
|
|||
|
|
glcm = get_glcm.calcu_glcm(img, mi, ma, nbit, slide_window, step, angle)
|
|||
|
|
|
|||
|
|
print('GLCM计算完成！')
|
|||
|
|
print('-----------------3. 提取纹理特征-------------------')
|
|||
|
|
|
|||
|
|
# 提取纹理特征用于显示
|
|||
|
|
# 只计算最后一个step和angle的特征用于显示
|
|||
|
|
i, j = len(step)-1, len(angle)-1 # 使用最后一个配置
|
|||
|
|
glcm_cut = glcm[:, :, i, j, :, :]
|
|||
|
|
mean = get_glcm.calcu_glcm_mean(glcm_cut, nbit)
|
|||
|
|
variance = get_glcm.calcu_glcm_variance(glcm_cut, nbit)
|
|||
|
|
homogeneity = get_glcm.calcu_glcm_homogeneity(glcm_cut, nbit)
|
|||
|
|
contrast = get_glcm.calcu_glcm_contrast(glcm_cut, nbit)
|
|||
|
|
dissimilarity = get_glcm.calcu_glcm_dissimilarity(glcm_cut, nbit)
|
|||
|
|
entropy = get_glcm.calcu_glcm_entropy(glcm_cut, nbit)
|
|||
|
|
energy = get_glcm.calcu_glcm_energy(glcm_cut, nbit)
|
|||
|
|
correlation = get_glcm.calcu_glcm_correlation(glcm_cut, nbit)
|
|||
|
|
Auto_correlation = get_glcm.calcu_glcm_Auto_correlation(glcm_cut, nbit)
|
|||
|
|
print('---------------4. Display and Result----------------')
|
|||
|
|
print(f'使用的波段: {band_index}')
|
|||
|
|
print(f'GLCM参数: distances={step}, angles={[f"{a*180/np.pi:.0f}°" for a in angle]}')
|
|||
|
|
|
|||
|
|
plt.figure(figsize=(10, 4.5))
|
|||
|
|
font = {'family' : 'Times New Roman',
|
|||
|
|
'weight' : 'normal',
|
|||
|
|
'size' : 12,
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
plt.subplot(2,5,1)
|
|||
|
|
plt.tick_params(labelbottom=False, labelleft=False)
|
|||
|
|
plt.axis('off')
|
|||
|
|
plt.imshow(img, cmap ='gray')
|
|||
|
|
plt.title(f'Band {band_index}', font)
|
|||
|
|
|
|||
|
|
plt.subplot(2,5,2)
|
|||
|
|
plt.tick_params(labelbottom=False, labelleft=False)
|
|||
|
|
plt.axis('off')
|
|||
|
|
plt.imshow(mean, cmap ='gray')
|
|||
|
|
plt.title('Mean', font)
|
|||
|
|
|
|||
|
|
plt.subplot(2,5,3)
|
|||
|
|
plt.tick_params(labelbottom=False, labelleft=False)
|
|||
|
|
plt.axis('off')
|
|||
|
|
plt.imshow(variance, cmap ='gray')
|
|||
|
|
plt.title('Variance', font)
|
|||
|
|
|
|||
|
|
plt.subplot(2,5,4)
|
|||
|
|
plt.tick_params(labelbottom=False, labelleft=False)
|
|||
|
|
plt.axis('off')
|
|||
|
|
plt.imshow(homogeneity, cmap ='gray')
|
|||
|
|
plt.title('Homogeneity', font)
|
|||
|
|
|
|||
|
|
plt.subplot(2,5,5)
|
|||
|
|
plt.tick_params(labelbottom=False, labelleft=False)
|
|||
|
|
plt.axis('off')
|
|||
|
|
plt.imshow(contrast, cmap ='gray')
|
|||
|
|
plt.title('Contrast', font)
|
|||
|
|
|
|||
|
|
plt.subplot(2,5,6)
|
|||
|
|
plt.tick_params(labelbottom=False, labelleft=False)
|
|||
|
|
plt.axis('off')
|
|||
|
|
plt.imshow(dissimilarity, cmap ='gray')
|
|||
|
|
plt.title('Dissimilarity', font)
|
|||
|
|
|
|||
|
|
plt.subplot(2,5,7)
|
|||
|
|
plt.tick_params(labelbottom=False, labelleft=False)
|
|||
|
|
plt.axis('off')
|
|||
|
|
plt.imshow(entropy, cmap ='gray')
|
|||
|
|
plt.title('Entropy', font)
|
|||
|
|
|
|||
|
|
plt.subplot(2,5,8)
|
|||
|
|
plt.tick_params(labelbottom=False, labelleft=False)
|
|||
|
|
plt.axis('off')
|
|||
|
|
plt.imshow(energy, cmap ='gray')
|
|||
|
|
plt.title('Energy', font)
|
|||
|
|
|
|||
|
|
plt.subplot(2,5,9)
|
|||
|
|
plt.tick_params(labelbottom=False, labelleft=False)
|
|||
|
|
plt.axis('off')
|
|||
|
|
plt.imshow(correlation, cmap ='gray')
|
|||
|
|
plt.title('Correlation', font)
|
|||
|
|
|
|||
|
|
plt.subplot(2,5,10)
|
|||
|
|
plt.tick_params(labelbottom=False, labelleft=False)
|
|||
|
|
plt.axis('off')
|
|||
|
|
plt.imshow(Auto_correlation, cmap ='gray')
|
|||
|
|
plt.title('Auto Correlation', font)
|
|||
|
|
|
|||
|
|
plt.tight_layout(pad=0.5)
|
|||
|
|
|
|||
|
|
# 创建输出文件名，包含波段信息
|
|||
|
|
output_filename = f'GLCM_Features_Band{band_index}.png'
|
|||
|
|
plt.savefig(output_filename
|
|||
|
|
, format='png'
|
|||
|
|
, bbox_inches = 'tight'
|
|||
|
|
, pad_inches = 0
|
|||
|
|
, dpi=300)
|
|||
|
|
print(f'结果图像已保存到: {output_filename}')
|
|||
|
|
|
|||
|
|
|
|||
|
|
end = time.time()
|
|||
|
|
print('Code run time:', end - start)
|