Nas últimas semanas, tenho trabalhado na implementação do algoritmo da Fortune em C ++. Esse algoritmo usa muitos pontos 2D e constrói um diagrama de Voronoi a partir deles. Se você não sabe o que é um diagrama de Voronoi, dê uma olhada na figura:


Para cada ponto de entrada, chamado de "site", precisamos encontrar muitos pontos mais próximos deste local do que de todos os outros. Esses conjuntos de pontos formam as células mostradas na imagem acima.

É notável no algoritmo da Fortune que ele construa esses diagramas com o tempo $$$O (n \ log n)$(ideal para um algoritmo de comparação), em que $$$n$É o número de lugares.

Estou escrevendo este artigo porque considero a implementação desse algoritmo uma tarefa muito difícil. No momento, esse é o algoritmo mais complicado que eu tive que implementar. Portanto, quero compartilhar os problemas que encontrei e como resolvê-los.

Como de costume, o código é publicado no github e todos os materiais de referência que usei estão listados no final do artigo.

Descrição do algoritmo da sorte

Não vou explicar como o algoritmo funciona, porque outras pessoas já o fizeram bem. Posso recomendar o estudo desses dois artigos: aqui e aqui . O segundo é muito interessante - o autor escreveu uma demonstração interativa em Javascript, que é útil para entender a operação do algoritmo. Se você precisar de uma abordagem mais formal com todas as evidências, recomendo a leitura do Capítulo 7 da Geometria Computacional, 3ª edição .

Além disso, prefiro lidar com detalhes de implementação que não estão bem documentados. E são eles que tornam a implementação do algoritmo tão complexa. Em particular, vou me concentrar nas estruturas de dados usadas.

Acabei de escrever um pseudo-código do algoritmo para que você tenha uma idéia da estrutura global:

  adicione um evento de local à fila de eventos para cada local
 até que a fila de eventos esteja vazia
     recuperar o evento principal
     se o evento for um evento de local
         insira um novo arco no litoral
         verifique se há novos eventos de círculo
     caso contrário
         crie um vértice no diagrama
         removemos do litoral um arco apertado
         excluir eventos inválidos
         verifique se há novos eventos de círculo 

Estrutura de dados do gráfico

O primeiro problema que encontrei foi escolher a maneira de armazenar o diagrama Voronoi.

Decidi usar uma estrutura de dados amplamente usada em geometria computacional chamada lista de arestas duplamente conectadas (DCEL).

