\( \newcommand{\E}{\mathrm{E}} \) \( \newcommand{\A}{\mathrm{A}} \) \( \newcommand{\R}{\mathrm{R}} \) \( \newcommand{\N}{\mathrm{N}} \) \( \newcommand{\Q}{\mathrm{Q}} \) \( \newcommand{\Z}{\mathrm{Z}} \) \( \def\ccSum #1#2#3{ \sum_{#1}^{#2}{#3} } \def\ccProd #1#2#3{ \sum_{#1}^{#2}{#3} }\)
CGAL - 3D Alpha Shapes
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Groups Pages
User Manual

Authors
Tran Kai Frank Da, Sébastien Loriot, and Mariette Yvinec
alphashape.png

Assume we are given a set \( S\) of points in 2D or 3D and we would like to have something like "the shape formed by these points". This is quite a vague notion and there are probably many possible interpretations, the alpha shape being one of them. Alpha shapes can be used for shape reconstruction from a dense unorganized set of data points. Indeed, an alpha shape is demarcated by a frontier, which is a linear approximation of the original shape [1].

As mentioned in Edelsbrunner's and Mücke's paper [2], one can intuitively think of an alpha shape as the following. Imagine a huge mass of ice-cream making up the space \( \mathbb{R}^3\) and containing the points as "hard" chocolate pieces. Using one of those sphere-formed ice-cream spoons, we carve out all parts of the ice-cream block we can reach without bumping into chocolate pieces, thereby even carving out holes in the inside (e.g. parts not reachable by simply moving the spoon from the outside). We will eventually end up with a (not necessarily convex) object bounded by caps, arcs and points. If we now straighten all "round" faces to triangles and line segments, we have an intuitive description of what is called the alpha shape of \( S\). The drawing above provides an example of this process in 2D (where our ice-cream spoon is simply a circle).

Alpha shapes depend on a parameter \( \alpha\) after which they are named. In the ice-cream analogy above, \( \alpha\) is the squared radius of the carving spoon. A very small value will allow us to eat up all of the ice-cream except the chocolate points themselves. Thus we already see that the alpha shape degenerates to the point-set \( S\) for \( \alpha \rightarrow 0\). On the other hand, a huge value of \( \alpha\) will prevent us even from moving the spoon between two points since it is too large and we will never spoon up the ice-cream lying in the inside of the convex hull of \( S\). Hence, the alpha shape becomes the convex hull of \( S\) as \( \alpha \rightarrow \infty\).

CGAL offers 2D and 3D alpha shapes. The GUDHI library offers a dD Alpha complex.

Definitions

We distinguish two versions of alpha shapes. Basic alpha shapes are based on the Delaunay triangulation. Weighted alpha shapes are based on its generalization, the regular triangulation (cf. Section Regular Triangulations), replacing the euclidean distance by the power to weighted points.

Let us consider the basic case with a Delaunay triangulation. We first define the alpha complex of the set of points \( S\). The alpha complex is a subcomplex of the Delaunay triangulation. For a given value of \( \alpha\), the alpha complex includes all the simplices in the Delaunay triangulation which have an empty circumscribing sphere with squared radius equal or smaller than \( \alpha\). Here "empty" means that the open sphere does not include any points of \( S\). The alpha shape is then simply the domain covered by the simplices of the alpha complex (see [2]).

In general, an alpha complex is a disconnected and non-pure complex, meaning in particular that the alpha complex may have singular faces. For \( 0 \leq k \leq d-1\), a \( k\)-simplex of the alpha complex is said to be singular if it is not a facet of a \( (k+1)\)-simplex of the complex. CGAL provides two versions of alpha shapes. In the general mode, the alpha shapes correspond strictly to the above definition. The regularized mode provides a regularized version of the alpha shapes. It corresponds to the domain covered by a regularized version of the alpha complex where singular faces are removed (See Figure 46.1 for an example).

gen-reg-ex.png
Figure 46.1 Comparison of general and regularized alpha-shape. Left: Some points are taken on the surface of a torus, three points being taken relatively far from the surface of the torus; Middle: The general alpha-shape (for a large enough alpha value) contains the singular triangle facet of the three isolated points; Right: The regularized version (for the same value of alpha) does not contains any singular facet.

The alpha shapes of a set of points \( S\) form a discrete family, even though they are defined for all real numbers \( \alpha\). The entire family of alpha shapes can be represented through the underlying triangulation of \( S\). In this representation each \( k\)-simplex of the underlying triangulation is associated with an interval that specifies for which values of \( \alpha\) the \( k\)-simplex belongs to the alpha complex. Relying on this fact, the family of alpha shapes can be computed efficiently and relatively easily. Furthermore, we can select the optimal value of \( \alpha\) to get an alpha shape including all data points and having less than a given number of connected components. Also, the alpha-values allows to define a filtration on the faces of the triangulation of a set of points. In this filtration, the faces of the triangulation are output in increasing order of the alpha value for which they appear in the alpha complex. In case of equal alpha values, lower dimensional faces are output first.

The definition is analog in the case of weighted alpha shapes. The input set is now a set of weighted points (which can be regarded as spheres) and the underlying triangulation is the regular triangulation of this set. Two spheres, or two weighted points, with centers \( C_1, C_2\) and radii \( r_1, r_2 \) are said to be orthogonal iff \( C_1C_2 ^2 = r_1^2 + r_2^2\) and suborthogonal iff \( C_1C_2 ^2 < r_1^2 + r_2^2\). For a given value of \( \alpha\), the weighted alpha complex is formed with the simplices of the regular triangulation triangulation such that there is a sphere orthogonal to the weighted points associated with the vertices of the simplex and suborthogonal to all the other input weighted points. Once again the alpha shape is then defined as the domain covered by a the alpha complex and comes in general and regularized versions.

Functionality

Family of Alpha Shapes

The class Alpha_shape_3<Dt,ExactAlphaComparisonTag> represents the whole family of alpha shapes for a given set of points. The class includes the underlying triangulation Dt of the set, and associates to each \( k\)-face of this triangulation an interval specifying for which values of \( \alpha\) the face belongs to the alpha complex. The second template parameter, ExactAlphaComparisonTag, is a tag that, when set to CGAL::Tag_true, triggers exact comparisons between alpha values.

The class provides functions to set and get the current \( \alpha\)-value, as well as an iterator that enumerates the \( \alpha\) values where the alpha shape changes.

Additionally, the class has a filtration member function that, given an output iterator with Object as value type, outputs the faces of the triangulation according to the order of apparition in the alpha complex when alpha increases.

Finally, it provides a function to determine the smallest value \( \alpha\) such that the alpha shape satisfies the following two properties:

  1. All data points are either on the boundary or in the interior of the regularized version of the alpha shape (no singular faces).
  2. The number of components is equal or less than a given number.

The current implementation is static, that is after its construction points cannot be inserted or removed.

Alpha Shape for a Fixed Alpha

Given a value of alpha, the class Fixed_alpha_shape_3<Dt> represents one alpha shape for a given set of points. The class includes the underlying triangulation Dt of the set, and associates to each \( k\)-face of this triangulation a classification type. This class is dynamic, that is after its construction points can be inserted or removed.

Classification and Iterators

Both classes provide member functions to classify for a (given) value of \( \alpha\) the different faces of the triangulation as EXTERIOR, SINGULAR, REGULAR or INTERIOR with respect to the alpha shape. A \( k\)-face on the boundary of the alpha complex is said to be: REGULAR if it is a subface of the alpha-complex which is a subface of a \( (k+1)\)-face of the alpha complex, and SINGULAR otherwise. A \( k\)-face of the alpha complex which is not on the boundary of the alpha complex is said to be INTERIOR. See Figure 46.2 for a 2D illustration.

classif.png
Figure 46.2 Classification of simplices, a 2D example. Left: The 2D Delaunay triangulation of a set of points; Right: Classification of simplices for a given alpha (the squared radius of the red circle). INTERIOR, REGULAR and SINGULAR simplices are depicted in black, green and blue respectively. EXTERIOR simplices are not depicted. The vertex s and the edge tu are SINGULAR since all higher dimension simplices they are incident to are EXTERIOR. The facet pqr is EXTERIOR because the squared radius of its circumscribed circle is larger than alpha.

The classes provide also output iterators to get for a given alpha value the vertices, edges, facets and cells of the different types (EXTERIOR, SINGULAR, REGULAR or INTERIOR).

Concepts and Models

We currently do not specify concepts for the underlying triangulation type. Models that work for a family of alpha-shapes are the instantiations of the classes Delaunay_triangulation_3 and Periodic_3_Delaunay_triangulation_3 (see example Example for Periodic Alpha Shapes). A model that works for a fixed alpha-shape are the instantiations of the class Delaunay_triangulation_3. Models that work for a weighted alpha-shape are the instantiations of the classes Regular_triangulation_3 and Periodic_3_regular_triangulation_3. The triangulation needs a geometric traits class and a triangulation data structure as template parameters.

Alpha Shapes

For the class Alpha_shape_3<Dt,ExactAlphaComparisonTag>, the requirements of the traits class are described in the concepts AlphaShapeTraits_3 in the non-weighted case and WeightedAlphaShapeTraits_3 in the weighted case. All CGAL kernels are models of both concepts.

The triangulation data structure of the triangulation has to be a model of the concept TriangulationDataStructure_3, and it must be parameterized with vertex and cell classes, which are model of the concepts AlphaShapeVertex_3 and AlphaShapeCell_3. The package provides by default the classes Alpha_shape_vertex_base_3<Gt> and Alpha_shape_cell_base_3<Gt>.

Fixed Alpha Shapes

For the class Fixed_alpha_shape_3<Dt>, the requirements of the traits class are described in the concepts FixedAlphaShapeTraits_3 in the non-weighted case and FixedWeightedAlphaShapeTraits_3 in the weighted case. All CGAL kernels are models of both concepts. The triangulation data structure of the triangulation has to be a model of the concept TriangulationDataStructure_3, and it must be parameterized with vertex and cell classes, which are model of the concepts FixedAlphaShapeVertex_3 and FixedAlphaShapeCell_3. The package provides models Fixed_alpha_shape_vertex_base_3<Gt> and Fixed_alpha_shape_cell_base_3<Gt>, respectively.

Triangulation data structure

Additionally requirements are put when using weighted or periodic triangulations as underlying triangulation:

Alpha_shape_3 vs. Fixed_alpha_shape_3

The class Alpha_shape_3<Dt,ExactAlphaComparisonTag> represents the whole family of alpha shapes for a given set of points while the class Fixed_alpha_shape_3<Dt> represents only one alpha shape (for a fixed alpha). When using the same kernel, Fixed_alpha_shape_3<Dt> is a lighter version. It is thus naturally much more efficient when the alpha-shape is needed for a single given value of alpha. In addition, note that the class Alpha_shape_3<Dt,ExactAlphaComparisonTag> requires constructions (squared radius of simplices) while the class Fixed_alpha_shape_3<Dt> uses only predicates. This implies that a certified construction of one (several) alpha-shape, using the Alpha_shape_3<Dt,ExactAlphaComparisonTag> requires a kernel with exact predicates and exact constructions (or setting ExactAlphaComparisonTag to Tag_true) while using a kernel with exact predicates is sufficient for the class Fixed_alpha_shape_3<Dt>. This makes the class Fixed_alpha_shape_3<Dt> even more efficient in this setting. In addition, note that the Fixed version is the only of the two that supports incremental insertion and removal of points.

We give the time spent while computing the alpha shape of a protein (considered as a set of weighted points) featuring 4251 atoms (using gcc 4.3 under Linux with -O3 and -DNDEBUG flags, on a 2.27GHz Intel(R) Xeon(R) E5520 CPU): Using Exact_predicates_inexact_constructions_kernel, building the regular triangulation requires 0.09s, then the class Fixed_alpha_shape_3<Dt> required 0.05s while the class Alpha_shape_3<Dt,ExactAlphaComparisonTag> requires 0.35s if ExactAlphaComparisonTag is Tag_false (and 0.70s with Tag_true). Using Exact_predicates_exact_constructions_kernel, building the regular triangulation requires 0.19s and then the class Alpha_shape_3<Dt,ExactAlphaComparisonTag> requires 0.90s.

Examples

Example for Basic Alpha-Shapes

This example builds a basic alpha shape using a Delaunay triangulation as underlying triangulation.


File Alpha_shapes_3/ex_alpha_shapes_3.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Alpha_shape_3.h>
#include <fstream>
#include <list>
#include <cassert>
typedef CGAL::Triangulation_data_structure_3<Vb,Fb> Tds;
typedef Gt::Point_3 Point;
typedef Alpha_shape_3::Alpha_iterator Alpha_iterator;
int main()
{
std::list<Point> lp;
//read input
std::ifstream is("./data/bunny_1000");
int n;
is >> n;
std::cout << "Reading " << n << " points " << std::endl;
Point p;
for( ; n>0 ; n--) {
is >> p;
lp.push_back(p);
}
// compute alpha shape
Alpha_shape_3 as(lp.begin(),lp.end());
std::cout << "Alpha shape computed in REGULARIZED mode by default"
<< std::endl;
// find optimal alpha value
Alpha_iterator opt = as.find_optimal_alpha(1);
std::cout << "Optimal alpha value to get one connected component is "
<< *opt << std::endl;
as.set_alpha(*opt);
assert(as.number_of_solid_components() == 1);
return 0;
}

Building Basic Alpha Shapes for Many Points

When many points are input in the alpha shape, say more than 10 000, it may pay off to use a Delaunay triangulation with the Fast_location policy as underlying triangulation in order to speed up point location queries (cf. Section The Location Policy Parameter).


File Alpha_shapes_3/ex_alpha_shapes_with_fast_location_3.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Alpha_shape_3.h>
#include <fstream>
#include <list>
#include <cassert>
typedef CGAL::Triangulation_data_structure_3<Vb,Fb> Tds;
typedef CGAL::Alpha_shape_3<Delaunay> Alpha_shape_3;
typedef K::Point_3 Point;
typedef Alpha_shape_3::Alpha_iterator Alpha_iterator;
typedef Alpha_shape_3::NT NT;
int main()
{
Delaunay dt;
std::ifstream is("./data/bunny_1000");
int n;
is >> n;
Point p;
std::cout << n << " points read" << std::endl;
for( ; n>0 ; n--) {
is >> p;
dt.insert(p);
}
std::cout << "Delaunay computed." << std::endl;
// compute alpha shape
Alpha_shape_3 as(dt);
std::cout << "Alpha shape computed in REGULARIZED mode by defaut."
<< std::endl;
// find optimal alpha values
Alpha_shape_3::NT alpha_solid = as.find_alpha_solid();
Alpha_iterator opt = as.find_optimal_alpha(1);
std::cout << "Smallest alpha value to get a solid through data points is "
<< alpha_solid << std::endl;
std::cout << "Optimal alpha value to get one connected component is "
<< *opt << std::endl;
as.set_alpha(*opt);
assert(as.number_of_solid_components() == 1);
return 0;
}

Example for Weighted Alpha-Shapes

The following examples build a weighted alpha shape requiring a regular triangulation as underlying triangulation. The alpha shape is built in GENERAL mode.


File Alpha_shapes_3/ex_weighted_alpha_shapes_3.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_3.h>
#include <CGAL/Alpha_shape_3.h>
#include <list>
typedef CGAL::Triangulation_data_structure_3<Vb,Cb> Tds;
typedef CGAL::Regular_triangulation_3<K,Tds> Triangulation_3;
typedef Alpha_shape_3::Cell_handle Cell_handle;
typedef Alpha_shape_3::Vertex_handle Vertex_handle;
typedef Alpha_shape_3::Facet Facet;
typedef Alpha_shape_3::Edge Edge;
typedef Triangulation_3::Weighted_point Weighted_point;
typedef Triangulation_3::Bare_point Bare_point;
int main()
{
std::list<Weighted_point> lwp;
//input : a small molecule
lwp.push_back(Weighted_point(Bare_point( 1, -1, -1), 4));
lwp.push_back(Weighted_point(Bare_point(-1, 1, -1), 4));
lwp.push_back(Weighted_point(Bare_point(-1, -1, 1), 4));
lwp.push_back(Weighted_point(Bare_point( 1, 1, 1), 4));
lwp.push_back(Weighted_point(Bare_point( 2, 2, 2), 1));
//build alpha_shape in GENERAL mode and set alpha=0
Alpha_shape_3 as(lwp.begin(), lwp.end(), 0, Alpha_shape_3::GENERAL);
//explore the 0-shape - It is dual to the boundary of the union.
std::list<Cell_handle> cells;
std::list<Facet> facets;
std::list<Edge> edges;
as.get_alpha_shape_cells(std::back_inserter(cells),
as.get_alpha_shape_facets(std::back_inserter(facets),
as.get_alpha_shape_facets(std::back_inserter(facets),
as.get_alpha_shape_edges(std::back_inserter(edges),
std::cout << " The 0-shape has : " << std::endl;
std::cout << cells.size() << " interior tetrahedra" << std::endl;
std::cout << facets.size() << " boundary facets" << std::endl;
std::cout << edges.size() << " singular edges" << std::endl;
return 0;
}

Example for Fixed Weighted Alpha-Shapes

Same example as previous but using a fixed value of alpha.


File Alpha_shapes_3/ex_fixed_weighted_alpha_shapes_3.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_3.h>
#include <CGAL/Fixed_alpha_shape_3.h>
#include <CGAL/Fixed_alpha_shape_vertex_base_3.h>
#include <CGAL/Fixed_alpha_shape_cell_base_3.h>
#include <list>
typedef CGAL::Triangulation_data_structure_3<Vb,Cb> Tds;
typedef CGAL::Regular_triangulation_3<K,Tds> Triangulation_3;
typedef CGAL::Fixed_alpha_shape_3<Triangulation_3> Fixed_alpha_shape_3;
typedef Fixed_alpha_shape_3::Cell_handle Cell_handle;
typedef Fixed_alpha_shape_3::Vertex_handle Vertex_handle;
typedef Fixed_alpha_shape_3::Facet Facet;
typedef Fixed_alpha_shape_3::Edge Edge;
typedef Triangulation_3::Weighted_point Weighted_point;
typedef Triangulation_3::Bare_point Bare_point;
int main()
{
std::list<Weighted_point> lwp;
//input : a small molecule
lwp.push_back(Weighted_point(Bare_point( 1, -1, -1), 4));
lwp.push_back(Weighted_point(Bare_point(-1, 1, -1), 4));
lwp.push_back(Weighted_point(Bare_point(-1, -1, 1), 4));
lwp.push_back(Weighted_point(Bare_point( 1, 1, 1), 4));
lwp.push_back(Weighted_point(Bare_point( 2, 2, 2), 1));
//build one alpha_shape with alpha=0
Fixed_alpha_shape_3 as(lwp.begin(), lwp.end(), 0);
//explore the 0-shape - It is dual to the boundary of the union.
std::list<Cell_handle> cells;
std::list<Facet> facets;
std::list<Edge> edges;
as.get_alpha_shape_cells(std::back_inserter(cells),
as.get_alpha_shape_facets(std::back_inserter(facets),
as.get_alpha_shape_facets(std::back_inserter(facets),
as.get_alpha_shape_edges(std::back_inserter(edges),
std::cout << " The 0-shape has : " << std::endl;
std::cout << cells.size() << " interior tetrahedra" << std::endl;
std::cout << facets.size() << " boundary facets" << std::endl;
std::cout << edges.size() << " singular edges" << std::endl;
return 0;
}

Building an Alpha Shapes with Exact Comparisons of Alpha Values

On some platforms, the alpha shapes of the set of points of this example cannot be computed when using a traits with inexact constructions. To be able to compute them with a traits with inexact constructions, the tag ExactAlphaComparisonTag is set to Tag_true.
File Alpha_shapes_3/ex_alpha_shapes_exact_alpha.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Alpha_shape_3.h>
#include <fstream>
#include <list>
#include <cassert>
typedef CGAL::Tag_true Alpha_cmp_tag;
//We use CGAL::Default to skip the complete declaration of base classes
typedef CGAL::Triangulation_data_structure_3<Vb,Fb> Tds;
//Alpha shape with ExactAlphaComparisonTag set to true (note that the tag is also
//set to true for Vb and Fb)
typedef Gt::Point_3 Point;
int main()
{
//Set of points for which the alpha shapes cannot be computed with
//a floating point alpha value (on certain platforms)
std::list<Point> lp;
lp.push_back(Point(49.2559,29.1767,6.7723));
lp.push_back(Point(49.3696,31.4775,5.33777));
lp.push_back(Point(49.4264,32.6279,4.6205));
lp.push_back(Point(49.3127,30.3271,6.05503));
// compute alpha shape
Alpha_shape_3 as(lp.begin(),lp.end(),0,Alpha_shape_3::GENERAL);
return 0;
}

Example for Periodic Alpha Shapes

The following example shows how to use a periodic Delaunay triangulation (Chapter 3D Periodic Triangulations) as underlying triangulation for the alpha shape computation. Usage of a weighted Delaunay periodic triangulation is presented in the example: ex_weighted_periodic_alpha_shapes_3.cpp.

In order to define the original domain and to benefit from the built-in heuristic optimizations of the periodic triangulation computation, it is recommended to first construct the triangulation and then construct the alpha shape from it. The alpha shape constructor that takes a point range can be used as well but in this case the original domain cannot be specified and the default unit cube will be chosen and no optimizations will be used.

It is also recommended to switch the triangulation to 1-sheeted covering if possible. Note that a periodic triangulation in 27-sheeted covering space is degenerate. In this case an exact constructions kernel needs to be used to compute the alpha shapes. Otherwise the results will suffer from round-off problems.


File Alpha_shapes_3/ex_periodic_alpha_shapes_3.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Random.h>
#include <CGAL/point_generators_3.h>
// Traits
// Vertex type
// Cell type
typedef CGAL::Triangulation_data_structure_3<AsVb,AsCb> Tds;
typedef CGAL::Alpha_shape_3<P3DT3> Alpha_shape_3;
typedef PK::Point_3 Point;
int main()
{
CGAL::Random random(7);
CGAL::Random_points_in_cube_3<Point, Creator> in_cube(1, random);
std::vector<Point> pts;
// Generating 1000 random points
for (int i=0 ; i < 1000 ; i++) {
Point p = *in_cube++;
pts.push_back(p);
}
// Define the periodic cube
P3DT3 pdt(PK::Iso_cuboid_3(-1,-1,-1,1,1,1));
// Heuristic for inserting large point sets (if pts is reasonably large)
pdt.insert(pts.begin(), pts.end(), true);
// As pdt won't be modified anymore switch to 1-sheeted cover if possible
if (pdt.is_triangulation_in_1_sheet()) pdt.convert_to_1_sheeted_covering();
std::cout << "Periodic Delaunay computed." << std::endl;
// compute alpha shape
Alpha_shape_3 as(pdt);
std::cout << "Alpha shape computed in REGULARIZED mode by default."
<< std::endl;
// find optimal alpha values
Alpha_shape_3::NT alpha_solid = as.find_alpha_solid();
Alpha_shape_3::Alpha_iterator opt = as.find_optimal_alpha(1);
std::cout << "Smallest alpha value to get a solid through data points is "
<< alpha_solid << std::endl;
std::cout << "Optimal alpha value to get one connected component is "
<< *opt << std::endl;
as.set_alpha(*opt);
assert(as.number_of_solid_components() == 1);
return 0;
}