Geometria computazionale in Python: dalla teoria all'applicazione

Quando le persone pensano alla geometria computazionale, secondo la mia esperienza, in genere pensano a una di queste due cose:

1. Wow, sembra complicato.
2. Oh sì, scafo convesso.

In questo post, vorrei fare un po' di luce sulla geometria computazionale, iniziando con una breve panoramica dell'argomento prima di passare a qualche consiglio pratico basato sulle mie esperienze (salta avanti se hai una buona padronanza dell'argomento).

Di cosa si tratta?

Mentre gli algoritmi di geometria computazionale dello scafo convesso sono tipicamente inclusi in un corso introduttivo di algoritmi, la geometria computazionale è un argomento molto più ricco che raramente riceve sufficiente attenzione dallo sviluppatore/scienziato informatico medio (a meno che tu non stia creando giochi o qualcosa del genere).

Teoricamente intrigante...

Da un punto di vista teorico, le questioni di geometria computazionale sono spesso estremamente interessanti; le risposte, convincenti; e le vie per cui sono raggiunte, variegate. Queste qualità da sole ne fanno un campo che vale la pena studiare, secondo me.

Ad esempio, considera il problema della galleria d'arte: possediamo una galleria d'arte e desideriamo installare telecamere di sicurezza per proteggere le nostre opere d'arte. Ma abbiamo un budget limitato, quindi vogliamo utilizzare il minor numero di fotocamere possibile. Di quante telecamere abbiamo bisogno?

Quando traduciamo questo in notazione geometrica computazionale, la "pianta" della galleria è solo un semplice poligono. E con un po' di olio di gomito, possiamo dimostrare che n/3 telecamere sono sempre sufficienti per un poligono su n vertici, non importa quanto sia disordinato. La dimostrazione stessa utilizza grafi doppi, teoria dei grafi, triangolazioni e altro.

Qui vediamo una tecnica di prova intelligente e un risultato abbastanza curioso da essere apprezzato da solo. Ma se la rilevanza teorica non ti basta...

E importante nella pratica

Come accennato in precedenza, lo sviluppo del gioco si basa fortemente sull'applicazione della geometria computazionale (ad esempio, il rilevamento delle collisioni spesso si basa sul calcolo dello scafo convesso di un insieme di oggetti); così come i sistemi di informazione geografica (GIS), che vengono utilizzati per memorizzare ed eseguire calcoli su dati geografici; e anche robotica (ad es. per problemi di visibilità e pianificazione).

Perché è così difficile?

Prendiamo un problema di geometria computazionale abbastanza semplice: dati un punto e un poligono, il punto si trova all'interno del poligono? (Questo è chiamato problema del punto nel poligono o PIP.)

PIP fa un ottimo lavoro nel dimostrare perché la geometria computazionale può essere (ingannevolmente) difficile. Per l'occhio umano, questa non è una domanda difficile. Vediamo il diagramma seguente ed è subito evidente che il punto è nel poligono:

Anche per poligoni relativamente complicati, la risposta non ci sfugge per più di un secondo o due. Ma quando inviamo questo problema a un computer, potrebbe vedere quanto segue:

 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)

Ciò che è intuitivo per il cervello umano non si traduce così facilmente nel linguaggio del computer.

Più astrattamente (e ignorando la necessità di rappresentare queste cose in codice), i problemi che vediamo in questa disciplina sono molto difficili da rigorare ("rendere rigorosi") in un algoritmo di geometria computazionale. Come descriveremmo lo scenario del punto nel poligono senza usare un linguaggio tautologico come "Un punto è all'interno di un poligono se è all'interno del poligono"? Molte di queste proprietà sono così fondamentali e così basilari che è difficile definirle concretamente.

Difficile, ma non impossibile. Ad esempio, potresti rigorizzare un punto nel poligono con le seguenti definizioni:

• Un punto è all'interno di un poligono se un raggio infinito che inizia nel punto si interseca con un numero dispari di bordi del poligono (nota come regola pari-dispari).
• Un punto è all'interno di un poligono se ha un numero di avvolgimento diverso da zero (definito come il numero di volte che la curva che definisce il poligono percorre il punto).

