Ограничение домена вершин Вороного, сгенерированного Qhull

Вкратце мой вопрос: возможно ли ограничить домен вершин Вороного, сгенерированных Qhull? Если да, то как это сделать?

Мой вопрос с контекстом: я работаю над визуализацией данных, в которой у меня есть точки в 2D-поле. Точки немного перекрываются, поэтому я их слегка «дрожу», чтобы они не перекрывались.

Мой текущий подход к этой задаче заключается в использовании алгоритма Ллойда для перемещения точек. Алгоритм Ллойда по существу берет начальные положения точек, вычисляет карту Вороного и перемещает каждую точку в центр ее области Вороного во время каждой итерации алгоритма.

Вот пример на 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, вверху) перекрываются все меньше и меньше, и это здорово!

введите здесь описание изображения

Проблема с этим подходом заключается в том, что я в настоящее время удаляю с графика те точки, границы области вороного которых лежат за пределами домена исходного набора данных. (Если я этого не сделаю, самые дальние точки быстро уйдут в гиперпространство и уйдут очень далеко от остальных точек.) Это в конечном итоге означает, что я в конечном итоге отбрасываю точки, что нехорошо.

Я думал, что смогу решить эту проблему, ограничив область Qhull Voronoi так, чтобы создавать вершины Voronoi только в пределах исходной области данных.

Можно ли таким образом ограничить Кхалла? Любая помощь, которую могут предложить другие, будет принята с благодарностью!


Обновлять

Получив отличный ответ @tfinniga, я составил сообщение в блоге с подробным описанием итерации Ллойда в ее ограниченном и неограниченные формы. Я также собрал небольшой пакет lloyd, чтобы упростить запуск ограниченных итераций Ллойда в наборе данных. Я хотел поделиться этими ресурсами на случай, если они помогут другим в проведении соответствующего анализа.


person duhaime    schedule 25.09.2018    source источник
comment
Что вы считаете исходной предметной областью? Выпуклая оболочка?   -  person Ripi2    schedule 26.09.2018
comment
@Ripi2, да, это было бы здорово. Я полагаю, что предполагал исходную область данных просто как min_x, max_x, min_y, max_y, но если можно ограничиться корпусом, это даже лучше   -  person duhaime    schedule 26.09.2018
comment
Если у вас уже есть граница, просто не допускайте создания/перемещения точки за пределы этой границы.   -  person Ripi2    schedule 26.09.2018
comment
@Ripi2 Да, это мой план, но я подумал, что было бы проще/чище/быстрее передать работу Qhull (т.е. попросить Qhull не создавать вершины voronoi за пределами указанного домена)   -  person duhaime    schedule 26.09.2018


Ответы (2)


Основная проблема, с которой вы сталкиваетесь, заключается в том, что алгоритм Ллойда ничем не ограничивается, поэтому точки взрываются. Два способа исправить это:

  1. Возьмите полученную диаграмму Вороного и вручную прикрепите ее к ограничивающему прямоугольнику перед вычислением центроида. Это даст вам правильное решение — вы можете проверить его на примере, который вы указали в статье в Википедии, и убедиться, что он соответствует.

  2. Добавьте искусственную границу точек. Например, вы можете добавить точки в 4 углах вашей ограничивающей рамки или несколько вдоль каждой стороны, которые вы не перемещаете. Это предотвратит взрыв внутренних точек. Это не даст вам «правильного» результата алгоритма Ллойда, но может дать вам вывод, полезный для вашей визуализации, и его проще реализовать.

person tfinniga    schedule 26.09.2018
comment
абсолютно идеально. Отсечение вершин Вороного оставило много точек, сгруппированных по краям (похоже, площадь этих ячеек приближается к 0), но добавление ограничивающих точек сработало отлично. Вы случайно не знаете какие-нибудь хорошие книги, в которых обсуждаются геометрические алгоритмы, такие как итерация Ллойда? Я узнал об этом только из случайного твита от разработчика WebGL, и я хотел бы узнать больше о других интересных алгоритмах, связанных с графикой и геометрией. Любые рекомендации будут очень признательны! - person duhaime; 30.09.2018
comment
@duhaime К сожалению, я не знаю хорошего ресурса, на который можно было бы указать вам - я узнал об алгоритме Ллойда только из вашего вопроса. - person tfinniga; 10.10.2018
comment
Кстати, если вы довольны результатами, которые вы получаете, я бы не стал с этим связываться, но есть несколько способов исправить свернутые точки - вместо того, чтобы устанавливать каждую точку в центр тяжести, вы можете переместить каждую точку на некоторый процент путь к центроиду. В качестве альтернативы вы могли бы иметь более сложный термин для движения, который включает поле расстояния со знаком, так что точки отталкиваются от краев ограничивающего прямоугольника. - person tfinniga; 10.10.2018
comment
аминь спасибо за ваше замечание. Я думал о каком-то скаляре, который бы контролировал притяжение к центроиду, но несколько итераций перемещения точек к точному центроиду в данный момент работают для моей конкретной проблемы. Это ловкий трюк! - person duhaime; 10.10.2018

К сожалению, ни одно из предложений @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()
person Paul Brodersen    schedule 11.05.2021