PHP: triangolazioni di Delaunay in 2D

Questa particolare tecnica di Geometria computazionale, che prende il nome dal russo Boris Nikolaevich Delone, francesizzato in Delaunay, consente di definire una griglia di triangoli in una superficie in 2D o in 3D, dato un insieme di punti P.
La griglia viene costruita in modo che, per ogni circonferenza circoscritta ad un triangolo, nessun punto di P (oltre a quelli che formano il triangolo stesso) giace all’interno della circonferenza. Questa condizione viene anche detta di cironferenza libera.

La libreria di funzioni, che propongo è la traduzione “letterale” in PHP delle sole funzioni per superfici in 2D, scritte in codice Visual Basic 6 da Franco Languasco, disponibili qui. Il pacchetto di Languasco è davvero molto ben fatto e contempla anche le funzioni per le triangolazioni in 3D, nonchè la visualizzazione grafica del risultato e la possibilità di ruotare i 3 assi. Il codice iniziale risale in realtà alle librerie GEOMPACK3 di Barry Joe originariamente scritte nel 1991 in linguaggio Fortran.

E’ bene precisare che avendo fatto una traduzione “letterale”, il mio codice PHP non risulta molto elegante, e soprattutto non sfrutta appieno le caratteristiche del linguaggio, sebbene già dai primi test, sembra che le prestazioni siano soddisfacenti!

Esistono diversi algoritmi per le triangolazioni di Delaunay, i più importanti sono:

Incrementale con complessità O(N²)
Dividi et impera con complessità O(NlogN)
Convex Hull con complessità O(NlogN)

L’algoritmo incrementale è quello utilizzato in questa libreria, sebbene non sia chiaramente il più efficiente, va più che bene per un insieme non troppo grande di punti.

Funzionamento:

L’algoritmo incrementale descritto da Dani Lischinski in questo documento, inizia con la costruzione di un primo triangolo, largo a sufficienza per contenere tutti i punti dati. I punti vengono aggiunti uno alla volta verificando che la condizione di circonferenza libera venga rispettata.
La prima operazione è quella di connettere i vertici del triangolo su cui giace il nuovo punto con il punto stesso.

algoritmo incrementale passo uno

Si vengono così a formare 3 nuovi triangoli. A questo punto i triangoli vengono verificati per la condizione di circonferenza libera.

algoritmo incrementale passo due

Se è rispettata, si procede con un nuovo punto altrimenti viene eseguita un’operazione che prende il nome di edge flip.

algoritmo incrementale passo tre

Questa operazione si basa sul concetto che la triangolazione di Delaunay massimizza l’angolo minimo.
Poiché in due triangoli che hanno un lato in comune, la somma degli angoli opposti a tale lato è sempre minore di 180°, sarà sufficiente scambiare le diagonali del quadrilatero formato dai due triangoli come mostrato nella figura sotto.

algoritmo incrementale passo tre

Nel codice della libreria di funzioni, ho volutamente lasciato tutti i commenti originali e ho commentato le istruzioni originali Visual Basic senza cancellarle, per poter controllare meglio eventuali errori di traduzione.
~~L’intera libreria è visibile qui~~. Questo, invece, è il codice del file di esempio:

<?php 
require_once("delaunay2D.php");
/*
================== PROVA =======================================================
001:   (4.111   0.668)
002:   (-4.209   -3.961)
003:   (-9.720   5.214)
004:   (4.181   -9.093)
005:   (7.252   5.810)
006:   (9.239   7.429)
007:   (8.991   -2.720)
Risultato:
001:   (2   1   3)
002:   (1   2   4)
003:   (1   4   7)
004:   (3   1   5)
005:   (5   1   7)
006:   (3   5   6)
007:   (6   5   7)
*/
//dati
 
  $labels = array("001","002","003","004","005","006","007");
  $XI =  array(4.111,-4.209,-9.720,4.181,7.252,9.239,8.991);
  $YI =  array(0.668,-3.961,5.214,-9.093,5.810,7.429,-2.720); 
 
if (count($XI) <> count($YI)){
  die("le coordinate delle X non sono in numero uguale alle Y");
}
else {
  $NP = count($XI); //N° di punti.  
}
 
echo "Numero punti: ".$NP."<br />";
 
$NT = 0; // N° di triangoli.
 
$IERR = 0;
// Inizializza vettori
$Vcl = array();
$Ind = array();
$TIL = array();
$TNBR = array();
 
    echo ('&lt;pre&gt;');
    for ($I = 1; $I <= $NP; $I++) {
        // Organizza le coordinate dei vertici come
        // richiesto dalle routines di delaunay2D:
        $Vcl[1][$I] = $XI[$I-1]; //Vettore con indice iniziale 1
        $Vcl[2][$I] = $YI[$I-1]; //Vettore con indice iniziale 1
        $Ind[$I] = $I;  // Triangolazione di tutti i vertici contenuti in VCL().
        echo ("Punto $I ".$labels[$I-1]." (".
              $Vcl[1][$I].",".$Vcl[2][$I].")<br />");
    }
    echo ('&lt;/pre&gt;');
 
    // Costruisce le triangolazioni
    $res = DTRIW2($NP, $Vcl, $Ind, $NT, $TIL, $TNBR, $IERR);
 
    echo ("<h2>N. triangoli: ".$NT."</h2>");
 
    if ($IERR == 0) {
      echo ('&lt;pre&gt;');
          for ($I = 1; $I <= $NT; $I++) {
              $Triang[1][$I] = $TIL[1][$I];
              $Triang[2][$I] = $TIL[2][$I];
              $Triang[3][$I] = $TIL[3][$I];
              echo ("Triangolo $I: (".$Triang[1][$I].",".$Triang[2][$I].
                    ",".$Triang[3][$I].")<br />");            
          }
      echo ('&lt;/pre&gt;');        
    }
    else {
      echo "Errore n.: $IERR  - ".zCodiciErrore($IERR)."<br />";
    }
?>

L’array $XI contiene le ascisse dei punti, mentre $YI le ordinate. Tutti i vettori devono avere indice iniziale 1. Le coordinate dei punti vengono passate mediante la matrice $Vcl e il vettore degli indici $Ind alla funzione DTRIW2(). L’output viene fornito passando per riferimento la matrice $TIL così come il vettore per i codici di errore $IERR.

Ho preparato anche questa piccola demo dove è anche possibile visualizzare graficamente i punti su una gif, e al passaggio del mouse, i triangoli che compongono la griglia. Per questo effetto ho utilizzato una map area che viene generata dallo script php e JQuery maphilight uno script che consente di evidenziare le forme geometriche definite nelle map area.

Download

~~In questo pacchetto è disponibile la libreria e il codice per gli esempi~~

Ringraziamenti

Ringrazio Franco Languasco per la sua ottima versione Visual Basic 6 senza la quale sarebbe stata dura! Federico Villa, Project Director di Target Tobrukper avermi dato lo spunto a scrivere questo post, e Alessandro Scoscia per le sue dritte geniali!

Conclusioni

Non ho trovato molto sulla generazione di griglie di superfici 2D utilizzando il metodo di Delaunay in PHP, questo potrebbe essere un inizio. Spero, in futuro di migliorare il “porting”, al momento piuttosto grezzo, di questa libreria. Infine, una raccomandazione: non ho potuto testare in modo esaustivo il codice, quindi prima di utilizzarla è bene fare un po’ di prove. Fatemi sapere, in caso di errori!

Riferimenti ed approfondimenti: