Geometria Computacional em Python: Da Teoria à Aplicação

Quando as pessoas pensam em geometria computacional, na minha experiência, elas normalmente pensam em uma de duas coisas:

1. Uau, isso parece complicado.
2. Ah sim, casco convexo.

Neste post, gostaria de lançar alguma luz sobre geometria computacional, começando com uma breve visão geral do assunto antes de passar para alguns conselhos práticos com base em minhas próprias experiências (pule adiante se você tiver um bom domínio sobre o assunto).

Qual é o alvoroço?

Embora os algoritmos de geometria computacional de casco convexo sejam normalmente incluídos em um curso introdutório de algoritmos, a geometria computacional é um assunto muito mais rico que raramente recebe atenção suficiente do desenvolvedor/cientista da computação médio (a menos que você esteja fazendo jogos ou algo assim).

Teoricamente intrigante…

Do ponto de vista teórico, as questões em geometria computacional são muitas vezes extremamente interessantes; as respostas, convincentes; e os caminhos pelos quais eles são alcançados, variados. Essas qualidades por si só tornam um campo que vale a pena estudar, na minha opinião.

Por exemplo, considere o problema da galeria de arte: possuímos uma galeria de arte e queremos instalar câmeras de segurança para proteger nossas obras de arte. Mas estamos com um orçamento apertado, então queremos usar o menor número de câmeras possível. De quantas câmeras precisamos?

Quando traduzimos isso para notação geométrica computacional, a 'planta' da galeria é apenas um simples polígono. E com alguma graxa de cotovelo, podemos provar que n/3 câmeras são sempre suficientes para um polígono em n vértices, não importa o quão confuso seja. A prova em si usa gráficos duais, alguma teoria dos grafos, triangulações e muito mais.

Aqui, vemos uma técnica de prova inteligente e um resultado curioso o suficiente para ser apreciado por si só. Mas se a relevância teórica não for suficiente para você…

E importante na prática

Como mencionei anteriormente, o desenvolvimento de jogos depende muito da aplicação da geometria computacional (por exemplo, a detecção de colisões geralmente depende do cálculo do casco convexo de um conjunto de objetos); como os sistemas de informação geográfica (GIS), que são usados ​​para armazenar e realizar cálculos sobre dados geográficos; e robótica também (por exemplo, para problemas de visibilidade e planejamento).

Por que é tão difícil?

Vamos pegar um problema de geometria computacional bastante simples: dado um ponto e um polígono, o ponto está dentro do polígono? (Isso é chamado de problema do ponto no polígono ou PIP.)

O PIP faz um ótimo trabalho ao demonstrar por que a geometria computacional pode ser (enganosamente) difícil. Para o olho humano, esta não é uma pergunta difícil. Vemos o diagrama a seguir e é imediatamente óbvio para nós que o ponto está no polígono:

Mesmo para polígonos relativamente complicados, a resposta não nos escapa por mais de um ou dois segundos. Mas quando alimentamos esse problema em um computador, ele pode ver o seguinte:

 poly = Polygon([Point(0, 5), Point(1, 1), Point(3, 0), Point(7, 2), Point(7, 6), Point(2, 7)]) point = Point(5.5, 2.5) poly.contains(point)

O que é intuitivo para o cérebro humano não se traduz tão facilmente em linguagem de computador.

De forma mais abstrata (e ignorando a necessidade de representar essas coisas em código), os problemas que vemos nesta disciplina são muito difíceis de serem rigorosos ('tornar rigoroso') em um algoritmo de geometria computacional. Como descreveríamos o cenário do ponto no polígono sem usar uma linguagem tautológica como 'Um ponto está dentro de um polígono se estiver dentro do polígono'? Muitas dessas propriedades são tão fundamentais e tão básicas que é difícil defini-las concretamente.

Difícil, mas não impossível. Por exemplo, você pode rigorizar o ponto-em-polígono com as seguintes definições:

• Um ponto está dentro de um polígono se qualquer raio infinito que começa no ponto cruza com um número ímpar de arestas do polígono (conhecida como regra par-ímpar).
• Um ponto está dentro de um polígono se tiver um número de enrolamento diferente de zero (definido como o número de vezes que a curva que define o polígono percorre o ponto).

