Cómo generar diagramas de Voronoi con el algoritmo de Fortune

• Fortune’s Algorithm puede crear un diagrama de Voronoi en tiempo O(n log n), pero es difícil de implementar, así que salvo que se generen repetidamente diagramas grandes, una implementación O(n²) o una librería suele ser más práctica
• Un diagrama de Voronoi divide el plano en regiones según el site más cercano, y las fronteras se forman con puntos que están a la misma distancia de dos site
• El algoritmo mantiene una sweep line que avanza de izquierda a derecha y una beachline formada por arcos parabólicos, y solo procesa site event y circle event
• Un site event inserta un nuevo arco, divide un arco existente y crea una incomplete edge, mientras que un circle event elimina el arco del medio y completa un Voronoi Vertex y una half edge
• En una implementación práctica hay que manejar en conjunto la cola de eventos, la beachline, el mapa de incomplete edge y el DCEL, y la eliminación de circle event inválidos junto con la limpieza de las edge restantes aumenta mucho la complejidad

---

Dificultad de implementación y alcance de uso

• Fortune’s Algorithm es un algoritmo que genera diagramas de Voronoi en O(n log n)
• Si el objetivo es usarlo en la práctica, conviene evaluar primero la escala requerida antes de implementarlo directamente
  • A menos que haya que generar varios diagramas grandes por segundo, una implementación O(n²) puede ser una opción más sencilla
  • Una alternativa más realista es usar una librería existente
• El resultado del algoritmo es visualmente interesante, pero su implementación es complicada y bastante frustrante

Conceptos básicos del diagrama de Voronoi

• Un diagrama de Voronoi es una forma de dividir el plano en varias regiones, y se usa con frecuencia en generación procedural de mapas
• Los puntos elegidos como entrada se llaman site o seed
• La cell correspondiente a cada site es el conjunto de puntos del plano que están más cerca de ese site que de cualquier otro
• Las fronteras de las celdas están formadas por puntos a la misma distancia de dos site
• El Voronoi Vertex, donde se encuentran las esquinas de las celdas, es un punto a la misma distancia de tres site

sweep line, beachline y event

• Fortune’s Algorithm usa una sweep line, una línea vertical que se mueve de izquierda a derecha
• Cuando la sweep line encuentra un site, se crea un arco parabólico con ese site como foco, y el arco crece a medida que la sweep line se aleja
• El punto donde se cruzan dos arcos de site distintos está a la misma distancia de ambos site, así que se convierte en la frontera entre sus celdas
• Cuando dos fronteras se encuentran, se genera un vértice del diagrama
• El frente formado por los arcos activos se llama beachline
• En una implementación real no se mueve la sweep line píxel por píxel, sino que solo se procesan event en puntos específicos que sí pueden calcularse
  • site event: se define con coordenadas de site conocidas de antemano, y al procesarlo se agrega un nuevo arco a la beachline
  • circle event: se define por tres arcos de la beachline, y al procesarlo se elimina un arco, se genera un Voronoi Vertex y se crea una half edge

Encontrar fronteras con parábolas

• En este algoritmo las parábolas no se manejan con la forma habitual y = ax^2 + bx + c, sino con su locus definition
• Una parábola se define por un focus point y una directrix
  • El focus point corresponde al site
  • La directrix corresponde a la sweep line
• La intersección de dos parábolas que usan la misma sweep line como directrix está a la misma distancia de los dos site
• Por eso, si se encuentra la intersección de esas dos parábolas, se obtiene la equiedge entre ambos site
• Se usa pseudocódigo para calcular la coordenada x de una parábola, junto con un ejemplo de cómo la intersección entre dos parábolas se desplaza por la frontera cuando cambia la posición de la sweep line

Representación de la beachline y manejo de site event

• Cada arco de la beachline puede representarse solo con las coordenadas de su site correspondiente
  • La sweep line se aplica en común a todos los arcos
  • En la implementación, los arcos se manejan como coordenadas 2D, no como objetos separados
• La beachline puede representarse como un orden simple de puntos
  • Ejemplo: [arc1, arc2], [arc1, arc2, arc3]
  • Un arco del mismo site puede aparecer varias veces en la beachline
  • Ejemplo: [arc1, arc3, arc1, arc2]
• Cuando ocurre un site event, se busca el arco de la beachline que sería alcanzado al trazar una línea hacia la izquierda desde el nuevo site, y el nuevo arco divide ese arco
• Si el nuevo site L divide a j en una beachline existente [.., i, j, k, ..], la estructura pasa a ser [.., i, j, L, j, k, ..]
• Los site se colocan en la cola según su coordenada x, y cada vez que se procesa uno se actualizan la beachline y los candidatos a event

