I have 3D clustered data:

data data

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

  • $\begingroup$ @belisarius, I guess tetgenconvexhull calculates polygondata from data point itself, see updated question. hope I was able to explain the question $\endgroup$
    – user1362
    Commented Jul 26, 2012 at 19:53
  • 1
    $\begingroup$ I think what you need is a 3D alpha shape. See here for example $\endgroup$ Commented Jul 26, 2012 at 21:15
  • 4
    $\begingroup$ 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... $\endgroup$ Commented Jul 26, 2012 at 21:16
  • 6
    $\begingroup$ In 2D you just need a child with a crayon ... $\endgroup$ Commented Jul 26, 2012 at 21:36
  • 1
    $\begingroup$ 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? $\endgroup$ Commented Jul 27, 2012 at 6:29

3 Answers 3


Here's a possible approach.

First use TetGen to tetrahedralize the data:

{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]

enter image description here

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

enter image description here


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]

A "Concave Hull"

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

  • $\begingroup$ 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.. $\endgroup$ Commented Jul 27, 2012 at 17:54
  • $\begingroup$ @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. $\endgroup$
    – Matt W-D
    Commented Jul 27, 2012 at 18:14
  • $\begingroup$ @PlatoManiac: I don't see any spikes; what version of Mathematica are you running? (I'm on v8.0.0 myself) $\endgroup$
    – Matt W-D
    Commented Jul 27, 2012 at 18:16
  • $\begingroup$ @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... $\endgroup$
    – Matt W-D
    Commented Jul 27, 2012 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}];
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}},
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).

  • $\begingroup$ How did you define pts? Thank You $\endgroup$
    – Fortsaint
    Commented Dec 14, 2016 at 1:22
  • $\begingroup$ @Fortsaint pts is the original data points $\endgroup$
    – Taiki
    Commented Dec 14, 2016 at 2:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.