Центроидальная вороной тесселяция

Question

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

import matplotlib.pyplot as pl
import numpy as np
import scipy as sp
import scipy.spatial
import sys

np.random.seed(1)
eps = sys.float_info.epsilon

n_robots = 10
robots = np.random.rand(n_robots, 2)
#print(robots)
bounding_box = np.array([0., 1., 0., 1.]) 

def in_box(robots, bounding_box):
    return np.logical_and(np.logical_and(bounding_box[0] <= robots[:, 0],
                                         robots[:, 0] <= bounding_box[1]),
                          np.logical_and(bounding_box[2] <= robots[:, 1],
                                         robots[:, 1] <= bounding_box[3]))


def voronoi(robots, bounding_box):
    i = in_box(robots, bounding_box)
    points_center = robots[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 = []
    ind = np.arange(points.shape[0])
    ind = np.expand_dims(ind,axis= 1)


    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 vor

def centroid_region(vertices):

    A = 0

    C_x = 0

    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]])

def plot(r,index):
    vor = voronoi(r, bounding_box)

    fig = pl.figure()
    ax = fig.gca()
# Plot initial points
    ax.plot(vor.filtered_points[:, 0], vor.filtered_points[:, 1], 'b.')
    print("initial",vor.filtered_points)
# Plot ridges points
    for region in vor.filtered_regions:
        vertices = vor.vertices[region, :]
        ax.plot(vertices[:, 0], vertices[:, 1], 'go')
# Plot ridges
    for region in vor.filtered_regions:
        vertices = vor.vertices[region + [region[0]], :]
        ax.plot(vertices[:, 0], vertices[:, 1], 'k-')
# Compute and plot centroids
    centroids = []
    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.')
    centroids = np.asarray(centroids)
    rob = np.copy(vor.filtered_points)
    # the below code is for the plotting purpose the update happens in the update function
    interim_x = np.asarray(centroids[:,0] - rob[:,0])
    interim_y = np.asarray(centroids[:,1] - rob[:,1])
    magn = [np.linalg.norm(centroids[i,:] - rob[i,:]) for i in range(rob.shape[0])]
    x = np.copy(interim_x)
    x = np.asarray([interim_x[i]/magn[i] for i in range(interim_x.shape[0])])
    y = np.copy(interim_y)
    y = np.asarray([interim_y[i]/magn[i] for i in range(interim_y.shape[0])])
    nor = np.copy(rob)
    for i in range(x.shape[0]):
        nor[i,0] = x[i]
        nor[i,1] = y[i]
    temp = np.copy(rob)
    temp[:,0] = [rob[i,0] + 0.5*interim_x[i] for i in range(rob.shape[0])]
    temp[:,1] = [rob[i,1] + 0.5*interim_y[i] for i in range(rob.shape[0])]
    ax.plot(temp[:,0] ,temp[:,1], 'y.' )
    ax.set_xlim([-0.1, 1.1])
    ax.set_ylim([-0.1, 1.1])
    pl.savefig("voronoi" + str(index) + ".png")
    return centroids

def update(rob,centroids):

  interim_x = np.asarray(centroids[:,0] - rob[:,0])
  interim_y = np.asarray(centroids[:,1] - rob[:,1])
  magn = [np.linalg.norm(centroids[i,:] - rob[i,:]) for i in range(rob.shape[0])]
  x = np.copy(interim_x)
  x = np.asarray([interim_x[i]/magn[i] for i in range(interim_x.shape[0])])
  y = np.copy(interim_y)
  y = np.asarray([interim_y[i]/magn[i] for i in range(interim_y.shape[0])])
  nor = [np.linalg.norm([x[i],y[i]]) for i in range(x.shape[0])]
  temp = np.copy(rob)
  temp[:,0] = [rob[i,0] + 0.5*interim_x[i] for i in range(rob.shape[0])]
  temp[:,1] = [rob[i,1] + 0.5*interim_y[i] for i in range(rob.shape[0])]
  return np.asarray(temp)

if __name__ == '__main__':
    for i in range(1):
        centroids = plot(robots,i)
        robots = update(robots,centroids)

Также вот изображение того, что делает код. Синие точки - это точки генератора, красные - центроиды, а желтые - промежуточные точки между синей и красной точками. Но, как вы можете видеть, желтые точки не находятся между синей и красной точками.

Answer 1

Проблема в том, что ваш points когда кормят Voronoi получить инфляцию во время построения тесселяции, и когда вы позже отфильтруете их, точки будут в неправильном порядке. Следовательно, когда вы установите vor.filtered_points = points_center в voronoi(), точки перемешиваются по сравнению с порядком регионов. Поэтому, пока вы правильно вычисляете средние точки, вы используете неправильные пары точек.

Я обвел две правильные пары зеленым и неправильную красным: Как видно из красного круга, базовая точка в граничной ячейке соединяется с центроидом соседней ячейки.

Решение простое: когда вы фильтруете регионы и находите регион для хранения, вам нужно собрать точку, которая попадает в соответствующий регион. Вы можете сделать это, сопоставив vor.points в vor.point_region и найти соответствующий регион, для которого вам нужно enumerate ваш regions:

# Compute Voronoi
vor = sp.spatial.Voronoi(points)
# Filter regions and select corresponding points
regions = []
points_to_filter = [] # we'll need to gather points too
ind = np.arange(points.shape[0])
ind = np.expand_dims(ind,axis= 1)

for i,region in enumerate(vor.regions): # enumerate the regions
    if not region: # nicer to skip the empty region altogether
        continue

    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 flag:
        regions.append(region)

        # find the point which lies inside
        points_to_filter.append(vor.points[vor.point_region == i][0,:])

vor.filtered_points = np.array(points_to_filter)
vor.filtered_regions = regions

С этими модификациями усреднение прекрасно работает: