mercoledì 16 agosto 2017

Convex Hull in 2D

Problema: Dato un insieme di $n$ punti nel piano, trovare il più piccolo involucro (o inviluppo) convesso che contiene gli $n$ punti. Trovare cioè i punti, tra gli $n$ dell'insieme dato, che si trovano sui lati del più piccolo poligono convesso che contiene tutti gli altri.

Soluzione: Graham Scan
Complessità: $O(n\log n)$



La prima cosa da fare (metodo find_lowest; si veda il codice a corredo dell'articolo) è quella di trovare il punto con ordinata più piccola e metterlo all'inizio dell'array che contiene le coordinate dei punti (in caso di più punti con stessa ordinata più piccola scegliere quello con ascissa più grande).

La seconda (metodi qsort, della libreria standard del C, e comp_vertices) è quella di ordinare in ordine crescente gli altri punti dell'array in base all'angolo formato dall'asse delle $x$ e il segmento $V_{low} - V_i$, con $low$ indice del punto con ordinata più piccola trovato in precedenza e $0 \le i < n$ con $i \ne low$.




Per questo passaggio si rivela utile quanto visto in [3] e [4]. Infatti, se l'area formata da tre vertici $A, B, C$ è positiva allora $A-C$ (e quindi il vertice $C$) avrà angolo maggiore rispetto a $A-B$ (e quindi al vertice $B$). Se invece l'area è negativa $A-C$ (e quindi il vertice $C$) avrà angolo minore rispetto a $A-B$ (e quindi al vertice $B$). In questo modo si è definito un metodo di comparazione tra i punti dell'insieme dato in input.




Per i punti che giacciono su una stessa retta non si può stabilire una precedenza ma si sfrutterà questa situazione per eliminare il punto più vicino a $V_{low}$ (vedi $V_7$, $V_4$, $V_8$ e $V_2$ della prima immagine; solo $V_4$ e $V_2$ vengono mantenuti. $V_7$ e $V_8$ vengono marcati per essere eliminati. Punti come $V_0$ e $V_1$, invece, non possono essere marcati allo stesso modo in questa prima fase dell'algoritmo).

Dopo queste due operazioni preliminari si è certi che i primi due punti dell'array si troveranno sui bordi dell'involucro convesso (vedi $V_3$ e $V_4$). Prima di passare alla sua costruzione definitiva, però, è necessario ripulire l'array di tutti quei punti che sono stati marcati per essere eliminati (metodo clean_up).
A questo punto si inseriscono i primi due punti dell'array su uno stack e si itera sui restanti punti (metodo graham) inserendo sullo stack solo quelli che si trovano a sinistra dei primi due sullo stack. Sfruttando quanto visto in [5], infatti, si può controllare la costruzione dell'involucro convesso in maniera agevole scartando tutti quei punti che si trovano all'interno o sul bordo dello stesso. Ad esempio, alla prima iterazione, sullo stack ci saranno i vertici $V_3$ e $V_4$. $V_0$ si trova a sinistra del segmento $V_3 - V_4$ e quindi verrà messo in cima allo stack. $V_1$, invece, si trova a destra del segmento $V_4 - V_0$, segnalandoci che il vertice $V_0$ era, in realtà, all'interno dell'involucro convesso. Si può quindi procedere alla rimozione di $V_0$ dalla cima dello stack e procedere con la prossima iterazione. $V_1$ si trova a sinistra del segmento $V_3 - V_4$ (i primi due vertici rimasti sullo stack dopo la rimozione di $V_0$) e viene messo in cima allo stack. $V_6$ giace sulla stessa retta di $V_4 - V_1$, segnalandoci che il vertice $V_1$ era, di fatto, sul bordo dell'involucro convesso. $V_1$ viene rimosso dalla cima dello stack e l'iterazione successiva dice che $V_6$ si trova a sinistra di $V_3 - V_4$ e quindi viene messo sulla cima dello stack. $V_2$ è a sinistra di $V_4 - V_6$ e viene messo in cima allo stack. $V_5$ è a sinistra di $V_6 - V_2$ e viene messo sulla cima dello stack. Infine, $V_3$ è a sinistra di $V_2 - V_5$ e viene messo in cima allo stack.

/*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 <cmath>
#include <deque>
 
typedef float Vertex_Coord[2];
 
struct Vertex
{
 int vnum;
 Vertex_Coord pos;
 bool remove;
};
 
Vertex *vlist;
 
void print_stack(const std::deque<Vertex*>&);
void print_list(const Vertex*, const int);
void find_lowest(const int);
int comp_vertices(const void*, const void*);
int clean_up(const int);
void graham(std::deque<Vertex*>&, const int);
 
bool left(const Vertex_Coordconst Vertex_Coordconst Vertex_Coord);
float area2(const Vertex_Coordconst Vertex_Coordconst Vertex_Coord);
 
inline float max_3(float afloat bfloat c)
{
 return ((a > b) ? ((a > c) ? a : c) : (b > c) ? b : c);
}
 
int main(int argccharargv[], charenv[])
{
 Vertex_Coord poly_convex_coord[] = {
 {3.0, 3.0}, {3.0, 5.0}, {0.0, 1.0}, {2.0, 5.0}, {-2.0, 2.0},
 {-3.0, 2.0}, {6.0, 5.0}, {-3.0, 4.0}, {-5.0, 2.0}, {-5.0, -1.0},
 {1.0, -2.0}, {-3.0, -2.0}, {4.0, 2.0}, {5.0, 1.0}, {-5.0, 1.0},
 {3.0, -2.0}, {0.0, 5.0}, {0.0, 0.0}, {7.0, 4.0} };
 
 int numv = sizeof(poly_convex_coord) / (sizeof(float) * 2);
 vlist = new Vertex[numv];
 
 for (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].remove = false;
 }
 
 find_lowest(numv);
 qsort(&vlist[1], numv - 1, sizeof(Vertex), comp_vertices);
 numv = clean_up(numv);
 std::deque<Vertex*> vstack;
 graham(vstack, numv);
 print_stack(vstack);
 
 return 0;
}
 
void graham(std::deque<Vertex*> &vstackconst int n)
{
 int i = 2;
 
 vstack.push_front(&vlist[0]);
 vstack.push_front(&vlist[1]);
 
 while (i < n)
 {
  if (left(vstack[1]->pos, vstack[0]->pos, vlist[i].pos))
  {
   vstack.push_front(&vlist[i]);
   i++;
  }
  else
   vstack.pop_front();
 }
}
 
int clean_up(const int n)
{
 int j = 0;
 
 for (int i = 0; i < n; i++)
 {
  if (!(vlist[i].remove))
  {
   vlist[j] = vlist[i];
   j++;
  }
  i++;
 }
 return j;
}
 
int comp_vertices(const void *lhsconst void *rhs)
{
 Vertex *p1 = (Vertex*)lhs;
 Vertex *p2 = (Vertex*)rhs;
 int a = area2(vlist[0].pos, p1->pos, p2->pos);
 
 if (a > 0)
  return -1;
 else if (a < 0)
  return 1;
 else // Collineari
 {
  // Se vertice p2 è più lontano |p2 - vlist[0]| > |p1 - vlist[0]|
  // quindi proiezione della lunghezza di p2 su entrambi gli assi cartesiani è 
  // maggiore di quella di p1
  int x = fabs(p1->pos[0] - vlist[0].pos[0]) - fabs(p2->pos[0] - vlist[0].pos[0]);
  int y = fabs(p1->pos[1] - vlist[0].pos[1]) - fabs(p2->pos[1] - vlist[0].pos[1]);
 
  // Se p2 è più lontano risultato sopra sarà negativo sia per x che per y
  if ((x < 0) || (y < 0))
  {
   p1->remove = true;
   return -1;
  }
  else if ((x > 0) || (y > 0)) // altrimenti p1 è più lontano
  {
   p2->remove = true;
   return 1;
  }
  else // se no, p1 e p2 sono coincidenti; eliminato quello con indice maggiore
  {
   if (p1->vnum > p2->vnum)
    p2->remove = true;
   else
    p1->remove = true;
 
   return 0;
  }
 }
}
 
void find_lowest(const int vnum)
{
 int m = 0;
 
 for (int i = 0; i < vnum; ++i)
  if (vlist[i].pos[1] < vlist[m].pos[1] || (vlist[i].pos[1] == vlist[m].pos[1] && vlist[i].pos[0] > vlist[m].pos[0]))
   m = i;
 
 Vertex t = vlist[0];
 vlist[0] = vlist[m];
 vlist[m] = t;
}
 
void print_stack(const std::deque<Vertex*> &s)
{
 size_t i = 0;
 
 for (std::deque<Vertex*>::const_iterator ci = s.begin(); ci != s.end(); ci++, i++)
  std::cout << (*ci)->vnum << ((i < (s.size() - 1)) ? ", " : "");
}
 
void print_list(const Vertex *vlistconst int n)
{
 for (int i = 0; i < n; i++)
 {
  std::cout << vlist[i].vnum << "(" << vlist[i].pos[0] << ", " << vlist[i].pos[1] << ")" << std::endl;
 }
}
 
bool left(const Vertex_Coord aconst Vertex_Coord bconst Vertex_Coord c)
{
 return area2(abc) > 0.0f;
}
 
float area2(const Vertex_Coord aconst Vertex_Coord bconst 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] Elementi di Geometria Computazionale
[3] Area di un poligono
[4] Collinearità di tre punti in 2D
[5] Verificare l'intersezione di due segmenti in 2D

Nessun commento:

Posta un commento