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)
|