ボロノイテッセレーションとそれに続くスムージングを実装しています。スムージングのためにロイドリラクゼーションを行うつもりでしたが、問題が発生しました。
ボロノイ辺の計算に次のモジュールを使用しています。
https://bitbucket.org/mozman/geoalg/src/5bbd46fa2270/geoalg/voronoi.py
スムージングのために、中心を計算できるように各ポリゴンのエッジを知る必要がありますが、残念ながらこのコードでは提供されていません。
私がアクセスできる情報は次のとおりです。
誰かがこれを計算するための比較的簡単な方法を見ることができますか？
重心を見つけるには、ウィキペディアで説明されている式を使用できます。
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です。
ボロノイ図のノードとエッジのコレクションに基づいて、個々のボロノイセルの重心を見つけるための単純なコードを次に示します。ノードに属するエッジを見つける方法を紹介し、前の図心と凸包のコードに依存します。
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'])
以下は、スクリプトの結果の画像です。青い線はエッジです。黒い四角はノードです。赤い四角は青い線が由来する頂点です。頂点とノードは任意に選択されました。赤い十字は図心です。実際のボロノイ分割ではありませんが、重心を取得するために使用される方法は、凸状のセルで構成される分割にも当てはまるはずです。
画像をレンダリングするための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ポイントを使用して、「最後の1ポイント」をより自然な方法で最初のポイントとして扱うことができます。
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 このリポジトリは都市のリストを受け取り、それらをポリゴンに変換します。