图像矢量化

python将图像目标区域转换为LineString、Polygon数据格式。

import os
import numpy as np
import cv2
from PIL import Image
import matplotlib.pyplot as plt
from tqdm import tqdm
import sknw
from skimage.morphology import skeletonize, medial_axis
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point, LineString, Polygon
from scipy.signal import savgol_filter
from rdp import rdp
from shapely.wkt import loads
from skimage.measure import find_contours
import scipy.signal as signal
from multiprocessing.dummy import Pool as ThreadPool
import time
from datetime import datetime
import sys
import  warnings
from coord_convert.transform import wgs2gcj, wgs2bd, gcj2wgs, gcj2bd, bd2wgs, bd2gcj
from pyproj import Transformer
import math
from shapely.geometry import Polygon,LineString
import shapely

# warnings.filterwarnings("ignore")

# 图像缩放：https://www.cnblogs.com/lfri/p/10596530.html
def douglas_peucker(l, epsilon=0.00007):
    """ 道格拉斯-普克算法
    """
    try:
        if (type(l) == LineString):
            in_lng_list, in_lat_list = l.xy[0], l.xy[1]
            ins = [[i, j] for i, j in zip(in_lng_list, in_lat_list)]
            outs = rdp(ins, epsilon)

            out_lng_list = [i[0] for i in outs]
            out_lat_list = [i[1] for i in outs]
            new_l = LineString([Point(x, y) for x, y in zip(out_lng_list, out_lat_list)])
            return new_l
        if (type(l) == Polygon):
            in_lng_list, in_lat_list = l.exterior.xy[0], l.exterior.xy[1]
            ins = [[i, j] for i, j in zip(in_lng_list, in_lat_list)]
            outs = rdp(ins, epsilon)

            out_lng_list = [i[0] for i in outs]
            out_lat_list = [i[1] for i in outs]
            #             print(len(out_lng_list))
            if (out_lng_list[0] != out_lng_list[-1] and out_lat_list[0] != out_lat_list[-1]):
                print('处理LinearRing问题。')
                out_lng_list.append(out_lng_list[0])
                out_lat_list.append(out_lat_list[0])
            print('------', out_lng_list, out_lat_list)
            new_l = Polygon([Point(x, y) for x, y in zip(out_lng_list, out_lat_list)])
            return new_l
    except:
        print('douglas_peucker error--->', str(l))

    return None


def sg_filter(l, m=5, n=2):
    """ sg滤波
    """
    # try:
    if (type(l) == LineString):
        x_list, y_list = l.xy[0], l.xy[1]
        if (len(x_list) < m):
            return l
        new_x_list = list(savgol_filter(np.array(x_list), m, n))
        new_y_list = list(savgol_filter(np.array(y_list), m, n))
        new_l = LineString([Point(x, y) for x, y in zip(new_x_list, new_y_list)])
        return new_l
    if (type(l) == Polygon):
        x_list, y_list = l.exterior.xy[0], l.exterior.xy[1]
        if (len(x_list) < m):
            return l
        new_x_list = list(savgol_filter(np.array(x_list), m, n))
        new_y_list = list(savgol_filter(np.array(y_list), m, n))
        new_l = Polygon(list(zip(new_x_list, new_y_list)))
        return new_l
    # except:
    # print(type(l) == Polygon)
    # print('sg_filter error--->', str(l))

    # return None


# === 函数定义 ===

# 时间戳格式转换
def trans_time(timeStamp):
    if (timeStamp == timeStamp):
        timeStamp = float(timeStamp)
        return datetime.fromtimestamp(int(timeStamp))
    else:
        return '-1'


# 计算时间戳秒差
def diff_seconds(timeStamp1, timeStamp2):
    return abs((trans_time(timeStamp1) - trans_time(timeStamp2)).total_seconds())


def b_image(mask):
    """ 二值化彩色图像。
    """
    gray = cv2.cvtColor(mask, cv2.COLOR_RGB2GRAY)
    ret, binary = cv2.threshold(gray, 0, 255, cv2.THRESH_OTSU)

    return binary


def add_black_edg(arr):
    """
    find_counters 边缘提取时无法封闭图像边界，需要加一层黑边。
    :param arr:
    :return:
    """
    arr[:, 0] = 0
    arr[:, -1] = 0
    arr[0, :] = 0
    arr[-1, :] = 0
    return arr