circle event y circumcircle

• En tres arcos de una beachline [.., i, j, k, ..], puede aparecer una situación en la que se encuentren dos fronteras y desaparezca el arco central j
• En ese momento existe un circumcircle que pasa por los tres site, y el centro del círculo está a la misma distancia de ellos
• El centro del circumcircle se convierte en un Voronoi Vertex
• El circle event se coloca en la cola de eventos tomando como referencia el circle point, es decir, el extremo derecho del círculo
• Si antes de llegar al circle point aparece un nuevo site dentro del círculo, el circle event anterior queda invalidado
  • Eso ocurre porque el nuevo site divide antes el arco central y la combinación de tres arcos deja de mantenerse
  • El triple original i, j, k desaparece y hay que revisar triples nuevos como i, j, L y L, j, k

incomplete edge y half edge

• Una incomplete edge es una línea cuyo extremo de un lado ya está fijado, pero cuyo otro extremo queda definido por la intersección de dos arcos parabólicos
• Cuando un site event inserta un nuevo arco, se crean dos incomplete edge
  • El punto fijo es la coordenada donde el nuevo arco se encontró con la beachline existente
  • Si el nuevo arco j divide al arco existente i, se crean edge correspondientes a las intersecciones [i, j] y [j, i]
• Cuando dos incomplete edge chocan en un circle event, ese punto de colisión se convierte en un Voronoi Vertex
• La incomplete edge existente se completa ahí como half edge, y entre los dos arcos que ahora quedan adyacentes se crea una nueva incomplete edge

Solo los círculos en sentido antihorario se convierten en circle event

• Si en la beachline aparece [i, j, k, j, i], tanto ijk como kji pueden formar un círculo, pero no ambos son circle event válidos
• El caso donde desaparece el arco central es solo el lado en que las fronteras realmente convergen
• En el programa, la orientation de tres puntos se determina con un determinante
  • Si el determinante es negativo, el orden es antihorario y sí hay circle event
  • Si el determinante es positivo, el orden es horario y no hay circle event
  • Si el determinante es 0, los tres puntos están alineados y no existe círculo

Flujo completo del algoritmo

• Los site de entrada se ordenan por coordenada x y se colocan en la cola como site event
• Se extrae y procesa el siguiente event hasta que la cola quede vacía
• Procesamiento de site event:
  • Se eliminan los circle event futuros en los que el nuevo site caería dentro del círculo
  • Se busca el arco de la beachline que el nuevo site va a dividir
  • Se inserta el nuevo arco y se divide el arco existente
  • Se agregan dos incomplete edge
  • Se revisa si los nuevos triples pueden generar circle event
• Procesamiento de circle event:
  • Se agrega el centro del circumcircle como Voronoi Vertex
  • Se elimina de la beachline el arco central
  • Se eliminan los circle event futuros que quedan inválidos por la desaparición de ese arco
  • Se revisan los triples de los arcos que ahora quedaron adyacentes para agregar nuevos circle event
• Cuando la cola queda vacía, las incomplete edge restantes se extienden hasta el borde del diagrama, y en los puntos donde tocan ese borde se crean Voronoi Vertex

Estructuras de datos en la implementación en Odin

• La implementación de ejemplo está escrita en Odin, un lenguaje alternativo a C
• Todo el código está en el repositorio RedPenguin101/voronoi
• Tipos básicos:
  • V2: punto 2D con forma [2]int
  • PointPair: par de dos V2
  • Event: estructura {site: bool, a, b, c: V2}
• El significado de Event cambia según el tipo
  • En un site event, a es la coordenada del site y b, c no se usan
  • En un circle event, a, b, c son los tres arcos de la beachline que generaron el event
• La estructura Fortune administra el siguiente estado
  • beachline: arreglo de V2
  • queue: arreglo de Event
  • incomplete_edges: mapa PointPair -> V2
  • vd: un DCEL que guarda el Voronoi Diagram

Partes omitidas o simplificadas en la implementación

• La beachline se representó con un vector, pero para mejorar la eficiencia sería más apropiado usar un binary tree
• Aunque conceptualmente la event queue es una priority queue, en la implementación de ejemplo se maneja insertando elementos ordenados en un arreglo
• La invalidación de circle event se hace recorriendo eventos futuros para verificarlos, y hay un TODO que indica que hace falta un método más rápido
• clean_beachline_edges es el procedimiento que recorta arcos innecesarios en ambos extremos de la beachline
• La implementación incluye manejo de excepciones, como site con la misma coordenada x, casos en que el circle point coincide con un site y choques entre puntos de referencia
• La etapa final, que limpia las incomplete edge restantes, las half edge sin twin y los vertex cuando la cola ya está vacía, se trata solo como un procesamiento matemático simple

