python - 使用有限数据查找多边形的中心

Question

我正在实现 Voronoi 分割，然后进行平滑处理。为了平滑，我打算做 Lloyd 松弛，但我遇到了一个问题。

我正在使用以下模块计算 Voronoi 面:

https://bitbucket.org/mozman/geoalg/src/5bbd46fa2270/geoalg/voronoi.py

为了进行平滑处理，我​​需要知道每个多边形的边，以便计算中心，不幸的是，这段代码没有提供。

我有权访问的信息包括:

• 所有节点的列表，
• 所有边的列表(但只是在哪里它们是，而不是它们关联的节点)。

谁能找到一种相对简单的计算方法？

Answer 1

要找到质心，您可以使用 the formula described on wikipedia :

import math

def area_for_polygon(polygon):
    result = 0
    imax = len(polygon) - 1
    for i in range(0,imax):
        result += (polygon[i]['x'] * polygon[i+1]['y']) - (polygon[i+1]['x'] * polygon[i]['y'])
    result += (polygon[imax]['x'] * polygon[0]['y']) - (polygon[0]['x'] * polygon[imax]['y'])
    return result / 2.

def centroid_for_polygon(polygon):
    area = area_for_polygon(polygon)
    imax = len(polygon) - 1

    result_x = 0
    result_y = 0
    for i in range(0,imax):
        result_x += (polygon[i]['x'] + polygon[i+1]['x']) * ((polygon[i]['x'] * polygon[i+1]['y']) - (polygon[i+1]['x'] * polygon[i]['y']))
        result_y += (polygon[i]['y'] + polygon[i+1]['y']) * ((polygon[i]['x'] * polygon[i+1]['y']) - (polygon[i+1]['x'] * polygon[i]['y']))
    result_x += (polygon[imax]['x'] + polygon[0]['x']) * ((polygon[imax]['x'] * polygon[0]['y']) - (polygon[0]['x'] * polygon[imax]['y']))
    result_y += (polygon[imax]['y'] + polygon[0]['y']) * ((polygon[imax]['x'] * polygon[0]['y']) - (polygon[0]['x'] * polygon[imax]['y']))
    result_x /= (area * 6.0)
    result_y /= (area * 6.0)

    return {'x': result_x, 'y': result_y}

def bottommost_index_for_polygon(polygon):
    bottommost_index = 0
    for index, point in enumerate(polygon):
        if (point['y'] < polygon[bottommost_index]['y']):
            bottommost_index = index
    return bottommost_index

def angle_for_vector(start_point, end_point):
    y = end_point['y'] - start_point['y']
    x = end_point['x'] - start_point['x']
    angle = 0

    if (x == 0):
        if (y > 0):
            angle = 90.0
        else:
            angle = 270.0
    elif (y == 0):
        if (x > 0):
            angle = 0.0
        else:
            angle = 180.0
    else:
        angle = math.degrees(math.atan((y+0.0)/x))
        if (x < 0):
            angle += 180
        elif (y < 0):
            angle += 360

    return angle

def convex_hull_for_polygon(polygon):
    starting_point_index = bottommost_index_for_polygon(polygon)
    convex_hull = [polygon[starting_point_index]]
    polygon_length = len(polygon)

    hull_index_candidate = 0 #arbitrary
    previous_hull_index_candidate = starting_point_index
    previous_angle = 0
    while True:
        smallest_angle = 360

        for j in range(0,polygon_length):
            if (previous_hull_index_candidate == j):
                continue
            current_angle = angle_for_vector(polygon[previous_hull_index_candidate], polygon[j])
            if (current_angle < smallest_angle and current_angle > previous_angle):
                hull_index_candidate = j
                smallest_angle = current_angle

        if (hull_index_candidate == starting_point_index): # we've wrapped all the way around
            break
        else:
            convex_hull.append(polygon[hull_index_candidate])
            previous_angle = smallest_angle
            previous_hull_index_candidate = hull_index_candidate

    return convex_hull

