En las últimas semanas, he estado trabajando en la implementación
del algoritmo de Fortune en C ++. Este algoritmo toma muchos puntos 2D y construye
un diagrama de Voronoi a partir de ellos. Si no sabe qué es un diagrama de Voronoi, eche un vistazo a la figura:
Para cada punto de entrada, que se llama un "sitio", necesitamos encontrar muchos puntos que estén más cerca de este lugar que de todos los demás. Tales conjuntos de puntos forman las celdas que se muestran en la imagen de arriba.
Es notable en el algoritmo de Fortune que construya tales diagramas a tiempo
(que es óptimo para un algoritmo de comparación), donde
Es el número de lugares.
Estoy escribiendo este artículo porque considero que la implementación de este algoritmo es una tarea muy difícil. Por el momento, este es el más complicado de los algoritmos que tuve que implementar. Por lo tanto, quiero compartir los problemas que encontré y cómo resolverlos.
Como de costumbre, el código se publica en
github , y todos los materiales de referencia que utilicé se enumeran al final del artículo.
Descripción del algoritmo de fortuna
No explicaré cómo funciona el algoritmo, porque otras personas ya lo han hecho bien. Puedo recomendar estudiar estos dos artículos:
aquí y
aquí . El segundo es muy interesante: el autor escribió una demostración interactiva en Javascript, que es útil para comprender el funcionamiento del algoritmo. Si necesita un enfoque más formal con toda la evidencia, le recomiendo leer el Capítulo 7 de
Geometría computacional, 3a edición .
Además, prefiero tratar con detalles de implementación que no están bien documentados. Y son ellos quienes hacen que la implementación del algoritmo sea tan compleja. En particular, me centraré en las estructuras de datos utilizadas.
Acabo de escribir un pseudocódigo del algoritmo para que tengas una idea de su estructura global:
agregar un evento de lugar a la cola de eventos para cada lugar
hasta que la cola de eventos esté vacía
recuperar el evento principal
si el evento es un evento de lugar
inserte un nuevo arco en la costa
buscar nuevos eventos circulares
de lo contrario
crear un vértice en el diagrama
quitamos de la costa un arco apretado
eliminar eventos inválidos
buscar nuevos eventos circulares
Estructura de datos del gráfico
El primer problema que encontré fue elegir la forma de almacenar el diagrama de Voronoi.
Decidí usar una estructura de datos que se usa ampliamente en geometría computacional llamada lista de bordes doblemente conectados (DCEL).
Mi clase
VoronoiDiagram usa cuatro contenedores como campos:
class VoronoiDiagram { public:
Hablaré en detalle sobre cada uno de ellos.
La clase de
Site describe el punto de entrada. Cada lugar tiene un índice, que es útil para ponerlo en la cola, coordenadas y un puntero a la celda (
face ):
struct Site { std::size_t index; Vector2 point; Face* face; };
Los vértices de la celda están representados por la clase
Vertex , solo tienen un campo de coordenadas:
struct Vertex { Vector2 point; };
Aquí está la implementación de los medios bordes:
struct HalfEdge { Vertex* origin; Vertex* destination; HalfEdge* twin; Face* incidentFace; HalfEdge* prev; HalfEdge* next; };
Te preguntarás, ¿qué es una media costilla? Un borde en el diagrama de Voronoi es común a dos celdas vecinas. En la estructura de datos DCEL, dividimos estos bordes en dos medios bordes, uno para cada celda, y están unidos por un puntero
twin . Además, el medio borde tiene un punto inicial y final. El campo
incidentFace indica la cara a la que pertenece el medio borde. Las celdas en DCEL se implementan como una lista doblemente cíclica de medios bordes, en la que los medios bordes adyacentes están conectados entre sí. Por lo tanto, los campos
prev y
next indican los medios bordes anteriores y siguientes en la celda.
La siguiente figura muestra todos estos campos para el medio borde rojo:
Finalmente, la clase
Face define la celda. Simplemente contiene un puntero a su lugar y otro a una de sus medias costillas. No importa cuál de los medios bordes esté seleccionado, porque la celda es un polígono cerrado. Por lo tanto, tenemos acceso a todos los medios bordes al atravesar una lista vinculada cíclica.
struct Face { Site* site; HalfEdge* outerComponent; };
Cola de eventos
La forma estándar de implementar una cola de eventos es con una cola de prioridad. En el proceso de procesamiento de eventos de lugar y círculo, es posible que debamos eliminar los eventos de círculo de la cola porque ya no son válidos. Pero la mayoría de las implementaciones de cola de prioridad estándar no le permiten eliminar un elemento que no está en la parte superior. Esto se aplica en particular a
std::priority_queue .
Hay dos formas de resolver este problema. El primero, más simple, es agregar una bandera
valid a los eventos.
valid se establece inicialmente en
true . Luego, en lugar de eliminar el evento circular de la cola, simplemente podemos establecer su bandera en
false . Finalmente, al procesar todos los eventos en el bucle principal, verificamos si el valor de indicador
valid del evento es
false , y si es así, simplemente omítelo y procesa el siguiente.
El segundo método que apliqué fue no usar
std::priority_queue . En su lugar, implementé mi propia cola de prioridad, que admite la eliminación de cualquier elemento contenido en ella. La implementación de tal cola es bastante simple. Elegí este método porque hace que el código del algoritmo sea más limpio.
Costa
La estructura de datos de la costa es una parte compleja del algoritmo. En caso de implementación incorrecta, no hay garantías de que el algoritmo se ejecute en
. La clave para obtener esta complejidad de tiempo es usar un árbol de auto-equilibrio. ¡Pero es más fácil decirlo que hacerlo!
Se recomienda que la mayoría de los recursos que he estudiado (los dos artículos mencionados anteriormente y el libro de
Geometría computacional ) implementen la costa como un árbol en el que los nodos internos indican puntos de ruptura y las hojas indican arcos. Pero no dicen nada sobre cómo equilibrar un árbol. Creo que ese modelo no es el mejor, y he aquí por qué:
- hay información redundante: sabemos que hay un punto de ruptura entre dos arcos adyacentes, por lo que no es necesario representar estos puntos como nodos
- no es adecuado para el autoequilibrio: solo se puede equilibrar el subárbol formado por puntos de ruptura. Realmente no podemos equilibrar todo el árbol, porque de lo contrario los arcos pueden convertirse en nodos internos y hojas de puntos de ruptura. Escribir un algoritmo para equilibrar solo el subárbol formado por nodos internos me parece una pesadilla.
Por lo tanto, decidí presentar la costa de manera diferente. En mi implementación, la costa también es un árbol, pero todos los nodos son arcos. Tal modelo no tiene ninguna de las desventajas enumeradas.
Aquí está la definición de
Arc Arc en mi implementación:
struct Arc { enum class Color{RED, BLACK};
Los primeros tres campos se utilizan para estructurar el árbol. El campo
leftHalfEdge indica el medio borde dibujado por el punto más a la izquierda del arco. Y
rightHalfEdge está en el borde medio dibujado por el punto extremo derecho. Se utilizan dos punteros,
prev y
next para obtener acceso directo al arco anterior y siguiente de la costa. Además, también le permiten omitir la costa como una lista doblemente vinculada. Finalmente, cada arco tiene un color que se usa para equilibrar la costa.
Para equilibrar la costa, decidí usar un
esquema rojo-negro . Al escribir código, me inspiró el libro
Introducción a los algoritmos . El Capítulo 13 describe dos algoritmos interesantes,
insertFixup y
deleteFixup , que equilibran el árbol después de la inserción o eliminación.
Sin embargo, no puedo usar el método de
insert que se muestra en el libro, porque las teclas se usan para encontrar el punto de inserción del nodo. No hay claves en el algoritmo de Fortune, solo sabemos que necesitamos insertar un arco antes o después de otro en la costa. Para implementar esto, creé los
insertAfter insertBefore e
insertAfter :
void Beachline::insertBefore(Arc* x, Arc* y) {
Insertar
y antes de
x se realiza en tres pasos:
- Encuentre un lugar para insertar un nuevo nodo. Para hacer esto, utilicé la siguiente observación: el niño izquierdo
x o el niño derecho x->prev es Nil , y el que es Nil está antes de x y después de x->prev . - Dentro de la línea costera, mantenemos la estructura de una lista doblemente vinculada, por lo que debemos actualizar los punteros
prev y next de los elementos x->prev , y y x . - Finalmente, simplemente llamamos al método
insertFixup descrito en el libro para equilibrar el árbol.
insertAfter se implementa de manera similar.
El método de eliminación tomado del libro se puede implementar sin cambios.
Límite de la tabla
Aquí está la salida del algoritmo de Fortune descrito anteriormente:
Hay un pequeño problema con algunos bordes de las celdas en el borde de la imagen: no se dibujan porque son infinitos.
Peor aún, una célula puede no ser un solo fragmento. Por ejemplo, si tomamos tres puntos que están en la misma línea, entonces el punto medio tendrá dos medios bordes infinitos que no están conectados entre sí. Esto no nos conviene muy bien, porque no podremos acceder a uno de los medios bordes, porque la celda es una lista vinculada de bordes.
Para resolver estos problemas, limitaremos el diagrama. Con esto quiero decir que limitaremos cada celda del diagrama para que ya no tengan bordes infinitos y cada celda sea un polígono cerrado.
Afortunadamente, el algoritmo de Fortune nos permite encontrar rápidamente bordes interminables: corresponden a medios bordes aún en la costa al final del algoritmo.
Mi algoritmo de restricción recibe un cuadro como entrada y consta de tres pasos:
- Proporciona la colocación de cada vértice del diagrama dentro del rectángulo.
- Cortar cada borde infinito.
- Cierra celdas.
La etapa 1 es trivial: solo necesitamos expandir el rectángulo si no contiene un vértice.
La etapa 2 también es bastante simple: consiste en calcular las intersecciones entre los rayos y el rectángulo.
La etapa 3 tampoco es muy complicada, solo requiere atención. Lo realizo en dos etapas. Primero, agrego los puntos de esquina del rectángulo a las celdas en los vértices de las cuales deberían estar. Luego me aseguro de que todos los vértices de la celda estén conectados por medios bordes.
Le recomiendo que estudie el código y haga preguntas si necesita detalles sobre esta parte.
Aquí está el diagrama de salida del algoritmo delimitador:
Ahora vemos que todos los bordes están dibujados. Y si aleja la imagen, puede asegurarse de que todas las celdas estén cerradas:
Intersección con rectángulo
Genial Pero la primera imagen del comienzo del artículo es mejor, ¿verdad?
En muchas aplicaciones, es útil tener la intersección entre el diagrama de Voronoi y el rectángulo, como se muestra en la primera imagen.
Lo bueno es que después de restringir el gráfico, es mucho más fácil de hacer. La mala noticia es que, aunque el algoritmo no es muy complicado, debemos tener cuidado.
La idea es esta: vamos alrededor de su medio borde de cada celda y verificamos la intersección entre el medio borde y el rectángulo. Son posibles cinco casos:
- La media costilla está completamente dentro del rectángulo: guardamos una media costilla
- La media costilla está completamente fuera del rectángulo: descartamos dicha media costilla
- La media costilla sale del rectángulo: truncamos la media costilla y la guardamos como la última media costilla que sale .
- La media costilla va dentro del rectángulo: truncamos la media costilla para conectarla con la última media costilla que salió (la guardamos en el caso 3 o 5)
- La media costilla cruza el rectángulo dos veces: truncamos la media costilla y agregamos una media costilla para asociarla con la última media costilla que salió , y luego la guardamos como la última última media costilla que salió .
Sí, ha habido muchos casos. Creé una imagen para mostrarlos a todos:
El polígono naranja es la celda original, y el rojo es la celda truncada. Las medias costillas truncadas están marcadas en rojo. Se han agregado costillas verdes para conectar las medias costillas que ingresan al rectángulo con las medias costillas que salen.
Aplicando este algoritmo a un diagrama acotado, obtenemos el resultado esperado:
Conclusión
El artículo resultó ser bastante largo. Y estoy seguro de que muchos aspectos aún no están claros para ti. Sin embargo, espero que te sea útil. Examine el código para más detalles.
Para resumir y asegurarme de que no perdimos el tiempo en vano, medí en mi computadora portátil (barata) el tiempo para calcular el diagrama de Voronoi para un número diferente de lugares:
- : 2 ms
- : 33 ms
- : 450 ms
- : 6600 ms
No tengo nada con lo que comparar estos indicadores, ¡pero parece que es increíblemente rápido!
Referencias