Finding centre of a polygon using limited data

Question

I'm implementing Voronoi tesselation followed by smoothing. For the smoothing I was going to do Lloyd relaxation, but I've encountered a problem.

I'm using following module for calculation of Voronoi sides:

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

For the smoothing I need to know the edges of each polygon so I can calculate the centre, which unfortunately this code doesn't provide.

Information I have access to consists of:

A list of all nodes,
A list of all edges (but just where they are, not what nodes they're associated with).

Can anyone see a relatively simple way to calculate this?

mgamba · Accepted Answer · 2013-01-07T07:38:42.480

For finding a centroid, you can use 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

I used a gift-wrapping algorithm to find the outside points (a.k.a. convex hull). There are a bunch of ways to do this, but gift-wrapping is nice because of its conceptual and practical simplicity. Here's an animated gif explaining this particular implementation:

step-by-step animated gif for counter-clockwise gift-wrapping, starting at the bottommost node

Here's some naive code to find centroids of the individual voronoi cells based on a collection of nodes and edges for a voronoi diagram. It introduces a method to find edges belonging to a node and relies on the previous centroid and convex-hull code:

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'])

Below is an image of the script results. The blue lines are edges. The black squares are nodes. The red squares are vertices that the blue lines are derived from. The vertices and nodes were chosen arbitrarily. The red crosses are centroids. While not an actual voronoi tesselation, the method used to procure the centroids should hold for tessalations composed of convex cells:

triangulated point cloud with calculated centroids and arbitrarily-chosen approximate centers

Here's the html to render the image:

<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>

This will likely get the job done. A more robust algo for finding which edges belong to a cell would be to use an inverse gift-wrapping method where edges are linked end-to-end and path choice at a split would be determined by angle. That method would not have a susceptibility to concave polygons and it would have the added benefit of not relying on the nodes.

That looks pretty impressive, but wouldn't I need to know which vertices are specifically linked to which node? — djcmm476, Jan 02 '13 at 19:32
Sorry for the misunderstanding. The script now calculates a convex hull for the polygon, and then finds the centroid of the hull. — mgamba, Jan 03 '13 at 09:09
Hmm, looks good. But would that not only work if you only had the vertices for the one polygon in play at once (as opposed to every vertex for every polygon)? — djcmm476, Jan 03 '13 at 17:43
Ah, I see. Ok, I'll post some pseudo code now and update it later. — mgamba, Jan 04 '13 at 22:55

unutbu · Answer 2 · 2013-08-06T20:49:59.797

4

This is @mgamba's answer, rewritten in a bit more of a python style. In particular, it uses itertools.cycle on the points so that "one-plus-the-last-point" can be treated as the first point in a more natural way.

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()

edited Aug 06 '13 at 20:49

answered Jan 02 '13 at 20:50

unutbu

842,883
184
1,785
1,677

I may be reading this code incorrectly, but doesn't that need the list of points for each cell? My problem is that I don't have a list of which edges comprise which cell. – djcmm476 Jan 02 '13 at 21:03
My code, and mgamba's are equivalent, except that we use different data structures for how a polygon is defined. What mgamba called `polygon` I called `points`. I chose to use a list of 2-tuples, since this structure is more generally useful. When doing mathematical calculations, there is no advantage to using a dict to represent a coordinate. If you do have polygons defined as lists of dicts, then it is easy to convert from `polygon` to `points` like this: `points = [(d['x'], d['y']) for d in polygon]`. – unutbu Jan 02 '13 at 21:16
I'm sorry if I am not answering your question directly; I do not know much about voronoi tesselation. I understand (perhaps mistakenly?) that you have a list of points which define a polygon, and wish to find its centroid. If that is not the case, can you give a concrete example of what data you have and what you expect as the answer? – unutbu Jan 02 '13 at 21:21
I didn't mean that, sorry. And again, I may just be reading the code wrong, but my code has a list of all the edges that make up the voronoi cells, it doesn't, however, know which edges form which cells, so you can't figure out the centre that way. I think yours might assume I know which edges go with which cell. – djcmm476 Jan 02 '13 at 21:21
(In response to your second post) I currently have a lot of polygons, and the data I have is all the vertices, and which pairs of the vertices join together to make a line. I don't have any information though about which vertices are used to make the edges of which polygons (sorry for the confusion). – djcmm476 Jan 02 '13 at 21:23
Are you looking for a way to compute a [Delaunay triagulation](http://en.wikipedia.org/wiki/Delaunay_triangulation#Relationship_with_the_Voronoi_diagram) from a Voronoi diagram? In other words, are the vertices of the Delaunay triangulation the centroid of the Voronoi cells? – unutbu Jan 02 '13 at 21:45
I have the delauney triangle points actually, but they only give you the point used to create the voronoi polygon, which isn't usually the centre point. I was going to use the centre points to try and normalise the polygons a little. – djcmm476 Jan 02 '13 at 21:49
So if I understand correctly, our new task is to find the nodes which form the outward-facing edges, and then apply the centroid formula? – mgamba Jan 02 '13 at 23:40
Yes, I can probably handle the centroid calculation, it's just finding the nodes that form the edges using the data I have. – djcmm476 Jan 03 '13 at 00:23
if the area is zero there will be a compilation error in case that the area is 0 – Dejell Dec 15 '13 at 12:00
i didn't understand what IT.cycle does – Dejell Dec 15 '13 at 13:01
I'd like to use your functions in an open-source project but cannot do that because by default all SO code is under [CC BY-SA license](https://meta.stackexchange.com/a/12539/155787). Can you explicitly state the license of your code, preferably something BSD-like? – letmaik Oct 23 '14 at 15:23
@neo: You are welcome to use the `centroid_of_polygon` code under the BSD license. Note that the `area_of_polygon` is [Joe Kington's](http://stackoverflow.com/a/4682656/190597). – unutbu Oct 23 '14 at 16:15

score 0 · Answer 3 · answered Jun 02 '15 at 13:31

0

Maybe this can help you: https://github.com/Bennyelg/geo_polygon_finder This repository receive a list of cities and translate them into polygons.

answered Jun 02 '15 at 13:31

Benjaminel

1
1

Finding centre of a polygon using limited data

3 Answers3

Linked