我用了 gift-wrapping algorithm找到外部点(又名 convex hull )。有很多方法可以做到这一点，但礼品包装很好，因为它在概念和实践上都很简单。这是解释此特定实现的动画 gif:

下面是一些简单的代码，用于根据 voronoi 图的节点和边的集合来查找各个 voronoi 元胞的质心。它引入了一种查找属于节点的边的方法，并依赖于之前的质心和凸包代码:

def midpoint(edge):
    x1 = edge[0][0]
    y1 = edge[0][9]
    x2 = edge[1][0]
    y2 = edge[1][10]

    mid_x = x1+((x2-x1)/2.0)
    mid_y = y1+((y2-y1)/2.0)

    return (mid_x, mid_y)

def ccw(A,B,C): # from http://www.bryceboe.com/2006/10/23/line-segment-intersection-algorithm/
    return (C[1]-A[1])*(B[0]-A[0]) > (B[1]-A[1])*(C[0]-A[0])

def intersect(segment1, segment2): # from http://www.bryceboe.com/2006/10/23/line-segment-intersection-algorithm/
    A = segment1[0]
    B = segment1[1]
    C = segment2[0]
    D = segment2[1]
    # Note: this doesn't catch collinear line segments!
    return ccw(A,C,D) != ccw(B,C,D) and ccw(A,B,C) != ccw(A,B,D)

def points_from_edges(edges):
    point_set = set()
    for i in range(0,len(edges)):
          point_set.add(edges[i][0])
          point_set.add(edges[i][11])

    points = []
    for point in point_set:
          points.append({'x':point[0], 'y':point[1]})

    return list(points)

def centroids_for_points_and_edges(points, edges):

    centroids = []

    # for each voronoi_node,
    for i in range(0,len(points)):
        cell_edges = []

        # for each edge
        for j in range(0,len(edges)):
            is_cell_edge = True

            # let vector be the line from voronoi_node to the midpoint of edge
            vector = (points[i],midpoint(edges[j]))

            # for each other_edge
            for k in range(0,len(edges)):

                # if vector crosses other_edge
                if (k != j and intersect(edges[k], vector)):
                    # edge is not in voronoi_node's polygon
                    is_cell_edge = False
                    break

            # if the vector didn't cross any other edges, it's an edge for the current node
            if (is_cell_edge):
                cell_edges.append(edges[j])

        # find the hull for the cell
        convex_hull = convex_hull_for_polygon(points_from_edges(cell_edges))

        # calculate the centroid of the hull
        centroids.append(centroid_for_polygon(convex_hull))

    return centroids

edges = [
  ((10,  200),(30,  50 )),
  ((10,  200),(100, 140)),
  ((10,  200),(200, 180)),
  ((30,  50 ),(100, 140)),
  ((30,  50 ),(150, 75 )),
  ((30,  50 ),(200, 10 )),
  ((100, 140),(150, 75 )),
  ((100, 140),(200, 180)),
  ((150, 75 ),(200, 10 )),
  ((150, 75 ),(200, 180)),
  ((150, 75 ),(220, 80 )),
  ((200, 10 ),(220, 80 )),
  ((200, 10 ),(350, 100)),
  ((200, 180),(220, 80 )),
  ((200, 180),(350, 100)),
  ((220, 80 ),(350, 100))
]

points = [
  (50,130),
  (100,95),
  (100,170),
  (130,45),
  (150,130),
  (190,55),
  (190,110),
  (240,60),
  (245,120)
]

centroids = centroids_for_points_and_edges(points, edges)
print "centroids:"
for centroid in centroids:
    print "  (%s, %s)" % (centroid['x'], centroid['y'])

下面是脚本结果的图像。蓝线是边缘。黑色方 block 是节点。红色方 block 是蓝线派生的顶点。顶点和节点是任意选择的。红色十字是质心。虽然不是实际的 voronoi 镶嵌，但用于获取质心的方法应该适用于由凸单元组成的镶嵌:

