Typical problem in computational geometry: given a set of points, what is its shape? This may often be intuitive for a human being, but try to define this concept mathematically, then implement it numerically.

One of the proposals is called the “α shape”, and is closely related to the Delaunay triangulation. More details below.

See also alpha shapes, for a recent application of this concept to molecular simulation. Check out the movie too!

An intuitive definition of the procedure is as follows. One may think of the whole of space as filled out with ice cream, with our points being chips in it. A spherical scoop carves out balls of ice cream, without removing the chips (notice we may also carve in between the points, reaching perhaps from a 4th dimension). The result of the process will be a region of space which may have a complicated shape, which is by bounded caps, arcs and points (see picture above). If we straighten all caps to triangles and all arcs to line segments, we end up with an α shape. This object has a well-defined border: the set of bounding points – these are the points that have been reached by the scoop but are still part of the α shape (i.e. the spherical scoop has been able to remove some of the ice cream around them, but not *all* of it).

What is α? It is related to the radius of the scoop, *R*. Historically, it is either *R* squared, or (confusingly) its inverse, 1/*R*. A sign is sometimes introduced; we have just described the “negative α” shape.

Two limits will perhaps clarify the procedure. If α is very small, all of the points will be reached; the resulting shape is, therefore, all of the points – but none of them will belong to the “border”. The other limit is perhaps more interesting: for values of α very large, the scoop takes half-spaces out, and can only reach the points that protrude from the surface. In this limit, the resulting α shape has as its border the very well-known *convex hull*.

The CGAL project (see previous post) includes routines to calculate these shapes. The α shape is evaluated once a Delaunay triangulation is computed: the *empty ball condition* permits to associate a given radius to each triangle, below which the scoop will be able to “get in”.

Let’s get a bit technical. One would need a preamble like this one:

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Delaunay_triangulation_3.h> #include <CGAL/Alpha_shape_3.h> #include <list> struct K : CGAL::Exact_predicates_inexact_constructions_kernel {}; typedef CGAL::Alpha_shape_cell_base_3<K> Fb; typedef CGAL::Triangulation_data_structure_3<Fb> Tds; typedef CGAL::Delaunay_triangulation_3<K,Tds> Triangulation_3; typedef CGAL::Alpha_shape_3<Triangulation_3> Alpha_shape_3; typedef Alpha_shape_3::Alpha_iterator Alpha_iterator; typedef Alpha_shape_3::Vertex_handle Vertex_handle; typedef std::list<Vertex_handle>::const_iterator Vl_it;

One can build a triangulation T, then its α shape:

Alpha_shape_3 as(T);

Then, select a value for α (or a close one)

Alpha_iterator aa = as.alpha_upper_bound(1.44); as.set_alpha(*aa);

Then, get its boundary thus:

std::list<Vertex_handle> al_vs; as.get_alpha_shape_vertices(back_inserter(al_vs), Alpha_shape_3::REGULAR);

(In 3D, it seems different in 2D!)

OK, I did warn that CGAL makes strong use of the STL. But, at the end of the day you have your list of external points, which you may traverse and print out:

for(Vl_it it=al_vs.begin();it!=al_vs.end(); it++) { double z=(*it)->point()[2]; if( z>0 ) { size++; sphere(image, (*it)->point(), 0, 0.5001 ); } }

The “sphere” function, incidentally, creates a chunk of script to be fed to POV-Ray in order to create a nice picture. More about this on another post!

Pingback: Natural coordinates « Daniel Duque Campayo

Could you leave a link to some code that uses the CGAL library to take a list of points as an input, and output a list of points as an output? Basically the code that you wrote, but in one file that has been tested.

Sure thing. The following code has been compiled against the latest version of the CGAL libs (3.8), and looks to be working fine.

https://docs.google.com/leaf?id=0B3zGtZjOS9UzMTlhODRlZDAtZDVkYy00NzQ0LTgyMDUtMTBhMzVkYWQwYTVk&hl=en_US

Amazing post. I am working with DEM and I need to find the boundary of a set of points. However, the points are in 2D. I tried to modify your code, but I couldn’t. Can you help me? Thanks.

You’d have to use 2D alpha shapes then, obviously. Other that, I am really busy these days. But I have some ideas: first, have a look at the manual,

http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Alpha_shapes_2/Chapter_main.html . The CGAL email list (cgal-discuss) is also great, try to look for info there: https://sympa.inria.fr/sympa/arc/cgal-discuss . If you still have problems, you can email the list for help, they usually are really kind.