A menos que você tenha alguma experiência com geometria computacional, essas definições provavelmente não farão parte do seu vocabulário existente. E talvez isso seja emblemático de como a geometria computacional pode levá-lo a pensar de forma diferente .

Apresentando o CCW

Agora que entendemos a importância e a dificuldade dos problemas de geometria computacional, é hora de começarmos a trabalhar.

Na espinha dorsal do assunto está uma operação primitiva enganosamente poderosa: no sentido anti-horário, ou 'CCW' para abreviar. (Vou avisá-lo agora: CCW aparecerá de novo e de novo.)

CCW toma três pontos A, B e C como argumentos e pergunta: esses três pontos compõem um giro no sentido anti-horário (vs. um giro no sentido horário)? Em outras palavras, A -> B -> C é um ângulo anti-horário?

Por exemplo, os pontos verdes são CCW, enquanto os pontos vermelhos não são:

Por que CCW é importante

CCW nos dá uma operação primitiva na qual podemos construir. Isso nos dá um lugar para começar a rigorizar e resolver problemas de geometria computacional.

Para lhe dar uma noção de seu poder, vamos considerar dois exemplos.

Determinando a Convexidade

A primeira: dado um polígono, você pode determinar se é convexo? A convexidade é uma propriedade inestimável: saber que seus polígonos são convexos geralmente permite melhorar o desempenho em ordens de magnitude. Como um exemplo concreto: há um algoritmo PIP bastante simples que é executado em tempo Log(n) para polígonos convexos, mas falha para muitos polígonos côncavos.

Intuitivamente, essa lacuna faz sentido: formas convexas são 'legais', enquanto formas côncavas podem ter bordas afiadas que se projetam para dentro e para fora - elas simplesmente não seguem as mesmas regras.

Um algoritmo de geometria computacional simples (mas não óbvio) para determinar a convexidade é verificar se cada trio de vértices consecutivos é CCW. Isso leva apenas algumas linhas de código de geometria Python (supondo que os points sejam fornecidos no sentido anti-horário - se os points estiverem no sentido horário, você desejará que todos os trigêmeos sejam no sentido horário):

 class Polygon(object): ... def isConvex(self): for i in range(self.n): # Check every triplet of points A = self.points[i % self.n] B = self.points[(i + 1) % self.n] C = self.points[(i + 2) % self.n] if not ccw(A, B, C): return False return True

Tente isso no papel com alguns exemplos. Você pode até usar esse resultado para definir a convexidade. (Para tornar as coisas mais intuitivas, observe que uma curva CCW de A -> B -> C corresponde a um ângulo menor que 180, que é uma maneira amplamente ensinada de definir convexidade.)

Interseção de Linha

Como segundo exemplo, considere a interseção do segmento de linha, que também pode ser resolvida usando apenas o CCW:

 def intersect(a1, b1, a2, b2): """Returns True if line segments a1b1 and a2b2 intersect.""" return ccw(a1, b1, a2) != ccw(a1, b1, b2) and ccw(a2, b2, a1) != ccw(a2, b2, b1)

Por que este é o caso? A interseção do segmento de linha também pode ser formulada como: dado um segmento com extremidades A e B, as extremidades C e D de outro segmento estão do mesmo lado de AB? Em outras palavras, se as curvas de A -> B -> C e A -> B -> D estiverem na mesma direção, os segmentos não podem se cruzar. Quando usamos esse tipo de linguagem, fica claro que tal problema é o pão com manteiga da CCW.

Uma definição rigorosa

Agora que temos um gostinho da importância do CCW, vamos ver como ele é calculado. Dados os pontos A, B e C:

 def ccw(A, B, C): """Tests whether the turn formed by A, B, and C is ccw""" return (Bx - Ax) * (Cy - Ay) > (By - Ay) * (Cx - Ax)

Para entender de onde vem essa definição, considere os vetores AB e BC. Se tomarmos seu produto vetorial, AB x BC, este será um vetor ao longo do eixo z. Mas em que direção (ou seja, +z ou -z)? Como se vê, se o produto vetorial for positivo, o giro será no sentido anti-horário; caso contrário, é no sentido horário.

