CGAL, Python, Octave

Having a good library for computational geometry is very important for me. A year ago I found CGAL and since then I am enjoying it. Sometimes the burden of C++ can be annoying (for rapid prototyping, for example), therefore lately I was trying the Python bindings and trying to get some functions working on GNU Octave ... English Version

Para mi es fundamental tener una buena bilioteca de geometría computacional a mano. Hace un año, laburando en un programita para seguir objetos en videos, me topé con CGAL y desde entonces es mi primera referencia a la hora de programar cuestiones geométricas (mallas, manipulacion de polígonos, etc...). CGAL está escrita en C++ y a veces programar es un poco tedioso debido al caracter declarativo del lenguaje. Por esta razón, últimamente estuve probando la versión para Python y estoy tratando de hacer funciones para GNU Octave (todavía sin éxito) ... versión en Español

Uno de los algoritmos que más me gustan (y que más he usado) de CGAL son las «formas alfa» (alpha shapes). Dada una nube de puntos ¿Cuál es el polígono (con vértices en puntos de la nube) que describe el borde (o contorno) de la nube? La Figura 1 ilustra el problema y una posible respuesta intuitiva (con el famoso método de unir los puntos con líneas).
Este algoritmo también se puede aplicar a nubes de puntos en 3D y CGAL trae una lista impresionante de funciones de manipulación de formas alfa. Sin embargo en 2D, obtener el polígono suele ser una pregunta fecuente en blogs y listas de mail. El siguiente código muestra una solución al problema en C++ (no pretendo ser eficiente! se agradecen correciones y mejorías):

//CGAL
#include < CGAL/basic.h >
#include < CGAL/Cartesian.h >
#include < CGAL/squared_distance_2.h >
#include < CGAL/Delaunay_triangulation_2.h >
#include < CGAL/Alpha_shape_2.h >
#include < CGAL/IO/io.h >
#include < CGAL/Regular_triangulation_2.h >
#include < CGAL/Regular_triangulation_vertex_base_2.h >
#include < CGAL/Regular_triangulation_face_base_2.h >
#include < CGAL/Triangulation_2.h >
#include < CGAL/Triangulation_face_base_2.h >
#include < CGAL/Triangulation_euclidean_traits_2.h >
#include < CGAL/Alpha_shape_vertex_base_2.h > 
#include < CGAL/Alpha_shape_face_base_2.h>
#include < CGAL/point_generators_2.h >
#include < CGAL/Kernel/function_objects.h >

// STL
#include < fstream >
#include < string >
#include < vector >

// Typedefstypedef double Coord_type;
typedef CGAL::Cartesian< coord_type > Rep;
typedef Rep::Point_2                Point_2;
typedef Rep::Vector_2               Vector_2;
typedef Rep::Segment_2              Segment;
typedef Rep::Line_2                 Line;
typedef Rep::Triangle_2             Triangle;
typedef Rep::Less_xy_2              Point_compare;
typedef CGAL::Triangulation_2< rep >  Triangulation;
typedef std::list< point_2 >          CGALPointlist;

//Delaunay triangulationtypedef Rep                                                Gt;
typedef CGAL::Alpha_shape_vertex_base_2< gt >                Vb;
typedef CGAL::Triangulation_face_base_2< gt >                Df;
typedef CGAL::Alpha_shape_face_base_2< Gt, Df >              Fb;
typedef CGAL::Triangulation_default_data_structure_2< Gt,Vb,Fb >  Tds;
typedef CGAL::Delaunay_triangulation_2< Gt,Tds >             Delaunay;

//Alpha shape and the types
typedef CGAL::Alpha_shape_2< delaunay >                      Alpha_shape;
typedef Alpha_shape::Face                                  Face;
typedef Alpha_shape::Vertex                                Vertex;
typedef Alpha_shape::Edge                                  Edge;
typedef Alpha_shape::Face_handle                           Face_handle;
typedef Alpha_shape::Vertex_handle                         Vertex_handle;
typedef Alpha_shape::Face_circulator                       Face_circulator;
typedef Alpha_shape::Vertex_circulator                     Vertex_circulator;
typedef Alpha_shape::Locate_type                           Locate_type;
typedef Alpha_shape::Face_iterator                         Face_iterator;
typedef Alpha_shape::Vertex_iterator                       Vertex_iterator;
typedef Alpha_shape::Edge_iterator                         Edge_iterator;
typedef Alpha_shape::Edge_circulator                       Edge_circulator;
typedef Alpha_shape::Alpha_iterator                        Alpha_iterator;