Guardar el diagrama de Voronoi con DCEL

• El diagrama de Voronoi suele almacenarse como Doubly Connected Edge List (DCEL)
• DCEL es una estructura de datos que permite representar y manipular con facilidad un complejo de celdas formado por vertex y edge
• Aunque la representación está centrada en edge, también guarda información de vertex y face
• Una edge normal no tiene dirección, pero en DCEL cada edge se almacena como dos half edge en direcciones opuestas
• En un diagrama de Voronoi, los vertex que se guardan en el DCEL no son los site, sino los Voronoi Vertex
• El destino de una edge E se obtiene con E.twin.origin, y la face derecha con E.twin.left

[GN⁺]

Opiniones en Hacker News

• Hace tiempo hice una implementación en ClojureScript que mostraba con una animación cómo avanza el algoritmo de Fortune: https://voronoi.ajwerner.net/#/app-diagrams
  Es un algoritmo realmente hermoso.
  Sin embargo, después de ese proyecto le agarré algo de antipatía al algoritmo de Fortune, porque no tiene buena estabilidad numérica con punto flotante.
  Puede romperse si los puntos son colineales o casi colineales según el criterio de punto flotante.
  Si mal no recuerdo, delaunator es mejor en ese aspecto: https://github.com/mapbox/delaunator

  • La animación es la mejor que he visto hasta ahora.
    En la página de referencia se ve un enlace a una implementación “old”; me pregunto si existe la posibilidad de que la versión animada actual también se publique como open source.

• Hace unos años hice esta visualización 3D: https://x.com/KangarooPhysics/status/1253336959755251716

• Hay una implementación en JavaScript de Raymond Hill, famoso por uBlock Origin: https://github.com/gorhill/Javascript-Voronoi
  Aquí la modifiqué un poco para que se moviera: https://animations.adgent.com/voronoi.html

  • Esta animación me recuerda al estilo de A Scanner Darkly (2006).
    Me pregunto si se podría tomar un video como entrada y pasarlo por un algoritmo que lo muestre al estilo Voronoi.
    Llegado a ese punto quizá, estrictamente hablando, ya no sería un diagrama de Voronoi, pero se vería bastante genial.

• D3.js tiene una implementación nueva: https://github.com/d3/d3-delaunay
  En la parte inferior de esa página hay una explicación del algoritmo de barrido usado y una lista de implementaciones en otros lenguajes además de JavaScript.
  El antiguo d3-voronoi está previsto para quedar obsoleto, pero se puede ver aquí: https://github.com/d3/d3-voronoi

• Si no te interesan las aristas y solo quieres pintar cada punto con un color distinto, puedes usar una variante de flood fill que parte de los puntos semilla.
  Basta con poner un píxel en la pila solo cuando la distancia de ese color sea menor que la del color que ya está pintado en ese píxel.

  • Puedes crear una escena 3D con conos circulares rectos de distintos colores, con sus vértices en cada punto del plano 2D, y poner el eje perpendicular al plano.
    Si renderizas con una proyección ortográfica 2D desde arriba de los vértices, el z-buffer conserva el píxel del vértice más cercano.
    Seguramente también se pueda hacer con shaders, pero la demo clásica de conos 3D es muy fácil de entender e implementar.

• Es interesante que D3 haya pasado del algoritmo de Fortune a https://mapbox.github.io/delaunator/
  La razón es que “al crear una triangulación de Delaunay o un diagrama de Voronoi, es 5 a 10 veces más rápido que d3-voronoi, es numéricamente más robusto, trae renderizado en Canvas incorporado y ofrece recorrido del grafo de Delaunay y varias mejoras”.

  • Si D3 considera que delaunator es la mejor opción para este tipo de efectos, ya no tengo excusa para no agregarlo a mi biblioteca de canvas, salvo mi costumbre natural de procrastinar.
    El código actual para calcular tiles es dolorosamente ingenuo.
    Nueva discusión: https://github.com/KaliedaRik/Scrawl-canvas/discussions/120

• Este artículo me hizo buscar dónde anda Steve últimamente.
  Lo conocí hace varias décadas.

• Lectura relacionada que vale la pena: https://news.ycombinator.com/item?id=37998923 - Crear diagramas de Voronoi y triangulaciones de Delaunay en O(n log n) con el algoritmo de Fortune (2020)
  En el artículo y la discusión anteriores también hay breves resúmenes de otros algoritmos.
  Personalmente, sigue gustándome más el Jump Flooding Algorithm: https://en.wikipedia.org/wiki/Jump_flooding_algorithm