Essa definição parecerá pouco intuitiva, a menos que você tenha um bom entendimento de álgebra linear, regra da mão direita, etc. Mas é por isso que temos abstração – quando você pensa em CCW, pense apenas em sua definição intuitiva em vez de em sua computação. O valor será imediatamente limpo.

Meu mergulho em geometria computacional e programação usando Python

No último mês, tenho trabalhado na implementação de vários algoritmos de geometria computacional em Python. Como vou desenhá-los nas próximas seções, vou dedicar um segundo para descrever meus aplicativos de geometria computacional, que podem ser encontrados no GitHub.

Nota: Minha experiência é reconhecidamente limitada. Como tenho trabalhado nessas coisas por meses em vez de anos, siga meu conselho com um grão de sal. Dito isso, aprendi muito nesses poucos meses, então espero que essas dicas sejam úteis.

Algoritmo de Kirkpatrick

No centro do meu trabalho estava uma implementação do Algoritmo de Kirkpatrick para localização de pontos. A declaração do problema seria algo como: dada uma subdivisão planar (um monte de polígonos não sobrepostos no plano) e um ponto P, qual polígono contém P? Pense em esteróides point-in-polygon – em vez de um único polígono, você tem um avião cheio deles.

Como um caso de uso, considere uma página da web. Quando um usuário clica no mouse, a página da Web precisa descobrir o que o usuário clicou o mais rápido possível. Foi o botão A? Seria o Link B? A página da web é composta de polígonos não sobrepostos, então o Algoritmo de Kirkpatrick estaria bem posicionado para ajudar.

Embora eu não discuta o algoritmo em profundidade, você pode aprender mais aqui.

Triângulo Delimitador Mínimo

Como uma subtarefa, também implementei o algoritmo de O'Rourke para calcular um triângulo delimitador/limitador mínimo (ou seja, encontrar o menor triângulo que encerra um conjunto de pontos convexos) em tempo linear.

Nota: O cálculo do triângulo delimitador mínimo não ajuda nem prejudica o desempenho assintótico do Algoritmo de Kirkpatrick, pois o cálculo em si é de tempo linear, mas é útil para fins estéticos.

Conselhos práticos, aplicações e preocupações

As seções anteriores focaram em por que a geometria computacional pode ser difícil de raciocinar com rigor.

Na prática, temos que lidar com toda uma nova série de preocupações.

Lembra do CCW? Como uma boa continuação, vamos ver mais uma de suas grandes qualidades: ela nos protege contra os perigos dos erros de ponto flutuante.

Erros de ponto flutuante: por que o CCW é rei

No meu curso de geometria computacional, Bernard Chazelle, um estimado professor que publicou mais artigos do que posso contar, estabeleceu como regra que não poderíamos mencionar ângulos ao tentar descrever um algoritmo ou uma solução.

Por quê? Os ângulos são confusos. Os ângulos são “sujos”. Quando você tem que calcular um ângulo, você precisa dividir, ou usar alguma aproximação (qualquer coisa envolvendo Pi, por exemplo) ou alguma função trigonométrica.

Quando você tiver que calcular um ângulo no código , quase sempre estará aproximando. Você estará errado por um pequeno grau de precisão de ponto flutuante - o que importa quando você está testando a igualdade. Você pode resolver para algum ponto no plano através de dois métodos diferentes e, é claro, esperar que p1.x == p2.x and p1.y == p2.y . Mas, na realidade, essa verificação falhará com frequência . Além disso (e obviamente), esses pontos terão hashes diferentes.

Para piorar as coisas, seu grau de erro aumentará à medida que suas pequenas diferenças se propagam através de seus cálculos. (Para alguns exemplos mais científicos, este artigo analisa o que pode dar errado ao calcular o casco convexo ou a triangulação de Delaunay.)

Então, o que nós podemos fazer sobre isso?

quase igual

Parte do problema de geometria computacional do Python é que estamos exigindo exatidão em um mundo onde as coisas raramente são exatas. Isso se tornará um problema com mais frequência do que ao lidar com ângulos. Considere o seguinte:

 # Define two random points p1 = RandomPoint() p2 = RandomPoint() # Take the line through them l1 = Line(p1, p2) # Shift both points up by sqrt(2) p1.y += sqrt(2) p2.y += sqrt(2) l2 = Line(p1, p2) # Slope 'should' be the same? if abs(l1.slope - l2.slope) > 0: print "Error!" # Error!

Na verdade, este código imprimirá “Erro!” aproximadamente 70% do tempo (empiricamente). Podemos abordar essa preocupação sendo um pouco mais brandos com nossa definição de igualdade; isto é, sacrificando um grau de precisão.

Uma abordagem que usei (e vi, por exemplo, em alguns módulos OpenCV) é definir dois números como iguais se eles diferirem apenas por algum pequeno valor epsilon. Em Python, você pode ter:

 def almostEqual(x, y, EPSILON=1e-5): return abs(x - y) < EPSILON class Point(object): ... def __eq__(self, that): return (almostEqual(self.x, that.x) and almostEqual(self.y, that.y))

Na prática, isso é muito útil. Raramente, se é que alguma vez, você calcularia dois pontos que diferem por menos de 1e-5 que na verdade deveriam ser pontos diferentes. Eu recomendo implementar esse tipo de substituição. Métodos semelhantes podem ser usados ​​para linhas, por exemplo:

 class Line(object): ... def __eq__(self, that): return (almostEqual(self.slope, that.slope) and almostEqual(self.intercept, that.intercept))

Soluções mais avançadas foram propostas, é claro. Por exemplo, a escola de pensamento 'computação geométrica exata' (descrita neste artigo) visa fazer com que todos os caminhos de decisão em um programa dependam apenas do sinal de algum cálculo, em vez de seu valor numérico exato, eliminando muitas das preocupações relacionadas para cálculos de ponto flutuante. Nossa abordagem de quase igualdade apenas arranha a superfície, mas muitas vezes será suficiente na prática.

CCW é rei

Em um nível mais alto, é (discutivelmente) problemático que definimos nossas soluções em termos de quantidades computacionais exatas como ângulos ou coordenadas de pontos. Em vez de abordar apenas os sintomas (ou seja, eliminar erros de ponto flutuante com almostEqual ), por que não abordar a causa? A solução: em vez de pensar em termos de ângulos, pense em termos de CCW , o que ajudará a abstrair as preocupações associadas à computação de ponto flutuante.

Aqui está um exemplo concreto: digamos que você tenha algum polígono convexo P , um vértice v e algum ponto u fora do polígono. Como você pode descobrir se a linha uv intercepta P acima ou abaixo de v , ou não cruza em tempo constante?

A solução de força bruta (além de ser linear em tempo, em vez de constante) seria problemática, pois você teria que calcular alguns pontos de interseção de linha exatos.

Uma abordagem de tempo constante que vi envolve:

• Calculando alguns ângulos usando arctan2 .
• Convertendo esses ângulos em graus multiplicando por 180/Pi.
• Examinando as relações entre esses vários ângulos.

Felizmente, o autor usou a técnica almostEqual acima para suavizar os erros de ponto flutuante.

Na minha opinião, seria melhor evitar completamente o problema de erros de ponto flutuante. Se você levar alguns minutos para analisar o problema no papel, poderá obter uma solução totalmente baseada em CCW. A intuição: se os vértices adjacentes a v estão do mesmo lado de uv , então a linha não se cruza; caso contrário, veja se u e v estão do mesmo lado da linha entre os vértices adjacentes e, dependendo do resultado, compare suas alturas.

Aqui está o código Python para testar a interseção acima de v (a interseção abaixo apenas inverte a direção das comparações):

 def intersectsAbove(verts, v, u): """ Returns True if uv intersects the polygon defined by 'verts' above v. Assumes v is the index of a vertex in 'verts', and u is outside of the polygon. """ n = len(verts) # Test if two adjacent vertices are on same side of line (implies # tangency) if ccw(u, verts[v], verts[(v - 1) % n]) == ccw(u, verts[v], verts[(v + 1) % n]): return False # Test if u and v are on same side of line from adjacent # vertices if ccw(verts[(v - 1) % n], verts[(v + 1) % n], u) == ccw(verts[(v - 1) % n], verts[(v + 1) % n], verts[v]): return uy > verts[v].y else: return uy < verts[v].y

