Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Como determinar um Diagrama de Voronoi usando o método da força bruta?

+1 voto
62 visitas
perguntada Abr 27 em Ciência da Computação por Felipe Amaral (16 pontos)  

Determinar as células do Diagrama de Voronoi em um plano euclidiano para um conjunto arbitrário de pontos usando um algoritmo computacional baseado diretamente em sua definição. Apresentar uma implementação deste algoritmo em python e dar exemplos práticos da importância deste diagrama para a solução de problemas do mundo real.

Compartilhe

1 Resposta

+1 voto
respondida Abr 27 por Felipe Amaral (16 pontos)  

Definição do Diagrama de Voronoi

Um Diagrama de Voronoi é uma subdivisão do espaço em células. Cada célula possui um ponto em seu centro e é definida de tal forma que a distância de qualquer ponto arbitrário localizado em seu interior está mais próximo de seu ponto central do que de que quaisquer outros pontos centrais das outras células.

Aplicações práticas

Os Diagramas de Voronoi possuem várias aplicações práticas em diversas áreas do conhecimento, como as ciências naturais, biológicas e sociais. Exemplos de aplicações podem ser vistos aqui.

Um problema clássico que ilustra a utilidade dos Diagramas de Voronoi é o da determinação das regiões de influência das agências de correios espalhadas em uma cidade, supondo que cada cliente irá optar por escolher a agência mais próxima geograficamente de sua residência.

Um caso famoso em relação à aplicação dos Diagramas de Voronoi se deu na Inglaterra em 1854, onde um médico usou tal diagrama para ilustrar como a maioria das pessoas que morreram no surto de cólera viviam mais perto de uma bomba de água infectada do que qualquer outra bomba.

Um algoritmo simples (e pouco eficiente) para a determinação do Diagrama de Voronoi

Dado um conjunto de pontos aleatórios em um plano bidimensional e valendo-se da definição do Diagrama de Voronoi podemos implementar um algoritmo para a determinação de cada uma de suas células de acordo com a seguinte lógica:

  1. Para o ponto central da célula p(xi,yi), determinar a reta perpendicular que passa pelo ponto médio e a reta que liga os pontos p(xi,yi) e p(xj,yj) para cada j diferente de i.
  2. Calcular todos os pontos de interseção entre todas as retas calculadas no passo anterior e também entre as retas que passam pelos pontos do retângulo que contém todos os pontos iniciais, isto é, p(xmin,ymin), p(xmax,ymin), p(xmax,ymax), p(xmin,ymax).
  3. Eliminar todos os pontos de interseção que se encontram do lado oposto entre as retas e o ponto central em análise p(xi,yi). Os pontos remanescentes constituirão os vértices de um polígono convexo.
  4. Ordenar todos os pontos do polígono em uma lista de tal forma que as linhas que ligam os pontos consecutivos da lista não cruzem com nenhuma outra linha. Tal ordenação pode ser alcançada usando um algoritmo para encontrar um fecho convexo (convex hull algorithm).

Assim, para determinar todo o Diagrama é só realizar os passos 1 a 4 para cada um dos pontos centrais.

Cabe observar que existem outros algoritmos mais eficientes para o cálculo do Diagrama de Voronoi, como o Algoritmo de Fortune. Este site apresenta uma boa visualização da evolução dos passos intermediários deste algoritmo até a construção final do Diagrama de Voronoi.

Uma implementação do algoritmo em Python

O código abaixo escrito em Python pode ser utilizado para calcular o Diagrama de Voronoi para um número arbitrário de pontos:

from scipy.spatial import ConvexHull
import matplotlib.pyplot as plt
from random import *

#=================================================================================
#Funções auxiliares
#=================================================================================

def RetaEntreDoisPontos(p1,p2):
    #Retorna os coeficientes a,b da inreta y=ax+b de igual distância entre dois pontos
    x1=p1[0]
    y1=p1[1]
    x2=p2[0]
    y2=p2[1]
    am=(x2-x1)/(y1-y2)
    bm=((y1+y2)*(y1-y2)+(x1+x2)*(x1-x2))/(2*(y1-y2))
    ans=(am,bm)
    return ans

