Voronoi泰森多边形应用案例-python

Voronoi泰森多边形应用案例-python

下面代码为Voronoi泰森多边形的应用案例。在下面代码中，我们会对输入的数据集，计算泰森多边形，之后在这个基础上，再对泰森多边形进行分割，之后我们可以通过书输入点，然后直接计算出其在我们分割的那个框里，再次基础上，结合泰森多边形的边界线知识，就可以很快的判断输入点的坐标距离输入数据集那个点最近。

import random 
import math 
from scipy.spatial import Voronoi, voronoi_plot_2d
import os

from scipy.spatial import cKDTree
import numpy as np
import matplotlib .pyplot as plt
import cv2
from sklearn.datasets import make_classification

x_min = -100
x_max =100
a=-1
b=1
num=10
y_min = -100
y_max =100
divide_num=5

x_lmin = -2
x_lmax =2
y_lmin = -2
y_lmax =2
#判断点是否在框中
def judge_frame(point):

    x, y = point[0], point[1]
    flag = (x >= x_min and x <= x_max and y >= y_min and y <= y_max)
    return flag

# 判断两线段相交


def judge_inter(A, B, C, D):

    def vector(piont_1, point_2):   # 输入numpy格式
        return point_2 - piont_1

    def vector_product(vec_1, vec_2):   # 计算向量
        return vec_1[0] * vec_2[1] - vec_2[0] * vec_1[1]

    AC = vector(A, C)
    AD = vector(A, D)
    BC = vector(B, C)
    BD = vector(B, D)
    CA = vector(C, A)
    CB = vector(C, B)
    DA = vector(D, A)
    DB = vector(D, B)

    # 相交返回true
    return (vector_product(AC, AD) * vector_product(BC, BD) <= 0 ) \
           and (vector_product(CA, CB) * vector_product(DA, DB) <= 0)


# 判断得出两线段的交点
def intersection(point_1, point_2):
    point_bound = [[[x_min, y_min], [x_max, y_min]], [[x_max, y_min], [x_max, y_max]],
                   [[x_min, y_max], [x_max, y_max]], [[x_min, y_max], [x_min, y_min]]]

    for line in point_bound:
        if judge_inter(point_1, point_2, line[0], line[1]):
            x1, x2, x3, x4 = point_1[0], point_2[0], line[0][0], line[1][0]
            y1, y2, y3, y4 = point_1[1], point_2[1], line[0][1], line[1][1]
            px = ((x1 * y2 - y1 * x2) * (x3 - x4) - (x1 - x2) * (x3 * y4 - y3 * x4)) / ((x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4))
            py = ((x1 * y2 - y1 * x2) * (y3 - y4) - (y1 - y2) * (x3 * y4 - y3 * x4)) / ((x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4))
            return [px, py]

# 求射线方程
def ray(barr, piont):

    # 计算障碍的中垂线
    def medLine(x1, y1, x2, y2):
        A = 2 * (x2 - x1)
        B = 2 * (y2 - y1)
        C = x1 ** 2 - x2 ** 2 + y1 ** 2 - y2 ** 2
        return A, B, C
    A, B, C = medLine(barr[0][0], barr[0][1], barr[1][0], barr[1][1])
    return [A, B, C]

#判断两线段是否相交
def Intersect(l1, l2):
    v1 = (l1[0] - l2[0], l1[1] - l2[1])
    v2 = (l1[0] - l2[2], l1[1] - l2[3])
    v0 = (l1[0] - l1[2], l1[1] - l1[3])
    a = v0[0] * v1[1] - v0[1] * v1[0]
    b = v0[0] * v2[1] - v0[1] * v2[0]

    temp = l1
    l1 = l2
    l2 = temp
    v1 = (l1[0] - l2[0], l1[1] - l2[1])
    v2 = (l1[0] - l2[2], l1[1] - l2[3])
    v0 = (l1[0] - l1[2], l1[1] - l1[3])
    c = v0[0] * v1[1] - v0[1] * v1[0]
    d = v0[0] * v2[1] - v0[1] * v2[0]

    if a*b < 0 and c*d < 0:
        return True
    else:
        return False

# 判断直线与线段相交并求出交点
def ray_inter_point(ray, point, barr):
    point_bound = [[[x_min, y_min], [x_max, y_min]], [[x_max, y_min], [x_max, y_max]],
                   [[x_min, y_max], [x_max, y_max]], [[x_min, y_max], [x_min, y_min]]]
    list_inter = []
    for bound in point_bound:
        sPoint, ePoint = bound[0], bound[1]
        line = [sPoint[1] - ePoint[1], ePoint[0] - sPoint[0], sPoint[0] * ePoint[1] - ePoint[0] * sPoint[1]]

        a0, b0, c0 = ray[0], ray[1], ray[2]
        a1, b1, c1 = line[0], line[1], line[2]

        D = a0 * b1 - a1 * b0
        x = (b0 * c1 - b1 * c0) / D
        y = (a1 * c0 - a0 * c1) / D

        list_inter.append([x, y])
    result=[]
    for i in list_inter:
            result.append(i)
    return(result)

path="C:/Users/gaoxi/Desktop/data/input.png"
# 使用 make_classification 创建数据

def save_point(self,vor):
		with open('C:/Users/gaoxi/Desktop/data/vor/input_vor.txt','w',encoding='utf-8') as file:
			 #vor.vertices
			 for i in  vor.vertices:
				 file.write(str(i[0])+"  "+str(i[1])+'\n')




#points=[]
#for i in range(num):
    
#        c=a + (b-a)*np.random.random()
#        d=a + (b-a)*np.random.random()
#        points.append([c,d])
#points = np.array(points)


points = np.load('C:/Users/gaoxi/Desktop/data/data/pose_det_results_list.npy')
barr_center = points


# 计算Voronoi图


vor = Voronoi(barr_center)




vertices, ridge_vertices = [], []

my_dict={}
new_dict = []

for k,v in vor.ridge_dict.items():
    my_dict[str(v[0])+str(v[1])]=k

    if -1 in v:
        new_dict.append(list(k))



print("my_dict",my_dict)
vers=[]

pots=[]

###   l1 [xa, ya, xb, yb]   l2 [xa, ya, xb, yb]

#求解voronio线段
for ver_rela in vor.ridge_vertices:
    print("ver_rela:",ver_rela)

    if -1 in ver_rela:
        continue
    else:
        # 两个都在
      
            if vor.vertices[ver_rela[0]].tolist() not in vertices:
                vertices.append(vor.vertices[ver_rela[0]].tolist())
            if vor.vertices[ver_rela[1]].tolist() not in vertices:
                vertices.append(vor.vertices[ver_rela[1]].tolist())
            pots.append(my_dict[str(ver_rela[0])+str(ver_rela[1])])
            vers.append([vor.vertices[ver_rela[0]].tolist(),vor.vertices[ver_rela[1]].tolist()])
        # 一个都不在
#求解voronio射线
for ver_rela in vor.ridge_vertices:
  #  print("ver_rela:",ver_rela)
    if -1 in ver_rela:
        pots.append(my_dict[str(ver_rela[0])+str(ver_rela[1])])
        temp = [0, 1]
        temp.pop(ver_rela.index(-1))
        if judge_frame(vor.vertices[ver_rela[temp[0]]])  or 1:
            # 与边界的交点
            barr_ind = new_dict.pop(0)
            points1=vor.vertices[ver_rela[temp[0]]].tolist()
            barr = [vor.points[barr_ind[0]], vor.points[barr_ind[1]]]
            f_ray = ray(barr, vor.vertices[ver_rela[temp[0]]])  # 改为直线方程
          #  print(f_ray)
            inter = ray_inter_point(f_ray, vor.vertices[ver_rela[temp[0]]], barr)
            points2_list=inter

            for i in points2_list:
                l1=[points1[0],points1[1],i[0],i[1]]
                flag=0
                #plt.plot([i[0], points1[0]],
                #               [i[1], points1[1]],
                #              color = 'r')
                for pointsz in vers:
                    l2=[pointsz[0][0],pointsz[0][1],pointsz[1][0],pointsz[1][1]]
                    if  Intersect(l1,l2) :
                        flag=1 
                        break                   
                if flag==0 :
                     # print("true")
                      vertices.append(i)
                      vers.append([points1,i])
                      break


            #print("points1:",points1)
            #print("points2:",points2_list)

       

vertices = np.array(vertices)


ridge_vertices = np.array(ridge_vertices)
#绘图
def test1():
    
    for i in range(num):
  
        plt.annotate(str(i), (vor.points[i][0], vor.points[i][1]), size = 10)

    for i in range(len(vor.vertices)):
  
        plt.annotate(str(i), (vor.vertices[i][0], vor.vertices[i][1]), size = 10,color='r')


    plt.xlim((-4, 4))
    plt.ylim((-4, 4))
    plt.scatter(vertices[:, 0], vertices[:, 1])
    plt.scatter(points[:, 0], points[:, 1])

    for pointsz in vers:
        plt.plot([pointsz[0][0], pointsz[1][0]],
                 [pointsz[0][1], pointsz[1][1]],
                 color = 'k')

  #  voronoi_plot_2d(vor)
   # print(vers)
    #print(pots)