A solução não é imediatamente óbvia a olho nu, mas está na linguagem de um algoritmo de geometria computacional: 'mesmo lado da linha' é um elemento clássico desse algoritmo confiável.

Feito é melhor que perfeito

Na literatura geométrica computacional, muitas vezes há uma quantidade razoável de magia envolvida em operações aparentemente simples. Isso lhe dá uma escolha: você pode fazer as coisas da maneira mais difícil, seguindo algum artigo que define uma solução incrivelmente avançada para um problema não tão avançado – ou você pode fazer as coisas da maneira mais fácil com um pouco de força bruta.

Novamente, usarei um exemplo: amostragem de um ponto interno aleatório de um polígono arbitrário. Em outras palavras, eu lhe dou um polígono simples e você me dá um ponto aleatório dentro dele (distribuído uniformemente pelo polígono).

Muitas vezes, os pontos internos são necessários para o teste. Nesse caso, você não tem nenhum requisito de tempo de execução específico no algoritmo de geometria computacional que os produz (dentro do razoável). A solução rápida e suja, que leva ~2 minutos para ser implementada, seria escolher um ponto aleatório dentro de uma caixa contendo o polígono e ver se o próprio ponto está dentro do polígono.

Por exemplo, podemos errar duas vezes e encontrar uma amostra válida apenas no terceiro ponto:

Aqui está o código:

 class Polygon(object): ... def interiorPoint(self): """Returns a random point interior point""" min_x = min([px for p in self.points]) max_x = max([px for p in self.points]) min_y = min([py for p in self.points]) max_y = max([py for p in self.points]) def x(): return min_x + random() * (max_x - min_x) def y(): return min_y + random() * (max_y - min_y) p = Point(x(), y()) while not self.contains(p): p = Point(x(), y()) return p def contains(self, p): for i in range(self.n): p1 = self.points[i] p2 = self.points[(i + 1) % self.n] p3 = self.points[(i + 2) % self.n] if not ccw(p1, p2, p3): return False return True

Isso é conhecido como amostragem de rejeição : pegue pontos aleatórios até que um satisfaça seus critérios. Embora possa exigir várias amostras para encontrar um ponto que atenda aos seus critérios, na prática, a diferença será insignificante para seu conjunto de testes. Então, por que trabalhar mais? Em resumo: não tenha medo de seguir o caminho sujo quando a ocasião exigir.

A propósito: se você quer um algoritmo exato para amostragem aleatória, há um inteligente aqui que eu implementei abaixo. A essência disso:

1. Triangule seu polígono (ou seja, divida-o em triângulos).
2. Escolha um triângulo com probabilidade proporcional à sua área.
3. Pegue um ponto aleatório de dentro do triângulo escolhido (uma operação de tempo constante).

Observe que esse algoritmo exige que você triangule seu polígono, o que impõe imediatamente um limite de tempo de execução diferente ao algoritmo, bem como a necessidade de ter uma biblioteca para triangular polígonos arbitrários (usei poly2tri com ligações Python).

 from p2t import CDT class Triangle(object): ... def area(self): return abs((Bx * Ay - Ax * By) + (Cx * By - Bx * Cy) + (Ax * Cy - Cx * Ay)) / 2 def interiorPoint(self): r1 = random() r2 = random() # From http://www.cs.princeton.edu/~funk/tog02.pdf return (1 - sqrt(r1)) * A + sqrt(r1) * (1 - r2) * B + r2 * sqrt(r1) * C class Polygon(object): ... def triangulate(self): # Triangulate poly with hole cdt = CDT(poly.points) triangles = cdt.triangulate() def convert(t): A = Point(tax, tay) B = Point(tbx, tby) C = Point(tcx, tcy) return Triangle(A, B, C) return map(convert, triangles) def interiorPoint(self): # Triangulate polygon triangles = self.triangulate() areas = [t.area() for t in triangles] total = sum(areas) # Calculate normalized areas probabilities = [area / total for area in areas] weighted_triangles = zip(triangles, probabilities) # Sample triangles according to area r = random() count = 0 for (triangle, prob) in weighted_triangles: count += prob # Take random point from chosen triangle if count > r: return triangle.interiorPoint()

Felizmente, o esforço extra é evidente no código. Lembre-se: como dizem no Facebook, “feito é melhor que perfeito”. O mesmo vale para problemas de geometria computacional.

Testes visuais e automatizados

Como muitos dos problemas com os quais você trabalha em geometria computacional são definidos em termos de qualidades ou quantidades facilmente visualizáveis, o teste visual é particularmente importante - embora insuficiente por si só. O conjunto de testes ideal terá uma combinação de testes automatizados visuais e aleatórios.

Mais uma vez, procedemos pelo exemplo. Considere testar nossa implementação do Algoritmo de Kirkpatrick. Em uma etapa, o algoritmo precisa limitar o polígono fornecido por um triângulo e triangular a região entre o polígono e o triângulo externo. Aqui está um exemplo visual, onde a linha verde sólida define o polígono inicial e as linhas tracejadas definem a região triangulada:

Confirmar que esta triangulação foi executada corretamente é muito difícil de verificar através do código, mas é imediatamente evidente ao olho humano. Nota: Eu sugiro usar o Matplotlib para ajudar em seus testes visuais - há um bom guia aqui.

Mais tarde, vamos querer verificar se o algoritmo localiza os pontos corretamente. Uma abordagem aleatória e automatizada seria gerar vários pontos internos para cada polígono e garantir que retornamos o polígono desejado. Em código:

 class TestLocator(unittest.TestCase): ... def runLocator(self, polygons): # Pre-process regions l = Locator(polygons) # Ensure correctness for polygon in polygons: # Test 100 random interior points per region for k in range(100): target = polygon.interiorPoint() target_polygon = l.locate(target) self.assertEqual(polygon, target_polygon) self.assertTrue(target_polygon.contains(target))

Poderíamos então usar o método runLocator em diferentes conjuntos de polígonos, dando-nos um conjunto de testes bem diversificado.

Soluções de código aberto

A geometria computacional tem um bom conjunto de bibliotecas e soluções de código aberto disponíveis, independentemente da linguagem de programação de sua escolha (embora as bibliotecas C++ pareçam surgir em uma quantidade desproporcional).

Os benefícios de usar soluções de código aberto existentes (como a computação científica em Python) são bem conhecidos e foram discutidos extensivamente, então não vou continuar aqui. Mas pensei em mencionar alguns recursos centrados em Python que achei úteis:

• poly2tri: uma ótima biblioteca para triangulações rápidas de polígonos. Também suporta (e isso muitas vezes é crucial) polígonos com buracos neles. Escrito em C++, poly2tri também tem ligações Python e foi bem fácil de colocar em funcionamento. Veja meu método triangulate acima para um gosto pelas chamadas de função.
• scipy.spatial: inclui funções para calcular cascos convexos, Triangulações de Delaunay e muito mais. Rápido (como sempre), confiável, etc. Nota: Achei útil usar meu próprio tipo de dados Point com um método toNumpy : def np(self): return [self.x, self.y] . Então, eu poderia facilmente chamar métodos scipy.spatial, por exemplo: scipy.spatial.ConvexHull(np.array(map(lambda p: p.np()), points)) .
• OpenCV: a biblioteca de visão computacional de código aberto tem alguns bons módulos de geometria computacional autônomos. Em particular, usei sua função de triângulo de fechamento mínimo por um tempo antes de implementá-la.

Conclusão

Espero que este post tenha lhe dado um gostinho da beleza da geometria computacional como desenvolvedor Python, um assunto rico em problemas fascinantes e aplicações igualmente fascinantes.

Na prática, as implementações geométricas computacionais apresentam desafios únicos que o levarão a exercitar novas e empolgantes habilidades de resolução de problemas.

Se você estiver interessado em saber mais ou tiver alguma dúvida para mim, posso ser contatado em [email protected].