# coding: utf-8 # Shape Feature Analysis Module # 用于分析二值化图像的形状特征 import numpy as np import cv2 import pandas as pd from skimage import measure, morphology from skimage.segmentation import watershed from skimage.feature import peak_local_max from scipy import ndimage import spectral import os import warnings from shapely.geometry import shape import geopandas as gpd import json warnings.filterwarnings('ignore') class ShapeFeatureConfig: """ 形状特征分析配置类 """ def __init__(self): # 输入数据设置 self.input_type = 'dat' # 'dat' 或 'shp' self.dat_file_path = None # dat文件路径 self.hdr_file_path = None # 对应的hdr文件路径 self.shp_file_path = None # shp文件路径 self.band_index = 0 # 要处理的波段索引（对于dat文件） # 连通域处理参数 self.connectivity = 8 # 连通性：4或8 self.min_area = 10 # 最小区域面积阈值 self.use_watershed = False # 是否使用分水岭算法分离相邻物体 # 分水岭算法参数 self.watershed_min_distance = 10 # 局部最大值的最小距离 # 输出设置 self.output_dir = "output" self.output_csv = True # 是否输出CSV文件 self.save_labeled_image = False # 是否保存标记后的图像 def validate(self): """验证配置参数的有效性""" if self.input_type not in ['dat', 'shp']: raise ValueError("input_type必须是'dat'或'shp'") if self.input_type == 'dat': if not self.dat_file_path or not os.path.exists(self.dat_file_path): raise FileNotFoundError(f"dat文件不存在: {self.dat_file_path}") if self.hdr_file_path and not os.path.exists(self.hdr_file_path): raise FileNotFoundError(f"hdr文件不存在: {self.hdr_file_path}") elif self.input_type == 'shp': if not self.shp_file_path or not os.path.exists(self.shp_file_path): raise FileNotFoundError(f"shp文件不存在: {self.shp_file_path}") if self.connectivity not in [4, 8]: raise ValueError("connectivity必须是4或8") if self.min_area < 1: raise ValueError("min_area必须大于等于1") return True def load_binary_image_from_dat(dat_path, hdr_path=None, band_index=0): """ 从dat文件加载二值化图像数据 参数: - dat_path: dat文件路径 - hdr_path: 对应的hdr文件路径（可选） - band_index: 要读取的波段索引 返回: - binary_image: 二值化图像数组 - metadata: 元数据字典 """ try: # 如果提供了hdr文件路径，使用spectral库读取ENVI格式 if hdr_path and os.path.exists(hdr_path): print(f"使用hdr文件: {hdr_path}") img_obj = spectral.open_image(hdr_path) image_data = img_obj.load() print(f"通过spectral库读取成功，形状: {image_data.shape}") # 如果没有hdr文件或读取失败，直接读取dat文件 else: print(f"直接读取dat文件: {dat_path}") # 从dat文件路径推断hdr文件路径 if not hdr_path: hdr_path = dat_path.replace('.dat', '.hdr') if os.path.exists(hdr_path): # 读取hdr文件获取元数据 with open(hdr_path, 'r') as f: hdr_content = f.read() # 解析hdr文件 lines = hdr_content.split('\n') metadata = {} for line in lines: line = line.strip() if '=' in line: key, value = line.split('=', 1) key = key.strip() value = value.strip() if key == 'samples': metadata['samples'] = int(value) elif key == 'lines': metadata['lines'] = int(value) elif key == 'bands': metadata['bands'] = int(value) elif key == 'data type': metadata['data_type'] = int(value) print(f"从hdr文件解析的元数据: {metadata}") # 读取dat文件 with open(dat_path, 'rb') as f: if metadata['data_type'] == 1: # 8-bit unsigned integer image_data = np.frombuffer(f.read(), dtype=np.uint8) else: raise ValueError(f"不支持的数据类型: {metadata['data_type']}") # 重塑为正确的形状 (lines, samples, bands) image_data = image_data.reshape(metadata['lines'], metadata['samples'], metadata['bands']) print(f"重塑后的数据形状: {image_data.shape}") else: # 如果没有hdr文件，使用spectral库尝试读取 try: img_obj = spectral.open_image(dat_path) image_data = img_obj.load() print(f"通过spectral库读取成功，形状: {image_data.shape}") except: raise FileNotFoundError(f"找不到对应的hdr文件: {hdr_path}") # 处理波段选择 if len(image_data.shape) == 3: if band_index >= image_data.shape[2]: raise ValueError(f"波段索引 {band_index} 超出范围 [0, {image_data.shape[2]-1}]") binary_image = image_data[:, :, band_index] else: binary_image = image_data # 确保是二值图像 binary_image = binary_image.astype(np.uint8) # 对于分割结果，通常已经是二值的，但确保值范围正确 unique_values = np.unique(binary_image) print(f"图像中的唯一值: {unique_values}") # 如果图像值范围不是0-1，转换为二值（假设非零值是前景） if binary_image.max() > 1: binary_image = (binary_image > 0).astype(np.uint8) print("已转换为二值图像") metadata = { 'shape': binary_image.shape, 'dtype': binary_image.dtype, 'source': 'dat', 'band_index': band_index, 'unique_values': unique_values.tolist() } print(f"最终图像形状: {binary_image.shape}, 数据类型: {binary_image.dtype}") return binary_image, metadata except Exception as e: print(f"读取dat文件失败: {e}") import traceback traceback.print_exc() raise def load_regions_from_shp(shp_path, image_shape=None): """ 从shp文件加载区域信息并转换为二值图像 参数: - shp_path: shp文件路径 - image_shape: 目标图像形状 (height, width)，用于创建二值图像 返回: - binary_image: 二值化图像数组 - metadata: 元数据字典 """ try: # 读取shp文件 gdf = gpd.read_file(shp_path) print(f"成功读取shp文件，包含 {len(gdf)} 个几何对象") if image_shape is None: # 如果没有指定图像形状，尝试从bounds推断 bounds = gdf.total_bounds # [minx, miny, maxx, maxy] # 这里需要根据实际坐标系进行转换，暂时使用简单的假设 height = int(bounds[3] - bounds[1]) width = int(bounds[2] - bounds[0]) if height <= 0 or width <= 0: raise ValueError("无法从shp文件推断图像尺寸，请提供image_shape参数") image_shape = (height, width) # 创建二值图像 binary_image = np.zeros(image_shape, dtype=np.uint8) # 将每个几何对象转换为像素坐标并填充 for idx, geom in enumerate(gdf.geometry): if geom is not None: # 这里需要坐标转换逻辑，暂时使用简化的实现 # 实际实现需要根据shp文件的坐标系进行准确转换 try: # 简化的几何转换（实际使用时需要更复杂的坐标变换） coords = np.array(geom.exterior.coords) if hasattr(geom, 'exterior') else np.array(geom.coords) # 归一化坐标到图像尺寸（这只是示例，实际需要正确的地理坐标转换） coords_normalized = coords.copy() if len(coords) > 0: # 创建标签图像，每个区域用不同的值标记 label_value = idx + 1 # 这里应该实现正确的几何到像素的转换 # 暂时用占位符实现 pass except Exception as e: print(f"处理几何对象 {idx} 时出错: {e}") continue # 如果没有有效的几何转换，创建一个示例图像 if binary_image.max() == 0: print("警告：未能从shp文件转换为有效的二值图像，使用示例数据") binary_image[100:200, 100:200] = 1 # 示例矩形区域 metadata = { 'shape': binary_image.shape, 'dtype': binary_image.dtype, 'source': 'shp', 'num_features': len(gdf) } return binary_image, metadata except Exception as e: print(f"读取shp文件失败: {e}") raise def apply_watershed_segmentation(binary_image, config): """ 使用分水岭算法分离相邻物体 参数: - binary_image: 二值化图像 - config: 配置对象 返回: - labeled_image: 分水岭分割后的标记图像 - num_objects: 物体数量 """ try: # 计算距离变换 distance = ndimage.distance_transform_edt(binary_image) # 找到局部最大值作为种子点 local_maxi = peak_local_max(distance, indices=False, footprint=np.ones((3, 3)), labels=binary_image, min_distance=config.watershed_min_distance) # 标记种子点 markers = ndimage.label(local_maxi)[0] # 应用分水岭算法 labels = watershed(-distance, markers, mask=binary_image) # 重新标记以确保连续的标签 unique_labels = np.unique(labels) unique_labels = unique_labels[unique_labels > 0] # 排除背景 relabeled_image = np.zeros_like(labels) for new_label, old_label in enumerate(unique_labels, 1): relabeled_image[labels == old_label] = new_label num_objects = len(unique_labels) print(f"分水岭分割完成，识别出 {num_objects} 个物体") return relabeled_image.astype(np.int32), num_objects except Exception as e: print(f"分水岭分割失败: {e}") raise def label_connected_components(binary_image, config): """ 标记连通域，可选择使用分水岭算法 参数: - binary_image: 二值化图像 - config: 配置对象 返回: - labeled_image: 标记后的图像 - num_objects: 物体数量 """ try: if config.use_watershed: # 使用分水岭算法 labeled_image, num_objects = apply_watershed_segmentation(binary_image, config) else: # 使用标准连通域标记 if config.connectivity == 8: structure = np.ones((3, 3), dtype=int) else: # connectivity == 4 structure = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]], dtype=int) labeled_image, num_objects = ndimage.label(binary_image, structure=structure) # 过滤小区域 if config.min_area > 1: component_sizes = np.bincount(labeled_image.ravel()) too_small = component_sizes < config.min_area too_small_mask = too_small[labeled_image] labeled_image[too_small_mask] = 0 # 重新标记以保持连续性 unique_labels = np.unique(labeled_image) unique_labels = unique_labels[unique_labels > 0] relabeled_image = np.zeros_like(labeled_image) for new_label, old_label in enumerate(unique_labels, 1): relabeled_image[labeled_image == old_label] = new_label labeled_image = relabeled_image num_objects = len(unique_labels) print(f"连通域标记完成，找到 {num_objects} 个有效区域") return labeled_image, num_objects except Exception as e: print(f"连通域标记失败: {e}") raise def calculate_shape_features(labeled_image, config): """ 计算所有形状特征指标 参数: - labeled_image: 标记后的图像 - config: 配置对象 返回: - features_df: 包含所有特征的DataFrame - contour_coords_dict: 轮廓坐标字典 """ try: # 获取区域属性 regions = measure.regionprops(labeled_image) if len(regions) == 0: print("警告：没有找到有效的区域") return pd.DataFrame(), {} features_list = [] contour_coords_dict = {} for region in regions: try: # 基本属性 label = region.label area = region.area perimeter = region.perimeter # 边界框和最小外接矩形 min_row, min_col, max_row, max_col = region.bbox bbox_width = max_col - min_col bbox_height = max_row - min_row bbox_area = bbox_width * bbox_height # 获取轮廓坐标 coords = region.coords # 像素坐标 # 计算各种形状特征 feature_dict = { 'label': label, 'area': area, 'perimeter': perimeter, 'bbox_width': bbox_width, 'bbox_height': bbox_height, 'bbox_area': bbox_area } # 1. 圆形度 (Circularity) = (4 * π * 面积) / (周长²) if perimeter > 0: circularity = (4 * np.pi * area) / (perimeter ** 2) else: circularity = 0 feature_dict['circularity'] = circularity # 2. 矩形度 (Rectangularity) = 面积 / 外接矩形面积 rectangularity = area / bbox_area if bbox_area > 0 else 0 feature_dict['rectangularity'] = rectangularity # 3. 长宽比 (Aspect Ratio) = 宽度 / 高度 aspect_ratio = bbox_width / bbox_height if bbox_height > 0 else 0 feature_dict['aspect_ratio'] = aspect_ratio # 4. 紧密度 (Compactness) = (周长²) / 面积 compactness = (perimeter ** 2) / area if area > 0 else 0 feature_dict['compactness'] = compactness # 5. 偏心率 (Eccentricity) - 使用regionprops的eccentricity eccentricity = region.eccentricity feature_dict['eccentricity'] = eccentricity # 6. 形状矩 (Shape Moments) - Hu矩 if hasattr(region, 'moments_hu'): hu_moments = region.moments_hu for i, hu_moment in enumerate(hu_moments): feature_dict[f'hu_moment_{i+1}'] = hu_moment else: # 如果没有Hu矩，计算为0 for i in range(7): feature_dict[f'hu_moment_{i+1}'] = 0 # 7. 凸性 (Convexity) = 面积 / 凸包面积 convex_area = region.convex_area convexity = area / convex_area if convex_area > 0 else 0 feature_dict['convexity'] = convexity # 8. 固体度 (Solidity) = 面积 / 凸包面积 (与凸性相同) solidity = convexity # 在scikit-image中，solidity就是这个定义 feature_dict['solidity'] = solidity # 9. 轮廓矩 (Contour Moments) - 使用中心矩 if hasattr(region, 'moments_central'): central_moments = region.moments_central feature_dict['central_moment_02'] = central_moments[0, 2] # μ02 feature_dict['central_moment_20'] = central_moments[2, 0] # μ20 feature_dict['central_moment_11'] = central_moments[1, 1] # μ11 else: feature_dict['central_moment_02'] = 0 feature_dict['central_moment_20'] = 0 feature_dict['central_moment_11'] = 0 # 10. 质心坐标 centroid_row, centroid_col = region.centroid feature_dict['centroid_row'] = centroid_row feature_dict['centroid_col'] = centroid_col # 11. 主要轴长度和角度 if hasattr(region, 'major_axis_length'): feature_dict['major_axis_length'] = region.major_axis_length feature_dict['minor_axis_length'] = region.minor_axis_length feature_dict['orientation'] = region.orientation else: feature_dict['major_axis_length'] = 0 feature_dict['minor_axis_length'] = 0 feature_dict['orientation'] = 0 # 12. 边界坐标序列（将保存到单独的JSON文件中） coords_list = coords.tolist() coord_key = f"region_{region.label}_coords" feature_dict['contour_coords'] = coord_key # 引用标识符 # 将轮廓坐标添加到字典中 contour_coords_dict[coord_key] = coords_list # 13. 其他有用的特征 feature_dict['equivalent_diameter'] = region.equivalent_diameter feature_dict['extent'] = region.extent # 面积与边界框面积之比 features_list.append(feature_dict) except Exception as e: print(f"计算区域 {region.label} 的特征时出错: {e}") continue # 创建DataFrame features_df = pd.DataFrame(features_list) print(f"成功计算 {len(features_df)} 个区域的形状特征") return features_df, contour_coords_dict except Exception as e: print(f"计算形状特征失败: {e}") raise def save_features_to_csv(features_df, output_path, config, contour_coords_dict=None): """ 将形状特征保存为CSV文件 参数: - features_df: 特征DataFrame - output_path: 输出文件路径（不含扩展名） - config: 配置对象 - contour_coords_dict: 轮廓坐标字典（可选） """ try: # 确保输出目录存在 output_dir = os.path.dirname(output_path) if output_dir and not os.path.exists(output_dir): os.makedirs(output_dir) # 如果有轮廓坐标数据，保存到单独的JSON文件 if contour_coords_dict is not None and len(contour_coords_dict) > 0: coords_json_path = f"{output_path}_contour_coords.json" with open(coords_json_path, 'w', encoding='utf-8') as f: json.dump(contour_coords_dict, f, indent=2, ensure_ascii=False) print(f"轮廓坐标已保存到: {coords_json_path}") # 保存为CSV文件 csv_path = f"{output_path}.csv" features_df.to_csv(csv_path, index=False, encoding='utf-8') print(f"形状特征已保存到: {csv_path}") print(f"共保存 {len(features_df)} 个区域的特征数据") # 打印特征统计信息 if len(features_df) > 0: numeric_columns = features_df.select_dtypes(include=[np.number]).columns print(f"数值特征列数: {len(numeric_columns)}") print(f"特征列: {list(features_df.columns)}") return csv_path except Exception as e: print(f"保存CSV文件失败: {e}") raise def save_labeled_image(labeled_image, output_path, config): """ 保存标记后的图像（可选） 参数: - labeled_image: 标记后的图像 - output_path: 输出文件路径（不含扩展名） - config: 配置对象 """ try: # 保存为PNG图像用于可视化 png_path = f"{output_path}_labeled.png" # 归一化标签图像以便可视化 if labeled_image.max() > 0: normalized_image = (labeled_image / labeled_image.max() * 255).astype(np.uint8) else: normalized_image = labeled_image.astype(np.uint8) # 创建彩色标签图像 colored_image = cv2.applyColorMap(normalized_image, cv2.COLORMAP_JET) # 保存图像 cv2.imwrite(png_path, colored_image) print(f"标记图像已保存到: {png_path}") # 同时保存为dat文件（ENVI格式） dat_path = f"{output_path}_labeled.dat" with open(dat_path, 'wb') as f: labeled_image.astype(np.int32).tofile(f) # 保存对应的hdr文件 hdr_path = f"{output_path}_labeled.hdr" header_content = f"""ENVI description = {{Labeled Image from Shape Analysis [Watershed: {config.use_watershed}]}} samples = {labeled_image.shape[1]} lines = {labeled_image.shape[0]} bands = 1 header offset = 0 file type = ENVI Standard data type = 3 interleave = bsq byte order = 0 band names = {{Labeled_Result}} classes = {labeled_image.max() + 1} class names = {{Background{' , '.join([f'Region_{i}' for i in range(1, labeled_image.max() + 1)])}}} """ with open(hdr_path, 'w', encoding='utf-8') as f: f.write(header_content) print(f"ENVI格式标记数据已保存到: {dat_path}") except Exception as e: print(f"保存标记图像失败: {e}") def analyze_shape_features(config): """ 执行形状特征分析的主函数 参数: - config: 配置对象 返回: - features_df: 特征DataFrame - labeled_image: 标记后的图像 """ try: print("=" * 60) print("开始形状特征分析") print("=" * 60) # 1. 加载数据 print("1. 加载输入数据...") if config.input_type == 'dat': binary_image, metadata = load_binary_image_from_dat( config.dat_file_path, config.hdr_file_path, config.band_index ) elif config.input_type == 'shp': binary_image, metadata = load_regions_from_shp( config.shp_file_path, image_shape=None ) else: raise ValueError(f"不支持的输入类型: {config.input_type}") print(f"数据加载完成，图像形状: {binary_image.shape}") binary_image = binary_image.squeeze() # 2. 连通域标记 print("2. 执行连通域标记...") labeled_image, num_objects = label_connected_components(binary_image, config) print(f"找到 {num_objects} 个连通域") if num_objects == 0: print("警告：没有找到有效的连通域") return pd.DataFrame(), labeled_image # 3. 计算形状特征 print("3. 计算形状特征...") features_df, contour_coords_dict = calculate_shape_features(labeled_image, config) if len(features_df) == 0: print("警告：未能计算出有效的形状特征") return features_df, labeled_image # 4. 保存结果 print("4. 保存分析结果...") # 生成输出文件名 timestamp = pd.Timestamp.now().strftime("%Y%m%d_%H%M%S") base_name = f"shape_features_{config.input_type}" if config.input_type == 'dat': base_name += f"_band{config.band_index}" if config.use_watershed: base_name += "_watershed" base_name += f"_{timestamp}" output_path = os.path.join(config.output_dir, base_name) # 保存CSV文件 if config.output_csv: csv_path = save_features_to_csv(features_df, output_path, config, contour_coords_dict) # 保存标记图像 if config.save_labeled_image: save_labeled_image(labeled_image, output_path, config) print("=" * 60) print("形状特征分析完成！") print(f"分析了 {len(features_df)} 个区域") print(f"计算了 {len(features_df.columns)} 个特征指标") if config.output_csv: print(f"结果已保存到: {csv_path}") print("=" * 60) return features_df, labeled_image except Exception as e: print(f"形状特征分析失败: {e}") raise def main(): """主函数：命令行接口""" import argparse parser = argparse.ArgumentParser( description='形状特征分析工具 - 从分割结果提取形状特征', formatter_class=argparse.RawDescriptionHelpFormatter, epilog=""" 使用示例: 1. 从dat文件分析形状特征: python shape_feature.py input.dat -t dat -b 0 -o shape_features 2. 从shp文件分析形状特征: python shape_feature.py input.shp -t shp -o shape_features 3. 启用分水岭算法分离物体: python shape_feature.py input.dat -t dat --use-watershed --watershed-min-distance 15 -o features 4. 自定义连通性和面积阈值: python shape_feature.py input.dat -t dat -c 4 --min-area 100 --save-labeled-image -o results 输入格式: - dat: ENVI格式的分割结果文件 (.dat + .hdr) - shp: Shapefile格式的矢量数据 形状特征包括: - 面积、周长、紧密度 - 圆度、矩形度、伸长度 - 质心坐标、边界框 - 方向角、主轴长度 - Hu矩不变量等几何特征 """ ) parser.add_argument('input_path', help='输入文件路径 (.dat分割结果 或 .shp矢量文件)') parser.add_argument('-t', '--input-type', choices=['dat', 'shp'], default='dat', help='输入文件类型 (默认: dat)') parser.add_argument('-o', '--output', required=True, help='输出文件路径前缀') # dat文件参数 parser.add_argument('-b', '--band-index', type=int, default=0, help='dat文件的波段索引 (默认: 0)') parser.add_argument('--hdr-path', help='dat文件对应的hdr文件路径 (默认自动查找)') # 连通域参数 parser.add_argument('-c', '--connectivity', type=int, choices=[4, 8], default=8, help='连通性，4或8 (默认: 8)') parser.add_argument('--min-area', type=int, default=10, help='最小区域面积阈值 (默认: 10)') # 分水岭参数 parser.add_argument('--use-watershed', action='store_true', help='启用分水岭算法分离相邻物体') parser.add_argument('--watershed-min-distance', type=int, default=10, help='分水岭局部最大值的最小距离 (默认: 10)') # 输出选项 parser.add_argument('--output-dir', default='output', help='输出目录 (默认: output)') parser.add_argument('--save-labeled-image', action='store_true', help='保存标记后的图像文件') args = parser.parse_args() try: print("=" * 70) print("形状特征分析工具") print("=" * 70) print(f"输入文件: {args.input_path}") print(f"输入类型: {args.input_type}") print(f"输出路径: {args.output}") print(f"输出目录: {args.output_dir}") if args.input_type == 'dat': print(f"波段索引: {args.band_index}") print(f"HDR文件: {args.hdr_path or '自动查找'}") print(f"连通性: {args.connectivity}") print(f"最小面积: {args.min_area}") print(f"使用分水岭: {args.use_watershed}") if args.use_watershed: print(f"分水岭最小距离: {args.watershed_min_distance}") print() # 创建配置 config = ShapeFeatureConfig() config.input_type = args.input_type config.output_dir = args.output_dir config.connectivity = args.connectivity config.min_area = args.min_area config.use_watershed = args.use_watershed config.watershed_min_distance = args.watershed_min_distance config.save_labeled_image = args.save_labeled_image # 设置输入文件路径 if args.input_type == 'dat': config.dat_file_path = args.input_path config.hdr_file_path = args.hdr_path config.band_index = args.band_index elif args.input_type == 'shp': config.shp_file_path = args.input_path # 验证配置 try: config.validate() print("配置验证通过") except Exception as e: print(f"✗ 配置验证失败: {e}") return 1 # 执行形状特征分析 print("\n开始形状特征分析...") features_df, labeled_image = analyze_shape_features(config) # 显示结果统计 if len(features_df) > 0: print("\n✓ 分析完成!") print(f"检测到的物体数量: {len(features_df)}") print(f"提取的特征数量: {len(features_df.columns)}") # 显示特征名称 feature_names = [col for col in features_df.columns if col not in ['label', 'centroid_row', 'centroid_col']] print(f"形状特征: {', '.join(feature_names[:10])}{'...' if len(feature_names) > 10 else ''}") # 显示统计信息 print("\n特征统计:") numeric_cols = features_df.select_dtypes(include=[np.number]).columns if len(numeric_cols) > 0: stats = features_df[numeric_cols].describe() print(f" 面积范围: {features_df['area'].min():.0f} - {features_df['area'].max():.0f} 像素") print(f" 平均面积: {features_df['area'].mean():.1f} 像素") print(f" 周长范围: {features_df['perimeter'].min():.1f} - {features_df['perimeter'].max():.1f}") print(f" 圆度范围: {features_df['circularity'].min():.3f} - {features_df['circularity'].max():.3f}") else: print("\n⚠ 未检测到任何物体，请检查输入数据或调整参数") print("\n" + "=" * 70) print("形状特征分析完成!") print("=" * 70) except Exception as e: print(f"✗ 处理失败: {e}") import traceback traceback.print_exc() return 1 return 0 if __name__ == '__main__': main()