Centroidal Voronoi tessellation

I am trying to build a bounded Voronoi diagram using the scipy package and in each iteration I compute the centroids of the Voronoi cells and move a bit say some delta towards the centroid and recompute the Voronoi diagram by updating the generator points. When I try to plot the updated points I get a weird error as in the point I plot is not where it is expected to be. Here's the code

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)

Also here is an image of what the code does. The blue points are the generator points, red are the centroids and yellow are supposed to be the midway points between the blue and red points. But as you can see the yellow points are not in between the blue and red points.

Solution

The problem is that your points when fed to Voronoi get inflated during the construction of the tessellation, and when you later filter them out the points are in the wrong order. Consequently when you set vor.filtered_points = points_center in voronoi(), the points are shuffled compared to the order of regions. So while you're computing the midpoints correctly, you're using the wrong pairs of points.

I circled two correct pairings in green and an incorrect one in red here: As you can see from the red circle, the basis point in an edge cell is paired with the centroid of an adjacent cell.

The solution is simple: when you're filtering the regions and find a region to keep, you need to gather the point which falls inside the corresponding region. You can do this by matching vor.points to vor.point_region and finding the corresponding region, for which you'll need to enumerate your 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

With these modifications the averaging works fine: