文章目录
- 应用背景
- 代码
- 结果
- 进阶版
- 获得断点后获取该断点所在线段的方向向量
- 获得路口点
- 稀疏化
应用背景
下面介绍一种方法，对骨架线的断点进行提取。
代码
import numpy as np
import cv2 as cv
import matplotlib.pyplot as plt
from skimage import morphology
import os
def detect_breakpoints(img):
r, c = img.shape
new_image = np.zeros((r, c))
detection_operator = np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]])
break_points = []
for i in range(r-2):
for j in range(c-2):
if img[i+1, j+1] == 0:
new_image[i + 1, j + 1] = 0
else:
if np.sum(img[i:i + 3, j:j + 3] * detection_operator) == 2:
new_image[i+1, j+1] = 1
break_points.append([i+1, j+1])
else:
new_image[i+1, j+1] = 0
return np.uint8(new_image), break_points结果
进阶版
获得断点后获取该断点所在线段的方向向量
def detectCenterDirection(img, centerPoint):
detection_operator = np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]])
x, y = centerPoint
dx, dy = x, y
whilestopSignal = False
count = 0
linePoints = []
linePoints.append((x, y))
while True:
forStopSiganl = False
for i in [-1, 0, 1]:
for j in [-1, 0, 1]:
start_center_x, start_center_y = x + i, y + j
# print(f'***********{linePoints}')
if (start_center_x, start_center_y) in linePoints:
# print(f'***********{linePoints}')
continue
subImg = img[start_center_x - 1: start_center_x + 2, start_center_y - 1: start_center_y + 2]
if np.sum(subImg * detection_operator) > 3:
dx, dy = start_center_x, start_center_y
linePoints.append((start_center_x, start_center_y))
whilestopSignal = True
forStopSiganl = True
break
if np.sum(subImg * detection_operator) == 3 and subImg[1, 1] == 1:
x, y = start_center_x, start_center_y
dx, dy = x, y
linePoints.append((start_center_x, start_center_y))
forStopSiganl = True
break
if forStopSiganl == True:
break
if whilestopSignal == True:
break
count += 1
if count > 8:
whilestopSignal = True
print(f'count:{count}, linePoints:{linePoints}')
return [centerPoint[0], centerPoint[1], dx, dy]获得路口点
def detect_intersection(img):
r, c = img.shape
new_image = np.zeros((r, c))
detection_operator = np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]])
break_points = []
for i in range(r-2):
for j in range(c-2):
if img[i+1, j+1] == 0:
new_image[i + 1, j + 1] = 0
else:
if np.sum(img[i:i + 3, j:j + 3] * detection_operator) > 3 and img[i+1, j+1] == 1:
new_image[i+1, j+1] = 1
break_points.append([i+1, j+1])
else:
new_image[i+1, j+1] = 0
return np.uint8(new_image), break_points稀疏化
利用欧式距离进行点的稀疏化操作。
v1
def 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 keep
else:
return 0
def 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 pointKeepv2
网格化后再稀疏
import numpy as np
import os
import gdalTools
import time
import sys
def detect_breakpoints2(img):
ignore_border = 4000
r, c = img.shape
new_image = np.zeros((r, c))
detection_operator = np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]])
break_points = []
for i in range(ignore_border, r-ignore_border):
for j in range(ignore_border, c-ignore_border):
if img[i+1, j+1] == 0:
new_image[i + 1, j + 1] = 0
else:
if np.sum(img[i:i + 3, j:j + 3] * detection_operator) == 2:
new_image[i+1, j+1] = 1
break_points.append([i+1, j+1])
else:
new_image[i+1, j+1] = 0
return break_points
def detect_breakpoints(img, interval=400):
ignore_border = 4000
width, height = img.shape
col = np.ceil(width / interval).astype(np.uint16)
row = np.ceil(height / interval).astype(np.uint16)
print(f"col :{col}, row: {row}")
length = col * row
new_image = np.zeros((width, height))
detection_operator = np.array([[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):
if img[i+1, j+1] == 0:
new_image[i + 1, j + 1] = 0
else:
val = np.sum(img[i:i + 3, j:j + 3] * detection_operator)
if val == 2 or val == 5:
sub_col = np.floor(i / interval).astype(np.uint16)
sub_row = np.floor(j / interval).astype(np.uint16)
new_image[i+1, j+1] = 1
index = sub_col * row + sub_row
break_points[index].append([i + 1, j + 1])
else:
new_image[i+1, j+1] = 0
return break_points, new_image
def 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 keep
else:
return 0
def 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 pointKeep
def point_sparse2(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 pointKeep
if __name__ == '__main__':
CoarseSegPath = r'D:\MyWorkSpace\MyUtilsCode\segmentation\CroplandPredict\eval8\skeleton.tif'
outPath = CoarseSegPath.replace("skeleton", "points")
im_proj, im_geotrans, im_width, im_height, coarseImg = gdalTools.read_img(CoarseSegPath)
coarseImg = np.where(coarseImg > 0, 1, 0)
break_points, mask = detect_breakpoints(coarseImg, 200)
gdalTools.write_img(outPath, im_proj, im_geotrans, mask)
# print(break_points)
break_points2 = []
for k, v in break_points.items():
# print(f"key:{k}, val:{v}")
break_points2 += v
print(f"detect {len(break_points2)}")
s1 = time.time()
sparsePoints = []
for k, v in break_points.items():
sub_break_points = point_sparse(v, 30)
sparsePoints += sub_break_points
print(f"process {time.time() - s1}, length:{len(sparsePoints)}")
s2 = time.time()
break_points2 = point_sparse(break_points2, 30)
print(f"process {time.time() - s2}, length:{len(break_points2)}")
速度要快不少
