MUI_IIMA_1314_Geometria_Computacional_Voronoi

1451 days ago by trinidad

Diagramas de Voronoi

Algoritmo lento

Comenzaremos con el algoritmo sencillo visto en las transparencias. Tenemos que considerar todas las combinaciones de tres vértices, generar la circunferencia implícita y comprobar si algún otro punto está dentro o no (si no hay puntos dentro, será un vértice del diagrama de Voronoi).

En primer lugar, creamos una función auxiliar que recibe tres puntos del plano y devuelve las posiciones x,y del centro de la circunferencia que pasa por los tres puntos y el radio al cuadrado de dicha circunferencia.

def circunferencia_circunscrita(p1,p2,p3): # Recibe tres puntos y calcula el centro de la circunferencia circunscrita d=2*(p1[0]*(p2[1]-p3[1])+p2[0]*(p3[1]-p1[1])+p3[0]*(p1[1]-p2[1])) xc=((p1[1]^2+p1[0]^2)*(p2[1]-p3[1])+(p2[1]^2+p2[0]^2)*(p3[1]-p1[1])+(p3[1]^2+p3[0]^2)*(p1[1]-p2[1]))/d yc=-((p1[1]^2+p1[0]^2)*(p2[0]-p3[0])+(p2[1]^2+p2[0]^2)*(p3[0]-p1[0])+(p3[1]^2+p3[0]^2)*(p1[0]-p2[0]))/d r2=(p1[0]-xc)^2+(p1[1]-yc)^2 return xc,yc,r2 
       

A continuación, crearemos una función que calcule el diagrama de voronoi de una lista de puntos. La función recibirá una lista de puntos y devolverá una lista tal que en cada posición hay un vértice del diagrama de voronoi y los índices de los puntos que equidistan de dicho vértice.

def voronoi_lento(p): diagrama=[] n=len(p) for i in [0..n-3]: for j in [i+1..n-2]: for k in [j+1..n-1]: # Obtenemos el centro y el radio de la circunferencia circunscrita xc,yc,r2=circunferencia_circunscrita(p[i],p[j],p[k]) esVertice=True # Comprobamos que ningún otro punto esté for l in [0..n-1]: if l!=i and l!=j and l!=k and r2>((p[l][0]-xc)^2+(p[l][1]-yc)^2): esVertice=False if esVertice: diagrama.append([(xc,yc),i,j,k]) return diagrama 
       

Vamos a probarlo con un ejemplo.

puntos = [(random(),random()) for _ in range(10)] vd=voronoi_lento(puntos);vd 
       

Dibujamos los puntos originales (en azul) y los del diagrama de Voronoi en verde.

p1=points(puntos,size=50,faceted=True) pv=points([k[0] for k in vd],color='green',size=50,faceted=True) p1+pv 
       

Representación gráfica

Vamos a dibujar las regiones del diagrama de Voronoi. En primer lugar vamos a calcular para cada punto, los vértices del diagrama de voronoi que rodean al punto y los representamos mediante la envolvente convexa.

def regiones_voronoi_lento(pts): vd=voronoi_lento(pts) rV=[[v[0] for v in vd if v[1]==k or v[2]==k or v[3]==k] for k in range(len(pts))] return sum([Polyhedron(vertices=rV[k]).projection().render_outline_2d(color='black') for k in range(len(pts))]) 
       

Vamos a probar la función anterior.

dvd=regiones_voronoi_lento(puntos) (dvd+p1+pv).show(axes=False) 
       

Podemos observar que sólo funciona con las regiones acotadas. Podríamos modificar el método anterior para considerar las regiones no acotadas, sin embargo un método más sencillo es símplemente añadir algunos puntos externos al diagrama suficientemente lejos del borde para que el diagrama que veamos sea el correcto.

dvd=regiones_voronoi_lento(puntos+[(-10,-10),(-10,10),(10,10),(10,-10)]) (p1+pv+dvd).show(axes=False,xmin=0,xmax=1,ymin=0,ymax=1) 
       

Diagrama de Voronoi (Delaunay) con scipy

Vamos a calcular el diagrama de Voronoi utilizando el paquete scipy. En esta versión no está disponible el cálculo directo del diagrama de Voronoi, por lo que vamos a calcularla a partir de la triangulación de Delaunay.

Recordemos que la triangulación de Delaunay consiste en unir los puntos con sus vecinos de Thiessen.

from scipy.spatial import Delaunay tri = Delaunay(puntos) 
       

Podemos obtener los vértices de la triangulación mediante el atributo vertices. Nótese que se corresponden con lo obtenido mediante la función voronoi_lento, sálvo que faltan los vértices del diagrama de voronoi.

tri.vertices 
       

Añadimos los vértices al triángulo usando la función para calcular la circunferencia cincunscrita que creamos en el primer método.

def voronoi(pts): tri = Delaunay(pts) return [[circunferencia_circunscrita(pts[t[0]],pts[t[1]],pts[t[2]])[0:2],t[0],t[1],t[2]] for t in tri.vertices] 
       

Veamos que coincide con el método anterior.

vd2=voronoi(puntos);vd2 
       
p1=points(puntos,size=50,faceted=True) pv2=points([k[0] for k in vd2],color='yellow',size=50,faceted=True) p1+pv2 
       

Volvemos a definir la función para colorear regiones, ahora empleando el nuevo método.

def regiones_voronoi(pts): vd3=voronoi(pts) rV=[[v[0] for v in vd3 if v[1]==k or v[2]==k or v[3]==k] for k in range(len(pts))] return sum([Polyhedron(vertices=rV[k]).projection().render_outline_2d(color='black',linestyle='--') for k in range(len(pts))]) 
       
dvd2=regiones_voronoi(puntos+[(-10,-10),(-10,10),(10,10),(10,-10)]) (p1+pv2+dvd2).show(axes=False,xmin=0,xmax=1,ymin=0,ymax=1) 
       

Enlaces

A continuación se incluyen algunos enlaces interesantes sobre diagramas de Voronoi.

Voronoi en Wikipedia

Algoritmo de Fortune paso a paso en Sage

Mosaicos "Voronoi" a partir de imágenes

Libreria para diagramas Voronoi con Python (disponible en Sage a partir de la versión 5.10)

 

Ejercicios

1. Genera y dibuja el diagrama de Voronoi de los puntos $(1,0),(0,1),(-1,0),(0,-1),(0,0)$ mediante los dos métodos.

2. Determina en qué región está el punto $(0.7,0.3)$ (sin utilizar el dibujo).

3. Calcula cuántos puntos puede procesar cada algoritmo en un segundo. Usa el comando %time para calcular el tiempo y ve aumentando el número de puntos hasta que el tiempo sea aproximadamente un segundo.

 
       

Carga el archivo que está en la dirección http://matematicas.unex.es/~trinidad/mui/municipios.csv como archivo de datos de esta hoja de Sage. Para relacionar un municipio con su código, puedes utilizar la base de datos del INE.

4. Filtra los datos para quedarte con los correspondientes a la provincia en la que naciste. Calcula el diagrama de Voronoi y dibújalo. Es posible que tengas que cambiar la función Polyhedron por la función que utilizamos para calcular la envolvente convexa, creando además una función auxiliar:

def envolvente(puntos):
  return sum([line(k) for k in convex_hull(puntos)])

5. Localiza el municipio en el que naciste. Colorea la región en la que está dicho municipio. Encuentra los vecinos de Thiessen del municipio en el que naciste. Estima la altura del municipio como la media de las alturas de sus vecinos.

6. Dibuja la triangulación de Delaunay de la provincia eliminando el municipio en el que naciste. Determina en qué triángulo estaría el municipio en el que naciste. Calcula la altura del municipio por triangulación con los tres vértices de dicho triángulo.

 
       

Funciones auxiliares

def cuenta_derecha(puntos,i,j): # Devuelve el número de puntos de la lista que están "a la derecha" de la recta que pasa por los puntos i y j, orientada en sentido i,j. n=len(puntos) ux=puntos[j][0]-puntos[i][0] uy=puntos[j][1]-puntos[i][1] contador=0 for k in [0..n-1]: vx=puntos[k][0]-puntos[i][0] vy=puntos[k][1]-puntos[i][1] if k!=i and k!=j and ux*vy-uy*vx<0: contador=contador+1 return contador def angulo(v): # Calcula el ángulo del vector v=(x,y) respecto de la horizontal if v[0]==0: if v[1]>0: alpha=pi/2 else: alpha=-pi/2 else: if v[0]>0: alpha=arctan(v[1]/v[0]) else: alpha=arctan(v[1]/v[0])+pi return alpha def convex_hull(puntos): # Calcula la envolvente convexa de una lista de puntos ch=[] n=len(puntos) for i in [0..n-1]: for j in [0..n-1]: if i!=j and cuenta_derecha(puntos,i,j)==0: alpha=angulo([puntos[j][0]-puntos[i][0],puntos[j][1]-puntos[i][1]]) ch=ch+[(puntos[i],puntos[j],alpha)] ch=sorted(ch,key=lambda k:k[2]) # Los ordenamos sólo para poder pintarlos como polígono. return [(k[0],k[1]) for k in ch] def envolvente(puntos): return sum([line(k) for k in convex_hull(puntos)]) def regiones_voronoi(pts): vd3=voronoi(pts) rV=[[v[0] for v in vd3 if v[1]==k or v[2]==k or v[3]==k] for k in range(len(pts))] return sum([envolvente(rV[k]) for k in range(len(pts))])