要找到一个质心,您可以使用维基百科上描述的公式:
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
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):
break
else:
convex_hull.append(polygon[hull_index_candidate])
previous_angle = smallest_angle
previous_hull_index_candidate = hull_index_candidate
return convex_hull
我使用了一种礼物包装算法来找到外部点(也称为凸包)。有很多方法可以做到这一点,但是礼物包装算法因其概念上和实际上的简单性而受欢迎。下面是一个动画GIF,解释了这个特定实现:
![step-by-step animated gif for counter-clockwise gift-wrapping, starting at the bottommost node](https://istack.dev59.com/cpadV.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):
return (C[1]-A[1])*(B[0]-A[0]) > (B[1]-A[1])*(C[0]-A[0])
def intersect(segment1, segment2):
A = segment1[0]
B = segment1[1]
C = segment2[0]
D = segment2[1]
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 i in range(0,len(points)):
cell_edges = []
for j in range(0,len(edges)):
is_cell_edge = True
vector = (points[i],midpoint(edges[j]))
for k in range(0,len(edges)):
if (k != j and intersect(edges[k], vector)):
is_cell_edge = False
break
if (is_cell_edge):
cell_edges.append(edges[j])
convex_hull = convex_hull_for_polygon(points_from_edges(cell_edges))
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分割,但用于获取质心的方法应适用于由凸多边形组成的分割:
![triangulated point cloud with calculated centroids and arbitrarily-chosen approximate centers](https://istack.dev59.com/uP7f9.webp)
以下是呈现图像的HTML代码:
这里是呈现图像的HTML代码:
<html>
<head>
<script>
window.onload = draw;
function draw() {
var canvas = document.getElementById('canvas').getContext('2d');
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);
}
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();
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);
}
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>
这可能能够完成工作。更强大的算法用于查找哪些边属于一个单元格将是使用反向礼品包装方法,其中边缘端对端链接,并且在拆分处的路径选择将由角度确定。该方法不会受到凹多边形的影响,并且具有不依赖节点的额外好处。