Soluzione: Ear clipping
Complessità: $O(n^2)$
Per triangolare un poligono la soluzione più semplice, ma anche più costosa, è quella di prendere tutti i segmenti che uniscono due vertici del poligono (ci sono $O(n^2)$ possibili candidati) e, per ogni segmento, controllare che sia interno al poligono e che non ci siano intersezioni con i lati del poligono (costo $O(n)$). Queste operazioni hanno quindi costo $O(n^3)$ e devono essere ripetute sino a quando non si trovano tutte le diagonali (che sono $n-3$). Il costo totale dunque arriva a $O(n^4)$. L'algoritmo è molto semplice ma con un piccolo sforzo in più si può ridurre il costo computazionale a $O(n^2)$ (in realtà si potrebbe arrivare fino a $O(n)$ ma il livello di complessità inizierebbe a farsi piuttosto elevato; per metodi più "semplici" con costo $O(n\log n)$ si vedano invece [2] e [3]).
L'idea iniziale è quella di ridurre il numero di possibili candidati a diagonali da $O(n^2)$ a $O(n)$, portando così il costo totale da $O(n^4)$ a $O(n^3)$. A questo scopo si può sfruttare il fatto che ogni poligono con almeno 4 vertici ha almeno 2 orecchie che non si sovrappongono. Per orecchie si intendono tre vertici consecutivi $a, b, c$ in cui $ac$ è la diagonale e $b$ è la punta dell'orecchio. La figura sopra rappresenta un pentagono con due orecchie che non si sovrappongono $v_4, v_0, v_1$ e $v_1, v_2, v_3$. Usando la ricorsione si divide il poligono in due parti: una già triangolata (l'orecchio) e l'altra da triangolare ricorsivamente allo stesso modo (trovando cioè l'orecchio del nuovo poligono diviso nella precedente iterazione).
Essendoci solo $O(n)$ diagonali prese in questo modo (che sono cioè basi di un orecchio) il costo totale passa, come già accennato, da $O(n^4)$ a $O(n^3)$. Prendendo come esempio la prima immagine (un pentagono), trovata la prima diagonale attraverso l'orecchio $v_4, v_0, v_1$ quello che resta è un quadrilatero con un'altra sola diagonale da trovare (ricorsivamente individuando un altro orecchio).
Per portare il costo a $O(n^2)$ è utile preliminarmente notare che l'eliminazione di un orecchio del poligono può incidere solo sullo stato dei vertici precedente e successivo alla punta dell'orecchio eliminato nell'iterazione corrente (in pratica i vertici della base dell'orecchio). Per stato del vertice si intende la proprietà che indica se il vertice è punta dell'orecchio o meno. Effettuando una operazione preliminare che verifica lo stato di tutti i vertici del poligono (effettuata una sola volta all'inizio e con costo $O(n^2)$) da lì in avanti, ad ogni iterazione, non resta che verificare lo stato dei soli vertici che formano la base/diagonale dell'orecchio eliminato (costo $O(n)$). Quest'ultima operazione naturalmente avviene solo se il vertice corrente è segnalato come punta dell'orecchio (non avviene per tutti i vertici) riducendo di fatto il costo totale portandolo a $O(n^2)$.
Infine, non resta che vedere la parte dell'algoritmo relativa alla verifica dell'intersezione delle diagonali. Bisogna prima di tutto prestare attenzione al fatto che una diagonale deve essere interna al poligono (l'algoritmo presentato finora, di fatto, non si è occupato di questo aspetto). E' utile a questo scopo dividere il problema in due parti: una che controlla l'intersezione di una diagonale con i lati del poligono che non si intersecano nei vertici della diagonale (sfruttando esclusivamente quanto visto in [4] per l'intersezione di due segmenti) e l'altra che prende in considerazione solo i lati che si intersecano nei vertici della diagonale (rigettando eventuali parallelismi tra la diagonale e questi lati). Per quanto riguarda quest'ultima parte, nella figura sotto si vede uno dei due possibili casi (il terzo caso, che gestisce tre vertici che giacciono sulla stessa retta, rientra nel primo). Sul lato sinistro della figura i vertici $v_{i-1}, v_i, v_{i+1}$ formano un angolo convesso. Per verificarlo basta notare che il vertice che precede la punta dell'orecchio ($v_{i-1}$ in questo caso) deve stare a sinistra di o sul segmento formato dagli altri due vertici ($v_i - v_{i+1}$ in questo caso). A questo punto, per verificare se la diagonale è all'interno del poligono, si può sfruttare il fatto che i vertici $v_{i-1}$ e $v_{i+1}$ devono stare sempre alla sinistra della diagonale se questa è presa prima in una direzione (ad esempio, da $v_i$ a $v_x$) e poi nell'altra (da $v_x$ a $v_i$). Nell'altro caso, lato destro della figura, i vertici formano un angolo concavo. Per verificare se la diagonale è all'interno del poligono si può sfruttare il fatto che i vertici $v_{i-1}$ e $v_{i+1}$ non devono stare entrambi alla sinistra della diagonale se questa è presa prima in una direzione e poi nell'altra (come descritto prima). In pratica è quasi la negazione della condizione precedente (quasi perché la semplice negazione del caso precedente permette anche di stare sulla diagonale).
La decisione di dividere quest'ultima parte in due operazioni distinte ha anche un beneficio dal punto di vista prestazionale. Entrambe le operazioni devono dare esito positivo perché il test sull'intersezione possa considerarsi superato. L'operazione che prende in considerazione solo i lati che si intersecano nei vertici della diagonale ha costo $O(1)$ mentre l'altra ha costo $O(n)$. In una espressione logica che verifica la veridicità di entrambe, mettere prima quella con costo $O(1)$ ha indubbi benefici, dato che l'altra non verrà valutata affatto nel caso la prima risulti falsa.
Nell'implementazione in codice C/C++ seguente il fulcro si trova nel metodo triangulate. In essa si possono notare due cicli innestati con il ciclo più esterno che si ferma quando il numero di vertici arriva a 3 (nel qual caso non si sono più orecchie) e con il ciclo più interno che chiama diagonal (costo $O(n)$) per i vertici precedente e successivo. Questo potrebbe far credere che il costo di triangulate sia $O(n^2)$ così da far salire il costo totale dell'intero algoritmo. In realtà, come già accennato, il flusso del codice entra nel ciclo più interno solo quando il vertice è marcato come orecchio, non per tutti i vertici indistintamente.
/*Per facilitare la stesura, la lettura e la comprensione del codice si è volutamente deciso di utilizzare un mix di C/C++ cercando di restare il più possibile all'interno del sottoinsieme C del C++ e di utilizzare solo quei costrutti e funzionalità di base del C++ che permettono una maggiore semplificazione.*/ #include <iostream> #include <cfloat> #include <cmath> typedef float Vertex_Coord[2]; struct vertex_node { int vnum; Vertex_Coord pos; bool ear; vertex_node *prev; vertex_node *next; }; void list_insert(vertex_node*&, vertex_node&); void list_delete(vertex_node*&, vertex_node&); void print_list(vertex_node*); void print_diagonal(vertex_node*, vertex_node*); bool diagonalie(vertex_node*, vertex_node*, vertex_node *); bool incone(vertex_node*, vertex_node*); bool diagonal(vertex_node*, vertex_node*, vertex_node *); void triangulate(vertex_node*&, unsigned int); bool left(const Vertex_Coord, const Vertex_Coord, const Vertex_Coord); bool left_on(const Vertex_Coord, const Vertex_Coord, const Vertex_Coord); bool collinear(const Vertex_Coord, const Vertex_Coord, const Vertex_Coord); bool AlmostEqual(float, float, float); bool proper_intersect(const Vertex_Coord, const Vertex_Coord, const Vertex_Coord, const Vertex_Coord); bool between(const Vertex_Coord, const Vertex_Coord, const Vertex_Coord); bool intersect(const Vertex_Coord, const Vertex_Coord, const Vertex_Coord, const Vertex_Coord); float area2(const Vertex_Coord, const Vertex_Coord, const Vertex_Coord); inline float max_3(float a, float b, float c) { return ((a > b) ? ((a > c) ? a : c) : (b > c) ? b : c); } int main(int argc, char* argv[], char* env[]) { Vertex_Coord poly_convex_coord[] = { {0.0, 0.0}, {10.0, 7.0}, {12.0, 3.0}, {20.0, 8.0}, {13.0, 17.0}, {10.0, 12.0}, {12.0, 14.0}, {14.0, 9.0}, {8.0, 10.0}, {6.0, 14.0}, {10.0, 15.0}, {7.0, 18.0}, {0.0, 16.0}, {1.0, 13.0}, {3.0, 15.0}, {5.0, 8.0}, {-2.0, 9.0}, {5.0, 5.0} }; int numv = sizeof(poly_convex_coord) / (sizeof(float) * 2); vertex_node *vlist = new vertex_node[numv]; // Matiene indirizzo originale di vlist per passarlo a delete[] dato che viene modificato // in triangulate e in metodi per cancellare nodi della lista vertex_node *tmp = vlist; for (unsigned int i = 0; i < numv; i++) { vlist[i].vnum = i; vlist[i].pos[0] = poly_convex_coord[i][0]; vlist[i].pos[1] = poly_convex_coord[i][1]; vlist[i].ear = false; vlist[i].next = nullptr; vlist[i].prev = nullptr; list_insert(vlist, vlist[i]); } //print_list(vlist); triangulate(vlist, numv); delete[] tmp; return 0; } void list_insert(vertex_node *&head, vertex_node &p) { if (head && head->next != nullptr) { p.next = head; p.prev = head->prev; head->prev = &p; p.prev->next = &p; } else { head = &p; head->next = head->prev = &p; } } void list_delete(vertex_node *&head, vertex_node &p) { if (head) { if (head == head->next) head = nullptr; else if (&p == head) head = head->next; p.next->prev = p.prev; p.prev->next = p.next; } } void print_list(vertex_node *head) { vertex_node *vp = head; int i; do { std::cout << "(" << vp->pos[0] << ", " << vp->pos[1] << ")" << std::endl; vp = vp->next; } while (vp != head); } void print_diagonal(vertex_node *a, vertex_node *b) { std::cout << "(" << a->vnum << ", " << b->vnum << ")" << std::endl; } bool diagonalie(vertex_node *a, vertex_node *b, vertex_node *vlist) { vertex_node *c, *c1; c = vlist; do { c1 = c->next; // controlla che non ci siano intersezioni tra diagonale ab e bordi del poligono; salta bordi incidenti a vertici a o b if ((c != a) && (c1 != a) && (c != b) && (c1 != b) && intersect(a->pos, b->pos, c->pos, c1->pos)) return false; c = c->next; } while (c != vlist); // vero solo se diagonale ab non è incidente con bordi del poligono; esclusi quelli incidenti a vertici a o b return true; } bool incone(vertex_node *a, vertex_node *b) { vertex_node *a0, *a1; a1 = a->next; a0 = a->prev; // se a è un vertice convesso if (left_on(a->pos, a1->pos, a0->pos)) return left(a->pos, b->pos, a0->pos) && left(b->pos, a->pos, a1->pos); // se a è un vertice concavo return !(left_on(a->pos, b->pos, a1->pos) && left_on(b->pos, a->pos, a0->pos)); } bool diagonal(vertex_node *a, vertex_node *b, vertex_node *vlist) { return incone(a, b) && incone(b, a) && diagonalie(a, b, vlist); } void ear_init(vertex_node *vlist) { vertex_node *v0, *v1, *v2; // inizializza ear per tutti i vertici v1 = vlist; do { v2 = v1->next; v0 = v1->prev; v1->ear = diagonal(v0, v2, vlist); v1 = v1->next; } while (v1 != vlist); } void triangulate(vertex_node *&vlist, unsigned int num_v) { vertex_node *v0, *v1, *v2, *v3, *v4; ear_init(vlist); while (num_v > 3) { v2 = vlist; do { if (v2->ear) { v3 = v2->next; v4 = v3->next; v1 = v2->prev; v0 = v1->prev; print_diagonal(v1, v3); v1->ear = diagonal(v0, v3, vlist); v3->ear = diagonal(v1, v4, vlist); v1->next = v3; v3->prev = v1; vlist = v3; num_v--; break; } v2 = v2->next; } while (v2 != vlist); } } bool intersect(const Vertex_Coord a, const Vertex_Coord b, const Vertex_Coord c, const Vertex_Coord d) { if (proper_intersect(a, b, c, d)) return true; else if (between(a, b, c) || between(a, b, d) || between(c, d, a) || between(c, d, b)) return true; else return false; return (left(a, b, c) xor left(a, b, d)) && (left(c, d, a) xor left(c, d, b)); } bool between(const Vertex_Coord a, const Vertex_Coord b, const Vertex_Coord c) { if (!collinear(a, b, c)) return false; if (a[0] != b[0]) return ((a[0] <= c[0]) && (c[0] <= b[0])) || ((a[0] >= c[0]) && (c[0] >= b[0])); else return ((a[1] <= c[1]) && (c[1] <= b[1])) || ((a[1] >= c[1]) && (c[1] >= b[1])); } bool proper_intersect(const Vertex_Coord a, const Vertex_Coord b, const Vertex_Coord c, const Vertex_Coord d) { if (collinear(a, b, c) || collinear(a, b, d) || collinear(c, d, a) || collinear(c, d, b)) return false; return (left(a, b, c) xor left(a, b, d)) && (left(c, d, a) xor left(c, d, b)); } bool AlmostEqual(float a, float b, float epsilon = FLT_EPSILON) { if (fabs(a - b) <= epsilon * max_3(1.0f, fabs(a), fabs(b))) return true; return false; } bool left(const Vertex_Coord a, const Vertex_Coord b, const Vertex_Coord c) { return area2(a, b, c) > 0.0f; } bool left_on(const Vertex_Coord a, const Vertex_Coord b, const Vertex_Coord c) { return area2(a, b, c) >= 0.0f; } bool collinear(const Vertex_Coord a, const Vertex_Coord b, const Vertex_Coord c) { return AlmostEqual(area2(a, b, c), 0.0f); } float area2(const Vertex_Coord a, const Vertex_Coord b, const Vertex_Coord c) { return (b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1]); }
Riferimenti
[1] Computational Geometry in C - O'Rourke
[2] Geometric Tools for Computer Graphics - Schneider, Eberly
[3] Computational Geometry, Algorithms and Applications - de Berg, et al.
[4] Verificare l'intersezione di due segmenti
Nessun commento:
Posta un commento