algorithm - hacer - poligonos de voronoi qgis
Algoritmo para calcular un diagrama de Voronoi en una esfera? (11)
Estoy buscando un algoritmo simple (si existe) para encontrar el diagrama de Voronoi para un conjunto de puntos en la superficie de una esfera. El código fuente sería genial. Soy un hombre Delphi (sí, lo sé ...), pero también como C-code.
Actualización en julio de 2016:
Gracias a una serie de voluntarios (especialmente Nikolai Nowaczyk y yo), ahora hay un código mucho más robusto / correcto para manejar diagramas de Voronoi en la superficie de una esfera en Python. Esto está disponible oficialmente como scipy.spatial.SphericalVoronoi
de la versión 0.18
de scipy en adelante. Hay un ejemplo de uso y trazado en los docs oficiales.
El algoritmo sigue la complejidad del tiempo cuadrático. Mientras que loglineal es el óptimo teórico para los diagramas de Voronoi en las superficies de las esferas, este es actualmente el mejor que hemos podido implementar. Si desea obtener más información y ayudar con el esfuerzo de desarrollo, existen algunos problemas abiertos relacionados con la mejora de la forma en que Python maneja los diagramas esféricos de Voronoi y las estructuras de datos relacionadas:
- Esfuerzo para mejorar el trazado de polígonos esféricos en matplotlib
- Esfuerzo para mejorar el manejo de los cálculos esféricos del área superficial del polígono en scipy
Para obtener más antecedentes sobre la teoría / desarrollo / desafíos relacionados con este código de Python y los esfuerzos de geometría computacional relacionados, también puede consultar algunas charlas de Nikolai y yo:
- Charla de Nikolai PyData London 2016
- Charla de Tyler PyData London 2015
- Tutorial de Tyler PyCon 2016 de Geometría Computacional
Respuesta Original:
De hecho, recientemente escribí código de Python de código abierto para diagramas de Voronoi en la superficie de una esfera: https://github.com/tylerjereddy/py_sphere_Voronoi
El uso, el algoritmo y las limitaciones están documentados en readthedocs ( http://py-sphere-voronoi.readthedocs.org/en/latest/voronoi_utility.html ). Hay algunos ejemplos detallados allí, pero ubicaré uno o dos a continuación también. El módulo también maneja el cálculo de las áreas superficiales de la región de Voronoi, aunque con algunos puntos débiles numéricos en la versión de desarrollo actual.
No he visto muchas implementaciones de código abierto bien documentadas para diagramas esféricos de Voronoi, pero ha habido un poco de rumores sobre la implementación de JavaScript en el sitio web de Jason Davies ( http://www.jasondavies.com/maps/voronoi/ ) . Sin embargo, no creo que su código esté abierto. También vi una publicación en el blog sobre el uso de Python para tratar parte del problema ( http://jellymatter.com/2014/01/29/voronoi-tessellation-on-the-surface-of-a-sphere-python-code/ ). Muchas de las fuentes de literatura primaria citadas en las publicaciones anteriores parecían muy difíciles de implementar (probé algunas de ellas) pero tal vez algunas personas encuentren útil mi implementación o incluso sugieran formas de mejorarla.
Ejemplos:
1) Produzca un diagrama de Voronoi para un conjunto pseudoaleatorio de puntos en la esfera de la unidad:
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.colors as colors
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import numpy as np
import scipy as sp
import voronoi_utility
#pin down the pseudo random number generator (prng) object to avoid certain pathological generator sets
prng = np.random.RandomState(117) #otherwise, would need to filter the random data to ensure Voronoi diagram is possible
#produce 1000 random points on the unit sphere using the above seed
random_coordinate_array = voronoi_utility.generate_random_array_spherical_generators(1000,1.0,prng)
#produce the Voronoi diagram data
voronoi_instance = voronoi_utility.Voronoi_Sphere_Surface(random_coordinate_array,1.0)
dictionary_voronoi_polygon_vertices = voronoi_instance.voronoi_region_vertices_spherical_surface()
#plot the Voronoi diagram
fig = plt.figure()
fig.set_size_inches(2,2)
ax = fig.add_subplot(111, projection=''3d'')
for generator_index, voronoi_region in dictionary_voronoi_polygon_vertices.iteritems():
random_color = colors.rgb2hex(sp.rand(3))
#fill in the Voronoi region (polygon) that contains the generator:
polygon = Poly3DCollection([voronoi_region],alpha=1.0)
polygon.set_color(random_color)
ax.add_collection3d(polygon)
ax.set_xlim(-1,1);ax.set_ylim(-1,1);ax.set_zlim(-1,1);
ax.set_xticks([-1,1]);ax.set_yticks([-1,1]);ax.set_zticks([-1,1]);
plt.tick_params(axis=''both'', which=''major'', labelsize=6)
2) Calcule las áreas de superficie de los polígonos de la región de Voronoi y verifique que la superficie reconstituida sea sensible:
import math
dictionary_voronoi_polygon_surface_areas = voronoi_instance.voronoi_region_surface_areas_spherical_surface()
theoretical_surface_area_unit_sphere = 4 * math.pi
reconstituted_surface_area_Voronoi_regions = sum(dictionary_voronoi_polygon_surface_areas.itervalues())
percent_area_recovery = round((reconstituted_surface_area_Voronoi_regions / theoretical_surface_area_unit_sphere) * 100., 5)
print percent_area_recovery
97.87551 #that seems reasonable for now
Aquí hay un buen programa de ejemplo del diagrama de Voronoi (incluido el código fuente de Delphi 5/6).
Creo que "puntos en la superficie de una esfera" significa que primero tienes que reasignarlos a coordenadas 2D, crear el diagrama de Voronoi y luego reasignarlos a las coordenadas de la superficie de la esfera. ¿Las dos fórmulas del artículo de mapeo UV de Wikipedia funcionan aquí?
Observe también que el diagrama de Voronoi tendrá una topología incorrecta (está dentro de un rectángulo y no se "ajusta"), aquí podría ayudar copiar todos los puntos de (0,0) - (x, y) al vecino regiones arriba (0, -y * 2) - (x, 0), debajo (0, y) - (x, y * 2), izquierda (-x, 0) - (0, y) y derecha (x, 0) - (x * 2, y). Espero que sepas a qué me refiero, no dudes en preguntar :)
Aquí hay un documento sobre diagramas esféricos de Voronoi .
O si comprendes Fortran (bleah!) Hay este sitio .
Citando de esta referencia: http://www.qhull.org/html/qdelaun.htm
Para calcular la triangulación de Delaunay de puntos en una esfera, calcule su casco convexo. Si la esfera es la esfera de la unidad en el origen, las normales facetas son los vértices Voronoi de la entrada.
Creo que el plano Voronoi para cada punto se puede construir utilizando una geometría no euclidiana. Lo que normalmente era una línea en un plano 2d, ahora es un "gran círculo" en la esfera (ver Wikipedia: geometría elíptica ). Es fácil encontrar qué puntos están en el lado equivocado de cualquier gran círculo entre dos puntos, simplemente girando la esfera de modo que el gran círculo divisorio sea el ecuador, y luego sean todos los puntos del otro hemisferio que el punto en el que se encuentre construyendo el avión Voronoi para.
Esta no es la respuesta completa, pero aquí es donde comenzaría ...
En resumen, prueba cssgrid desde NCAR Graphics . Escribí una respuesta más larga para una pregunta similar en codereview.stackexchange.com .
Ha pasado un tiempo desde que la pregunta fue respondida, pero encontré dos documentos que implementan el algoritmo de Fortune (eficiencia O (N lg N), memoria O (N)) sobre la superficie de la esfera. Tal vez un futuro espectador encuentre útil esta información.
- "Barriendo la Esfera" de Dinis y Mamede, publicado en el Simposio Internacional de 2010 sobre Diagramas de Voronoi en Ciencia e Ingeniería. Se puede comprar en http://dx.doi.org/10.1109/ISVD.2010.32
- "Algoritmo de barrido de plano para la teselación de la esfera de Voronoi" por Zheng et al. No estoy seguro de que se haya publicado debido a la primera, pero está fechada el 13 de diciembre de 2011. Está disponible de forma gratuita en e-lc.org/tmp/Xiaoyu__Zheng_2011_12_05_14_35_11.pdf
Estoy trabajando en ellos en este momento, así que no puedo explicarlo bien. La idea básica es que el algoritmo de Fortune funciona en la superficie de la esfera siempre que calcules las parábolas de delimitación de los puntos correctamente. Debido a que la superficie de la esfera se envuelve, también puede usar una lista circular para contener la línea de playa y no preocuparse por manejar las celdas en el borde del espacio rectangular. Con eso, puedes barrer desde el polo norte de la esfera hacia el sur y hacer una copia de seguridad otra vez, saltando a sitios que introducen nuevos puntos a la línea de la playa (agregando una parábola a la línea de la playa) o la introducción de vértices celulares (eliminando un parábola de la línea de la playa).
Ambos documentos esperan un alto nivel de comodidad con álgebra lineal para comprender los conceptos, y ambos me siguen perdiendo en el momento en que comienzan a explicar el algoritmo en sí. Tampoco proporcionamos el código fuente, desafortunadamente.
Hay un documento de INRIA sobre la Triangulación de Delaunay (DT) de los puntos que se encuentran en una esfera: CAROLI, Manuel, et al. Triangulaciones de Delaunay robustas y eficientes de puntos en o cerca de una esfera. 2009. donde hablan de una implementación en CGAL .
El documento se refiere a varias implementaciones disponibles de algoritmos DT.
Citando del periódico:
Una respuesta fácil y estándar consiste en calcular el casco convexo tridimensional de los puntos, que es notoriamente equivalente.
para calcular el casco convexo, el artículo sugiere:
- Hull, un programa para cascos convexos.
- Qhull .
- Cascos convexos tridimensionales. en FORTRAN. Cáscaras convexas tridimensionales.
- STRIPACK en FORTRAN.
La clase DT C ++ de CGAL tiene el método dual
para obtener el diagrama de Voronoi.
De acuerdo con este post de Monique Teillaud (una de las autoras del mencionado artículo), me parece que en noviembre de 2012 la implementación aún no estaba lista.
Si tus puntos están dentro de un hemisferio, podrías hacer una proyección gnomónica desde coordenadas esféricas a planas, y luego triangular, ya que los círculos grandes se convierten en líneas rectas de la distancia más corta.
Tenga en cuenta que la triangulación de Delaunay en una esfera es solo el casco convexo. Por lo tanto, puede calcular el casco convexo tridimensional (por ejemplo, utilizando CGAL) y tomar el doble.
CGAL está trabajando en el paquete "kernel esférico", que permitiría calcular exactamente este tipo de cosas. Desafortunadamente, aún no se ha lanzado , pero tal vez sea en su próximo lanzamiento, ya que ya lo mencionaron en una charla técnica de google en marzo.