A meno che tu non abbia una certa esperienza con la geometria computazionale, queste definizioni probabilmente non faranno parte del tuo vocabolario esistente. E forse questo è emblematico di come la geometria computazionale possa spingerti a pensare in modo diverso .

Presentazione in senso antiorario

Ora che abbiamo un'idea dell'importanza e della difficoltà dei problemi di geometria computazionale, è tempo di bagnarci le mani.

Alla spina dorsale del soggetto c'è un'operazione primitiva ingannevolmente potente: in senso antiorario, o 'CCW' in breve. (Ti avverto ora: CCW apparirà ancora e ancora.)

CCW prende tre punti A, B e C come argomenti e chiede: questi tre punti compongono un giro in senso antiorario (contro un giro in senso orario)? In altre parole, A -> B -> C è un angolo antiorario?

Ad esempio, i punti verdi sono CCW, mentre i punti rossi non sono:

Perché CCW è importante

CCW ci fornisce un'operazione primitiva su cui possiamo costruire. Ci dà un punto di partenza per rigorizzare e risolvere problemi di geometria computazionale.

Per darti un'idea del suo potere, consideriamo due esempi.

Determinazione della convessità

Il primo: dato un poligono, puoi determinare se è convesso? La convessità è una proprietà inestimabile: sapere che i tuoi poligoni sono convessi spesso ti consente di migliorare le prestazioni di ordini di grandezza. Come esempio concreto: esiste un algoritmo PIP abbastanza semplice che viene eseguito in tempo Log(n) per poligoni convessi, ma fallisce per molti poligoni concavi.

Intuitivamente, questo divario ha senso: le forme convesse sono "belle", mentre le forme concave possono avere spigoli vivi che sporgono dentro e fuori, semplicemente non seguono le stesse regole.

Un semplice (ma non ovvio) algoritmo di geometria computazionale per determinare la convessità consiste nel verificare che ogni tripletta di vertici consecutivi sia CCW. Questo richiede solo poche righe di codice della geometria Python (supponendo che i points siano forniti in ordine antiorario, se i points sono in senso orario, vorrai che tutte le terzine siano in senso orario):

 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

Prova questo su carta con alcuni esempi. Puoi anche usare questo risultato per definire la convessità. (Per rendere le cose più intuitive, nota che una curva CCW da A -> B -> C corrisponde a un angolo inferiore a 180, che è un modo ampiamente insegnato per definire la convessità.)

Intersezione di linea

Come secondo esempio, considera l'intersezione del segmento di linea, che può anche essere risolta usando solo 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)

Perché è così? L'intersezione del segmento di linea può anche essere formulata come: dato un segmento con le estremità A e B, le estremità C e D di un altro segmento giacciono sullo stesso lato di AB? In altre parole, se le spire da A -> B -> C e A -> B -> D sono nella stessa direzione, i segmenti non possono intersecarsi. Quando usiamo questo tipo di linguaggio, diventa chiaro che un tale problema è il pane quotidiano di CCW.

Una definizione rigorosa

Ora che abbiamo un assaggio dell'importanza del CCW, vediamo come viene calcolato. Dati i punti 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)

Per capire da dove deriva questa definizione, consideriamo i vettori AB e BC. Se prendiamo il loro prodotto incrociato, AB x BC, questo sarà un vettore lungo l'asse z. Ma in quale direzione (cioè, +z o -z)? A quanto pare, se il prodotto incrociato è positivo, il giro è in senso antiorario; altrimenti è in senso orario.

Questa definizione sembrerà non intuitiva a meno che tu non abbia una buona comprensione dell'algebra lineare, della regola della mano destra, ecc. Ma ecco perché abbiamo l'astrazione: quando pensi in senso antiorario, pensa solo alla sua definizione intuitiva piuttosto che al suo calcolo. Il valore sarà immediatamente chiaro.

Il mio tuffo nella geometria computazionale e nella programmazione usando Python

Nell'ultimo mese, ho lavorato all'implementazione di diversi algoritmi di geometria computazionale in Python. Dato che li attingerò nelle prossime sezioni, mi prenderò un secondo per descrivere le mie applicazioni di geometria computazionale, che possono essere trovate su GitHub.

Nota: la mia esperienza è certamente limitata. Dato che lavoro su questa roba da mesi piuttosto che da anni, prendi il mio consiglio con le pinze. Detto questo, ho imparato molto in quei pochi mesi, quindi spero che questi suggerimenti si rivelino utili.

