I'm trying to identify the shape and the boundaries of objects. The data was generated with a laser scanner and represent the surface of an area in which many boulders lay on the ground and making up heaps.

the data is available here

enter image description here

In the center of the plot you can see a bigger boulder an i want to seperate it from the other stuff around. For this goal I tried different methods which I want to present you and hope that you could enhance these idea(s)

But at first i will reduce the number of samples in the data, because for the most ways the data is to big and need enormous computation time:

data3D = Import["data-a.dat"];
data3Ds = Round[data3D, 0.05] // DeleteDuplicates;

My first intention was to work with ImageProcessing like in this question: Image segmentation and object separation in 3D using Mathematica, because i worked a lot with 2D image processing. At first sight it was a nice idea to create an 3D image out of the data and worked more or less. But the result was not satisfying. The main problem is that all filters and transformations detect the edges between the data points and the NULL-values and doesn't (what i want to do) segment between the data. The presented DistanceTransfrom was created for volumetric data, in my example it is a surface in a three dimensional space. I would like to show some examples and the results, but the .nb-script crashed and I have no access to the results anymore...

As second try I used the idea of the documentation of VertexNormals. The idea was to get the orientation of smaller parts of the surface. For this an interpolation of the data is needed. The documentation suggest to use ListSurfacePlot3D.

g = ListSurfacePlot3D[data3Ds, Mesh -> None, MaxPlotPoints -> 50]

enter image description here

{pl, vl} = 
    GraphicsComplex[pl_, prims_, VertexNormals -> vl_, 
      opts___?OptionQ] :> {pl, vl}, Infinity];
vl1 = Table[Norm[vl[[i]] - {0, 0, 1}], {i, Length[vl]}];
nf = Nearest[pl -> vl1];
vl2 = ParallelMap[StandardDeviation[nf[pl[[#]], 50]] &, 
ListPointPlot3D[List /@ pl, 
 PlotStyle -> ({PointSize[Medium], 
      Blend[{{.1, Darker[Green]}, {0.45, Yellow}, {0.7, 
         Red}}, #1]} & /@ vl2)]

vl is a list of vectors representing the orientation of a point pts and with the length of 1. vl1 is the difference of the orientation between all vl and the vertical. At the end vl2 represent the deviation of the orientation of 50 surrounding data points at one point.

enter image description here

The result is not that bad, but at the one hand I have no idea how to seperate the different parts. At the other hand this version is very slow and only for the reduced data points. Additionally the ListSurfacePlot3D creates some artefacts like in this question.

At the end I tried another interpolation method and worked with tetrahedral meshs of the TetGenLink Package. The simple interpolation function doesn't work good because the data is mainly unstructured (see this link and the related questions) and have more than on z-value for some x-y-pairs.

Maybe you have better ideas for this segmentation problem or could give me some advice to think into other directions.

I would apppreciate no solution for MMA version 10 because I'm just working with version 9.

Thank you!

Edit 1 (Delaunay)

I've tried already the TetGenLink-Package with the this help.

<< TetGenLink`
{pts, tetrahedra} = TetGenDelaunay[data3Ds];
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]])]
radii = csr[pts[[#]]] & /@ 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]
polys = externalfaces@faces@alphashape[0.1];
Graphics3D[GraphicsComplex[pts, Polygon@polys], Boxed -> False]

enter image description here

The result is not that bad and works quite good with version 9. Just I don't have any idea to work further with this Graphic.

  • $\begingroup$ Do you have a set of criteria that is commonly used for this type of segmentation. Features that might be unique to, or indicative of, boulders ? $\endgroup$ Aug 27, 2014 at 13:56
  • $\begingroup$ @image_doctor no I don't have any criteria. Maybe that's way I'm stucked. Do you have any suggestions? The only idea was that the edge of a boulder is defined by the different orientation of the adjacent faces. $\endgroup$
    – zeckra
    Aug 27, 2014 at 14:09

2 Answers 2


You can fit your data to a plane surface and use ComponentMeasurements to detect the bumps on the surface.

data3Ds = Round[data3D, 0.02] // DeleteDuplicates;
plane = LinearModelFit[data3D, {x, y}, {x, y}]
data3Dp = 
  ParallelMap[{#[[1]], #[[2]], #[[3]] - plane[#[[1]], #[[2]]]} &, 
projection = Rasterize[
  ListPlot3D[Reverse[SortBy[data3Dp, Last]], InterpolationOrder -> 0, 
   ColorFunction -> GrayLevel, Mesh -> False, PlotRange -> All, 
   ViewPoint -> {0, 0, Infinity}, Boxed -> False, Axes -> False, 
   BoundaryStyle -> None, MaxPlotPoints -> Infinity, 
   PlotRangePadding -> 0]]
bin = DeleteSmallComponents[MorphologicalBinarize[projection, .99]]
boulders = 
 ComponentMeasurements[bin, {"Centroid", "EquivalentDiskRadius"}]
{1 -> {{177.991, 202.474}, 75.1919}}
 Graphics[{Red, Thick, 
   Circle[#[[1]], #[[2]]] & /@ boulders[[All, 2]]}]]

Mathematica graphics

The only thing that's left is to transform the boulders coordinates back to the initial data.


This will not work on version 9...

Because there are too much points (about 300295) my laptop can not handle the whole data. The main idea is to use Delaunay triangulation DelaunayMesh in Mathematica.

dat = data3D[[1 ;; 3000, All]];
mesh = DelaunayMesh[dat];
HighlightMesh[mesh, {Style[0, Directive[PointSize[Medium], Black]], 
  Style[2, Opacity[0.1]]}]

You get


You can further reduce the points by using interpolation or some other techniques. Alternatively, you divide the meshes and combine them later. Once you get 3d result you can segment the 3d image based on other techniques explained in related questions.

  • $\begingroup$ Could the V9 TetGenDelaunay be used instead of the V10 DelaunayMesh? (I'm not familiar with TetGenLink.) $\endgroup$
    – Michael E2
    Aug 27, 2014 at 13:13
  • 3
    $\begingroup$ @MichaelE2, yes, because V10 DelaunayMesh uses TetGenLink which is available since V8 $\endgroup$
    – user21
    Aug 27, 2014 at 13:16
  • $\begingroup$ @s.s.o that is the reason why I created data3DS which contains only 6063 data points $\endgroup$
    – zeckra
    Aug 27, 2014 at 13:23
  • $\begingroup$ @MichaelE2 yes it works with V9. I've eddited my question and added some code which works fine $\endgroup$
    – zeckra
    Aug 27, 2014 at 13:24
  • 2
    $\begingroup$ @zeckra Instead of dealing with geometric models how about defining problem as optimization task and fitting biggest blob or sphere? If I am not mistaken MichaelE2 has related nice solution for 2d circle. $\endgroup$
    – s.s.o
    Aug 27, 2014 at 13:29

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.