// Here define your favorite vector
typedef struct Vector{
    double x;
    double y;
    double z;
} vect;

using namespace std;
void alpha_shape(vector< vect > cloud, vector< vect > & vert, int & count)
{ 
Delaunay          tr1; 
CGALPointlist     L; 
Alpha_shape       A;  

//this is the alpha value for 
//the alpha_shape 
double alpha_index;
CGAL::set_ascii_mode(std::cout);

// We load the points of the cloud into a CGAL Point List
int N = cloud.size();
for (int i = 0; i < N; i++)
{
  Point_2 p(cloud[i].x,cloud[i].y);
  L.push_back(p);
}

//We make the alpha-shape here
A.make_alpha_shape(L.begin(), L.end());  

// This find the minimum value for alpha 
// such that there is only one connected component 
// that can have holes 
Alpha_iterator opt = A.find_optimal_alpha(1);
      cout << "Optimal alpha value to get one connected component is " <<*opt<< endl;
alpha_index=*opt;  
A.set_mode(Alpha_shape::REGULARIZED); 
A.set_alpha(alpha_index); 
Alpha_shape::Alpha_shape_edges_iterator it = A.alpha_shape_edges_begin();  
int Anumedges=1;    

// ****    
//What next is building the boundary in order    
// Is there a more efficient way of doing this?    
// ****        
// Count the number of edges. 

for(;it!=A.alpha_shape_edges_end();it++)
  Anumedges++; 

//Order the vertices  
Alpha_shape::Alpha_shape_edges_iterator  et = A.alpha_shape_edges_begin();  
it = A.alpha_shape_edges_begin();  
CGALPointlist vertices;      
int acount=0; 
vertices.push_back((A.segment(*et)).vertex(0)); 
while(acount < Anumedges)
 {
  while(A.segment(*et).vertex(1) != A.segment(*it).vertex(0))
   it++;
  vertices.push_back((A.segment(*it)).vertex(0));
  acount++;
  et = it;
  it = A.alpha_shape_edges_begin();
  }
 count = acount;
 int i = 0;
 vect temp;
 for(CGALPointlist::iterator ip=vertices.begin(); ip!=vertices.end();ip++)
 {
  temp.x = (*ip).x();
  temp.y=(*ip).y();
  vert.push_back(temp);
 }
 return;
}

Ok, no creo que haga falta decir que el codigo es grande, no? Uno puede poner los «include» y demas cosas en un archivo aparte, pero sigue siendo bastante molesto. Pero bueno, funciona muy bien.
El código que sigue hace lo mismo pero con Python. Como estoy compartiendo el código con un colega, lo escribí un poco mas estructurado

''' Border of a set of points using Alpha shapes '''
from __future__ import division
from CGAL.Kernel import *

from CGAL.Triangulations_2 import *
from CGAL.Alpha_shapes_2 import *

import matplotlib.pyplot as plt
from numpy import mod

def get_border(L, alpha=None,interactive=False):

  print 'Number of points in list ', len(L)
  A = Alpha_shape_2()
  A.make_alpha_shape(L)
  op = A.find_optimal_alpha(1)
  if alpha == None:
    alpha = op.next()
    print 'Optimal alpha value to get one connected component is ', \
    alpha
  else:
    print 'Using provided alpha value ',alpha
    print 'Optimal alpha value to get one connected component is ', \
    op.next()

  m = A.set_mode()
  A.set_mode(m.REGULARIZED)
  A.set_alpha(alpha)

  # Get the vertices ordered by edges incidence
  alpha_shape_vertices = []
  for it in A.alpha_shape_vertices:
    alpha_shape_vertices.append(it)
  
  print 'Vertices in alpha-shape ', len(alpha_shape_vertices)

  if interactive:
   while True:
      plt.ion()
      plt.figure(1)
      plt.clf()
      xo = [p.x() for p in L]
      yo = [p.y() for p in L]
      x = [p.point().x() for p in alpha_shape_vertices]
      y = [p.point().y() for p in alpha_shape_vertices]

      plt.plot(xo,yo,'or',x,y,'ob')

      plt.show()
      try:
        alpha = float(raw_input('New alpha value [Enter = finish]: '))
      except ValueError:
        break
          
      A.set_alpha(alpha)

      # Get the vertices ordered by edges incidence
      alpha_shape_vertices = []
      for it in A.alpha_shape_vertices:
          alpha_shape_vertices.append(it)
             
              
      print 'Vertices in alpha-shape ', len(alpha_shape_vertices)
      
  return A,alpha

