# Finding a Concave Hull

I have a 3d clustered data:

Is there any other way to get concavehull of 3D data points?

-
@belisarius, I am using Mma 8.0.4.0 – mathew Jul 26 '12 at 19:34
I think what you need is a 3D alpha shape. See here for example – Simon Woods Jul 26 '12 at 21:15
It seems to be challenging just to define what the "concave hull" should be mathematically. Perhaps minimum surface area enclosing polyhedron? Then the computation is NP-hard... – Joseph O'Rourke Jul 26 '12 at 21:16
In 2D you just need a child with a crayon ... – image_doctor Jul 26 '12 at 21:36
This is a interesting question for MMA fans. One reference is arxiv.org/pdf/math.CO/9410208.pdf. However alpha shapes are implemented in cgal.org/Manual/3.2/doc_html/cgal_manual/Alpha_shapes_3/… question is can we use the LibraryFunctionLoad capability of MMA here with out going into much detail of CGAL? – PlatoManiac Jul 27 '12 at 6:29

Here's a possible approach.

First use TetGen to tetrahedralize the data:

Needs["TetGenLink"]
{pts, tetrahedra} = TetGenDelaunay[data3D];


Next define a function to compute the radius of the circumsphere of a tetrahedron (formula from Wikipedia)

csr[{aa_, bb_, cc_, dd_}] :=
With[{a = aa - dd, b = bb - dd, c = cc - dd},
Norm[a.a Cross[b, c] + b.b Cross[c, a] + c.c Cross[a, b]]/(2 Norm[a.Cross[b, c]])]


Now calculate the radius of each tetrahedron in your data, and define a function to pick out only those tetrahedra with radius below a given threshold:

radii = csr[pts[[#]]] & /@ tetrahedra;
alphashape[rmax_] := Pick[tetrahedra, radii, r_ /; r < rmax]


This function converts a tetrahedron into a list of 4 triangular faces:

faces[tetras_] :=
Flatten[tetras /. {a_, b_, c_, d_} :> {{a, b, c}, {a, b, d}, {a, c, d}, {b, c, d}}, 1]


We also need to remove internal polygons. Internal polygons appear twice, since they belong to two tetrahedra, so we simply need to prune the list of faces to those which only appear once in the list.

externalfaces[faces_] := Cases[Tally[Sort /@ faces], {face_, 1} :> face]


Now we can generate the graphics for a particular value of rmax

polys = externalfaces@faces@alphashape[2.0];
Graphics3D[GraphicsComplex[pts, Polygon@polys], Boxed -> False]


Here's an animation showing the effect of different values of rmax

-

As pointed out in the comments, there's really no mathematical definition of a concave hull.

Of course, just because there's no mathematical definition does not preclude coming up with something that sort of works. I can think of two ways to do this:

Easy Way, Not General

Your data roughly has axial symmetry parallel to the x-axis. Moreover, all of your coordinates appear to be integers.

So for every integral value of x, you can make a slice of your data between x and x+1. Then, for each slice, compute the convex hull. Finally, stitch all of the results together, like so:

data3D = Import["http://dl.dropbox.com/u/68983831/process.vtk", "VertexData"];

xvals = (#1[[1]] & ) /@ data3D;

slices = Table[Select[data3D, x <= #1[[1]] <= x + 1 & ], {x, Min[xvals], Max[xvals] - 1}];

Show[(Graphics3D[GraphicsComplex[#1, Polygon[TetGenConvexHull[#1]]]] & )
/@ slices]


I'd call this a "duct tape" solution.

Addendum: This works on Mathematica v8.0.0, however it appears that v8.0.4.0 makes a change to the output of TetGenConvexHull. Here's a workaround for the later version:

Show[(Graphics3D[(GraphicsComplex[#1, Polygon[#2]] & ) @@ TetGenConvexHull[#1]] & )
/@ slices]


Complicated Way, More General (And Not Implemented ;-)

A complicated but perhaps more general solution to this problem would be to use Voronoi diagrams. Similar techniques have been proposed for removing data points from banks in k-Nearest Neighbor searches. Here's a paper on this if you're interested.

Sketch of the algorithm:

• Compute the Voronoi diagram for your points
• For each cell in the Voronoi diagram

• Check to see if all of the cell's faces are touching faces of other cells
• If the cell is surrounded by other cells, remove it

This would effectively keep only "exterior" points in your set.

This approach has drawbacks: the first is that it will sometimes discard points that, intuitively, a child with a crayon would probably not. The second is that it's slow. Finally, it violates the KISS design pattern, so it's a chore to write (hence, I didn't bother).

-
If I evaluate your above code I am getting some funny spiky graphics not the picture your answer accompany...can you cross check it once.. – PlatoManiac Jul 27 '12 at 17:54
@Jay: I'm running Mathematica v8.0.0, and I noted that TetGenConvexHull just returns a surface instead of {pts, surface} like your version does. I can't test this, but maybe you could try Show[(Graphics3D[(GraphicsComplex[#1, Polygon[#2]] & ) @@ TetGenConvexHull[#1]] & ) /@ slices] instead of my Show code. – Matt W-D Jul 27 '12 at 18:14
@PlatoManiac: I don't see any spikes; what version of Mathematica are you running? (I'm on v8.0.0 myself) – Matt W-D Jul 27 '12 at 18:16
@Jay: You know... this is really a separate question. Can you break this out and start a new question? I'll try to give it a shot... – Matt W-D Jul 27 '12 at 22:30

Let me update Simon's code for Mathematica 10. We no longer need to explicitly load TetGenLink.

tetrahedra = Level[MeshPrimitives[DelaunayMesh[data3D], 3], {-3}];
alphashape[rmax_] := Pick[tetrahedra, radii, r_ /; r < rmax]
faces[tetras_] := Flatten[
tetras /. {a_, b_, c_, d_} :> {{a, b, c}, {a, b, d}, {a, c, d}, {b, c, d}},
1
];
externalfaces[faces_] := Cases[Tally[Sort /@ faces], {face_, 1} :> face];
rmax = 500;
polys = externalfaces@faces@alphashape[rmax];
Graphics3D[GraphicsComplex[pts, Polygon@polys], Boxed -> False]


For demonstration, see #79313 (Approach #2).

-