约束 Qhull 生成的 Voronoi 顶点的域
Constraining the Domain of Voronoi Vertices Generated by Qhull
我的问题简而言之：是否可以约束 Qhull 生成的 Voronoi 顶点的域？如果可以，怎么办？
我的上下文问题：我正在研究数据可视化，其中我在 2D 字段中有点。这些点有点重叠，所以我稍微"jittering"它们，以免它们重叠。
我目前执行此任务的方法是使用 Lloyd's algorithm 移动点。 Lloyd 算法本质上采用初始点位置，计算 Voronoi 映射，并在算法的每次迭代期间将每个点移动到其 Voronoi 区域的中心。
这是 Python 中的示例：
from scipy.spatial import Voronoi
from scipy.spatial import voronoi_plot_2d
import matplotlib.pyplot as plt
import numpy as np
import sys
class Field():
'''
Create a Voronoi map that can be used to run Lloyd relaxation on an array of 2D points.
For background, see: https://en.wikipedia.org/wiki/Lloyd%27s_algorithm
'''
def __init__(self, arr):
'''
Store the points and bounding box of the points to which Lloyd relaxation will be applied
@param [arr] arr: a numpy array with shape n, 2, where n is number of points
'''
if not len(arr):
raise Exception('please provide a numpy array with shape n,2')
x = arr[:, 0]
y = arr[:, 0]
self.bounding_box = [min(x), max(x), min(y), max(y)]
self.points = arr
self.build_voronoi()
def build_voronoi(self):
'''
Build a Voronoi map from self.points. For background on self.voronoi attrs, see:
https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.spatial.Voronoi.html
'''
eps = sys.float_info.epsilon
self.voronoi = Voronoi(self.points)
self.filtered_regions = [] # list of regions with vertices inside Voronoi map
for region in self.voronoi.regions:
inside_map = True # is this region inside the Voronoi map?
for index in region: # index = the idx of a vertex in the current region
# check if index is inside Voronoi map (indices == -1 are outside map)
if index == -1:
inside_map = False
break
# check if the current coordinate is in the Voronoi map's bounding box
else:
coords = self.voronoi.vertices[index]
if not (self.bounding_box[0] - eps <= coords[0] and
self.bounding_box[1] + eps >= coords[0] and
self.bounding_box[2] - eps <= coords[1] and
self.bounding_box[3] + eps >= coords[1]):
inside_map = False
break
# store hte region if it has vertices and is inside Voronoi map
if region != [] and inside_map:
self.filtered_regions.append(region)
def find_centroid(self, vertices):
'''
Find the centroid of a Voroni region described by `vertices`, and return a
np array with the x and y coords of that centroid.
The equation for the method used here to find the centroid of a 2D polygon
is given here: https://en.wikipedia.org/wiki/Centroid#Of_a_polygon
@params: np.array `vertices` a numpy array with shape n,2
@returns np.array a numpy array that defines the x, y coords
of the centroid described by `vertices`
'''
area = 0
centroid_x = 0
centroid_y = 0
for i in range(len(vertices)-1):
step = (vertices[i, 0] * vertices[i+1, 1]) - (vertices[i+1, 0] * vertices[i, 1])
area += step
centroid_x += (vertices[i, 0] + vertices[i+1, 0]) * step
centroid_y += (vertices[i, 1] + vertices[i+1, 1]) * step
area /= 2
centroid_x = (1.0/(6.0*area)) * centroid_x
centroid_y = (1.0/(6.0*area)) * centroid_y
return np.array([centroid_x, centroid_y])
def relax(self):
'''
Moves each point to the centroid of its cell in the Voronoi map to "relax"
the points (i.e. jitter them so as to spread them out within the space).
'''
centroids = []
for region in self.filtered_regions:
vertices = self.voronoi.vertices[region + [region[0]], :]
centroid = self.find_centroid(vertices) # get the centroid of these verts
centroids.append(list(centroid))
self.points = centroids # store the centroids as the new point positions
self.build_voronoi() # rebuild the voronoi map given new point positions
##
# Visualize
##
# built a Voronoi diagram that we can use to run lloyd relaxation
field = Field(points)
# plot the points with no relaxation relaxation
plt = voronoi_plot_2d(field.voronoi, show_vertices=False, line_colors='orange', line_alpha=0.6, point_size=2)
plt.show()
# relax the points several times, and show the result of each relaxation
for i in range(6):
field.relax() # the .relax() method performs lloyd relaxation
plt = voronoi_plot_2d(field.voronoi, show_vertices=False, line_colors='orange', line_alpha=0.6, point_size=2)
plt.show()
正如我们所见，在每次迭代（第 2 帧和第 3 帧）期间，原始数据集中的点（第 1 帧，顶部）重叠越来越少，这太棒了！
这种方法的问题是我目前正在从绘图中删除那些 voronoi 区域边界超出初始数据集域的点。 （如果我不这样做，最外面的点会迅速射入超空间并远离其余点。）这最终意味着我最终会丢弃点，这是不好的。
我想我可以通过约束 Qhull Voronoi 域来解决这个问题，以便只在原始数据域内创建 Voronoi 顶点。
是否可以通过这种方式来约束Qhull？其他人可以提供的任何帮助将不胜感激！
更新
在收到@tfinniga 的出色回复后，我整理了一个 blog post detailing Lloyd iteration in its bounded and unbounded forms. I also put together a little package lloyd 以便更容易地 运行 对数据集进行有界劳埃德迭代。我想分享这些资源，以防他们帮助其他人进行相关分析。
您 运行 遇到的核心问题是没有任何限制 Lloyd 算法的东西，所以点数激增。解决此问题的两种方法：
获取您得到的 voronoi 图，并在计算质心之前将其手动剪辑到您的边界矩形。这将为您提供正确的解决方案 - 您可以根据您在维基百科文章中链接的示例对其进行测试，并确保它匹配。
添加点的人工边界。例如，您可以在边界框的 4 个角处添加点，或者在您不移动的每边添加一些点。这些将阻止内部点爆炸。这不会为您提供 Lloyd 算法的 'correct' 结果，但可能会为您提供对可视化有用的输出，并且更容易实现。
不幸的是，@tfinniga 的两个建议都没有给我满意的结果。
在我的手中，边界框角落的人工边界点似乎不会限制布局。边界点似乎只有在沿着边界框的边缘密集放置时才有效，这会大大减慢 Voronoi 曲面细分的计算速度。
将 Voronoi 顶点或计算出的 Voronoi 质心简单剪裁到边界框适用于仅在一维上超出边界框的点。否则，多个点最终可能会被分配到边界框的同一个角，这会导致未定义的行为，因为 Voronoi 细分仅针对一组唯一位置定义。
相反，在我下面的实现中，我根本不更新点，新位置在边界框之外。约束布局所需的额外计算非常小，而且——据我的测试所知——这种方法不会 运行 进入任何突破性的边缘情况。
#!/usr/bin/env python
"""
Implementation of a constrained Lloyds algorithm to remove node overlaps.
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi
def lloyds(positions, origin=(0,0), scale=(1,1), total_iterations=3):
positions = positions.copy()
for _ in range(total_iterations):
centroids = _get_voronoi_centroids(positions)
is_valid = _is_within_bbox(centroids, origin, scale)
positions[is_valid] = centroids[is_valid]
return positions
def _get_voronoi_centroids(positions):
voronoi = Voronoi(positions)
centroids = np.zeros_like(positions)
for ii, idx in enumerate(voronoi.point_region):
# ignore voronoi vertices at infinity
# TODO: compute regions clipped by bbox
region = [jj for jj in voronoi.regions[idx] if jj != -1]
centroids[ii] = _get_centroid(voronoi.vertices[region])
return centroids
def _get_centroid(polygon):
# TODO: formula may be incorrect; correct one here:
# https://en.wikipedia.org/wiki/Centroid#Of_a_polygon
return np.mean(polygon, axis=0)
def _is_within_bbox(points, origin, scale):
minima = np.array(origin)
maxima = minima + np.array(scale)
return np.all(np.logical_and(points >= minima, points <= maxima), axis=1)
if __name__ == '__main__':
positions = np.random.rand(200, 2)
adjusted = lloyds(positions)
fig, axes = plt.subplots(1, 2, sharex=True, sharey=True)
axes[0].plot(*positions.T, 'ko')
axes[1].plot(*adjusted.T, 'ro')
for ax in axes:
ax.set_aspect('equal')
fig.tight_layout()
plt.show()
我的问题简而言之：是否可以约束 Qhull 生成的 Voronoi 顶点的域？如果可以，怎么办？
我的上下文问题：我正在研究数据可视化，其中我在 2D 字段中有点。这些点有点重叠，所以我稍微"jittering"它们，以免它们重叠。
我目前执行此任务的方法是使用 Lloyd's algorithm 移动点。 Lloyd 算法本质上采用初始点位置，计算 Voronoi 映射，并在算法的每次迭代期间将每个点移动到其 Voronoi 区域的中心。
这是 Python 中的示例：
from scipy.spatial import Voronoi
from scipy.spatial import voronoi_plot_2d
import matplotlib.pyplot as plt
import numpy as np
import sys
class Field():
'''
Create a Voronoi map that can be used to run Lloyd relaxation on an array of 2D points.
For background, see: https://en.wikipedia.org/wiki/Lloyd%27s_algorithm
'''
def __init__(self, arr):
'''
Store the points and bounding box of the points to which Lloyd relaxation will be applied
@param [arr] arr: a numpy array with shape n, 2, where n is number of points
'''
if not len(arr):
raise Exception('please provide a numpy array with shape n,2')
x = arr[:, 0]
y = arr[:, 0]
self.bounding_box = [min(x), max(x), min(y), max(y)]
self.points = arr
self.build_voronoi()
def build_voronoi(self):
'''
Build a Voronoi map from self.points. For background on self.voronoi attrs, see:
https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.spatial.Voronoi.html
'''
eps = sys.float_info.epsilon
self.voronoi = Voronoi(self.points)
self.filtered_regions = [] # list of regions with vertices inside Voronoi map
for region in self.voronoi.regions:
inside_map = True # is this region inside the Voronoi map?
for index in region: # index = the idx of a vertex in the current region
# check if index is inside Voronoi map (indices == -1 are outside map)
if index == -1:
inside_map = False
break
# check if the current coordinate is in the Voronoi map's bounding box
else:
coords = self.voronoi.vertices[index]
if not (self.bounding_box[0] - eps <= coords[0] and
self.bounding_box[1] + eps >= coords[0] and
self.bounding_box[2] - eps <= coords[1] and
self.bounding_box[3] + eps >= coords[1]):
inside_map = False
break
# store hte region if it has vertices and is inside Voronoi map
if region != [] and inside_map:
self.filtered_regions.append(region)
def find_centroid(self, vertices):
'''
Find the centroid of a Voroni region described by `vertices`, and return a
np array with the x and y coords of that centroid.
The equation for the method used here to find the centroid of a 2D polygon
is given here: https://en.wikipedia.org/wiki/Centroid#Of_a_polygon
@params: np.array `vertices` a numpy array with shape n,2
@returns np.array a numpy array that defines the x, y coords
of the centroid described by `vertices`
'''
area = 0
centroid_x = 0
centroid_y = 0
for i in range(len(vertices)-1):
step = (vertices[i, 0] * vertices[i+1, 1]) - (vertices[i+1, 0] * vertices[i, 1])
area += step
centroid_x += (vertices[i, 0] + vertices[i+1, 0]) * step
centroid_y += (vertices[i, 1] + vertices[i+1, 1]) * step
area /= 2
centroid_x = (1.0/(6.0*area)) * centroid_x
centroid_y = (1.0/(6.0*area)) * centroid_y
return np.array([centroid_x, centroid_y])
def relax(self):
'''
Moves each point to the centroid of its cell in the Voronoi map to "relax"
the points (i.e. jitter them so as to spread them out within the space).
'''
centroids = []
for region in self.filtered_regions:
vertices = self.voronoi.vertices[region + [region[0]], :]
centroid = self.find_centroid(vertices) # get the centroid of these verts
centroids.append(list(centroid))
self.points = centroids # store the centroids as the new point positions
self.build_voronoi() # rebuild the voronoi map given new point positions
##
# Visualize
##
# built a Voronoi diagram that we can use to run lloyd relaxation
field = Field(points)
# plot the points with no relaxation relaxation
plt = voronoi_plot_2d(field.voronoi, show_vertices=False, line_colors='orange', line_alpha=0.6, point_size=2)
plt.show()
# relax the points several times, and show the result of each relaxation
for i in range(6):
field.relax() # the .relax() method performs lloyd relaxation
plt = voronoi_plot_2d(field.voronoi, show_vertices=False, line_colors='orange', line_alpha=0.6, point_size=2)
plt.show()
正如我们所见，在每次迭代（第 2 帧和第 3 帧）期间，原始数据集中的点（第 1 帧，顶部）重叠越来越少，这太棒了！
这种方法的问题是我目前正在从绘图中删除那些 voronoi 区域边界超出初始数据集域的点。 （如果我不这样做，最外面的点会迅速射入超空间并远离其余点。）这最终意味着我最终会丢弃点，这是不好的。
我想我可以通过约束 Qhull Voronoi 域来解决这个问题，以便只在原始数据域内创建 Voronoi 顶点。
是否可以通过这种方式来约束Qhull？其他人可以提供的任何帮助将不胜感激！
更新
在收到@tfinniga 的出色回复后，我整理了一个 blog post detailing Lloyd iteration in its bounded and unbounded forms. I also put together a little package lloyd 以便更容易地 运行 对数据集进行有界劳埃德迭代。我想分享这些资源，以防他们帮助其他人进行相关分析。
您 运行 遇到的核心问题是没有任何限制 Lloyd 算法的东西，所以点数激增。解决此问题的两种方法：
获取您得到的 voronoi 图，并在计算质心之前将其手动剪辑到您的边界矩形。这将为您提供正确的解决方案 - 您可以根据您在维基百科文章中链接的示例对其进行测试，并确保它匹配。
添加点的人工边界。例如，您可以在边界框的 4 个角处添加点，或者在您不移动的每边添加一些点。这些将阻止内部点爆炸。这不会为您提供 Lloyd 算法的 'correct' 结果，但可能会为您提供对可视化有用的输出，并且更容易实现。
不幸的是，@tfinniga 的两个建议都没有给我满意的结果。
在我的手中，边界框角落的人工边界点似乎不会限制布局。边界点似乎只有在沿着边界框的边缘密集放置时才有效，这会大大减慢 Voronoi 曲面细分的计算速度。
将 Voronoi 顶点或计算出的 Voronoi 质心简单剪裁到边界框适用于仅在一维上超出边界框的点。否则，多个点最终可能会被分配到边界框的同一个角，这会导致未定义的行为，因为 Voronoi 细分仅针对一组唯一位置定义。
相反，在我下面的实现中，我根本不更新点，新位置在边界框之外。约束布局所需的额外计算非常小，而且——据我的测试所知——这种方法不会 运行 进入任何突破性的边缘情况。
#!/usr/bin/env python
"""
Implementation of a constrained Lloyds algorithm to remove node overlaps.
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi
def lloyds(positions, origin=(0,0), scale=(1,1), total_iterations=3):
positions = positions.copy()
for _ in range(total_iterations):
centroids = _get_voronoi_centroids(positions)
is_valid = _is_within_bbox(centroids, origin, scale)
positions[is_valid] = centroids[is_valid]
return positions
def _get_voronoi_centroids(positions):
voronoi = Voronoi(positions)
centroids = np.zeros_like(positions)
for ii, idx in enumerate(voronoi.point_region):
# ignore voronoi vertices at infinity
# TODO: compute regions clipped by bbox
region = [jj for jj in voronoi.regions[idx] if jj != -1]
centroids[ii] = _get_centroid(voronoi.vertices[region])
return centroids
def _get_centroid(polygon):
# TODO: formula may be incorrect; correct one here:
# https://en.wikipedia.org/wiki/Centroid#Of_a_polygon
return np.mean(polygon, axis=0)
def _is_within_bbox(points, origin, scale):
minima = np.array(origin)
maxima = minima + np.array(scale)
return np.all(np.logical_and(points >= minima, points <= maxima), axis=1)
if __name__ == '__main__':
positions = np.random.rand(200, 2)
adjusted = lloyds(positions)
fig, axes = plt.subplots(1, 2, sharex=True, sharey=True)
axes[0].plot(*positions.T, 'ko')
axes[1].plot(*adjusted.T, 'ro')
for ax in axes:
ax.set_aspect('equal')
fig.tight_layout()
plt.show()