Triangolazione di Delaunay in 2D

Problema: Dato un insieme $\mathbf{P}$ di $n$ punti nel piano, trovare la triangolazione di Delaunay $\mathbf{D}(\mathbf{P})$ di $\mathbf{P}$. Trovare cioè quella triangolazione tale che nessun punto di $\mathbf{P}$ sia interno al cerchio circoscritto di ogni triangolo di $\mathbf{D}(\mathbf{P})$.

Soluzione: Convex Hull in 3D
Complessità: $O(n^2)$

La relazione tra il problema della triangolazione di Delaunay (in dimensione $d$) e quello dell'involucro convesso (in dimensione $d+1$) permette di sfruttare quanto visto in [3] per risolvere il problema con solo una decina di righe di codice.
L'idea principale, in sintesi, è quella di proiettare i punti di $\mathbf{P}$ sul paraboloide $z = x^2 + y^2$, completare l'involucro convesso di questi punti proiettati, scartare le facce superiori e prendere solo le facce inferiori dell'involucro. Queste ultime, di fatto, individuano e descrivono proprio i triangoli di $\mathbf{D}(\mathbf{P})$.

A tale scopo basta dimostrare che, per ogni triangolo inferiore dell'involucro convesso 3D, la sua proiezione sul piano $xy$ non contiene altri punti di $\mathbf{P}$. A livello intuitivo si immagini di prendere un piano che interseca 3 punti di un triangolo inferiore qualsiasi dell'involucro convesso 3D.

L'intersezione tra piano e paraboloide è un ellisse che verrà proiettata come una circonferenza sul piano $xy$. Si immagini ora di traslare il piano intersecante verticalmente verso il basso, fino renderlo tangente all'involucro convesso in un solo punto. Questo movimento individua così centro (la proiezione del punto tangente) e raggio (quantità traslata) della circonferenza proiettata sul piano $xy$ e che passa per i 3 punti di $\mathbf{P}$ corrispondenti alle rispettive proiezioni sul paraboloide. Considerando solo le facce inferiori, tutti gli altri punti di $\mathbf{P}$ proiettati sul paraboloide si troveranno al di sopra del piano intersecante i 3 punti e quindi fuori dal raggio della circonferenza che verrà proiettata sul piano $xy$ e che quindi non conterrà altri punti di $\mathbf{P}$. Per una dimostrazione formale e completa si veda [1] o [2].
Tuttavia è interessante notare che i vertici dei triangoli adesso sono in ordine inverso rispetto a quanto definito per l'algoritmo dell'involucro convesso in 3D. Questo perché $\mathbf{D}(\mathbf{P})$ è la proiezione sul piano $xy$ della parte inferiore dell'involucro convesso appena descritto. E' come se si tagliasse l'involucro in due parti (una superiore e l'altra inferiore) e si scartocciasse la parte inferiore applicandola sul piano $xy$, guardandola di fatto dall'interno e non più dall'esterno.

Per quanto riguarda l'implementazione in codice C/C++ le uniche modifiche al codice visto in [3] riguardano l''inizializzazione dei vertici (in cui la coordinata $z$ è ora definita come $z = x^2 + y^2$) e l'aggiunta dei metodi lower_faces e norm_z. Il primo si occupa di trovare tutte le facce inferiori dell'involucro convesso e lo fa attraverso il secondo metodo che si occupa di calcolare la coordinata $z$ della normale alla faccia.
Per quanto detto in [4], infatti, la coordinata $z$ di un vettore non è altro che la lunghezza, con segno, della sua proiezione sul versore $\mathbf{k}$. Si ha cioè che $z = \pm\mathbf{|proj_k v|} = \mathbf{|v|}cos\theta = \mathbf{v} \cdot \mathbf{k}$. Se $\frac{\pi}{2} < \theta < -\frac{\pi}{2}$ allora $z < 0$.

struct Face
{
 Edge *edge[3];
 Vertex *vertex[3];
 Face *prev;
 Face *next;
 bool visible;
 bool lower;
};
 
int main(int argc, char* argv[], char* env[])
{
 Vertex_Coord coord[] = {
 {31.0, -76.0}, {-13.0, 21.0}, {-63.0, -83.0}, {-5.0, -66.0}, {87.0, -94.0},
 {40.0, 71.0}, {23.0, -46.0}, {64.0, -80.0}, {0.0, -57.0}, {-14.0, 2.0} };
 int numv = sizeof(coord) / (sizeof(float) * 3);
 
 Vertex *vlist = nullptr;
 Edge *elist = nullptr;
 Face *flist = nullptr;
 
 for (int i = 0; i < numv; i++)
 {
  Vertex *tmpv = new Vertex;
  tmpv->vnum = i;
  tmpv->pos[0] = coord[i][0];
  tmpv->pos[1] = coord[i][1];
  float z = tmpv->pos[0] * tmpv->pos[0] + tmpv->pos[1] * tmpv->pos[1]; // z = x^2 + y^2
  tmpv->pos[2] = z;
  tmpv->duplicate = nullptr;
  tmpv->onhull = false;
  tmpv->mark = false;
  tmpv->next = nullptr;
  tmpv->prev = nullptr;
 
  list_insert<Vertex>(vlist, tmpv);
 }
 
 double_triangle(vlist, elist, flist);
 construct_hull(vlist, elist, flist);
 lower_faces(flist);
 print_faces(flist);
 
 return 0;
}
 
void lower_faces(Face* flist)
{
 Face *f = flist;
 
 do {
  if (norm_z(f) < 0.0f)
   f->lower = true;
  else
   f->lower = false;
 
  f = f->next;
 } while (f != flist);
}
 
float norm_z(Face *f)
{
 Vertex *a, *b, *c;
 
 a = f->vertex[0];
 b = f->vertex[1];
 c = f->vertex[2];
 
 return (b->pos[0] - a->pos[0]) * (c->pos[1] - a->pos[1]) - (b->pos[1] - a->pos[1]) * (c->pos[0] - a->pos[0]);
}
 
void print_faces(Face *flist)
{
 Face *f = flist;
 std::cout << "Faces of Delaunay Triangulation:\t(vertex indices)\n";
 std::cout << "v0\tv1\tv2\n";
 do {
  if (f->lower)
   std::cout << f->vertex[0]->vnum << "\t" << f->vertex[1]->vnum << "\t" << f->vertex[2]->vnum << std::endl;
 
  f = f->next;
 } while (f != flist);
}

Riferimenti
[1] Computational Geometry in C - O'Rourke
[2] http://www.cs.wustl.edu/~pless/506/l17.html
[3] Convex Hull in 3D
[4] Collinearità e Complanarità in 3D