djcmm476 - 1 year ago 85
Python Question

# Finding centre of a polygon using limited data

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.

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

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:

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

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:

Here's the html to render the image:

``````<html>
<script>
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>