Algoritmo di Kirkpatrick

Al centro del mio lavoro c'era l'implementazione dell'algoritmo di Kirkpatrick per la localizzazione dei punti. L'affermazione del problema sarebbe qualcosa del tipo: data una suddivisione planare (un gruppo di poligoni non sovrapposti nel piano) e un punto P, quale poligono contiene P? Pensa al punto nel poligono con gli steroidi: invece di un singolo poligono, ne hai un piano pieno.

Come caso d'uso, considera una pagina web. Quando un utente fa clic con il mouse, la pagina Web deve capire su cosa l'utente ha fatto clic il più rapidamente possibile. Era il pulsante A? Era il collegamento B? La pagina web è composta da poligoni non sovrapposti, quindi l'algoritmo di Kirkpatrick sarebbe ben posizionato per dare una mano.

Anche se non discuterò in modo approfondito l'algoritmo, puoi saperne di più qui.

Triangolo di delimitazione minimo

Come attività secondaria, ho anche implementato l'algoritmo di O'Rourke per calcolare un triangolo di chiusura/limite minimo (ovvero, trovare il triangolo più piccolo che racchiude un insieme convesso di punti) in tempo lineare.

Nota: il calcolo del triangolo di delimitazione minimo non aiuta o danneggia le prestazioni asintotiche dell'algoritmo di Kirkpatrick poiché il calcolo stesso è in tempo lineare, ma è utile per scopi estetici.

Consigli pratici, applicazioni e preoccupazioni

Le sezioni precedenti si sono concentrate sul motivo per cui può essere difficile ragionare rigorosamente sulla geometria computazionale.

In pratica, dobbiamo affrontare tutta una nuova serie di preoccupazioni.

Ricordi CCW? Come bel passaggio, vediamo ancora un'altra delle sue grandi qualità: ci protegge dai pericoli degli errori in virgola mobile.

Errori in virgola mobile: perché CCW è il re

Nel mio corso di geometria computazionale, Bernard Chazelle, uno stimato professore che ha pubblicato più articoli di quanti ne possa contare, ha stabilito che non si possono menzionare gli angoli quando si tenta di descrivere un algoritmo o una soluzione.

Come mai? Gli angoli sono disordinati. Gli angoli sono "sporchi". Quando devi calcolare un angolo, devi dividere o usare qualche approssimazione (qualsiasi cosa che coinvolga Pi, per esempio) o qualche funzione trigonometrica.

Quando devi calcolare un angolo in code , sarai quasi sempre approssimativo. Sarai fuori da un piccolo grado di precisione in virgola mobile, che conta quando stai testando l'uguaglianza. Puoi risolvere un punto del piano con due metodi diversi e, ovviamente, aspettarti che p1.x == p2.x and p1.y == p2.y . Ma, in realtà, questo controllo fallirà spesso . Inoltre (e ovviamente), questi punti avranno hash diversi.

A peggiorare le cose, il tuo grado di errore aumenterà man mano che le tue piccole differenze si propagano attraverso i tuoi calcoli. (Per alcuni esempi più scientifici, questo articolo esamina cosa può andare storto quando si calcola lo scafo convesso o la triangolazione di Delaunay.)

Allora, cosa possiamo fare al riguardo?

quasi uguale

Parte del problema della geometria computazionale di Python è che richiediamo esattezza in un mondo in cui le cose sono raramente esatte. Questo diventerà un problema più spesso rispetto a quando si maneggiano gli angoli. Considera quanto segue:

 # 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!

In effetti, questo codice stamperà "Errore!" circa il 70% delle volte (empiricamente). Possiamo affrontare questa preoccupazione essendo leggermente più indulgenti con la nostra definizione di uguaglianza; cioè sacrificando un certo grado di precisione.

Un approccio che ho usato (e visto, ad esempio, in alcuni moduli OpenCV) è definire due numeri uguali se differiscono solo per un piccolo valore epsilon. In Python, potresti avere:

 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))

In pratica, questo è molto utile. Raramente, se mai, calcoleresti due punti che differiscono di meno di 1e-5 che in realtà sono intesi come punti diversi. Consiglio vivamente di implementare questo tipo di override. Metodi simili possono essere utilizzati per le linee, ad esempio:

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