def orderAlphaShape(A,L):

  print 'Ordering border. This may take a while...'
  # Get the vertices ordered by edges incidence
  alpha_shape_edges = []
  alpha_shape_vertices = []
  for it in A.alpha_shape_edges:
    alpha_shape_edges.append(A.segment(it))
  for it in A.alpha_shape_vertices:
    alpha_shape_vertices.append(it.point())

  nver = len(alpha_shape_vertices)
  incidency = list()
  ind = 0
  for edge in alpha_shape_edges:
    vert1 = edge.vertex(0)
    vert2 = edge.vertex(1)
    incidency.append([alpha_shape_vertices.index(vert1),alpha_shape_vertices.index(vert2)])
              
  border = list()
  borders = list()
  border_index = list()
  indeces = list()
  vert1 = [ i[0] for i in incidency]
  vert2 = [ i[1] for i in incidency]
  
  i = 0
  while len(border) < nver:
    vert = alpha_shape_vertices[vert1[i]]
    if vert not in border:
      border.append(vert)
      border_index.append(L.index(vert))

    vert1.pop(i)
    vert = vert2[i]
    vert2.pop(i)
    try:
      i = vert1.index(vert)
    except ValueError:
      print "Found inital vertex!", vert
      border.append(alpha_shape_vertices[vert])
      borders.append(border)
      indeces.append(border_index)
      border = list()
      border_index = list()
      if len(vert1) == 0:
        break
      else:
        i = 0

  return borders, indeces

def main(infile, outfile, alpha=None,interactive=None):
  #read the data
  L = list()
  with open(infile,'r') as f:
    for line in f:
      s = line.split()
      arg = []
      for n in s:
        if n != ' ':
            arg.append(float(n))
            
      L.append(Point_2(arg[0],arg[1]))
  f.close()
  
  # Get vertices in alpha shape
  [A,alpha] = get_border(L,alpha,interactive)

  # Get the ordered border and the indeces of the vertices
  # in the border in the original list
  [borders, indeces]=orderAlphaShape(A,L)
  # Different borders (like holes ar multiple components)
  # are sublists of borders
  
  #write output
  nborder=0
  with open(outfile,'w') as f:
    for border,index in zip(borders,indeces):
    
      if nborder > 0:
        info = 'NaN NaN NaN\n'
        f.write(info)
      nborder += 1
      
      for p,i in zip(border,index):
        info = str(i) + ' ' + str(p.x()) + ' ' + \
        str(p.y()) + '\n'
        f.write(info)

  f.close()

  return borders

Si bien parece la misma cantidad de código, la parte que ejecuta la forma alfa es mucho mas corta y simple. Luego viene todo un código que fui inventando para obtener los vertices ordenados de tal forma que sea facil hacer un dibujo del borde. El resultado de usar este algoritmo con un valor para alfa de 0.01 lo muestro en la Figura 2. Usando el valor óptimo, dado por el algoritmo, se ve en la Figura 3.

