Finding a Concave Hull

Question

I have 3D clustered data:

Is there any other way to get the concave hull of 3D data points?

Answer 1

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

Answer 2

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}]; 

Needs["TetGenLink`"]; 

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).

Answer 3

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

tetrahedra = Level[MeshPrimitives[DelaunayMesh[data3D], 3], {-3}];
radius[p_] := Sqrt[Area[Circumsphere[p]]/(4 Pi)];
radii = radius /@ tetrahedra;
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).