Object separation from a 2.5D surface

Question

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

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]

{pl, vl} = 
  First@Cases[g, 
    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]] &, 
  Range[Length[pl]]]
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.

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]

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.

Answer 1

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]]]} &, 
   data3Ds];
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}}

Show[projection, 
 Graphics[{Red, Thick, 
   Circle[#[[1]], #[[2]]] & /@ boulders[[All, 2]]}]]

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

Answer 2

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.