限制由Qhull生成的Voronoi顶点的定义域

Question

限制由Qhull生成的Voronoi顶点的定义域

numpyscipygeometrycomputational-geometryqhull

4

我的问题简而言之：是否可以限制Qhull生成的Voronoi顶点的定义域？如果可以，如何做到？

我的问题和背景：我正在处理一个数据可视化项目，在2D场中有一些点。这些点有些重叠在一起，所以我会稍微“抖动”它们，以使它们不重叠。

我目前解决这个问题的方法是使用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 的出色回复后，我撰写了一篇博客文章，详细介绍了Lloyd迭代的有界和无界形式。我还制作了一个小软件包lloyd，使得在数据集上运行有界Lloyd迭代更加容易。我希望分享这些资源，以帮助其他人进行相关分析。

- duhaime

你认为什么是“原始数据域”？凸包吗？ - Ripi2

@Ripi2 好的，那太棒了。我想我一开始考虑的是将原始数据域仅限于 min_x、max_x、min_y 和 max_y，但如果可以通过凸壳加以限制，那就更好了。 - duhaime

如果您已经有了边界，则只需防止创建/移动点超出该边界即可。 - Ripi2

@Ripi2 是的，那是我的计划，但我认为将工作推给Qhull会更简单/更清洁/更快（即要求Qhull不要在指定域之外创建Voronoi顶点）。 - duhaime

2个回答

1

很遗憾，@tfinniga的建议都没有给我满意的结果。

在我的实现中，边界框角上的人工边界点似乎不能限制布局。只有当边界点密集地放置在边界框的边缘上时，边界点才能起作用，这会大大减慢沃罗诺伊图的计算速度。

对于仅在一个维度上超出边界框的点，对沃罗诺伊顶点或计算出的沃罗诺伊质心进行简单剪切可以起到作用。否则，多个点可能被分配到边界框的同一个角落，这将导致未定义行为，因为沃罗诺伊图仅针对一组唯一位置定义。

相反，在我的下面的实现中，我只是不更新新位置在边界框外的点。需要约束布局的额外计算量微不足道，并且--根据我的测试所示--这种方法不会遇到任何破坏性的边缘情况。

#!/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()

- Paul Brodersen

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- tfinniga · Accepted Answer

你遇到的核心问题是Lloyd算法没有边界限制，因此点会爆炸。解决方法有两种：

在计算重心之前，将你得到的Voronoi图手动剪裁到你的边界矩形。这将给你正确的解决方案-你可以将其与维基百科文章中链接的示例进行测试，并确保它匹配。
添加人工边界点。例如，您可以在边界框的4个角落或每条边上添加一些不移动的点。这将阻止内部点爆炸。这不会给你Lloyd算法的“正确”结果，但可能会为你的可视化输出提供有用的结果，并且更容易实现。