voronoi and lloyd relaxation using python/scipy

How to determine, using Qhull, which voronoi cells (by index) are "proper" (made of "existing vertices")

I am trying to perform a constrained relaxation using LLoyds algorithm and input generated by scipy.spatial Voronoi (which is a wrapper for Qhull).

In terms of code it looks like:

points = [n  for n in itertools.product(xrange(3),xrange(3))]
vor = Voronoi(points)

vor2 = lloyd(vor) # my relaxation function  - not relevant to the question

The output graph generated by the code looks ok (see below) but the data in the vor structure is not enough to perform the Lloyds relaxation. It is because I should only move the points that are inside the valid voronoi cells ( #4 in the image). The other should be left as they are. Qhull messes the order of points/regions so I can not estimate which region belongs to which point.

Here is the illustration of the problem:

print vor.vertices
#[[ 0.5  0.5]
# [ 1.5  0.5]
# [ 0.5  1.5]
# [ 1.5  1.5]]

print vor.regions  
# [[], [-1, 0], [-1, 1], [1, -1, 0], [3, -1, 2], [-1, 3], [-1, 2], [3, 2, 0, 1], [2, -1, 0], [3, -1, 1]]

print vor.points
# [[ 0.  0.]
#  [ 0.  1.]
#  [ 0.  2.]
#  [ 1.  0.]
#  [ 1.  1.]
#  [ 1.  2.]
#  [ 2.  0.]
#  [ 2.  1.]
#  [ 2.  2.]]
print vor.point_region
# [1 8 6 3 7 4 2 9 5]

Now i should find somehow that vor.regions[7] is a region that belongs to the point vor.points[4]. How to achieve this?

You have a region you know, a point you don't, and you know that vor.point_region[point] == region. For a single region, you can figure out the corresponding point as:

point = np.argwhere(vor.point_region == region)

You can also create a region_point indexing array to figure out multiple points from an array of regions as:

region_point = np.argsort(vor.point_region)
points = region_point[regions-1]