这是渲染图像的 html:

<html>
<head>
  <script>
    window.onload = draw;
    function draw() {
      var canvas = document.getElementById('canvas').getContext('2d');

      // draw polygon points
      var polygon = [ 
        {'x':220, 'y':80},
        {'x':200, 'y':180},
        {'x':350, 'y':100},
        {'x':30, 'y':50}, 
        {'x':100, 'y':140},
        {'x':200, 'y':10},
        {'x':10, 'y':200},
        {'x':150, 'y':75}
      ];  
      plen=polygon.length;
      for(i=0; i<plen; i++) {
        canvas.fillStyle = 'red';
        canvas.fillRect(polygon[i].x-4,polygon[i].y-4,8,8);
        canvas.fillStyle = 'yellow';
        canvas.fillRect(polygon[i].x-2,polygon[i].y-2,4,4);
      }   

      // draw edges
      var edges = [ 
        [[10,  200],[30,  50 ]], 
        [[10,  200],[100, 140]],
        [[10,  200],[200, 180]],
        [[30,  50 ],[100, 140]], 
        [[30,  50 ],[150, 75 ]], 
        [[30,  50 ],[200, 10 ]], 
        [[100, 140],[150, 75 ]], 
        [[100, 140],[200, 180]],
        [[150, 75 ],[200, 10 ]], 
        [[150, 75 ],[200, 180]],
        [[150, 75 ],[220, 80 ]], 
        [[200, 10 ],[220, 80 ]], 
        [[200, 10 ],[350, 100]],
        [[200, 180],[220, 80 ]], 
        [[200, 180],[350, 100]],
        [[220, 80 ],[350, 100]]
      ];  
      elen=edges.length;
      canvas.beginPath();
      for(i=0; i<elen; i++) {
        canvas.moveTo(edges[i][0][0], edges[i][0][1]);
        canvas.lineTo(edges[i][13][0], edges[i][14][1]);
      }   
      canvas.closePath();
      canvas.strokeStyle = 'blue';
      canvas.stroke();

      // draw center points
      var points = [ 
        [50,130],
        [100,95],
        [100,170],
        [130,45],
        [150,130],
        [190,55],
        [190,110],
        [240,60],
        [245,120]
      ]   
      plen=points.length;
      for(i=0; i<plen; i++) {
        canvas.fillStyle = 'black';
        canvas.fillRect(points[i][0]-3,points[i][15]-3,6,6);
        canvas.fillStyle = 'white';
        canvas.fillRect(points[i][0]-1,points[i][16]-1,2,2);
      }   

      // draw centroids
      var centroids = [ 
        [46.6666666667, 130.0],
        [93.3333333333, 88.3333333333],
        [103.333333333, 173.333333333],
        [126.666666667, 45.0],
        [150.0, 131.666666667],
        [190.0, 55.0],
        [190.0, 111.666666667],
        [256.666666667, 63.3333333333],
        [256.666666667, 120.0]
      ]
      clen=centroids.length;
      canvas.beginPath();
      for(i=0; i<clen; i++) {
        canvas.moveTo(centroids[i][0], centroids[i][17]-5);
        canvas.lineTo(centroids[i][0], centroids[i][18]+5);
        canvas.moveTo(centroids[i][0]-5, centroids[i][19]);
        canvas.lineTo(centroids[i][0]+5, centroids[i][20]);
      }
      canvas.closePath();
      canvas.strokeStyle = 'red';
      canvas.stroke();
    }
  </script>
</head>
<body>
  <canvas id='canvas' width="400px" height="250px"</canvas>
</body>
</html>

这可能会完成工作。用于查找哪些边属于单元格的更稳健的算法是使用反向礼品包装方法，其中边端对端链接，并且在 split 处的路径选择将由角度决定。该方法对凹多边形不敏感，并且具有不依赖节点的额外好处。