# coding=utf-8
#author:        4N
#createtime:    2022/3/8
#email:         nheweijun@sina.com


from osgeo import ogr
from osgeo.ogr import *
import math
import os
import copy
os.environ["CUDA_VISIBLE_DEVICES"] = "-1"

from test.zonghe.ShapeData import ShapeData
from test.zonghe.Delauney import DelauneyNode,DelauneyLine,DelauneyMinTree
from test.zonghe.RepresentPoint import RepresentPoint
from test.zonghe.Util import Util
import cv2

class Zonghe:

    data:ShapeData
    area_threshold = 100.0
    buffer_threshold = 100.0
    distance_buffer_threshold = 2.114691025217302e-04

    def __init__(self,data):
        self.data = data


    #只处理Polygon，MultiPolygon要拆开
    def get_polygon_reprecent_point(self,polygon:Geometry):

        polygon_env = polygon.GetEnvelope()

        area_thre = polygon.GetArea()/self.area_threshold

        line_length = math.sqrt(math.pow(polygon_env[1] - polygon_env[0],2) + math.pow(polygon_env[3] - polygon_env[2],2))


        short_length = min(polygon_env[1]-polygon_env[0],polygon_env[3]-polygon_env[2])



        center:Geometry = polygon.Centroid()
        # print(center)


        # 提前return
        if polygon.Contains(center):
            # print("\"" + center.ExportToWkt() + "\",")
            return RepresentPoint(center, polygon)


        last_state = None
        line = None
        calculate_time = 0

        direction = self.get_direction(polygon, center, line_length)


        in_hole = False
        in_hole_line = None


        #无孔Polygon
        if polygon.GetGeometryCount() == 1:
            polygon_line = ogr.ForceToLineString(polygon)

        # 有孔Polygon,要分2种情况
        else:
            polygon_line = ogr.ForceToLineString(polygon.GetGeometryRef(0))
            for i in range(1,polygon.GetGeometryCount()):
                hole_line = polygon.GetGeometryRef(i)
                hole_polygon : Geometry = ogr.ForceToPolygon(hole_line)
                if hole_polygon.Contains(center):
                    in_hole = True
                    direction = range(0,190,10)
                    in_hole_line = ogr.ForceToLineString(hole_line)
                    # print("in_hole")

        angle1,angle2=None,None

        for angle in direction:

            if in_hole:
                line = self.get_doubleline(center.GetX(),center.GetY(),angle,line_length)
            else:
                line = self.get_line(center.GetX(),center.GetY(),angle,line_length)

            oneside_area, otherside_area = self.jude(polygon, line, angle, line_length)


            if last_state is None:

                angle2 = angle

                if oneside_area == 0 and otherside_area==0:
                    last_state = 0
                else:
                    last_state = 1 if oneside_area < otherside_area else -1
            else:

                angle1 = angle2
                angle2 = angle

                if oneside_area == 0 and otherside_area==0:
                    now_state = 0
                else:
                    now_state = 1 if oneside_area < otherside_area else -1

                if last_state* now_state < 0:
                    # print(angle)
                    break
                else:
                    last_state = now_state

        line = self.division(polygon, center, in_hole, line_length, angle1, angle2, area_thre)

        # 提前return
        if in_hole:

            point_intersection: Geometry = line.Intersection(in_hole_line)
            plist = []
            for i in range(point_intersection.GetGeometryCount()):
                pi = point_intersection.GetGeometryRef(i)
                plist.append(pi)
            plist_d = [(pi,center.Distance(pi)) for pi in plist]
            plist_d.sort(key=lambda x:x[1])

            nearest_point= plist_d[0][0]
            firt_short_length = short_length/20.0
            radline,other_point = self.get_radline(center.GetX(),center.GetY(),nearest_point.GetX(),nearest_point.GetY(),firt_short_length)
            # print("\"" + other_point.ExportToWkt() + "\",")


            while not polygon.Contains(other_point):

                firt_short_length = firt_short_length / 2
                radline, other_point = self.get_radline(center.GetX(), center.GetY(), nearest_point.GetX(),
                                                        nearest_point.GetY(), firt_short_length)



            return RepresentPoint(other_point,polygon)


        point_intersection:Geometry = line.Intersection(polygon_line)
        xlist = []
        ylist = []
        plist = []

        for i in range(point_intersection.GetGeometryCount()):
            pi = point_intersection.GetGeometryRef(i)
            xlist.append(pi.GetX())
            ylist.append(pi.GetY())
            plist.append(pi)

        reprecent_point = [(min(xlist)+max(xlist))/2.0,(min(ylist)+max(ylist))/2.0]

        reprecent_point = self.creat_point(reprecent_point)

        if not polygon.Contains(reprecent_point):
            plist_d = [(pi,center.Distance(pi)) for pi in plist]
            plist_d.sort(key=lambda x:x[1])
            reprecent_point = self.creat_point(((plist_d[0][0].GetX()+ plist_d[1][0].GetX())/2.0, (plist_d[0][0].GetY()+plist_d[1][0].GetY())/2.0))

        # print("\""+reprecent_point.ExportToWkt() + "\",")

        rp = RepresentPoint(reprecent_point,polygon)

        return rp





    def get_line(self, x0, y0 , angle, r):

        line: Geometry = ogr.Geometry(ogr.wkbLineString)
        line.AddPoint(x0, y0)
        hudu = angle/360.0 * 2 * math.pi
        dx = math.sin(hudu) * r
        dy = math.cos(hudu) * r
        line.AddPoint(x0+dx, y0+dy)

        return line

    def get_doubleline(self, x0, y0 , angle, r):

        line: Geometry = ogr.Geometry(ogr.wkbLineString)

        hudu = angle/360.0 * 2 * math.pi
        dx = math.sin(hudu) * r
        dy = math.cos(hudu) * r

        hudu = (angle+180)/360.0 * 2 * math.pi
        dx2 = math.sin(hudu) * r
        dy2 = math.cos(hudu) * r

        line.AddPoint(x0 + dx2, y0 + dy2)
        line.AddPoint(x0+dx, y0+dy)

        return line


    def get_radline(self,x0,y0,x1,y1,len_size):
        a2 = math.pow((x1-x0),2)
        b2 = math.pow((y1-y0),2)
        len_size = len_size+ math.sqrt(a2+b2)

        line: Geometry = ogr.Geometry(ogr.wkbLineString)
        line.AddPoint(x0, y0)
        tanx = abs(y1-y0) / abs(x1-x0)
        cosx = 1 / math.sqrt(1 + tanx *tanx)
        dx = cosx* len_size
        dy = tanx * dx

        if x1 == x0:
            x2 = x0
            if y1 > y0 :
                y2 = y0 + abs(dy)
            else:
                y2 = y0 - abs(dy)
        else:
            if x1 < x0:
                x2 = x0 - dx
            else:
                x2 = x0 + dx

            if y1 < y0:
                y2 = y0 - dy
            else:
                y2 = y0 + dy

        line.AddPoint(x2, y2)

        other_point = ogr.Geometry(ogr.wkbPoint)
        other_point.AddPoint(x2, y2)

        return line , other_point

    def get_direction(self,polygon,center,line_length):
        '''
        判断旋转方向
        :param polygon:
        :param center:
        :param line_length:
        :return:
        '''

        line = self.get_line(center.GetX(),center.GetY(),0,line_length)

        oneside_area, otherside_area = self.jude(polygon,line,0,line_length)

        if oneside_area == 0 and otherside_area == 0 :
            return range(0, 390, 30)
        else:
            line = self.get_line(center.GetX(), center.GetY(), 10, line_length)
            oneside_area_contrast, otherside_area_contrast = self.jude(polygon, line, 30, line_length)

            if oneside_area_contrast == 0 and otherside_area_contrast == 0:
                return range(360, -30, -30)
            else:
                if abs(oneside_area - otherside_area) > abs(oneside_area_contrast - otherside_area_contrast):
                    return range(0, 390, 30)
                else:
                    return range(360, -30, -30)


    def jude(self,polygon,line,angle,line_length):

        line_buffer: Geometry = line.Buffer(line_length / self.buffer_threshold)
        clip_geom: Geometry = polygon.Difference(line_buffer)

        oneside_area = 0
        otherside_area = 0

        gc = clip_geom.GetGeometryCount()
        if gc <= 1:
            oneside_area = 0
            otherside_area = 0
        else:
            triangle1, triangle2 = self.get_side_triangle(line, angle)


            for gi in range(clip_geom.GetGeometryCount()):
                clip_geom_i: Geometry = clip_geom.GetGeometryRef(gi)
                it = clip_geom_i.Intersection(triangle1)

                if clip_geom_i.Intersect(triangle1) and it.Buffer(line_length / self.buffer_threshold).Intersect(
                        line_buffer):
                    oneside_area += clip_geom_i.GetArea()
                else:
                    otherside_area += clip_geom_i.GetArea()
        return oneside_area, otherside_area

    def division(self,polygon,center,in_hole,line_length,angle1,angle2,area_thre):

        mid_angle = (angle1 + angle2)/2.0
        if in_hole:
            line = self.get_doubleline(center.GetX(), center.GetY(), mid_angle, line_length)
        else:
            line = self.get_line(center.GetX(), center.GetY(), mid_angle, line_length)

        one,other = self.jude(polygon,line,mid_angle,line_length)

        while abs(one - other) > area_thre:
            if one > other:
                if angle1 > angle2:
                    return self.division(polygon, center, in_hole, line_length, angle2, mid_angle, area_thre)
                else:
                    return self.division(polygon, center, in_hole, line_length, mid_angle, angle1, area_thre)
            else:
                if angle1 > angle2:
                    return self.division(polygon, center, in_hole, line_length, mid_angle, angle1, area_thre)
                else:
                    return self.division(polygon, center, in_hole, line_length, angle2, mid_angle, area_thre)
        return line


    def get_side_triangle(self,line:Geometry,angle):
        '''
        获取方向三角形
        :param line:
        :param angle:
        :return:
        '''

        start = line.GetPoint(0)
        end = line.GetPoint(1)


        angle_t = angle-30
        hudu = angle_t/360.0 * 2 * math.pi

        r = line.Length() / (2*math.cos( math.pi/6))
        dx = math.sin(hudu) * r
        dy = math.cos(hudu) * r


        ring1 = ogr.Geometry(ogr.wkbLinearRing)
        ring1.AddPoint(start[0], start[1])
        ring1.AddPoint(start[0] + dx, start[1] + dy)
        ring1.AddPoint(end[0], end[1])
        ring1.AddPoint(start[0], start[1])
        triangle1 = ogr.Geometry(ogr.wkbPolygon)
        triangle1.AddGeometry(ring1)


        angle_t = angle+30
        hudu = angle_t/360.0 * 2 * math.pi

        r = line.Length() / (2*math.cos( math.pi/6))

        dx = math.sin(hudu) * r
        dy = math.cos(hudu) * r

        ring2 = ogr.Geometry(ogr.wkbLinearRing)
        ring2.AddPoint(start[0], start[1])
        ring2.AddPoint(start[0]+dx, start[1]+dy)
        ring2.AddPoint(end[0], end[1])
        ring2.AddPoint(start[0], start[1])
        triangle2 = ogr.Geometry(ogr.wkbPolygon)
        triangle2.AddGeometry(ring2)

        return triangle1,triangle2


    def creat_point(self,tuple):
        p: Geometry = ogr.Geometry(ogr.wkbPoint)
        p.AddPoint(tuple[0], tuple[1])
        return p




    def create_delauney(self,rps):
        multipoint: Geometry = ogr.Geometry(ogr.wkbMultiPoint)

        rp_dict = {}
        line_set = set()
        delauney_lines = []

        for rp in rps:
            multipoint.AddGeometry(rp.point)

            rp_dict[rp.point.GetPoint(0).__str__()] = rp

        delauney:Geometry =  multipoint.DelaunayTriangulation()

        for index in range(delauney.GetGeometryCount()):
            triangle = delauney.GetGeometryRef(index)

            triangle_line:Geometry = ogr.ForceToLineString(triangle)


            for pi in range(triangle_line.GetPointCount()-1):
                p = triangle_line.GetPoint(pi)
                p_next = triangle_line.GetPoint(pi+1)

                if p[0] < p_next[0]:
                    key = "{}{}".format(p.__str__(),p_next.__str__())
                elif p[0] > p_next[0]:
                    key = "{}{}".format(p_next.__str__(), p.__str__())
                else:
                    if p[1] < p_next[1]:
                        key = "{}{}".format(p.__str__(), p_next.__str__())
                    else:
                        key = "{}{}".format(p_next.__str__(), p.__str__())

                if not line_set.__contains__(key):
                    line_set.add(key)
                    rp1 = rp_dict.get(p.__str__())
                    rp2 = rp_dict.get(p_next.__str__())

                    delauney_line = DelauneyLine(rp1,rp2)
                    delauney_lines.append(delauney_line)

        return delauney_lines



    def create_min_delauney(self,delauney_lines,rps,polygons):
        '''
        使用prim算法计算最小生成树
        :param delauney_lines:
        :param polygons:
        :return:
        '''

        kv = {}
        rps_set = set(rps)

        for delauney_line in delauney_lines:
            if kv.get(delauney_line.rp1):

                kv[delauney_line.rp1].add(delauney_line)
            else:
                kv[delauney_line.rp1] = set()
                kv[delauney_line.rp1].add(delauney_line)

            if kv.get(delauney_line.rp2):

                kv[delauney_line.rp2].add(delauney_line)
            else:

                kv[delauney_line.rp2] = set()
                kv[delauney_line.rp2].add(delauney_line)


        #初始化树
        delauney_lines = sorted(delauney_lines,key=lambda x:x.distance)


        delauney_min_tree = DelauneyMinTree()
        delauney_min_tree.append(delauney_lines[0])

        rps_set.remove(delauney_lines[0].rp1)
        rps_set.remove(delauney_lines[0].rp2)
        kv[delauney_lines[0].rp1].remove(delauney_lines[0])
        kv[delauney_lines[0].rp2].remove(delauney_lines[0])


        while rps_set.__len__()>0:

            relate_delauney_line = delauney_min_tree.get_other_relate_delauney_lines(kv)

            relate_delauney_line_sort = sorted(relate_delauney_line,key=lambda x:x.distance)

            min_distance_line = relate_delauney_line_sort[0]

            #删除
            other_rp = delauney_min_tree.get_other_rp(min_distance_line)
            if other_rp is None:
                kv[min_distance_line.rp1].remove(min_distance_line)
                kv[min_distance_line.rp2].remove(min_distance_line)
                continue

            delauney_min_tree.append(min_distance_line)

            rps_set.remove(other_rp)

            kv[min_distance_line.rp1].remove(min_distance_line)
            kv[min_distance_line.rp2].remove(min_distance_line)



        for l in delauney_min_tree.delauney_lines_list:
            line = l.line
            # print("\"" + line.ExportToWkt() + "\",")

        # print(delauney_min_tree.delauney_lines_list.__len__())

        return delauney_min_tree

    def cut_delauney_min_tree(self,delauney_min_tree:DelauneyMinTree,threshold):

        trees = []
        kv = {}

        count = 0

        for delauney_line in delauney_min_tree.delauney_lines_list:

            dn1 = kv.get(delauney_line.rp1)
            dn2 = kv.get(delauney_line.rp2)


            if not dn1:
                dn1 = DelauneyNode(delauney_line.rp1)
                kv[delauney_line.rp1] = dn1

                if not dn2:

                    trees.append(dn1)

                    dn2 = DelauneyNode(delauney_line.rp2)
                    kv[delauney_line.rp2] = dn2

                    if delauney_line.distance < threshold:
                        dn1.add_children(dn2,delauney_line)
                        count +=1

                    else:
                        #断开
                        trees.append(dn2)
                else:

                    if delauney_line.distance < threshold:
                        dn2.add_children(dn1,delauney_line)
                        count += 1
                    else:
                        trees.append(dn1)
            else:

                dn2 = DelauneyNode(delauney_line.rp2)
                kv[delauney_line.rp2] = dn2

                if delauney_line.distance < threshold:
                    dn1.add_children(dn2,delauney_line)
                    count += 1
                else:
                    trees.append(dn2)


        return trees

    def merge(self,trees):

        polygons = []
        for tree in trees:

            merge_polygon_raw = list()

            dls = tree.get_delauney_lines()
            for dl in dls:
                p1 = dl.rp1.polygon
                p2 = dl.rp2.polygon

                merge_polygon_raw.append(p1)
                merge_polygon_raw.append(p2)
                line = dl.line

                p1_lines = self.get_polygon_lines(p1)
                p2_lines = self.get_polygon_lines(p2)

                p1_intersect_lines = [l for l in p1_lines if l.Intersect(line)]
                p2_intersect_lines = [l for l in p2_lines if l.Intersect(line)]

                if p1_intersect_lines.__len__()>1:
                    p1_intersect_lines = sorted(p1_intersect_lines,key=lambda x:p2.Distance(x))
                    p1_intersect_line_rep = p1_intersect_lines[0]
                elif p1_intersect_lines.__len__()==1:
                    p1_intersect_line_rep = p1_intersect_lines[0]
                else:
                    continue


                if p2_intersect_lines.__len__()>1:
                    p2_intersect_lines = sorted(p2_intersect_lines,key=lambda x:p1.Distance(x))
                    p2_intersect_line_rep = p2_intersect_lines[0]
                elif p2_intersect_lines.__len__()==1:
                    p2_intersect_line_rep = p2_intersect_lines[0]
                else:
                    continue


                # #两条线没有距离跳过
                # if p2_intersect_line_rep.Distance(p1_intersect_line_rep) == 0:
                #     continue
                # else:

                mid_polygon = self.create_mid_polygon(p1_intersect_line_rep,p2_intersect_line_rep,p1,p2)

                if mid_polygon.GetArea() > 0.0002116361000058077*0.0002116361000058077:
                    continue

                merge_polygon_raw.append(mid_polygon)

            merge_polygon = None
            for p in merge_polygon_raw:
                if merge_polygon is None:
                    merge_polygon = p
                else:
                    merge_polygon = merge_polygon.Union(p)

            polygons.append(merge_polygon)


        return polygons


    def merge2(self,trees):

        polygons = []




        for tree in trees:

            tree_polygons = tree.get_polygons()

            points = []
            for p in tree_polygons:

                if p.GetGeometryCount() == 1:
                    polygon_line: Geometry = ogr.ForceToLineString(p)

                    polygon_line.Con
                # 有孔Polygon
                else:
                    polygon_line: Geometry = ogr.ForceToLineString(p.GetGeometryRef(0))

                points.extend(polygon_line.GetPoints()[:-1])

            multipoint: Geometry = ogr.Geometry(ogr.wkbMultiPoint)

            for p in  points:
                poi = ogr.Geometry(ogr.wkbPoint)
                poi.AddPoint(p[0],p[1])
                multipoint.AddGeometry(poi)

            delauney: Geometry = multipoint.DelaunayTriangulation()
            merge_polygon = None

            for p in delauney:
                # if multi_polygon.Contains(p):
                #     continue

                pline:Geometry = ogr.ForceToLineString(p)
                pline :list = pline.GetPoints()
                pline.append(pline[1])


                pass_or_not = False
                for indx in range(3):
                    p1 = pline[indx]
                    p2 = pline[indx+1]

                    ppline:Geometry = ogr.Geometry(ogr.wkbLineString)
                    ppline.AddPoint(p1[0],p1[1])
                    ppline.AddPoint(p2[0], p2[1])



                    p3:Geometry = ogr.Geometry(ogr.wkbPoint)
                    p3.AddPoint(pline[indx+2][0], pline[indx+2][1])

                    if p3.Distance(ppline) > 0.0002116361000058077 or ppline.Length() > 0.0002116361000058077:
                        pass_or_not = True
                        break

                if pass_or_not:
                    continue

                if merge_polygon is None:
                    merge_polygon = p
                else:
                    merge_polygon = merge_polygon.Union(p)



            polygons.append(merge_polygon)


        return polygons

    def merge_circle(self):
        pass


    def create_mid_polygon(self,p1_intersect_line_rep,p2_intersect_line_rep,p1,p2):
        '''
        创建衔接polygon
        :param line:
        :return:
        '''

        if p1_intersect_line_rep.Intersect(p2_intersect_line_rep):

            p1_rep_p1 = p1_intersect_line_rep.GetPoints()[0]
            p1_rep_p2 = p1_intersect_line_rep.GetPoints()[1]

            p2_rep_p1 = p2_intersect_line_rep.GetPoints()[0]
            p2_rep_p2 = p2_intersect_line_rep.GetPoints()[1]

            polygon = self.create_polygon(p1_rep_p1, p2_rep_p2, p1_rep_p2, p2_rep_p1)

        else:

            p1_rep_p1 = p1_intersect_line_rep.GetPoints()[0]
            p1_rep_p2 = p1_intersect_line_rep.GetPoints()[1]

            p2_rep_p1 = p2_intersect_line_rep.GetPoints()[0]
            p2_rep_p2 = p2_intersect_line_rep.GetPoints()[1]

            line1 = ogr.Geometry(ogr.wkbLineString)


            line1.AddPoint(p1_rep_p1[0], p1_rep_p1[1])
            line1.AddPoint(p2_rep_p1[0], p2_rep_p1[1])

            line2 = ogr.Geometry(ogr.wkbLineString)
            line2.AddPoint(p1_rep_p2[0], p1_rep_p2[1])
            line2.AddPoint(p2_rep_p2[0], p2_rep_p2[1])

            if line1.Intersect(line2):
                polygon = self.create_polygon(p1_rep_p1, p1_rep_p2, p2_rep_p1, p2_rep_p2)
            else:

                polygon = self.create_polygon(p1_rep_p1, p1_rep_p2, p2_rep_p2, p2_rep_p1)
                if polygon.Intersect(p1) or polygon.Intersect(p2):
                    polygon = self.create_polygon(p1_rep_p1, p1_rep_p2, p2_rep_p1, p2_rep_p2)


        return polygon

    def create_polygon(self,p1,p2,p3,p4):

        ring = ogr.Geometry(ogr.wkbLinearRing)
        ring.AddPoint(p1[0], p1[1])
        ring.AddPoint(p2[0], p2[1])
        ring.AddPoint(p3[0], p3[1])
        ring.AddPoint(p4[0], p4[1])
        ring.AddPoint(p1[0], p1[1])

        poly = ogr.Geometry(ogr.wkbPolygon)
        poly.AddGeometry(ring)
        return poly


    def get_polygon_lines(self,polygon):
        #无孔Polygon
        if polygon.GetGeometryCount() == 1:
            polygon_line : Geometry = ogr.ForceToLineString(polygon)
        # 有孔Polygon
        else:
            polygon_line : Geometry = ogr.ForceToLineString(polygon.GetGeometryRef(0))

        points = polygon_line.GetPoints()

        lines = []
        for index in range(len(points)-1):

            point_start = points[index]
            point_end = points[index+1]
            line = ogr.Geometry(ogr.wkbLineString)
            line.AddPoint(point_start[0], point_start[1])
            line.AddPoint(point_end[0], point_end[1])

            lines.append(line)

        return lines

    def process(self):
        polygons = self.data.get_polygons()
        rps = []
        for poly in polygons:
            rp = self.get_polygon_reprecent_point(poly)
            rps.append(rp)

        delauney_lines = self.create_delauney(rps)

        min_delauney = self.create_min_delauney(delauney_lines,rps,polygons)

        trees = self.cut_delauney_min_tree(min_delauney,self.distance_buffer_threshold)


        polygons = self.merge(trees)

        for polygon in polygons:
            print("\""+polygon.ExportToWkt() + "\",")

if __name__ == '__main__':


    # print(sd.get_polygons()[0])
    # sd.close()

    create = True
    create = False

    if create:

        sd = ShapeData(r"J:\Data\制图综合\example.shp")

        ext =  sd.layer.GetExtent()
        threshold = (ext[1]-ext[0])/20
        print(threshold)


        zh = Zonghe(sd)
        zh.process()
    else:
        wkts = [
   ]

        result = r"J:\Data\制图综合\zhonghe3.shp"
        ShapeData.create_shp_fromwkts(result,"zh",wkts)