1 方法概述

首先基于语义分割模型 UNet 对目标区域进行粗提取，两分类（背景和鱼塘）模型获取前景区域，三分类（背景、鱼塘和边界）模型获取标记区域，针对困难边界，结合局部分割方法对边缘进行修复，最后利用分水岭算法对两个类别（养殖塘类和边界类）进行融合。养殖塘地块提取流程如下图所示。

2 标记区域提取

我们选用语义分割UNet作为本次试验的标记区域提取模型。为了获取地块级别的养殖塘结果，受文献启发，我们对二值图g进行标签增强操作来创造与边界相关的第三个类，这个操作通过形态学操作完成。在训练过程中，使用边界增强后的样本进行粗提取模型训练，得到的概率分布图将会有额外的类别来表示养殖塘边界的概率分布，使得粗提取模型能够更好区分养殖塘地块边缘和内部区域。

2 前景区域提取

我们选用语义分割UNet作为本次试验的前景区域提取模型。目标类别：背景和鱼塘。

4 困难边界修复

困难边界区域推荐。面对粗分割结果时，首先需要确定应该细化结果的哪一部分。本文设计了一种有效的滑窗算法，沿着养殖塘边界提取一系列困难边界。具体来说，我们设计了一个 7×7 大小的算子，其值全为 1，提取粗分割结果的养殖塘边界二值图，对边界二值图以步长为 2 进行滑窗，并计算对应位置的大小 v，当 v 大于 10 时，该位置视为困难边界。如图（b）所示，获取的推荐区域仍然存在很大的重叠和冗余，因此我们进一步应用区域中心点之间的欧氏距离来稀疏化推荐区域，其结果如图（c）所示。根据经验，重叠越大，分割性能可以提高，同时需要承受较大的计算成本。我们可以调整欧氏距离阈值来控制重叠的程度，以实现更好的速度和精度的平衡。根据稀疏化的推荐区域，在原始影像中生成对应位置的影像切片，并输入到下面的局部边缘修复网络中。

(a) 粗分割结果；（b）困难边界区域推荐；（c）稀疏化的推荐区域；（d）原始影像；（c）根据推荐区域生成的影像切片；（d）局部边缘修复网络的预测结果；（e）融合结果。

5 源码

https://github.com/SonwYang/DifficultBoundaryRepair

patch_predict.py


import os
os.environ["KMP_DUPLICATE_LIB_OK"]="TRUE"import torch
import gdalTools
import cv2
import numpy as np
from geoClip import *
import tqdmimport matplotlib.pyplot as pltdef point_sparse(pointlist, thresh=12):order = pointlist.copy()pointKeep = []while len(order) > 0:point = order[0]otherPoints = order[1:]keep = compute_euc(otherPoints, point, thresh)if keep == 0:pointKeep.append(point)order.pop(0)return pointKeepdef detect_hard_edges(img, step=1, ingore_size=32):r, c = img.shapenew_image = np.zeros((r, c))detection_operator = np.array([[1, 1, 1, 1, 1, 1, 1],[1, 1, 1, 1, 1, 1, 1],[1, 1, 1, 1, 1, 1, 1],[1, 1, 1, 1, 1, 1, 1],[1, 1, 1, 1, 1, 1, 1],[1, 1, 1, 1, 1, 1, 1],[1, 1, 1, 1, 1, 1, 1]])hard_edges_points = []for i in range(ingore_size, r-ingore_size, step):for j in range(ingore_size, c-ingore_size, step):val = np.sum(img[i:i + 7, j:j + 7] * detection_operator)pointVal = img[i+4, j+4]if val > 10:new_image[i+4, j+4] = 1hard_edges_points.append([i+4, j+4])else:new_image[i+1, j+1] = 0return np.uint8(new_image), hard_edges_pointsdef detect_hard_edges2(img, interval=512, ignore_border=32, value=10):width, height = img.shapecol = np.ceil(width / interval).astype(np.uint16)row = np.ceil(height / interval).astype(np.uint16)print(f"col :{col}, row: {row}")length = col * rowdetection_operator = np.array([[1, 1, 1, 1, 1, 1, 1],[1, 1, 1, 1, 1, 1, 1],[1, 1, 1, 1, 1, 1, 1],[1, 1, 1, 1, 1, 1, 1],[1, 1, 1, 1, 1, 1, 1],[1, 1, 1, 1, 1, 1, 1],[1, 1, 1, 1, 1, 1, 1]])break_points = {}for i in range(length):break_points[i] = []for i in range(ignore_border, width-ignore_border):for j in range(ignore_border, height-ignore_border):val = np.sum(img[i:i + 7, j:j + 7] * detection_operator)if val > value:sub_col = np.floor(i / interval).astype(np.uint16)sub_row = np.floor(j / interval).astype(np.uint16)index = sub_col * row + sub_rowbreak_points[index].append([i + 4, j + 4])return break_pointsdef compute_euc(pointlist, point, thresh=12):if len(pointlist) > 0:distances = np.sqrt(np.sum(np.asarray(np.array(point) - np.array(pointlist)) ** 2, axis=1))keep = np.where(distances < thresh, 1, 0).sum()return keepelse:return 0def normal_img(img):min = np.min(img)max = np.max(img)img = (img - min) / (max - min) * 255return img.astype(np.uint8)def inference(point, im_data, model):x, y = pointsubImg = im_data[:, x - half_width:x + half_width, y - half_width:y + half_width]vis = subImg.transpose((1, 2, 0))[..., :3]img = np.array(subImg, dtype=np.float32)img /= 255.tensor = torch.from_numpy(np.array([img])).cuda()outs = model(tensor).cpu().detach()prob = outs[0]_, maxprob = torch.max(prob, 0)maxprob = maxprob.numpy()[12:half_width * 2 - 12, 12:half_width * 2 - 12]return maxprobdef get_labels():"""Load the mapping that associates classes with label colorsReturns:np.ndarray with dimensions (13, 3)"""return np.asarray([[0, 0, 0],[255, 85, 0],[150, 250, 0],[150, 200, 150],[200, 0, 200],[150, 0, 250],[150, 150, 250],[250, 200, 0],[200, 200, 0],[100, 56, 250],[250, 150, 56],[200, 0, 0],[250, 0, 150],[200, 150, 150],[250, 150, 125],[0, 0, 200],[0, 100, 50],[0, 150, 200],[0, 200, 250],[255, 255, 255],# [255, 255, 0],[250, 128, 114],[255, 20, 147],[105, 139, 34],[139, 58, 58],[255,69,0],[255,165,0],[128,128,0],[64,224,208]])def decode_segmap(label_mask, n_classes):"""Decode segmentation class labels into a color imageArgs:label_mask (np.ndarray): an (M,N) array of integer values denotingthe class label at each spatial location.plot (bool, optional): whether to show the resulting color imagein a figure.Returns:(np.ndarray, optional): the resulting decoded color image."""label_colours = get_labels()color_length = len(label_colours)r = label_mask.copy()g = label_mask.copy()b = label_mask.copy()for ll in range(0, n_classes):if ll < color_length:r[label_mask == ll] = label_colours[ll, 0]g[label_mask == ll] = label_colours[ll, 1]b[label_mask == ll] = label_colours[ll, 2]else:lll = ll%(color_length-1)if lll == 0:lll = 13r[label_mask == ll] = label_colours[lll, 0]g[label_mask == ll] = label_colours[lll, 1]b[label_mask == ll] = label_colours[lll, 2]rgb = np.zeros((label_mask.shape[0], label_mask.shape[1], 3))rgb[:, :, 0] = rrgb[:, :, 1] = grgb[:, :, 2] = breturn rgb.astype(np.uint8)if __name__ == '__main__':import timefrom skimage.morphology import remove_small_objects, opening, squareimgPath = r'D:\BaiduNetdiskDownload\data_3groups\GF3_Validation_sub2\entropy_shannon_subset_Feature_sub2.tif'CoarseSegPath = r'D:\MyWorkSpace\paper\fishpond\data_evaluation\g2_otherNet\unet3\poly.tif'modelPath = r'D:\MyWorkSpace\paper\fishpond\PatchSeg\ckpts_cls3_res34_Unet\epoch0196_model.pth'outPath = 'poly.tif'euc_threshold = 8val_threshold = 10half_width = 32half_width2 = half_width - 12clip_size = 0model = torch.load(modelPath, map_location="cpu").cuda()im_proj, im_geotrans, _, _, im_data = gdalTools.read_img(imgPath)arr = [normal_img(im_data[idx]) for idx in range(len(im_data))]im_data = np.dstack(arr).transpose(2, 0, 1)#io.imsave("test.png", im_data.transpose(1, 2, 0))_, _, _, _, coarseImg = gdalTools.read_img(CoarseSegPath)coarseImg = np.where(coarseImg == 1, 1, 0).astype(np.uint8)## remove small areas# coarseImg = coarseImg.astype(bool)# coarseImg = remove_small_objects(coarseImg, 200)# coarseImg = coarseImg.astype(np.uint8)# gdalTools.write_img("remove_small_areas.tif", im_proj, im_geotrans, coarseImg)result = coarseImg.copy()visImg = np.array(decode_segmap(coarseImg, 2))visImg2 = visImg.copy()im_width, im_height = coarseImg.shapeimg_vis = np.zeros_like(coarseImg)PFMask = np.zeros_like(coarseImg)img_vis = [img_vis for _ in range(3)]img_vis = np.dstack(img_vis)contours, hierarchy = cv2.findContours(coarseImg, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)img_vis2 = cv2.drawContours(img_vis.copy(), contours, -1, (255, 255, 255), 1)mask = np.where(img_vis2[:, :, 0] > 0, 1, 0)hard_points = detect_hard_edges2(mask, ignore_border=half_width, value=val_threshold)break_points2 = []for k, v in hard_points.items():# print(f"key:{k}, val:{v}")break_points2 += vprint(f"detect {len(break_points2)}")sparsePoints = []for k, v in hard_points.items():sub_break_points = point_sparse(v, euc_threshold)sparsePoints += sub_break_pointsprint(f"The sparse op is finished, length:{len(sparsePoints)}")count = 0for i, point in tqdm.tqdm(enumerate(sparsePoints)):x, y = pointsubImg = im_data[:, x-half_width:x+half_width, y-half_width:y+half_width]#PFMask[x - half_width2: x + half_width2, y-half_width2: y + half_width2] = 1vis = subImg.transpose((1, 2, 0))[..., :3]img = np.array(subImg, dtype=np.float32)img /= 255.try:tensor = torch.from_numpy(np.array([img])).cuda()outs = model(tensor).cpu().detach()prob = outs[0]_, maxprob = torch.max(prob, 0)maxprob = maxprob.numpy()[clip_size:half_width * 2 - clip_size, clip_size:half_width * 2 -clip_size]result[x-half_width+clip_size:x+half_width-clip_size, y-half_width+clip_size:y+half_width-clip_size] = maxprobexcept:count += 1print(f"failure : {count}")continuegdalTools.write_img(outPath, im_proj, im_geotrans, result)

论文

https://www.mdpi.com/2072-4292/15/9/2246

基于SAR影像的鱼塘提取相关推荐

基于遥感影像的道路提取论文、开源代码和数据集汇总
文章目录前言 2017 DeepRoadMapper Topology Loss 2018 RoadTracer iterative-deep-learning 2019 Leveraging Cr ...
论文笔记（五）FWENet：基于SAR图像的洪水水体提取深度卷积神经网络（CVPR）
FWENet: a deep convolutional neural network for flood water body extraction based on SAR images 作者:J ...
[SARscape] 多时相SAR影像的应用 - 监督分类、提取水稻种植区 - 以Sentinel-1A数据为例
SARscape工具操作可查看:SARscape工具集 SAR影像分类 [基于SAR影像进行图像分类]根据特定物体后向散射系数的特征进行影像分类举例:水稻种植区提取技术(Sentinel-1A数据为 ...
sar偏移量追踪技术_论文推荐 | 吴文豪：基于几何配准的多模式SAR影像配准及其误差分析...
<测绘学报> 构建与学术的桥梁拉近与权威的距离基于几何配准的多模式SAR影像配准及其误差分析吴文豪1, 张磊2, 李陶3, 龙四春1, 段梦4, 周志伟5, 祝传广1, 蒋廷臣61. ...
基于高分辨率影像城市绿地信息提取_武汉大学眭海刚教授等：多时相遥感影像变化检测方法综述...
导语随着遥感对地观测技术的快速发展和应用,变化检测技术体系也在不断地发展和演化. 遥感影像变化检测到底是什么?为什么如此热门? 近年兴起的基于深度学习的变化检测新方法.高分辨率遥感影像场景变化分析. ...
基于高光谱影像的农作物检测应用简介
声明 1)本文仅供学术交流,非商用.具体引用的资料请看参考文献.如果某部分不小心侵犯了大家的利益,请联系博主删除. 2)本人才疏学浅,整理总结的时候难免出错,还望各位前辈不吝指正,谢谢. 联系方式:k ...
基于多阈值的形态提取遥感图像中的沿海线的特征方法（Qu Jishuang）
论文阅读笔记 Title:A multi-threshold based morphological approach for extracting coastal line feature in r ...
基于Sentinel-2的杞县大蒜提取试验
基于Sentinel-2的杞县大蒜提取试验一.前言杞县大蒜已有2000多年的栽培历史是全国农产品地理标志之一,2009年9月14日,中华人民共和国农业部正式批准对"杞县大蒜"实 ...
python实现dem输出三维模型_资源三号卫星影像立体像对如何提取DEM数据的方法
原标题:资源三号卫星影像立体像对如何提取DEM数据的方法 OrthoMapping是ArcGIS 10.5推出的基于无人机.大飞机.卫星拍摄的原始影像获取专业级别信息产品的生产能力.使用OrthoMa ...

基于SAR影像的鱼塘提取

文章目录

相关背景