利用有限的数据求多边形的中心

Question

我正在实现Voronoi tesselation，然后是平滑。为了平滑，我本来要做劳埃德放松，但我遇到了一个问题。

我使用以下模块计算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'])

下面是脚本结果的图像。蓝线是边缘。黑色方块是节点。红方格是由蓝线导出的顶点。顶点和节点的选择是任意的。红色的十字架是质心。虽然不是实际的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>

这很有可能完成任务。要找出哪些边缘属于一个单元格，一个更健壮的方法是使用反向礼品包装方法，其中边被连接到端到端，在分割处的路径选择将由角度来决定。这种方法不会有凹多边形的易感性，而且它还有不依赖于节点的额外好处。

Answer 2

这是@mgamba的答案，更像是python风格的重写。特别是，它在点上使用itertools.cycle，这样“一加最后一点”就可以以更自然的方式被看作是第一个点。

import itertools as IT

def area_of_polygon(x, y):
    """Calculates the signed area of an arbitrary polygon given its verticies
    http://stackoverflow.com/a/4682656/190597 (Joe Kington)
    http://softsurfer.com/Archive/algorithm_0101/algorithm_0101.htm#2D%20Polygons
    """
    area = 0.0
    for i in xrange(-1, len(x) - 1):
        area += x[i] * (y[i + 1] - y[i - 1])
    return area / 2.0

def centroid_of_polygon(points):
    """
    http://stackoverflow.com/a/14115494/190597 (mgamba)
    """
    area = area_of_polygon(*zip(*points))
    result_x = 0
    result_y = 0
    N = len(points)
    points = IT.cycle(points)
    x1, y1 = next(points)
    for i in range(N):
        x0, y0 = x1, y1
        x1, y1 = next(points)
        cross = (x0 * y1) - (x1 * y0)
        result_x += (x0 + x1) * cross
        result_y += (y0 + y1) * cross
    result_x /= (area * 6.0)
    result_y /= (area * 6.0)
    return (result_x, result_y)

def demo_centroid():
    points = [
        (30,50),
        (200,10),
        (250,50),
        (350,100),
        (200,180),
        (100,140),
        (10,200)
        ]
    cent = centroid_of_polygon(points)
    print(cent)
    # (159.2903828197946, 98.88888888888889)

demo_centroid()

Answer 3

也许这可以帮助您：finder这个存储库接收一个城市列表并将它们转换为多边形。