Minha classe VoronoiDiagram usa quatro contêineres como campos:

 class VoronoiDiagram { public: // ... private: std::vector<Site> mSites; std::vector<Face> mFaces; std::list<Vertex> mVertices; std::list<HalfEdge> mHalfEdges; } 

Vou falar em detalhes sobre cada um deles.

A classe Site descreve o ponto de entrada. Cada local possui um índice, útil para colocá-lo na fila, coordenadas e um ponteiro para a célula ( face ):

 struct Site { std::size_t index; Vector2 point; Face* face; }; 

Os vértices da célula são representados pela classe Vertex , eles possuem apenas um campo de coordenadas:

 struct Vertex { Vector2 point; }; 

Aqui está a implementação das meias arestas:

 struct HalfEdge { Vertex* origin; Vertex* destination; HalfEdge* twin; Face* incidentFace; HalfEdge* prev; HalfEdge* next; }; 

Você pode se perguntar, o que é meia costela? Uma aresta no diagrama de Voronoi é comum a duas células vizinhas. Na estrutura de dados do DCEL, dividimos essas arestas em duas meias arestas, uma para cada célula, e elas são vinculadas por um ponteiro twin . Além disso, a meia borda tem um ponto inicial e final. O campo incidentFace indica a face à qual a meia borda pertence. As células no DCEL são implementadas como uma lista cíclica de meias arestas duplamente ligada, na qual as meias arestas adjacentes são conectadas. Portanto, os campos prev e next indicam as meias-arestas anteriores e seguintes na célula.

A figura abaixo mostra todos esses campos para a meia borda vermelha:


Finalmente, a classe Face define a célula. Ele simplesmente contém um ponteiro para o seu lugar e outro para uma de suas meias costelas. Não importa qual das meias arestas está selecionada, porque a célula é um polígono fechado. Assim, temos acesso a todas as meias-arestas enquanto percorremos uma lista vinculada cíclica.

 struct Face { Site* site; HalfEdge* outerComponent; }; 

Fila de eventos

A maneira padrão de implementar uma fila de eventos é com uma fila prioritária. No processo de processamento de eventos de local e círculo, talvez seja necessário remover eventos de círculo da fila porque eles não são mais válidos. Mas a maioria das implementações de fila de prioridade padrão não permite excluir um item que não está no topo. Isso se aplica em particular ao std::priority_queue .

Existem duas maneiras de resolver esse problema. O primeiro, mais simples, é adicionar um sinalizador valid aos eventos. valid é inicialmente definido como true . Em seguida, em vez de remover o evento circle da fila, podemos simplesmente definir seu sinalizador como false . Finalmente, ao processar todos os eventos no loop principal, verificamos se o valor do sinalizador valid do evento é false e, nesse caso, simplesmente pule e processe o próximo.

O segundo método que apliquei não era usar std::priority_queue . Em vez disso, implementei minha própria fila de prioridades, que suporta a remoção de qualquer elemento contido nela. A implementação dessa fila é bastante simples. Eu escolhi esse método porque torna o código do algoritmo mais limpo.

Costa

A estrutura de dados da linha de costa é uma parte complexa do algoritmo. No caso de implementação incorreta, não há garantias de que o algoritmo será executado em $$$O (n \ log n)$. A chave para obter essa complexidade de tempo é usar uma árvore de auto-equilíbrio. Mas é mais fácil falar do que fazer!

A maioria dos recursos que estudei (os dois artigos mencionados acima e o livro Geometria Computacional ) são aconselhados a implementar o litoral como uma árvore na qual os nós internos indicam pontos de interrupção e as folhas indicam arcos. Mas eles não dizem nada sobre como equilibrar uma árvore. Eu acho que esse modelo não é o melhor, e aqui está o porquê:

• há informações redundantes: sabemos que existe um ponto de interrupção entre dois arcos adjacentes; portanto, não é necessário representar esses pontos como nós
• é inadequado para o auto-equilíbrio: somente a subárvore formada por pontos de interrupção pode ser equilibrada. Realmente não podemos equilibrar a árvore inteira, porque, caso contrário, os arcos podem se tornar nós internos e folhas de pontos de interrupção. Escrever um algoritmo para equilibrar apenas a subárvore formada por nós internos parece um pesadelo para mim.

Por isso, decidi apresentar o litoral de maneira diferente. Na minha implementação, o litoral também é uma árvore, mas todos os nós são arcos. Esse modelo não possui nenhuma das desvantagens listadas.

Aqui está a definição do Arc arco na minha implementação:

 struct Arc { enum class Color{RED, BLACK}; // Hierarchy Arc* parent; Arc* left; Arc* right; // Diagram VoronoiDiagram::Site* site; VoronoiDiagram::HalfEdge* leftHalfEdge; VoronoiDiagram::HalfEdge* rightHalfEdge; Event* event; // Optimizations Arc* prev; Arc* next; // Only for balancing Color color; }; 

Os três primeiros campos são usados ​​para estruturar a árvore. O campo leftHalfEdge indica a meia aresta desenhada pelo ponto mais à esquerda do arco. E rightHalfEdge está na meia-borda desenhada pelo ponto extremo direito. Dois ponteiros, prev e next são usados ​​para obter acesso direto ao arco anterior e próximo do litoral. Além disso, eles também permitem que você ignore o litoral como uma lista duplamente vinculada. Finalmente, cada arco tem uma cor usada para equilibrar a costa.

Para equilibrar a costa, decidi usar um esquema vermelho-preto . Ao escrever código, fui inspirado pelo livro Introdução aos algoritmos . O capítulo 13 descreve dois algoritmos interessantes, insertFixup e deleteFixup , que equilibram a árvore após a inserção ou exclusão.

No entanto, não posso usar o método de insert mostrado no livro, porque as chaves são usadas para encontrar o ponto de inserção do nó. Não há chaves no algoritmo da Fortune, apenas sabemos que precisamos inserir um arco antes ou depois de outro no litoral. Para implementar isso, criei os insertAfter insertBefore e insertAfter :

 void Beachline::insertBefore(Arc* x, Arc* y) { // Find the right place if (isNil(x->left)) { x->left = y; y->parent = x; } else { x->prev->right = y; y->parent = x->prev; } // Set the pointers y->prev = x->prev; if (!isNil(y->prev)) y->prev->next = y; y->next = x; x->prev = y; // Balance the tree insertFixup(y); } 

A inserção de y antes de x é realizada em três etapas:

1. Encontre um local para inserir um novo nó. Para fazer isso, usei a seguinte observação: o filho esquerdo x ou o filho direito x->prev é Nil , e o que é Nil é anterior a x e depois de x->prev .
2. Dentro do litoral, mantemos a estrutura de uma lista duplamente vinculada, portanto, devemos atualizar os ponteiros prev e next dos elementos x->prev , y x .
3. Finalmente, chamamos simplesmente o método insertFixup descrito no livro para equilibrar a árvore.

insertAfter é implementado da mesma forma.

O método de remoção retirado do livro pode ser implementado sem alterações.

Limite de gráfico

Aqui está a saída do algoritmo da Fortune descrito acima:


Há um pequeno problema com algumas bordas das células na borda da imagem: elas não são desenhadas porque são infinitas.

Pior, uma célula pode não ser um único fragmento. Por exemplo, se pegarmos três pontos que estão na mesma linha, o ponto médio terá duas meias arestas infinitas que não estão conectadas. Isso não nos serve muito, porque não conseguiremos acessar uma das meias arestas, porque a célula é uma lista vinculada de arestas.

Para resolver esses problemas, limitaremos o diagrama. Com isso, quero dizer que limitaremos cada célula do diagrama para que elas não tenham mais arestas infinitas e cada célula seja um polígono fechado.

Felizmente, o algoritmo da Fortune nos permite encontrar rapidamente arestas infinitas: elas correspondem a meias arestas ainda no litoral no final do algoritmo.

Meu algoritmo de restrição recebe uma caixa como entrada e consiste em três etapas:

1. Ele fornece o posicionamento de cada vértice do diagrama dentro do retângulo.
2. Corte todas as arestas infinitas.
3. Fecha células.

O estágio 1 é trivial - basta expandir o retângulo se ele não contiver um vértice.

O estágio 2 também é bastante simples - consiste em calcular as interseções entre os raios e o retângulo.

O estágio 3 também não é muito complicado, apenas requer atenção. Eu o faço em duas etapas. Primeiro, adiciono os pontos de canto do retângulo às células nos vértices dos quais deveriam estar. Então, asseguro-me de que todos os vértices da célula estejam conectados por meia-arestas.

Eu recomendo que você estude o código e faça perguntas se precisar de detalhes sobre esta parte.

Aqui está o diagrama de saída do algoritmo delimitador:


Agora vemos que todas as arestas são desenhadas. E se você diminuir o zoom, verifique se todas as células estão fechadas:

Interseção com retângulo

Ótimo! Mas a primeira imagem do começo do artigo é melhor, certo?

Em muitas aplicações, é útil ter a interseção entre o diagrama de Voronoi e o retângulo, como mostrado na primeira imagem.

O bom é que, depois de restringir o gráfico, é muito mais fácil fazer isso. A má notícia é que, embora o algoritmo não seja muito complicado, precisamos ter cuidado.

A idéia é a seguinte: contornamos a meia aresta de cada célula e verificamos a interseção entre a meia aresta e o retângulo. São possíveis cinco casos:

1. A meia costela está completamente dentro do retângulo: economizamos uma meia costela
2. A meia costela está completamente fora do retângulo: descartamos essa meia costela
3. A meia costela sai do retângulo: nós truncamos a meia costela e a salvamos como a última meia costela que sai .
4. A meia costela vai para dentro do retângulo: nós truncamos a meia costela para conectá-la à última meia costela que saiu (nós a salvamos no caso 3 ou 5)
5. A meia costela cruza o retângulo duas vezes: truncamos a meia costela e adicionamos uma meia costela para conectá-la à última meia costela que saiu e, em seguida, salvamos como a nova última meia costela que saiu .

Sim, houve muitos casos. Criei uma imagem para mostrar a todos:


O polígono laranja é a célula original e o vermelho é a célula truncada. As meias costelas truncadas estão marcadas em vermelho. Costelas verdes foram adicionadas para conectar as meias costelas que entram no retângulo com as meias costelas saindo.

Aplicando esse algoritmo a um diagrama delimitado, obtemos o resultado esperado:

Conclusão

O artigo acabou sendo bastante longo. E tenho certeza de que muitos aspectos ainda não estão claros para você. No entanto, espero que seja útil para você. Examine o código para obter detalhes.

Para resumir e garantir que não perdemos tempo em vão, medi no meu laptop (barato) o tempo para calcular o diagrama de Voronoi para um número diferente de locais:

• $$$n = 1000$: 2 ms
• $$$n = $ 10.00$: 33 ms
• $$$n = $ 100.00$: 450 ms
• $$$n = $ 1.000.00$: 6600 ms

Não tenho nada para comparar esses indicadores, mas parece que é incrivelmente rápido!

Referências

• Artigo original de Stephen Fortune
• Geometria Computacional, 3ª edição de Mark de Berg, Otfried Cheong, Marc van Kreveld e Mark Overmars
• Algoritmo da sorte: uma explicação intuitiva em jacquesheunis.com
• Algoritmo e implementação da Fortune no blog.ivank.net
• Introdução aos Algoritmos, 3ª edição de Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest e Clifford Stein