# coding: utf-8 # The code is written by Linghui import numpy as np from skimage import data from matplotlib import pyplot as plt import time from PIL import Image import spectral import os import warnings warnings.filterwarnings('ignore') # 导入GLCM计算模块 try: # 尝试相对导入 from . import get_glcm except ImportError: try: # 尝试绝对导入 import get_glcm except ImportError: # 如果都失败，尝试从当前目录导入 import sys current_dir = os.path.dirname(os.path.abspath(__file__)) if current_dir not in sys.path: sys.path.insert(0, current_dir) import get_glcm # 使用get_glcm模块进行GLCM计算 class SpatialFeatureConfig: """ 空间特征提取配置类 """ def __init__(self): # GLCM参数 self.nbit = 64 # gray levels self.mi = 0 # min gray self.ma = 255 # max gray self.slide_window = 7 # sliding window size # GLCM计算参数 self.step = [2] # step distances self.angle = [0] # angles in radians # 数据处理参数 self.band_index = 25 # 要处理的波段索引 # 文件路径 self.image_path = r"C:\Program Files\Spectronon3\_internal\examples\leaf_small.bip.hdr" # 输出设置 self.output_dir = "output" self.save_dat = True # 是否保存为dat文件 self.save_png = False # 是否保存为PNG文件 def validate(self): """验证配置参数的有效性""" if not os.path.exists(self.image_path): raise FileNotFoundError(f"图像文件不存在: {self.image_path}") if self.band_index < 0: raise ValueError("波段索引不能为负数") if self.nbit <= 0: raise ValueError("nbit必须大于0") if self.slide_window <= 0 or self.slide_window % 2 == 0: raise ValueError("slide_window必须为正奇数") return True def run_glcm_analysis(config): """ 执行GLCM纹理特征分析 参数: config: SpatialFeatureConfig对象 """ start = time.time() print('---------------0. Parameter Setting-----------------') # 验证配置 try: config.validate() print("配置验证通过") except Exception as e: print(f"配置验证失败: {e}") return print('-------------------1. Load Data---------------------') # 加载高光谱数据 img, h, w = load_hyperspectral_data(config.image_path, config.band_index) # 提取纹理特征 features_dict = extract_texture_features(img, config) print('---------------4. Display and Result----------------') print(f'使用的波段: {config.band_index}') print(f'GLCM参数: distances={config.step}, angles={[f"{a*180/np.pi:.0f}°" for a in config.angle]}') print(f'图像尺寸: {h}x{w}') print(f'处理时间: {time.time() - start:.2f} 秒') # 保存特征数据 if config.save_dat: output_filename = f"spatial_features_band{config.band_index}" output_path = os.path.join(config.output_dir, output_filename) save_features_to_dat(features_dict, output_path, config) print('GLCM纹理特征分析完成！') return features_dict 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 def save_features_to_dat(features_dict, output_path, config): """ 将纹理特征保存为dat文件和对应的hdr文件 参数: - features_dict: 包含各种纹理特征的字典 - output_path: 输出文件路径（不含扩展名） - config: 配置对象 """ try: # 确保输出目录存在 output_dir = os.path.dirname(output_path) if output_dir and not os.path.exists(output_dir): os.makedirs(output_dir) # 准备要保存的特征 feature_names = ['mean', 'variance', 'homogeneity', 'contrast', 'dissimilarity', 'entropy', 'energy', 'correlation', 'Auto_correlation'] # 将所有特征堆叠成多波段数据 feature_list = [] band_names = [] for feature_name in feature_names: if feature_name in features_dict: feature_data = features_dict[feature_name] feature_list.append(feature_data) band_names.append(feature_name) if not feature_list: print("没有找到要保存的特征数据") return # 堆叠所有特征为多波段数组 stacked_features = np.stack(feature_list, axis=-1) print(f"特征数据形状: {stacked_features.shape}") # 保存为dat文件（二进制格式） dat_path = f"{output_path}.dat" with open(dat_path, 'wb') as f: # 将数据转换为float32格式保存 stacked_features.astype(np.float32).tofile(f) print(f"特征数据已保存到: {dat_path}") # 保存对应的hdr文件 hdr_path = f"{output_path}.hdr" _save_hdr_file(hdr_path, stacked_features.shape, config, band_names) print(f"头文件已保存到: {hdr_path}") except Exception as e: print(f"保存dat文件时出错: {e}") def _save_hdr_file(hdr_path, data_shape, config, band_names): """ 保存ENVI头文件 参数: - hdr_path: 头文件路径 - data_shape: 数据形状 (height, width, bands) - config: 配置对象 - band_names: 波段名称列表 """ height, width, bands = data_shape # 计算数据范围（为每个波段计算） data_range_info = f"Data ranges vary per band (normalized 0-1)" header_content = f"""ENVI description = {{Spatial Texture Features [Band: {config.band_index}, Window: {config.slide_window}, Step: {config.step}, Angles: {[f"{a*180/np.pi:.0f}°" for a in config.angle]}]}} samples = {width} lines = {height} bands = {bands} header offset = 0 file type = ENVI Standard data type = 4 interleave = bip byte order = 0 band names = {{{', '.join(band_names)}}} data range = {{{data_range_info}}} GLCM_params = {{nbit: {config.nbit}, window: {config.slide_window}, step: {config.step}, angles: {config.angle}}} """ with open(hdr_path, 'w', encoding='utf-8') as f: f.write(header_content) def extract_texture_features(img, config): """ 提取纹理特征 参数: - img: 输入图像 - config: 配置对象 返回: - features_dict: 包含各种纹理特征的字典 """ print('------------------2. 计算GLCM特征---------------------') img = img.squeeze() # 使用get_glcm计算GLCM特征 glcm = get_glcm.calcu_glcm(img, config.mi, config.ma, config.nbit, config.slide_window, config.step, config.angle) print('GLCM计算完成！') print('-----------------3. 提取纹理特征-------------------') # 提取纹理特征 # 使用最后一个step和angle的特征 i, j = len(config.step)-1, len(config.angle)-1 glcm_cut = glcm[:, :, i, j, :, :] features_dict = { 'image_shape': img.shape, 'mean': get_glcm.calcu_glcm_mean(glcm_cut, config.nbit), 'variance': get_glcm.calcu_glcm_variance(glcm_cut, config.nbit), 'homogeneity': get_glcm.calcu_glcm_homogeneity(glcm_cut, config.nbit), 'contrast': get_glcm.calcu_glcm_contrast(glcm_cut, config.nbit), 'dissimilarity': get_glcm.calcu_glcm_dissimilarity(glcm_cut, config.nbit), 'entropy': get_glcm.calcu_glcm_entropy(glcm_cut, config.nbit), 'energy': get_glcm.calcu_glcm_energy(glcm_cut, config.nbit), 'correlation': get_glcm.calcu_glcm_correlation(glcm_cut, config.nbit), 'Auto_correlation': get_glcm.calcu_glcm_Auto_correlation(glcm_cut, config.nbit) } return features_dict def main(): """主函数：执行空间特征提取流程""" start = time.time() print('---------------0. Parameter Setting-----------------') # 创建配置对象 config = SpatialFeatureConfig() # 可以在这里修改配置参数 # config.band_index = 30 # 修改波段索引 # config.step = [2, 4, 8] # 修改步长 # config.angle = [0, np.pi/4, np.pi/2, np.pi*3/4] # 修改角度 # 验证配置 try: config.validate() print("配置验证通过") except Exception as e: print(f"配置验证失败: {e}") return print('-------------------1. Load Data---------------------') # 加载高光谱数据 img, h, w = load_hyperspectral_data(config.image_path, config.band_index) # 提取纹理特征 features_dict = extract_texture_features(img, config) print('---------------4. Display and Result----------------') print(f'使用的波段: {config.band_index}') print(f'GLCM参数: distances={config.step}, angles={[f"{a*180/np.pi:.0f}°" for a in config.angle]}') print(f'图像尺寸: {h}x{w}') print(f'处理时间: {time.time() - start:.2f} 秒') # 保存特征数据 if config.save_dat: output_filename = f"spatial_features_band{config.band_index}" output_path = os.path.join(config.output_dir, output_filename) save_features_to_dat(features_dict, output_path, config) print('处理完成！') if __name__ == '__main__': main()