Encontrar regiones voronoi que contienen una lista de coordenadas arbitrarias

Estoy trabajando con un algoritmo que, para cada iteración, necesita encontrar a qué región de un diagrama Voronoi pertenece un conjunto de coordinadores arbitrarios. es decir, en qué región se encuentra cada coordenada. (Podemos suponer que todas las coordenadas pertenecerán a una región, si eso marca alguna diferencia).

Todavía no tengo ningún código que funcione en Python, pero el pseudo código se parece a esto:

## we are in two dimensions and we have 0<x<1, 0<y<1.

for i in xrange(1000):
  XY = get_random_points_in_domain()
  XY_candidates = get_random_points_in_domain()
  vor = Voronoi(XY) # for instance scipy.spatial.Voronoi
  regions = get_regions_of_candidates(vor,XY_candidates) # this is the function i need

  ## use regions for something

Sé que el scipy.Delaunay tiene una función llamada find_simplex que hará lo que yo quiera por simplicidad en una triangulación de Delaunay, pero necesito el diagrama Voronoi, y construir ambos es algo que deseo evitar.

Preguntas:

1. ¿Hay alguna biblioteca de algún tipo que me permita hacer esto fácilmente?

2. Si no, ¿hay un buen algoritmo que pueda ver que me permita hacerlo de manera eficiente?

Actualizar

La solución de Jamie es exactamente lo que quería. Aunque me da un poco de vergüenza no haberlo pensado yo mismo ...

Respuestas a la pregunta(1)

Su respuesta a la pregunta