Naturalmente sono state proposte soluzioni più avanzate. Ad esempio, la scuola di pensiero del "calcolo geometrico esatto" (descritta in questo articolo) mira a fare in modo che tutti i percorsi decisionali in un programma dipendano esclusivamente dal segno di un calcolo, piuttosto che dal suo esatto valore numerico, eliminando molte delle preoccupazioni relative ai calcoli in virgola mobile. Il nostro approccio di quasi uguaglianza graffia solo la superficie, ma nella pratica sarà spesso sufficiente.

CCW è il re

A un livello superiore, è (probabilmente) problematico definire le nostre soluzioni in termini di quantità computazionali esatte come angoli o coordinate puntiformi. Piuttosto che affrontare i sintomi da soli (ad esempio, lavare gli errori in virgola mobile con almostEqual ), perché non affrontare la causa? La soluzione: invece di pensare in termini di angoli, pensare in termini di CCW , che aiuterà ad astrarre le preoccupazioni associate al calcolo in virgola mobile.

Ecco un esempio concreto: supponiamo di avere un poligono convesso P , un vertice v e un punto u al di fuori del poligono. Come puoi capire se la linea uv interseca P sopra o sotto v , o per niente, in tempo costante?

La soluzione della forza bruta (oltre ad essere un tempo lineare, piuttosto che costante) sarebbe problematica in quanto dovresti calcolare alcuni punti di intersezione esatti della linea.

Un approccio a tempo costante che ho visto comporta:

• Calcolo di alcuni angoli usando arctan2 .
• Conversione di questi angoli in gradi moltiplicando per 180/Pi.
• Esaminare le relazioni tra questi vari angoli.

Fortunatamente, l'autore ha utilizzato la tecnica almostEqual sopra per appianare gli errori in virgola mobile.

A mio parere, sarebbe meglio evitare del tutto il problema degli errori in virgola mobile. Se ti dedichi qualche minuto per esaminare il problema sulla carta, puoi ottenere una soluzione basata interamente su CCW. L'intuizione: se i vertici adiacenti a v sono dalla stessa parte di uv , allora la retta non si interseca; altrimenti, controlla se u e v sono sullo stesso lato della linea tra i vertici adiacenti e, a seconda del risultato, confronta le loro altezze.

Ecco il codice Python per testare l'intersezione sopra v (l'intersezione sotto inverte semplicemente la direzione dei confronti):

 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

La soluzione non è immediatamente ovvia ad occhio nudo, ma è nel linguaggio di un algoritmo di geometria computazionale: "lo stesso lato della linea" è un elemento classico di quel fidato algoritmo.

Fatto e 'meglio che perfetto

Nella letteratura geometrica computazionale, c'è spesso una discreta quantità di magia coinvolta in operazioni apparentemente semplici. Questo ti dà una scelta: puoi fare le cose nel modo più difficile, seguendo alcuni documenti che definiscono una soluzione incredibilmente avanzata a un problema non così avanzato, oppure puoi fare le cose nel modo più semplice con un po' di forza bruta.

Ancora una volta, userò un esempio: campionare un punto interno casuale da un poligono arbitrario. In altre parole, ti do un semplice poligono e tu mi dai un punto casuale al suo interno (distribuito uniformemente sul poligono).

Spesso sono necessari punti interni per il test. In tal caso, non hai requisiti di runtime specifici sull'algoritmo di geometria computazionale che li produce (entro limiti ragionevoli). La soluzione rapida e sporca, che richiede circa 2 minuti per l'implementazione, sarebbe quella di scegliere un punto casuale all'interno di una casella contenente il poligono e vedere se il punto stesso si trova all'interno del poligono.

Ad esempio, potremmo sbagliare due volte e trovare un campione valido solo sul terzo punto:

Ecco il codice:

 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

Questo è noto come campionamento del rifiuto : prendi punti casuali fino a quando uno soddisfa i tuoi criteri. Sebbene possano essere necessari diversi campioni per trovare un punto che soddisfi i tuoi criteri, in pratica la differenza sarà trascurabile per la tua suite di test. Allora perché lavorare di più? In sintesi: non abbiate paura di prendere la strada sporca quando l'occasione lo richiede.