#plt.show()
test1()################################################
x_list=np.linspace(x_lmin, x_lmax, num=divide_num)
y_list=np.linspace(y_lmin,  x_lmax, num=divide_num)
print("x_list",x_list)
print("y_list",y_list)

frame=[]


for i in range(len(x_list)-1):
    x1=x_list[i]
    x2=x_list[i+1]
    for j in range(len(y_list)-1):
        y1=y_list[j]
        y2=y_list[j+1]
        frame.append([[x1,y1],[x2,y1],[x2,y2],[x1,y2]])

print(frame)

def test2():
    for points in frame:
        for i in range(4):
            if i<3:
                 plt.plot([points[i][0], points[i+1][0]],
                         [points[i][1], points[i+1][1]],
                         color = 'k')
            else:
                  plt.plot([points[i][0], points[0][0]],
                         [points[i][1], points[0][1]],
                         color = 'k')
    plt.show()
###   l1 [xa, ya, xb, yb]   l2 [xa, ya, xb, yb]

index=[]


#检测框中通过那些射线或者线段
def  frame_points():
    for j in range(len(frame)):
        print(j)
        fra=frame[j]
        x_minz=fra[0][0]
        x_maxz=fra[1][0]
        y_minz=fra[0][1]
        y_maxz=fra[2][1]
        a=[]
        for i in range(len(vers)):
            #print(vers[i])
            posts=vers[i]
            l2=[posts[0][0],posts[0][1],posts[1][0],posts[1][1]]
            if  x_minz<posts[0][0] and posts[0][0]<x_maxz and x_minz<posts[1][0] and posts[1][0]<x_maxz:
                if y_minz<posts[0][1] and posts[0][1]<y_maxz and y_minz<posts[1][1] and posts[1][1]<y_maxz:
                   
                    a.append(pots[i][0])
                    a.append(pots[i][1])
                    continue
            for z in range(4):
                if z<3:
                     l1=[fra[z][0],fra[z][1],fra[z+1][0],fra[z+1][1]]
                     if Intersect(l1, l2):
                         a.append(pots[i][0])
                         a.append(pots[i][1])
                         break
                else:
                     l1=[fra[z][0],fra[z][1],fra[0][0],fra[0][1]]
                     if Intersect(l1, l2):
                         a.append(pots[i][0])
                         a.append(pots[i][1])
                         break
        if len(a)==0:
            x=(x_minz+x_maxz)/2
            y=(y_minz+y_maxz)/2
            min_distance=(points[0][0]-x)**2+(points[0][1]-y)**2
            index_min=0
            for k in range(1,len(points)):
                if  (points[k][0]-x)**2+(points[k][1]-y)**2<min_distance:
                    index_min=k
                    min_distance=(points[k][0]-x)**2+(points[k][1]-y)**2
            index.append([index_min])
        else:
              a=list(set(a))
              index.append(a)



frame_points()
po=[1.48,0.52]
#寻找查找点在那个框中
def find(po):
    position=0
    for i in range(len(x_list)-1):
        x1=x_list[i]
        x2=x_list[i+1]
        if x1<=po[0] and po[0]<=x2:

            for j in range(len(y_list)-1):
                y1=y_list[j]
                y2=y_list[j+1]
                print(y1,y2)
                if  y1<=po[1] and po[1]<=y2:
                    return position
                else: 
                    position=position+1
              #  print("position",position)
        else:
          #  print("sucess")
            position=position+len(x_list)-1
          #  print("position",position)
    
position=find(po)



print("position:",position)

print("index of position is:",index[position])


#寻找最近点得索引
def find_min_index(po,index_list):
            x=po[0]
            y=po[1]
            min_distance=(points[index_list[0]][0]-x)**2+(points[index_list[0]][1]-y)**2
            index_min=index_list[0]
            for k in range(1,len(index_list)):
                if  (points[index_list[k]][0]-x)**2+(points[index_list[k]][1]-y)**2<min_distance:
                    index_min=index_list[k]
                    min_distance=(points[index_list[k]][0]-x)**2+(points[index_list[k]][1]-y)**2
            return index_min


find_index=find_min_index(po,index[position])

print("find_index:",find_index)
test2()


os.system("pause")