# Concave polygon 'hull' finding

I implemented an algorithm to find the alpha shape of a set of points. The alpha shape is a concave hull for a set of points, whose shape depends on a parameter alpha deciding which points make up the hull.

I have resolved the set of points in concave hull. These points make up a concave polygon. I would like to order these points in a clockwise manner.

Ordering points in a clockwise manner is straightforward when it is a convex shape. With a concave thing, I really don't know what to do. What algorithm is behind? I looked at 'non crossing shortest path' algorithms, 'shortest non crossing Hamiltonian path' algorithms, but I am still not convinced of any of these solutions.

Is there a simple solution?

EDIT

• I've edited my answer with an example, but I should point out that in the graphic you have shown, if the edges in red are those that remain AFTER the alpha shape has done its work, then alpha was probably too small. See that some of the vertices from the original set have been completely dropped (lost) from the tessellation.
– user840
Sep 18, 2012 at 18:11
• @woodchips first, thanks a lot for the time you took to answer. your explanation is useful and raises a question: to obtain that all edges inside alpha shape 'crust' are detected as alpha neighbors (edges belonging to alpha shape), did you use this criterion : if alpha is between alpha_min and alpha_max, 'keep edge, with alpha_min and alpha_max the min and max distance between one of the edge's end vertex, and the voronoi diagram edge that is dual to this vertex' ? I used this criterion only to select edges in alpha shape, and even with higher alpha, I miss points 'inside' alpha shape. Sep 18, 2012 at 20:31
• @woodchips original paper by Edelsbrunner mentions this criterion, but maybe there is some explanation for your results compared to mine, which is a different criterion on your side. Sep 18, 2012 at 20:32
• @woodchips (my favorite one is Python!) Sep 18, 2012 at 20:33
• @woodchips in other words, my question is: why does it happen that edges inside alpha shape are duplicate in your alpha neighbors' list? what in your algorithm is cause for this? Thanks again for your very detailed explanations. Sep 19, 2012 at 6:58

I'm not sure what is your problem.

An alpha shape code starts with a delaunay tessellation (at least mine does.) It erodes any simplex that fails the alpha shape criterion, one simplex at a time. What remains is still a tessellation of the domain.

Next, I would find the outer boundary of that complex, by listing all facets (in 2-d, these are simply edges) of each simplex. These edges are defined by pairs of indexes into the list of vertices. Look for duplicate facets, and exclude all edges that appear twice in the list. (Make sure that you catch edge [1,2] as a duplicate of edge [2,1].) What remains in the list must be only those on the outer surface since they are unshared facets.

Thus in 2-d, we now have a list of edges that comprise the enclosing polygon, which may be concave depending on the value of alpha to generate that alpha shape. Connect the edges into a polygon, traversing one to the next in order. Again, since the list of edges is really a list of references into the list of vertices, that ordering will be simple. If your choice ends up the counter-clockwise polygon, flip them.

Edit:

Perhaps the only way to explain this is by an example. I'm using my own set of tools to compute the alpha shape, etc., but the ideas are the same no matter what tools you use. All computations were done in MATLAB.

I'll start with 121 lattice points in the plane, then exclude a few, by taking a bite out of one edge, thus excluding 34 of those points.

xy = lattice({0:.1:1,0:.1:1})
xy =
0            0
0          0.1
0          0.2
0          0.3
0          0.4
0          0.5
0          0.6
0          0.7
0          0.8
0          0.9
0            1
0.1            0
0.1          0.1

...

1          0.7
1          0.8
1          0.9
1            1
>> k = sqrt(sum(bsxfun(@minus,[.8 .5],xy).^2,2)) < .35;
>> sum(k)
ans =
34
>> xy(k,:) = [];
>> plot(xy(:,1),xy(:,2),'o')


Now, I can form a triangulation of that set. Again, I'll use my own tools, but in matlab, delaunay.m does exactly this.

>> tri = delaunayn(xy)
ans =
84    79    83
86    80    85
75    79    84
84    66    75
75    66    71
85    80    76
76    67    85
...


The tessellation itself is simply a set of references. Each row of that array is one triangle, so we have one of the triangles as [84 79 83], which are references into those points from our original point set. The triangulation is like this:

Note the long facets that span the hole along the right edge. Now, compute an alpha shape. The maximum nearest neighbor distance is 0.1, so here I'll use an alpha ball of radius 0.2 (twice the maximum nearest neighbor distance.)

>> sc = alphashape(xy,.2);
>> tria = sc.tessellation;


So the alpha shape is in yellow, with black edges for its triangles. I've overlaid the delaunay triangulation on top, with blue edges.

The alpha shape has 130 triangles in it, since each triangle has 3 edges, I will have a list of 130x3 = 390 edges, some of which are replicates. There really are only 216 edges in that list.

>> edgelist = [tria(:,[1 2]);tria(:,[1 3]);tria(:,[2 3])];
>> size(unique(sort(edgelist,2),'rows'))
ans =
216     2


What we wish to see however, are only the edges on the boundary of that alpha shape. These edges are those that appear in the list exactly once. Some careful sorting can resolve the boundary edges.

>> sortrows(sort(edgelist,2))
ans =
1     2
1    12
2     3
2    12
2    12
2    13
2    13
2    14
2    14
3     4
3    14
3    14
3    15
3    15
4     5
4    15
4    15
5     6
...


Thus we see that edge [1 2] appeared once in the list, but edge [2 12] shows up twice, as does edge [2 13]. These latter edges are interior edges that we wish to remove if we want to see the boundary. In MATLAB, I'd do it like this to remove those that appeared twice in the list:

>> k = find(all(diff(sortrows(sort(edgelist,2))) == 0,2));
>> edgelist(k,:) = []
edgelist =
79    83
80    85
75    79
76    80
50    51
60    61
49    50
67    68
58    59
78    82
1     2
2    12
10    21
10    20
...


42 edges remain in the 1-manifold representing the boundary polygon.

I still am not sure where you got stuck, but that should clear up some issues. Just use those 42 edges as your polygon, connecting the dots. How you should do this in your favorite language is nothing I can solve, since I don't even know what that is.

• Thanks! My alpha shape algo (Edelbrunner's) goes as you say: Delaunay triangulation and criterion. Like you I retrieve edges that are in the alpha shape. But I did not implement a rule to select the outermost edges. In fact, with a well chosen alpha, I obtained, in the particular case of my point sets, that all points in edges in the alpha shape are on the concave hull. But you are right, in general case, output of alpha shape can contain points inside the concave hull, and a decision rule must be implemented to discard edges inside. Sep 16, 2012 at 13:48
• You make me realize that I have to work with 'alpha neighbor' edges and not 'alpha extreme' points, so that there is no more problem with polygon ordering. Thanks a lot!! Which alpha shape algo did you implement? Sep 16, 2012 at 13:49
• It is simple as I described. (In 2-d) exclude those edges that appear twice in the set of all edges. This MUST return the enclosing polygon as a list of edges. Then just order those edges.
– user840
Sep 16, 2012 at 13:51
• yes I am coding it : ) in fact I was asking for your 'main' alpha shape algo (before alpha shape output). I implemented Edelbrunner's. I still struggle to gently automate choice for alpha. do you have a hint? Sep 16, 2012 at 13:57
• My (MATLAB) code was originally written almost 15 years ago to solve the 2-d and 3-d cases, but it remains proprietary to my former employer. I've since written my own n-d alpha shape code in MATLAB with a few options in how the erosion is done. It has not been posted online.
– user840
Sep 16, 2012 at 13:58