从Voronoi单元格获取有界多边形坐标

下一步该怎么做?只需根据边界框过滤点，边或面.然后根据众所周知的公式获取每个面孔的质心，以计算多边形的质心.这是结果的图像(质心为红色): 在显示代码之前一定很有趣太好了！看来行得通.如果经过一轮迭代后，我尝试在质心(红色)而不是初始点(蓝色)上重新运行算法，该怎么办?如果我一次又一次尝试怎么办? 第2步 步骤10 第25步 酷！ Voronoï细胞趋向于将其能量 ... 最小化这是代码import matplotlib.pyplot as plimport numpy as npimport scipy as spimport scipy.spatialimport syseps = sys.float_info.epsilonn_towers = 100towers = np.random.rand(n_towers, 2)bounding_box = np.array([0., 1., 0., 1.]) # [x_min, x_max, y_min, y_max]def in_box(towers, bounding_box): return np.logical_and(np.logical_and(bounding_box[0] <= towers[:, 0], towers[:, 0] <= bounding_box[1]), np.logical_and(bounding_box[2] <= towers[:, 1], towers[:, 1] <= bounding_box[3]))def voronoi(towers, bounding_box): # Select towers inside the bounding box i = in_box(towers, bounding_box) # Mirror points points_center = towers[i, :] points_left = np.copy(points_center) points_left[:, 0] = bounding_box[0] - (points_left[:, 0] - bounding_box[0]) points_right = np.copy(points_center) points_right[:, 0] = bounding_box[1] + (bounding_box[1] - points_right[:, 0]) points_down = np.copy(points_center) points_down[:, 1] = bounding_box[2] - (points_down[:, 1] - bounding_box[2]) points_up = np.copy(points_center) points_up[:, 1] = bounding_box[3] + (bounding_box[3] - points_up[:, 1]) points = np.append(points_center, np.append(np.append(points_left, points_right, axis=0), np.append(points_down, points_up, axis=0), axis=0), axis=0) # Compute Voronoi vor = sp.spatial.Voronoi(points) # Filter regions regions = [] for region in vor.regions: flag = True for index in region: if index == -1: flag = False break else: x = vor.vertices[index, 0] y = vor.vertices[index, 1] if not(bounding_box[0] - eps <= x and x <= bounding_box[1] + eps and bounding_box[2] - eps <= y and y <= bounding_box[3] + eps): flag = False break if region != [] and flag: regions.append(region) vor.filtered_points = points_center vor.filtered_regions = regions return vordef centroid_region(vertices): # Polygon's signed area A = 0 # Centroid's x C_x = 0 # Centroid's y C_y = 0 for i in range(0, len(vertices) - 1): s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1]) A = A + s C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s A = 0.5 * A C_x = (1.0 / (6.0 * A)) * C_x C_y = (1.0 / (6.0 * A)) * C_y return np.array([[C_x, C_y]])vor = voronoi(towers, bounding_box)fig = pl.figure()ax = fig.gca()# Plot initial pointsax.plot(vor.filtered_points[:, 0], vor.filtered_points[:, 1], 'b.')# Plot ridges pointsfor region in vor.filtered_regions: vertices = vor.vertices[region, :] ax.plot(vertices[:, 0], vertices[:, 1], 'go')# Plot ridgesfor region in vor.filtered_regions: vertices = vor.vertices[region + [region[0]], :] ax.plot(vertices[:, 0], vertices[:, 1], 'k-')# Compute and plot centroidscentroids = []for region in vor.filtered_regions: vertices = vor.vertices[region + [region[0]], :] centroid = centroid_region(vertices) centroids.append(list(centroid[0, :])) ax.plot(centroid[:, 0], centroid[:, 1], 'r.')ax.set_xlim([-0.1, 1.1])ax.set_ylim([-0.1, 1.1])pl.savefig("bounded_voronoi.png")sp.spatial.voronoi_plot_2d(vor)pl.savefig("voronoi.png")I have points (e.g., lat, lon pairs of cell tower locations) and I need to get the polygon of the Voronoi cells they form.from scipy.spatial import Voronoitower = [[ 24.686 , 46.7081], [ 24.686 , 46.7081], [ 24.686 , 46.7081]]c = Voronoi(towers)Now, I need to get the polygon boundaries in lat,lon coordinates for each cell (and what was the centroid this polygon surrounds). I need this Voronoi to be bounded as well. Meaning that the boundaries don't go to infinity, but rather within a bounding box. 解决方案 Given a rectangular bounding box, my first idea was to define a kind of intersection operation between this bounding box and the Voronoï diagram produce by scipy.spatial.Voronoi. An idea not necessarily great, since this requires to code a large number of basic functions of computational geometry.However, here is the second idea (hack?) that came to my mind: the algorithms to compute the Voronoï diagram of a set of n points in the plane have a time complexity of O(n ln(n)). What about adding points to constraint the Voronoï cells of the initial points to lie in the bounding box?Solution for a bounded Voronoï diagramA picture is worth a great speech:What I did here? That's pretty simple! The initial points (in blue) lie in [0.0, 1.0] x [0.0, 1.0]. Then I get the points (in blue) on the left (i.e. [-1.0, 0.0] x [0.0, 1.0]) by a reflection symmetry according to x = 0.0 (left edge of the bounding box). With reflection symmetries according to x = 1.0, y = 0.0 and y = 1.0 (other edges of the bounding box), I get all the points (in blue) I need to do the job.Then I run scipy.spatial.Voronoi. The previous image depicts the resulting Voronoï diagram (I use scipy.spatial.voronoi_plot_2d).What to do next? Just filter points, edges or faces according to the bounding box. And get the centroid of each face according to the well-known formula to compute centroid of polygon. Here is an image of the result (centroids are in red):Some fun before showing codeGreat! It seems to work. What if after one iteration I try to re-run the algorithm on the centroids (in red) rather than the initial points (in blue)? What if I try it again and again?Step 2Step 10Step 25Cool! Voronoï cells tend to minimize their energy...Here is the codeimport matplotlib.pyplot as plimport numpy as npimport scipy as spimport scipy.spatialimport syseps = sys.float_info.epsilonn_towers = 100towers = np.random.rand(n_towers, 2)bounding_box = np.array([0., 1., 0., 1.]) # [x_min, x_max, y_min, y_max]def in_box(towers, bounding_box): return np.logical_and(np.logical_and(bounding_box[0] <= towers[:, 0], towers[:, 0] <= bounding_box[1]), np.logical_and(bounding_box[2] <= towers[:, 1], towers[:, 1] <= bounding_box[3]))def voronoi(towers, bounding_box): # Select towers inside the bounding box i = in_box(towers, bounding_box) # Mirror points points_center = towers[i, :] points_left = np.copy(points_center) points_left[:, 0] = bounding_box[0] - (points_left[:, 0] - bounding_box[0]) points_right = np.copy(points_center) points_right[:, 0] = bounding_box[1] + (bounding_box[1] - points_right[:, 0]) points_down = np.copy(points_center) points_down[:, 1] = bounding_box[2] - (points_down[:, 1] - bounding_box[2]) points_up = np.copy(points_center) points_up[:, 1] = bounding_box[3] + (bounding_box[3] - points_up[:, 1]) points = np.append(points_center, np.append(np.append(points_left, points_right, axis=0), np.append(points_down, points_up, axis=0), axis=0), axis=0) # Compute Voronoi vor = sp.spatial.Voronoi(points) # Filter regions regions = [] for region in vor.regions: flag = True for index in region: if index == -1: flag = False break else: x = vor.vertices[index, 0] y = vor.vertices[index, 1] if not(bounding_box[0] - eps <= x and x <= bounding_box[1] + eps and bounding_box[2] - eps <= y and y <= bounding_box[3] + eps): flag = False break if region != [] and flag: regions.append(region) vor.filtered_points = points_center vor.filtered_regions = regions return vordef centroid_region(vertices): # Polygon's signed area A = 0 # Centroid's x C_x = 0 # Centroid's y C_y = 0 for i in range(0, len(vertices) - 1): s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1]) A = A + s C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s A = 0.5 * A C_x = (1.0 / (6.0 * A)) * C_x C_y = (1.0 / (6.0 * A)) * C_y return np.array([[C_x, C_y]])vor = voronoi(towers, bounding_box)fig = pl.figure()ax = fig.gca()# Plot initial pointsax.plot(vor.filtered_points[:, 0], vor.filtered_points[:, 1], 'b.')# Plot ridges pointsfor region in vor.filtered_regions: vertices = vor.vertices[region, :] ax.plot(vertices[:, 0], vertices[:, 1], 'go')# Plot ridgesfor region in vor.filtered_regions: vertices = vor.vertices[region + [region[0]], :] ax.plot(vertices[:, 0], vertices[:, 1], 'k-')# Compute and plot centroidscentroids = []for region in vor.filtered_regions: vertices = vor.vertices[region + [region[0]], :] centroid = centroid_region(vertices) centroids.append(list(centroid[0, :])) ax.plot(centroid[:, 0], centroid[:, 1], 'r.')ax.set_xlim([-0.1, 1.1])ax.set_ylim([-0.1, 1.1])pl.savefig("bounded_voronoi.png")sp.spatial.voronoi_plot_2d(vor)pl.savefig("voronoi.png") 这篇关于从Voronoi单元格获取有界多边形坐标的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！