def PontoIntersecaoDuasRetas(reta1,reta2):
    #Retorna o ponto p(x,y) originário da interseção de duas retas 
    #Os inputs das retas são do tipo (a,b)=> y=ax+b  ou (k)=> x=k
    ans=[]    
    if(isinstance(reta1, tuple) and isinstance(reta2, tuple)):
        a1=reta1[0]
        b1=reta1[1]
        a2=reta2[0]
        b2=reta2[1]
        if (a1-a2)!=0:
            x=(b2-b1)/(a1-a2)
            y=a1*x+b1
            ans=(x,y)
    if(isinstance(reta1, tuple) and not isinstance(reta2, tuple)):   
      a1=reta1[0]
      b1=reta1[1]
      x=reta2
      y=a1*x+b1
      ans=(x,y)
    if(not isinstance(reta1, tuple) and isinstance(reta2, tuple)):   
      a2=reta2[0]
      b2=reta2[1]
      x=reta1
      y=a2*x+b2
      ans=(x,y)    
    return ans

def RetasEntreVariosPontos(p,VetorDeOutrosPontos):
    #Retorna as retas que separam com igual distância o ponto p de cada um dos 
    #pontos no vetor VetorDeOutrosPontos
    ans=[]
    for i in VetorDeOutrosPontos:
        reta=RetaEntreDoisPontos(p,i)
        ans.append(reta)
    return ans

def PosicionamentoPontoEmRelacaoAReta(p,reta):
    #Retorna 0 se o ponto está na reta, 1 se o ponto p está a direita da reta y=-ax+b e -1 se está à esquerda
    ans=-1
    tol=10**(-6)
    if(isinstance(reta, tuple)):
        xp=p[0]
        yp=p[1]
        a=reta[0]
        b=reta[1]
        yr=a*xp+b
        if abs(yp-yr)<=tol:
            ans=0
        if abs(yp-yr)>tol and yp>yr:
            ans=1
    if(not isinstance(reta, tuple)):
        xp=p[0]
        xr=reta
        if abs(xp-xr)<=tol:
          ans=0
        if abs(xp-xr)>tol and xp>xr:
            ans=1
    return ans

def PontoIntersecaoEntreVariasRetas(ListaDeRetas):
    #Retorna uma lista com os pontos de interseção de várias retas
    ans=[]
    for i in range(0,len(ListaDeRetas)):
       for j in range(0,len(ListaDeRetas)):
            if i<j:
                ponto=PontoIntersecaoDuasRetas(ListaDeRetas[i],ListaDeRetas[j])
                if len(ponto) != 0:
                    ans.append(ponto)
    return ans

def EliminaPontosNoLadoOpostoDasRetas(p,ListaDeRetas,ListaDePontos):
    #Dado um ponto específico, elimina da ListaDePontos aqueles pontos
    #que estão no lado oposto da reta correspondente
    ans=ListaDePontos.copy()
    PontosEliminados=[]
    for i in range(0,len(ListaDeRetas)):
        for j in ans:
            Posicionamento1=PosicionamentoPontoEmRelacaoAReta(p,ListaDeRetas[i])
            Posicionamento2=PosicionamentoPontoEmRelacaoAReta(j,ListaDeRetas[i])
            if (Posicionamento1==1 and Posicionamento2==-1) or (Posicionamento1==-1 and Posicionamento2==1):
               PontosEliminados.append(j)
    for i in PontosEliminados:
        try:
           ans.remove(i)
        except ValueError:
            pass  # do nothing!

    return ans

def OrdenaListaDePontosDoPoligono(ListaDePontos):
    #Retorna uma lista de pontos ordenada do polígono em um sentido anti-horário
    if len(ListaDePontos)<3:
        ans=ListaDePontos
    if len(ListaDePontos)>2:
        ch=ConvexHull(ListaDePontos)
        ans=[]
        for i in range(0,len(ch.vertices)):
            ans.append(ListaDePontos[ch.vertices[i]])
        ans.append(ans[0])
    return ans

def CriaListaDePontos(NumeroDePontos,xmin,xmax,ymin,ymax):
    ans=[]
    for i in range(0,NumeroDePontos):
        x=xmin+(xmax-xmin)*random()
        y=xmin+(xmax-xmin)*random()
        ans.append((x,y))
    return ans

#=================================================================================
#Chamada principal
#=================================================================================
#Parâmetros iniciais
numerodepontos=5
xmin=0
xmax=10
ymin=0
ymax=10
ListaDePontos=CriaListaDePontos(numerodepontos,xmin,xmax,ymin,ymax)  