def split_gray_image(gray):
    """
    gray_level = [["CheDaoXianBaiShiXian",1],
                ["CheDaoXianBaiXuXian",2],
                ["CheDaoXianHuangShiXian",3],
                ["CheDaoXianHuangXuXian",4],
                ["RenXingHengDao",5],
                ["DaoLiuDai",6],
                ["JianTou",7],
                ["DaoLuMian",8]]
    """
    gray = (255 - gray) // 20
    CheDaoXianBaiShiXian = np.where(gray == 1, 1, 0)
    CheDaoXianBaiXuXian = np.where(gray == 2, 1, 0)
    CheDaoXianHuangShiXian = np.where(gray == 3, 1, 0)
    CheDaoXianHuangXuXian = np.where(gray == 4, 1, 0)
    DaoLuMian = add_black_edg(np.where(gray == 8, 1, 0))
    JianTou = np.where(gray == 7, 1, 0)
    RenXing = np.where(gray == 5, 1, 0)
    DaoLiuDai = np.where(gray == 6, 1, 0)
    Back = add_black_edg(np.where(gray == 0, 1, 0))

    LvHuaDai = add_black_edg(np.where(gray == 9, 1, 0))
    ZhaLan = add_black_edg(np.where(gray == 10, 1, 0))

    return CheDaoXianBaiShiXian, CheDaoXianBaiXuXian, CheDaoXianHuangShiXian, CheDaoXianHuangXuXian, JianTou, DaoLuMian, Back, RenXing, DaoLiuDai, LvHuaDai, ZhaLan


def extract_lane_by_sheet(seg_image, filename):
    #     print(filename)
    """ 提取图幅的线拓扑结构。
        input
            image_name 108689_108710-53301_53322.png
        output
            gdf_edges, gdf_nodes 线GeoDataFrame, 点GeoDataFrame

    """

    # 形态学处理
    # seg_image = cv2.medianBlur(seg_image, 3) # 中值滤波
    kernel = np.ones((3, 3), np.uint8)
    seg_image = cv2.dilate(np.uint8(seg_image), kernel, iterations=1)

    # 形态学处理后提取拓扑
    #     skeleton = skeletonize(cv2.cvtColor(seg_image, cv2.COLOR_BGR2GRAY)//255)
    skeleton = skeletonize(seg_image)
    graph = sknw.build_sknw(skeleton)

    # =============================================================================
    # 线，draw edges by pts
    # 包括LineString、少量Point，少量MultiLineString
    # =============================================================================

    gdf_edges_List = []
    # gdf_edges = gpd.GeoDataFrame()
    for (s, e) in graph.edges():
        ps = graph[s][e]['pts']
        edge_x, edge_y = ps[:, 1], ps[:, 0]
        if len(ps) > 1:  # trans to line
            lng_list = [x for x in edge_x]
            lat_list = [y for y in edge_y]
            lst = LineString([Point(x, y) for x, y in zip(lng_list, lat_list)])
            lst = sg_filter(lst)
            lst = douglas_peucker(lst, 5)  # 0.000005)
            gdf_edges_List.append(lst)

    gdf_edges = gpd.GeoDataFrame([], geometry=gdf_edges_List)
    gdf_edges['_id'] = [filename for i in range(len(gdf_edges))]

    if (len(gdf_edges) > 0):
        gdf_edges['length'] = [gdf_edges.iloc[i]['geometry'].length for i in range(len(gdf_edges))]

    return gdf_edges


def extract_arrow_by_sheet(seg_image, filename):
    """ 提取图幅的面拓扑结构。
        input
            image_name 108689_108710-53301_53322.png
        output
            gdf_polygons 面GeoDataFrame

    """

    # 形态学处理
    #     seg_image = cv2.medianBlur(seg_image, 3) # 中值滤波
    # seg_image = cv2.cvtColor(seg_image, cv2.COLOR_BGR2GRAY)
    # 形态学处理后提取拓扑
    contours = find_contours(seg_image, 0.8)
    # # Display the image and plot all contours found
    # fig, ax = plt.subplots()
    # ax.imshow(seg_image, cmap=plt.cm.gray)
    #
    # for n, contour in enumerate(contours):
    #     ax.plot(contour[:, 1], contour[:, 0], linewidth=2, color='red')
    #
    # ax.axis('image')
    # ax.set_xticks([])
    # ax.set_yticks([])
    # plt.show()

    # =============================================================================
    # 面，draw polygon
    # =============================================================================

    gdf_polygons_List = []
    # gdf_polygons = gpd.GeoDataFrame()
    for i in range(len(contours)):
        contour = contours[i]
        polygon_x, polygon_y = contour[:, 1], contour[:, 0]
        lng_list = [x for x in polygon_x]
        lat_list = [y for y in polygon_y]

        pol = Polygon(list(zip(lng_list, lat_list)))
        pol = sg_filter(pol)  # , 21, 3)
        gdf_polygons_List.append(pol)
    gdf_polygons = gpd.GeoDataFrame([], geometry=gdf_polygons_List)
    gdf_polygons['_id'] = [filename for i in range(len(gdf_polygons))]
    if (len(gdf_polygons) > 0):
        gdf_polygons['area'] = [gdf_polygons.iloc[i]['geometry'].area for i in range(len(gdf_polygons))]
    gdf_polygons = gdf_polygons[gdf_polygons["area"] > 20]
    return gdf_polygons


def main(path):
    "img 0/1两灰度值图像"
    img = cv2.imread(path,0)
    filename = "tmp"
    if np.sum(img):
       CD = extract_lane_by_sheet(img, filename)
       if len(CD):
          CD.to_file(filename, driver='GeoJSON')