import numpy as np import cv2 from osgeo import gdal from collections import Counter from typing import Optional, Union, Tuple gdal.UseExceptions() def laplacian_filter(image): """ 拉普拉斯算子纹理提取 使用二阶偏微分对图像进行卷积，获取纹理图像 公式: η²g(m,n) = g(m,n+1) + g(m,n-1) + g(m+1,n) + g(m-1,n) - 4g(m,n) 参数: image: 输入图像（全波段叠加后的灰度图像F） 返回: L: 拉普拉斯滤波后的纹理图像 """ # 拉普拉斯算子核（4邻域，中心5个像元） kernel = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]], dtype=np.float32) filtered_image = cv2.filter2D(image, -1, kernel) return filtered_image def apply_threshold(image, threshold): """ 二值化处理（分割耀光区域） 公式: - 对于L: W = {1, L>0; 0, L≤0} - 对于F: S = {1, F>N; 0, F≤N} 参数: image: 输入图像 threshold: 阈值 返回: binary_image: 二值化图像（0或1） """ _, binary_image = cv2.threshold(image, threshold, 1, cv2.THRESH_BINARY) return binary_image def morphological_dilation(image, iterations=1): """ 形态学处理（膨胀操作） 对耀光纹理区域进行扩展，使得检测到的耀光区域更加连贯 参数: image: 输入二值图像 iterations: 膨胀迭代次数 返回: dilated_image: 膨胀后的图像 """ kernel = np.ones((3, 3), np.uint8) dilated_image = cv2.dilate(image, kernel, iterations=iterations) return dilated_image def calculate_area_difference(W, S): """ 计算面积差值 公式: rq = area(W) - area(S) 其中 area(W) 表示拉普拉斯纹理提取区域的面积（用于定位耀光位置） area(S) 表示通过阈值分割全波段叠加图像得到的耀光区域面积 参数: W: 纹理提取区域（二值图像） S: 阈值分割区域（二值图像） 返回: rq: 面积差值 area_W: 纹理提取区域面积 area_S: 阈值分割区域面积 """ area_W = np.sum(W) area_S = np.sum(S) rq = area_W - area_S return rq, area_W, area_S def multi_band_weighted_sum(image_bands, water_mask): """ 全波段叠加：将所有波段的遥感反射率加权叠加形成一个灰度图像 公式: F = Σ(i=1 to n) R(λi) × G 其中 G 为水体二值化图像，F 为灰度图像 参数: image_bands: 多波段图像数组，形状为 (rows, cols, bands) water_mask: 水体二值化图像 G，形状为 (rows, cols)，值为0或1 返回: F: 全波段加权叠加后的灰度图像 """ # 确保water_mask是二值化的（0或1） if water_mask.dtype != np.float32 and water_mask.dtype != np.float64: water_mask = water_mask.astype(np.float32) # 对每个波段进行加权叠加：F = Σ R(λi) × G F = np.zeros((image_bands.shape[0], image_bands.shape[1]), dtype=np.float32) for band_idx in range(image_bands.shape[2]): F += image_bands[:, :, band_idx] * water_mask return F def find_optimal_threshold(F, u=0.1, q=50, r=20, max_iterations=None): """ 通过迭代找到最佳阈值 算法原理： 1. 使用拉普拉斯算子提取耀光纹理信息（用于定位耀光的位置） 2. 通过全波段叠加图像F进行阈值分割提取耀光面积 3. 当纹理提取区域的面积与通过阈值分割得到的耀光区域面积差最小时，确定最佳的阈值 步骤： 1. 对全波段叠加图像F进行拉普拉斯纹理提取，得到L 2. 对L进行二值化：W = {1, L>0; 0, L≤0}（用于定位耀光位置） 3. 对W和S进行形态学膨胀r次 4. 设定阈值N的初始值为F的最小值，表示为Nf 5. 对初始值叠加数值u更新阈值，迭代q次 6. 每次迭代记录W与S的面积差值到r行、q列的数组R 7. 寻找R中每列数组中最小值所在行数，形成新的数组aind 8. 统计aind的众数得到M，即为最佳的叠加次数 9. 最佳阈值：Nf = min(F) + u·M 参数: F: 全波段加权叠加后的灰度图像 u: 阈值更新步长，默认0.1（论文参数） q: 迭代次数，默认50（论文参数） r: 形态学膨胀次数，默认20（论文参数） max_iterations: 最大迭代次数（如果指定，则使用此值替代q） 返回: optimal_threshold: 最佳阈值Nf optimal_S: 最佳阈值对应的耀光区域S（通过阈值分割F得到） optimal_W: 最佳阈值对应的纹理区域W（通过拉普拉斯提取得到） M: 最佳迭代次数索引 area_differences: 面积差值矩阵R (r行, q列) thresholds: 每次迭代的阈值列表 """ # 如果指定了max_iterations，使用它替代q if max_iterations is not None: q = max_iterations # 初始化阈值为影像非0最小值 F_nonzero = F[F > 0] if len(F_nonzero) > 0: min_value = np.min(F_nonzero) else: # 如果所有值都为0，使用一个很小的正数作为最小值 min_value = np.finfo(np.float32).eps print("警告: F中所有值都为0，使用极小值作为最小值") # 步骤1: 对全波段叠加图像F进行拉普拉斯纹理提取（用于定位耀光位置） print("进行拉普拉斯纹理提取...") L = laplacian_filter(F) # 步骤2: 对L进行二值化，得到纹理区域W # W = {1, L>0; 0, L≤0} W = apply_threshold(L, 0.0) # 阈值为0，即L>0为1，L≤0为0 # 步骤3: 对W进行形态学膨胀r次 print(f"对纹理区域W进行形态学膨胀{r}次...") W_dilated = morphological_dilation(W, iterations=r) # 存储每次迭代的面积差值（r行，q列） # 注意：论文中r是膨胀次数，但这里R矩阵的r行应该对应不同的膨胀次数 # 根据论文描述，应该是迭代q次，每次记录面积差值 # 但论文提到"叠加q次"和"一共迭代r次"，这里理解为迭代q次，每次对W和S都膨胀r次 area_differences = [] # 存储每次迭代的面积差值 thresholds = [] W_masks = [] S_masks = [] print(f"开始迭代计算最佳阈值（迭代{q}次，步长u={u}）...") # 迭代更新阈值（迭代q次） for i in range(q): # 当前阈值：Nf = min(F) + u·(i+1) current_threshold = min_value + u * (i + 1) thresholds.append(current_threshold) # 步骤4: 对全波段叠加图像F进行阈值分割，得到耀光区域S（用于提取耀光面积） # S = {1, F>N; 0, F≤N} S = apply_threshold(F, current_threshold) # 步骤5: 对S进行形态学膨胀r次 S_dilated = morphological_dilation(S, iterations=r) # 步骤6: 计算面积差值：纹理提取区域面积 vs 阈值分割区域面积 rq, area_W, area_S = calculate_area_difference(W_dilated, S_dilated) area_differences.append(rq) W_masks.append(W_dilated.copy()) S_masks.append(S_dilated.copy()) if (i + 1) % 10 == 0: print(f" 迭代 {i+1}/{q}: 阈值={current_threshold:.4f}, 面积差值={rq:.2f}") # 步骤7: 寻找R中每列数组中最小值所在行数 # 注意：论文中提到"r行、q列的数组R"，但根据算法描述，应该是q次迭代 # 这里理解为：将area_differences重新组织成矩阵形式（如果需要） # 但根据论文描述，应该是直接找到最小面积差值对应的迭代次数 area_differences_array = np.array(area_differences) # 步骤8: 找到最小面积差值所在的行数（索引） # argmin_q(R(r,q))：找到最小面积差值所在的行数 min_indices = np.where(area_differences_array == np.min(area_differences_array))[0] # 步骤9: 通过众数统计找到最频繁出现的行数，确定最佳的叠加次数M # Mode(aind)：通过众数统计找到最频繁出现的行数 if len(min_indices) > 0: # 如果最小值出现多次，使用众数统计 counter = Counter(min_indices) most_common = counter.most_common(1)[0] M = most_common[0] # 最佳迭代次数索引（从0开始） else: M = 0 # 步骤10: 计算最终阈值 # Nf = min(F) + u·M # 注意：M是索引（从0开始），所以实际迭代次数是M+1 optimal_threshold = min_value + u * (M + 1) # 获取最佳阈值对应的掩膜 optimal_W = W_masks[M] optimal_S = S_masks[M] print(f"最佳迭代次数索引: M={M} (第{M+1}次迭代)") print(f"最佳阈值: Nf={optimal_threshold:.4f}") print(f"最小面积差值: {np.min(area_differences_array):.2f}") return optimal_threshold, optimal_S, optimal_W, M, area_differences, thresholds def generate_glint_mask(bsq_file, water_mask=None, u=0.1, q=50, r=20, max_iterations=None, output_file=None): """ 生成耀光掩膜 算法流程： 1. 全波段叠加：F = Σ(i=1 to n) R(λi) × G（G为水体二值化图像） 2. 拉普拉斯算子提取纹理信息（用于定位耀光位置） 3. 通过阈值分割全波段叠加图像提取耀光面积 4. 当纹理提取区域面积与阈值分割区域面积差最小时，确定最佳阈值 参数: bsq_file: 输入的BSQ文件路径 water_mask: 水体二值化图像G，可以是： - None: 自动生成（基于所有像素，即全为1的掩膜） - numpy数组: 直接使用数组作为掩膜，形状为 (rows, cols)，值为0或1 - 文件路径: 栅格文件路径(.tif/.dat)，将自动读取 u: 阈值更新步长，默认0.1（论文参数） q: 迭代次数，默认50（论文参数） r: 形态学膨胀次数，默认20（论文参数） max_iterations: 最大迭代次数（如果指定，则使用此值替代q） output_file: 输出文件路径，如果为None则自动生成 返回: tuple: (耀光掩膜文件路径, 纹理提取图像文件路径) - 耀光掩膜文件路径: 通过阈值分割全波段叠加图像得到的最终掩膜（S掩膜） - 纹理提取图像文件路径: 拉普拉斯纹理提取后的二值化图像（W掩膜） """ # 读取BSQ文件 bsq_dataset = gdal.Open(bsq_file) if bsq_dataset is None: raise ValueError(f"无法打开文件: {bsq_file}") # 获取影像数据 bands = bsq_dataset.RasterCount rows = bsq_dataset.RasterYSize cols = bsq_dataset.RasterXSize print(f"影像尺寸: {rows} x {cols}, 波段数: {bands}") # 读取所有波段 print("正在读取所有波段数据...") image_bands = np.zeros((rows, cols, bands), dtype=np.float32) for band in range(bands): image_bands[:, :, band] = bsq_dataset.GetRasterBand(band + 1).ReadAsArray().astype(np.float32) if (band + 1) % 20 == 0: print(f" 已读取 {band+1}/{bands} 个波段") # 处理水体掩膜G if water_mask is None: # 如果没有提供水体掩膜，使用全图（所有像素为1） print("未提供水体掩膜，使用全图进行处理") G = np.ones((rows, cols), dtype=np.float32) elif isinstance(water_mask, np.ndarray): # 如果直接提供了numpy数组 if water_mask.shape != (rows, cols): raise ValueError(f"水体掩膜尺寸 {water_mask.shape} 与影像尺寸 {(rows, cols)} 不匹配") G = water_mask.astype(np.float32) # 确保是二值化的（0或1） G = np.where(G > 0, 1.0, 0.0) elif isinstance(water_mask, str): # 如果是文件路径，读取文件 print(f"从文件读取水体掩膜: {water_mask}") water_dataset = gdal.Open(water_mask) if water_dataset is None: raise ValueError(f"无法打开水体掩膜文件: {water_mask}") if water_dataset.RasterXSize != cols or water_dataset.RasterYSize != rows: raise ValueError(f"水体掩膜尺寸与影像尺寸不匹配") G = water_dataset.GetRasterBand(1).ReadAsArray().astype(np.float32) water_dataset = None # 确保是二值化的（0或1） G = np.where(G > 0, 1.0, 0.0) else: raise ValueError(f"不支持的水体掩膜类型: {type(water_mask)}") print(f"水体掩膜统计: 水体像素数={np.sum(G)}, 总像素数={rows*cols}, 水体比例={np.sum(G)/(rows*cols)*100:.2f}%") # 步骤1: 全波段叠加 F = Σ(i=1 to n) R(λi) × G print("开始全波段叠加...") F = multi_band_weighted_sum(image_bands, G) print(f"全波段叠加完成，F值范围: [{np.min(F):.4f}, {np.max(F):.4f}]") print("开始计算最佳阈值...") # 找到最佳阈值 optimal_threshold, glint_mask, texture_mask, optimal_iteration, area_diffs, thresholds = find_optimal_threshold( F, u=u, q=q, r=r, max_iterations=max_iterations ) print(f"\n=== 最佳阈值计算结果 ===") print(f"最佳阈值: {optimal_threshold:.4f}") print(f"最佳迭代次数: {optimal_iteration + 1}") print(f"最小面积差值: {np.min(area_diffs):.4f}") print(f"纹理提取区域面积: {np.sum(texture_mask)}") print(f"阈值分割区域面积: {np.sum(glint_mask)}") # 设置输出文件路径 if output_file is None: output_file = 'glint_mask.tif' # 生成纹理提取图像输出路径 texture_output_file = output_file.replace('.tif', '_texture.tif') if texture_output_file == output_file: # 如果没有.tif扩展名 texture_output_file = output_file + '_texture.tif' # 保存纹理提取图像（W掩膜，拉普拉斯提取的纹理区域） print(f"\n保存输出文件...") driver = gdal.GetDriverByName('GTiff') texture_dataset = driver.Create(texture_output_file, cols, rows, 1, gdal.GDT_Byte) texture_dataset.SetGeoTransform(bsq_dataset.GetGeoTransform()) texture_dataset.SetProjection(bsq_dataset.GetProjection()) texture_mask_uint8 = (texture_mask * 255).astype(np.uint8) texture_dataset.GetRasterBand(1).WriteArray(texture_mask_uint8) texture_dataset = None print(f"纹理提取图像已保存至: {texture_output_file}") # 保存耀光掩膜（S掩膜，通过阈值分割全波段叠加图像得到） out_dataset = driver.Create(output_file, cols, rows, 1, gdal.GDT_Byte) # 设置地理变换和投影信息 out_dataset.SetGeoTransform(bsq_dataset.GetGeoTransform()) out_dataset.SetProjection(bsq_dataset.GetProjection()) # 写入掩膜数据（转换为0-255范围） glint_mask_uint8 = (glint_mask * 255).astype(np.uint8) out_dataset.GetRasterBand(1).WriteArray(glint_mask_uint8) # 关闭数据集 out_dataset = None bsq_dataset = None print(f"耀光掩膜已保存至: {output_file}") return output_file, texture_output_file # 使用示例 if __name__ == "__main__": bsq_file = r"D:\BaiduNetdiskDownload\yaobao\test_glint.bsq" # 输入的BSQ文件 output_file = r'D:\BaiduNetdiskDownload\yaobao\glint\lapulas_otsu_glint_mask.tif' # water_mask_file = r'path/to/water_mask.tif' # 可选：水体掩膜文件路径 # 示例1: 使用论文默认参数（q=50, r=20, u=0.1） mask_file, texture_file = generate_glint_mask( bsq_file, water_mask=None, u=0.1, q=50, r=20, output_file=output_file ) print(f"\n处理完成，耀光掩膜保存在: {mask_file}") print(f"纹理提取图像保存在: {texture_file}") # 示例2: 使用水体掩膜文件 # mask_file, texture_file = generate_glint_mask( # bsq_file, # water_mask=water_mask_file, # u=0.1, # q=50, # r=20, # output_file=output_file # ) # 示例3: 使用numpy数组作为水体掩膜 # import numpy as np # water_mask_array = np.ones((rows, cols), dtype=np.float32) # 示例：全为1的掩膜 # mask_file, texture_file = generate_glint_mask( # bsq_file, # water_mask=water_mask_array, # u=0.1, # q=50, # r=20, # output_file=output_file # )