One of the features I have used the most from CGAL is the 2D Alpha-shapes. Alpha-shapes are the answer to the problem of finding the boundary of a cloud of points. In Figure 1 I show a cloud of points randomly distributed on the unit disk and a boundary that may have been drawn by a person by linking points with segments (but it wasn't! :D). The following C++ code (partly taken from the CGAL examples) does the work, but is a little bit long (I wonder how bad is my way of ordering the boundary! help here?)

//CGAL
#include < CGAL/basic.h >
#include < CGAL/Cartesian.h >
#include < CGAL/squared_distance_2.h >
#include < CGAL/Delaunay_triangulation_2.h >
#include < CGAL/Alpha_shape_2.h >
#include < CGAL/IO/io.h >
#include < CGAL/Regular_triangulation_2.h >
#include < CGAL/Regular_triangulation_vertex_base_2.h >
#include < CGAL/Regular_triangulation_face_base_2.h >
#include < CGAL/Triangulation_2.h >
#include < CGAL/Triangulation_face_base_2.h >
#include < CGAL/Triangulation_euclidean_traits_2.h >
#include < CGAL/Alpha_shape_vertex_base_2.h > 
#include < CGAL/Alpha_shape_face_base_2.h>
#include < CGAL/point_generators_2.h >
#include < CGAL/Kernel/function_objects.h >

// STL
#include < fstream >
#include < string >
#include < vector >

// Typedefstypedef double Coord_type;
typedef CGAL::Cartesian< coord_type > Rep;
typedef Rep::Point_2                Point_2;
typedef Rep::Vector_2               Vector_2;
typedef Rep::Segment_2              Segment;
typedef Rep::Line_2                 Line;
typedef Rep::Triangle_2             Triangle;
typedef Rep::Less_xy_2              Point_compare;
typedef CGAL::Triangulation_2< rep >  Triangulation;
typedef std::list< point_2 >          CGALPointlist;

//Delaunay triangulationtypedef Rep                                                Gt;
typedef CGAL::Alpha_shape_vertex_base_2< gt >                Vb;
typedef CGAL::Triangulation_face_base_2< gt >                Df;
typedef CGAL::Alpha_shape_face_base_2< Gt, Df >              Fb;
typedef CGAL::Triangulation_default_data_structure_2< Gt,Vb,Fb >  Tds;
typedef CGAL::Delaunay_triangulation_2< Gt,Tds >             Delaunay;

//Alpha shape and the types
typedef CGAL::Alpha_shape_2< delaunay >                      Alpha_shape;
typedef Alpha_shape::Face                                  Face;
typedef Alpha_shape::Vertex                                Vertex;
typedef Alpha_shape::Edge                                  Edge;
typedef Alpha_shape::Face_handle                           Face_handle;
typedef Alpha_shape::Vertex_handle                         Vertex_handle;
typedef Alpha_shape::Face_circulator                       Face_circulator;
typedef Alpha_shape::Vertex_circulator                     Vertex_circulator;
typedef Alpha_shape::Locate_type                           Locate_type;
typedef Alpha_shape::Face_iterator                         Face_iterator;
typedef Alpha_shape::Vertex_iterator                       Vertex_iterator;
typedef Alpha_shape::Edge_iterator                         Edge_iterator;
typedef Alpha_shape::Edge_circulator                       Edge_circulator;
typedef Alpha_shape::Alpha_iterator                        Alpha_iterator;

// Here define your favorite vector
typedef struct Vector{
    double x;
    double y;
    double z;
} vect;

using namespace std;
void alpha_shape(vector< vect > cloud, vector< vect > & vert, int & count)
{ 
Delaunay          tr1; 
CGALPointlist     L; 
Alpha_shape       A;  

//this is the alpha value for 
//the alpha_shape 
double alpha_index;
CGAL::set_ascii_mode(std::cout);

// We load the points of the cloud into a CGAL Point List
int N = cloud.size();
for (int i = 0; i < N; i++)
{
  Point_2 p(cloud[i].x,cloud[i].y);
  L.push_back(p);
}

//We make the alpha-shape here
A.make_alpha_shape(L.begin(), L.end());  

// This find the minimum value for alpha 
// such that there is only one connected component 
// that can have holes 
Alpha_iterator opt = A.find_optimal_alpha(1);
      cout << "Optimal alpha value to get one connected component is " <<*opt<< endl;
alpha_index=*opt;  
A.set_mode(Alpha_shape::REGULARIZED); 
A.set_alpha(alpha_index); 
Alpha_shape::Alpha_shape_edges_iterator it = A.alpha_shape_edges_begin();  
int Anumedges=1;    

// ****    
//What next is building the boundary in order    
// Is there a more efficient way of doing this?    
// ****        
// Count the number of edges. 

for(;it!=A.alpha_shape_edges_end();it++)
  Anumedges++; 

//Order the vertices  
Alpha_shape::Alpha_shape_edges_iterator  et = A.alpha_shape_edges_begin();  
it = A.alpha_shape_edges_begin();  
CGALPointlist vertices;      
int acount=0; 
vertices.push_back((A.segment(*et)).vertex(0)); 
while(acount < Anumedges)
 {
  while(A.segment(*et).vertex(1) != A.segment(*it).vertex(0))
   it++;
  vertices.push_back((A.segment(*it)).vertex(0));
  acount++;
  et = it;
  it = A.alpha_shape_edges_begin();
  }
 count = acount;
 int i = 0;
 vect temp;
 for(CGALPointlist::iterator ip=vertices.begin(); ip!=vertices.end();ip++)
 {
  temp.x = (*ip).x();
  temp.y=(*ip).y();
  vert.push_back(temp);
 }
 return;
}

Quite long, uh? But it works nicely. Now the next Python snippet does a similar job, but since I am sharing this code with colleagues, i did it a little bit more structured.

''' Border of a set of points using Alpha shapes '''
from __future__ import division
from CGAL.Kernel import *

from CGAL.Triangulations_2 import *
from CGAL.Alpha_shapes_2 import *

import matplotlib.pyplot as plt
from numpy import mod

def get_border(L, alpha=None,interactive=False):

  print 'Number of points in list ', len(L)
  A = Alpha_shape_2()
  A.make_alpha_shape(L)
  op = A.find_optimal_alpha(1)
  if alpha == None:
    alpha = op.next()
    print 'Optimal alpha value to get one connected component is ', \
    alpha
  else:
    print 'Using provided alpha value ',alpha
    print 'Optimal alpha value to get one connected component is ', \
    op.next()

  m = A.set_mode()
  A.set_mode(m.REGULARIZED)
  A.set_alpha(alpha)

  # Get the vertices ordered by edges incidence
  alpha_shape_vertices = []
  for it in A.alpha_shape_vertices:
    alpha_shape_vertices.append(it)
  
  print 'Vertices in alpha-shape ', len(alpha_shape_vertices)

  if interactive:
   while True:
      plt.ion()
      plt.figure(1)
      plt.clf()
      xo = [p.x() for p in L]
      yo = [p.y() for p in L]
      x = [p.point().x() for p in alpha_shape_vertices]
      y = [p.point().y() for p in alpha_shape_vertices]

      plt.plot(xo,yo,'or',x,y,'ob')

      plt.show()
      try:
        alpha = float(raw_input('New alpha value [Enter = finish]: '))
      except ValueError:
        break
          
      A.set_alpha(alpha)

      # Get the vertices ordered by edges incidence
      alpha_shape_vertices = []
      for it in A.alpha_shape_vertices:
          alpha_shape_vertices.append(it)
             
              
      print 'Vertices in alpha-shape ', len(alpha_shape_vertices)
      
  return A,alpha

def orderAlphaShape(A,L):

  print 'Ordering border. This may take a while...'
  # Get the vertices ordered by edges incidence
  alpha_shape_edges = []
  alpha_shape_vertices = []
  for it in A.alpha_shape_edges:
    alpha_shape_edges.append(A.segment(it))
  for it in A.alpha_shape_vertices:
    alpha_shape_vertices.append(it.point())

  nver = len(alpha_shape_vertices)
  incidency = list()
  ind = 0
  for edge in alpha_shape_edges:
    vert1 = edge.vertex(0)
    vert2 = edge.vertex(1)
    incidency.append([alpha_shape_vertices.index(vert1),alpha_shape_vertices.index(vert2)])
              
  border = list()
  borders = list()
  border_index = list()
  indeces = list()
  vert1 = [ i[0] for i in incidency]
  vert2 = [ i[1] for i in incidency]
  
  i = 0
  while len(border) < nver:
    vert = alpha_shape_vertices[vert1[i]]
    if vert not in border:
      border.append(vert)
      border_index.append(L.index(vert))

    vert1.pop(i)
    vert = vert2[i]
    vert2.pop(i)
    try:
      i = vert1.index(vert)
    except ValueError:
      print "Found inital vertex!", vert
      border.append(alpha_shape_vertices[vert])
      borders.append(border)
      indeces.append(border_index)
      border = list()
      border_index = list()
      if len(vert1) == 0:
        break
      else:
        i = 0

  return borders, indeces

def main(infile, outfile, alpha=None,interactive=None):
  #read the data
  L = list()
  with open(infile,'r') as f:
    for line in f:
      s = line.split()
      arg = []
      for n in s:
        if n != ' ':
            arg.append(float(n))
            
      L.append(Point_2(arg[0],arg[1]))
  f.close()
  
  # Get vertices in alpha shape
  [A,alpha] = get_border(L,alpha,interactive)

  # Get the ordered border and the indeces of the vertices
  # in the border in the original list
  [borders, indeces]=orderAlphaShape(A,L)
  # Different borders (like holes ar multiple components)
  # are sublists of borders
  
  #write output
  nborder=0
  with open(outfile,'w') as f:
    for border,index in zip(borders,indeces):
    
      if nborder > 0:
        info = 'NaN NaN NaN\n'
        f.write(info)
      nborder += 1
      
      for p,i in zip(border,index):
        info = str(i) + ' ' + str(p.x()) + ' ' + \
        str(p.y()) + '\n'
        f.write(info)

  f.close()

  return borders

The code runs quite fast, except for my part of the code, for the ordering of the border. I have to work on that! I put the results on Figures 1 and 2. Figure 1 is made using a value for alpha of 0.01, and Figure 2 is using the optimal value. Do you see the diferences?

Fig 1. Alpha-shape for a set of 1000 random points uniformly distributed on the unit disk. The value of alpha is 0.01.

Fig 2. Alpha-shape for a set of 1000 random points uniformly distributed on the unit disk, using the optimal value for alpha.