A proposito: se vuoi un algoritmo esatto per il campionamento casuale, ce n'è uno intelligente qui che ho implementato di seguito. Il succo di esso:

1. Triangola il tuo poligono (cioè, spezzalo in triangoli).
2. Scegli un triangolo con probabilità proporzionale alla sua area.
3. Prendi un punto casuale all'interno del triangolo scelto (un'operazione a tempo costante).

Nota che questo algoritmo richiede di triangolare il tuo poligono, il che impone immediatamente un diverso limite di runtime all'algoritmo, così come la necessità di avere una libreria per triangolare poligoni arbitrari (ho usato poly2tri con i collegamenti 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()

Si spera che lo sforzo extra sia evidente dal codice. Ricorda: come si dice su Facebook, "fatto è meglio che perfetto". Lo stesso vale per i problemi di geometria computazionale.

Test visivi e automatizzati

Poiché molti dei problemi su cui lavori nella geometria computazionale sono definiti in termini di qualità o quantità facilmente visualizzabili, il test visivo è particolarmente importante, sebbene di per sé insufficiente. La suite di test ideale avrà una combinazione di test automatici visivi e randomizzati.

Ancora una volta, procediamo con l'esempio. Prendi in considerazione la possibilità di testare la nostra implementazione dell'algoritmo di Kirkpatrick. Ad un certo punto, l'algoritmo deve delimitare il poligono dato da un triangolo e triangolare la regione tra il poligono e il triangolo esterno. Ecco un esempio visivo, in cui la linea verde continua definisce il poligono iniziale e le linee tratteggiate definiscono la regione triangolare:

Confermare che questa triangolazione sia stata eseguita correttamente è molto difficile da verificare attraverso il codice, ma è immediatamente evidente all'occhio umano. Nota: consiglio vivamente di utilizzare Matplotlib per aiutare nei test visivi: qui c'è una bella guida.

Successivamente, vorremo verificare che l'algoritmo localizzi correttamente i punti. Un approccio randomizzato e automatizzato consiste nel generare un gruppo di punti interni per ogni poligono e assicurarsi di restituire il poligono desiderato. Nel codice:

 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))

Potremmo quindi utilizzare il metodo runLocator su diversi insiemi di poligoni, fornendoci una suite di test ben diversificata.

Soluzioni open source

La geometria computazionale ha una bella suite di librerie e soluzioni open source disponibili indipendentemente dal linguaggio di programmazione scelto (sebbene le librerie C++ sembrino crescere in modo sproporzionato).

I vantaggi dell'utilizzo di soluzioni open source esistenti (come con il calcolo scientifico in Python) sono ben noti e sono stati ampiamente discussi, quindi non ne parlerò qui. Ma ho pensato di menzionare alcune risorse incentrate su Python che ho trovato utili:

• poly2tri: una grande libreria per triangolazioni veloci di poligoni. Supporta anche (e questo è spesso cruciale) poligoni con fori al loro interno. Scritto in C++, poly2tri ha anche collegamenti Python ed è stato abbastanza facile da avviare e funzionare. Vedi il mio metodo triangulate sopra per un assaggio delle chiamate di funzione.
• scipy.spatial: include funzioni per il calcolo di scafi convessi, triangolazioni Delaunay e altro. Veloce (come sempre), affidabile, ecc. Nota: ho trovato utile usare il mio tipo di dati Point con un metodo toNumpy : def np(self): return [self.x, self.y] . Quindi, potrei facilmente chiamare i metodi scipy.spatial, ad esempio: scipy.spatial.ConvexHull(np.array(map(lambda p: p.np()), points)) .
• OpenCV: la libreria di computer vision open-source ha dei bei moduli di geometria computazionale autonomi. In particolare, ho usato per un po' la sua funzione di triangolo che racchiude il minimo prima di implementarla da solo.

Conclusione

Spero che questo post ti abbia dato un assaggio della bellezza della geometria computazionale come sviluppatore Python, un argomento ricco di problemi affascinanti e applicazioni altrettanto affascinanti.

In pratica, le implementazioni geometriche computazionali presentano sfide uniche che ti spingeranno a esercitare nuove ed entusiasmanti capacità di risoluzione dei problemi.

Se sei interessato a saperne di più o hai domande da farmi, posso essere contattato a [email protected].