DGtal
0.9.3beta

In this part of the manual, we describe how to compute wellknown geometric structures such as the convex hull or the alphashape of a 2D point set. The proposed implementations provide exact results as long as the underlying geometric predicates are robust (see Implementation of geometric predicates for more details about geometric predicates).
Related examples are exampleConvexHull2D.cpp and exampleAlphaShape.cpp.
Let us recall that a polygon is a weakly externally visible polygon (WEVP) if for every point on it, there exists an infinite halfline beginning at this point, which does not intersect the polygon at anywhere else.
The two basic procedures functions::Hull2D::buildHullWithStack (the object having a stack interface is passed by reference) and functions::Hull2D::buildHullWithAdaptedStack (the object having a stack interface is passed by copy) retrieve the extremal vertices of a WEVP in lineartime. This is an implementation of a wellknown method, called Sklansky's scan, Graham's scan or 3coins algorithm, which has been proven to correctly retrieve the convex hull of any WEVP [Toussaint and Avis, 1982: [70]].
Three convex hull algorithms use this basic procedure:
Before calling these algoritms, you must use the functions::Hull2D namespace defined in the file Hull2DHelpers.h.
Moreover, you need to create an orientation functor...
... and, if required, a predicate
See Implementation of geometric predicates for more details about geometric predicates.
The more versatile algorithm is Andrew's one [Andrew, 1979: [4]], also known to be the monotonechain algorithm. First, points are sorted along the horizontal axis. Then, the lower and upper convex hulls are computed by a simple Graham scan from left to right and from right to left. The resulting algorithm is correct because sorting points by increasing xcoordinates (and increasing ycoordinates for a same xcoordinate) produces a WEVP. The overall time complexity is \( O(n\log{(n)}) \) for a 2D point set of size \( n \).
In order to retrieve the extremal points of a range of points, you may call the procedure andrewConvexHullAlgorithm as follows:
Note that in this example the extremal points are pushed to a back insertion sequence called 'res'.
You can however change the predicate to accept aligned points and therefore all points lying on the boundary of the convex hull.
You can also change the predicate so that the resulting list of vertices is clockwiseoriented.
Graham's algorithm is another wellknown algorithm to compute the convex hull of a 2D point set of size \( n \) in \( O(n\log{(n)}) \) [Graham, 1972: [35]]. The method can be coarsely described as follows:
We have chosen the point of maximal xcoordinate (and maximal ycoordinate) to be the pole. Hence, the algorithm is correct since sorting the points by increasing angle about such a pole produces again a WEVP.
The main advantage of Graham's algorithm over Andrew's one is that only one scan is required instead of two to retrieve the whole convex hull.
You may call the procedure grahamConvexHullAlgorithm as follows:
The two previous algorithms are offline and compute the convex hull of any set of \( n \) 2D points in \( O(n\log{(n)}) \). Melkman's algorithm is online and computes the convex hull of a simple polygonal line, ie. without selfintersection, in lineartime (amortized constant time per update).
Note that Graham's scan, used by the two previous algorithms, runs in lineartime on weakly externally visible polygons, which is a subset of simple polygonal lines. As a consequence, Graham's scan may fail for simple polygonal lines that are not weakly externally visible.
In exampleConvexHull2D.cpp, you may find a minimal example:
On such example, Graham's scan fails because the polygonal line is not weakly externally visible. However, Melkman's algorithm produces the right output:
Melkman's algorithm works fine because it is based on a double ended queue and it assumes that the input points form a simple polygonal line. Due to this assumption, any new point cannot be located in the cone formed by the first and last edge of the current convex hull (in red in the figure below). As a consequence, it is enough to update the convex hull by a Graham's scan from the front and/or from the back of the deque and it is never required to remove/insert points in the middle of the deque.
You may use MelkmanConvexHull class as follows:
The class MelkmanConvexHull provides an online way of computing the convex hull of a sequence of points, which form a simple polygonal line: each time a point is added, the current convex hull is updated accordingly. However, if you are not interested in the online feature, you may call the offline procedure melkmanConvexHullAlgorithm as follows:
If you want to construction a convex hull by considering the front and back of a contour you will need to use the MelkmanConvexHull::reverse() method. For instance, to recover the convex hull of the point P0 of the following figure, you will need to do as follows:
Such sequence will produce the correct result given on the image (2) of the following figure.
From a convex hull it can be useful to compute its associated thickness. We propose the implementation of the well known rotating caliper algorithm allowing to compute all the antipodal pairs in linear time [Shamos , 1978: [66]].
Two definitions of thickness can be computed from an iterator on the convex hull. The first one is the vertical/horizontal thickness obtained by projection on the main xy axis from antipodal pairs. The second one is given from the normal projection on the edge of the same antipodal pairs. The two definitions are illustrated on the following figures:
To compute the thickness of a convex hull you have simply to call the computeHullThickness procedure:
To recover the antipodal pairs for which the thickness is minimal you can also use the same method with an antipodal pair:
And display the result:
By using the sequence of points of the first example, you will obtain a thickness of 12 with the following antipodal pair:
The alphashape of a finite set of points has been introduced in [Edelsbrunner et. al., 1983: [31]]. This geometric structure is based on the following concept of generalized disk. For an arbitrary real \( \alpha \), a generalized disk of radius \( 1/\alpha \) is defined as
We also define the direct (resp. indirect) generalized disk of radius \( 1/\alpha \) passing by two points \( P \) and \( Q \) as the unique generalized disk of radius \( 1/\alpha \) with \( P, Q \) on its boundary and such that \( P, Q \) and its center are counterclockwise oriented (resp. clockwise oriented).
The \(\alpha\)hull of a point set is the intersection of all closed generalized disks of radius \( 1/\alpha \) that contains the point set. A point \( P \) of this point set is \(\alpha\)extreme if there exists a closed generalized disk of radius \( 1/\alpha \) such that \( P \) lies on its boundary and which contains all the other points. The \(\alpha\)shape is the straightline graph whose vertices are the \(\alpha\)extreme points and whose edges connect any two \(\alpha\)extreme points such that there exists a closed generalized disk of radius \( 1/\alpha \) with both points on its boundary and which contains all the other points.
Any \(\alpha\)shape of a set of points is a subgraph of either the closest point or the farthest point Delaunay triangulation. Moreover, the set of \(\alpha_1\)extreme points is a subset of the set of \(\alpha_2\)extreme points if \(\alpha_1 \geq \alpha_2 \).
Various \(\alpha\)shapes of the digitization of a disk are depicted below for \( \alpha \in \{ 1, 1/\sqrt{5}, 1/5, 0, 1/9, 1/8 \} \).
You may observe that:
Even if the convex hull of a set of points is its \(0\)shape, the following properties are not true for general \(\alpha\)shapes, \( \alpha \neq 0 \):
As a consequence, the current implementation cannot compute the alphashape of an arbitrary set of points. Such an algorithm would have required to compute the closest point and farthest point Delaunay triangulation. Though, we provide procedures, based on the 3coins algorithm, which compute in lineartime the alphashape of specific polygons, for which the two above properties remain true: closedGrahamScanFromVertex (the starting point is expected to be extremal) and closedGrahamScanFromAnyPoint (no assumption about the starting point is required).
This class of polygons is the set of simple polygons such that:
Let us consider a connected and convex digital object such that its 4connected (resp. 8connected) boundary is simply 4connected (resp. 8connected). It is clear that its 4connected (resp. 8connected) boundary belongs to the above class of polygons provided that \( 1 \leq \alpha < 1/R_{min} \) (resp. \( \sqrt{2}/2 \leq \alpha < 1/R_{min} \)). We can therefore compute any alphashape of the digitization of a compact and convex set by a lineartime scan of the points of its boundary.
In such cases, you can perform the computation in two steps. First, you must define a predicate from an instance of InGeneralizedDiskOfGivenRadius as follows:
Then, you can call the procedure closedGrahamScanFromAnyPoint as follows:
In order to create an instance of InGeneralizedDiskOfGivenRadius, you must provide the value of \( 1/\alpha \) by three arguments:
For example, the following line set \( \alpha \) to \( 1/9 \).
The above illustrations are generated by the program written in exampleAlphaShape.cpp.