我正在实施 Voronoi 细分，然后进行平滑处理。对于平滑，我打算做劳埃德放松，但我遇到了一个问题。
我正在使用以下模块来计算 Voronoi 边：
https://bitbucket.org/mozman/geoalg/src/5bbd46fa2270/geoalg/voronoi.py
对于平滑，我需要知道每个多边形的边缘，以便我可以计算中心，不幸的是这段代码没有提供。
我可以访问的信息包括：
谁能看到一个相对简单的计算方法？
要找到质心，您可以使用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
我使用礼物包装算法来找到外部点（又名凸包）。有很多方法可以做到这一点，但礼品包装很好，因为它的概念和实用简单。这是一个解释这个特定实现的动画 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>
这可能会完成工作。用于查找哪些边缘属于单元格的更强大的算法是使用逆礼物包装方法，其中边缘端到端链接，并且分割处的路径选择将由角度确定。该方法对凹多边形不敏感，并且具有不依赖节点的额外好处。
这是@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()
也许这可以帮助你： https ://github.com/Bennyelg/geo_polygon_finder 这个存储库接收城市列表并将它们转换为多边形。