#Gráfico com os pontos iniciais
for i in range(0,len(ListaDePontos)):
    x=ListaDePontos[i][0]
    y=ListaDePontos[i][3]
    plt.scatter(x,y,color='black')
    plt.text(x+0.5, y+0.5, str(i))
axes = plt.gca()
axes.set_xlim([xmin,xmax])
axes.set_ylim([ymin,ymax])

#Cálculo da célula do ponto j (retas auxiliares em verde, célula de Voronoi em vermelho)
j=1
p=ListaDePontos[j]
VetorDeOutrosPontos =ListaDePontos.copy()
VetorDeOutrosPontos.pop(j)
ListaDeRetas=RetasEntreVariosPontos(p,VetorDeOutrosPontos)
for i in range(0,len(ListaDeRetas)):
    plt.plot([xmin,xmax],[xmin*ListaDeRetas[i][0]+ListaDeRetas[i][4],xmax*ListaDeRetas[i][0]+ListaDeRetas[i][5]],'g--')    
ListaDeRetas.append((0,ymin))
ListaDeRetas.append((0,ymax))
ListaDeRetas.append(xmin)
ListaDeRetas.append(xmax)
ListaDePontosdasIntersecoes=PontoIntersecaoEntreVariasRetas(ListaDeRetas)
ListaDePontosRemanecentes=EliminaPontosNoLadoOpostoDasRetas(p,ListaDeRetas,ListaDePontosdasIntersecoes)
ListaDePontosPoligono=OrdenaListaDePontosDoPoligono(ListaDePontosRemanecentes)
for i in range(0, len(ListaDePontosPoligono)-1):
        x1=ListaDePontosPoligono[i][0]
        y1=ListaDePontosPoligono[i][6]
        x2=ListaDePontosPoligono[i+1][0]
        y2=ListaDePontosPoligono[i+1][7]
        plt.plot([x1,x2], [y1,y2], 'r-')
axes = plt.gca()
axes.set_xlim([xmin,xmax])
axes.set_ylim([ymin,ymax])

#Cálculo do Diagrama Completo, isto é, todos os pontos j
for i in range(0,len(ListaDePontos)):
    x=ListaDePontos[i][0]
    y=ListaDePontos[i][8]
    plt.scatter(x,y,color='black')
    plt.text(x+0.5, y+0.5, str(i))
axes = plt.gca()
axes.set_xlim([xmin,xmax])
axes.set_ylim([ymin,ymax])
for j in range(0,len(ListaDePontos)):
    p=ListaDePontos[j]
    VetorDeOutrosPontos =ListaDePontos.copy()
    VetorDeOutrosPontos.pop(j)
    ListaDeRetas=RetasEntreVariosPontos(p,VetorDeOutrosPontos)
    #for i in range(0,len(ListaDeRetas)):
    #    plt.plot([xmin,xmax],[xmin*ListaDeRetas[i][0]+ListaDeRetas[i][9],xmax*ListaDeRetas[i][0]+ListaDeRetas[i][10]],'g--')    
    ListaDeRetas.append((0,ymin))
    ListaDeRetas.append((0,ymax))
    ListaDeRetas.append(xmin)
    ListaDeRetas.append(xmax)
    ListaDePontosdasIntersecoes=PontoIntersecaoEntreVariasRetas(ListaDeRetas)
    ListaDePontosRemanecentes=EliminaPontosNoLadoOpostoDasRetas(p,ListaDeRetas,ListaDePontosdasIntersecoes)
    ListaDePontosPoligono=OrdenaListaDePontosDoPoligono(ListaDePontosRemanecentes)
    for i in range(0, len(ListaDePontosPoligono)-1):
            x1=ListaDePontosPoligono[i][0]
            y1=ListaDePontosPoligono[i][11]
            x2=ListaDePontosPoligono[i+1][0]
            y2=ListaDePontosPoligono[i+1][12]
            plt.plot([x1,x2], [y1,y2], 'r-')
axes = plt.gca()
axes.set_xlim([xmin,xmax])
axes.set_ylim([ymin,ymax])

A figura abaixo foi gerada utilizando o código acima. As células do Diagrama de Voronoi correspondem aos polígonos dados pelos segmentos em vermelho. Um destaque é dado para a célula associada ao ponto 1, que foi formada pelas interseções das retas intermediárias (em verde), calculadas conforme o passo 1 descrito acima.

A imagem será apresentada aqui.

...