Extracting polygons from 3D contour plot surface

Question

Edit: This question turned out to be two parts so I am going to make this question about only the first part a kguler provided an excellent answer.

Here is a better representation. My actual data consists of many points which I bin onto a grid.

data = RandomVariate[NormalDistribution[1, 3], {10^5, 3}];
rounded = Round@data;
gathered = Gather[rounded];
offset = -Min /@ Transpose@rounded + 1;
counts = Length /@ gathered;
(*I offset the data so I can use SparseArray from 1*)
coords = offset + # & /@ (First /@ gathered); 
sparesearray = SparseArray[coords -> counts];

Each grid contains the count of the number of points that were in the grid cube (voxel). I now create a contour plot.

ListContourPlot3D[sparesearray, Contours -> {2.5}]

In this case its ugly since its random data. However, in my real data it has discrete volumes. I want to find the points in the original data (or the offset data).

I could simply use the bins with more than n counts but I was hoping to take advantage of the interpolation that ListContourPlot3D uses.

How do I extract the polygons rendered by ListContourPlot3D?

early version below for posterity:

I have a set of data points in 3D and I want to extract the subset of points that is contained within the surface generated by a contour plot.

data = Table[x^2 + y^2 + z^2 + RandomReal[0.1], {x, -2, 2, 0.2}, {y, -2, 2, 
  0.2}, {z, -2, 2, 0.2}]

The sample data here is on a regular grid but my real data is not on a regular grid. It is defined: (x, y, z, value).

plot = ListContourPlot3D[data, Contours -> {1}, Mesh -> None]

InputForm returns a GraphicsComplex object but I'm not really sure how to determine the polygons in the surface and then how to determine which of the points are inside the surface and which are outside.

Answer 1

From the Documentation Center Normal > Details

Normal[expr] converts GraphicsComplex objects into ordinary lists of graphics primitives and directives.

So, for

  data = Table[x^2 + y^2 + z^2 + RandomReal[0.1], {x, -2, 2, 0.2}, 
    {y, -2, 2,  0.2}, {z, -2, 2, 0.2}];
  plot = ListContourPlot3D[data, Contours -> {1}, Mesh -> None,  ImageSize -> 400];

you can get the ordinary list of polygons using

  Cases[Normal[plot[[1]]], Polygon[__], {0, Infinity}]

and

  polygonCoords = 
  Cases[Normal[plot[[1]]], Polygon[x_, ___] :> x, {0, Infinity}]

gives the coordinates of the polygons

grphcs = Graphics3D[{Opacity[.9], EdgeForm[Opacity[.3]], 
Polygon[#,  VertexColors -> Table[Hue[RandomReal[]], {Length[#]}]] & /@ 
 Cases[Normal[plot[[1]]], Polygon[x_, ___] :> x, {0, Infinity}]},
Lighting -> "Neutral", ImageSize -> 400];
Row[{plot, grphcs}]

Update: OP's new example:

 data = RandomVariate[NormalDistribution[1, 3], {10^5, 3}];
 rounded = Round@data; gathered = Gather[rounded];
 offset = -Min /@ Transpose@rounded + 1;
 counts = Length /@ gathered;
 coords = offset + # & /@ (First /@ gathered);
 sparesearray = SparseArray[coords -> counts];

  plot2 = ListContourPlot3D[sparesearray, Contours -> {2.5},  ImageSize -> 400];
  grphcs2 = Graphics3D[{EdgeForm[Opacity[.3]],
   Polygon[#, VertexColors -> Table[Hue[RandomReal[]], {Length[#]}]] & /@ 
  Cases[Normal[plot2[[1]]], Polygon[x_, ___] :> x, {0, Infinity}]},
   Lighting -> "Neutral", ImageSize -> 400];
  Row[{plot2, grphcs2}]

Note: still remains open is the second part of OP's question (determining the sets of points inside and outside the surface).

Answer 2

Here is a straightforward alternative in Version 10 that works for your purpose, since you just want to extract the polygons on the surface. The approach is to discretize the graphics object and collect the Polygons:

data = Table[x^2 + y^2 + z^2 + RandomReal[0.1], {x, -2, 2, 0.2}, {y, -2, 2, 
             0.2}, {z, -2, 2, 0.2}];

gr = ListContourPlot3D[data, Contours -> {1}, Mesh -> None, ImageSize -> 400];

Now we discretize it and extract the polygons:

bdg = DiscretizeGraphics[Normal[plot /. (Lighting -> _) :> Lighting -> Automatic]];

The replacement trick I used above is only because of one of many bugs in both DiscretizeGraphics and BoundaryDiscretizeGraphics e.g. see here

Getting the polygons is as easy as doing:

polys = MeshPrimitives[bdg, 2];

We can plot them to be sure:

polygr = Graphics3D[Transpose[{RandomColor[MeshCellCount[bdg, 2]], polys}]];

Row[{gr, polygr}]

For the second case, procedure is exactly the same:

data = RandomVariate[NormalDistribution[1, 3], {10^5, 3}]; 
rounded = Round@data; 
gathered = Gather[rounded]; 
offset = -Min /@ Transpose@rounded + 1; 
counts = Length /@ gathered; 
coords = offset + # & /@ (First /@ gathered);
sparesearray = SparseArray[coords -> counts];

gr2 = ListContourPlot3D[sparesearray, Contours -> {2.5}, ImageSize -> 400];

bdg2 = DiscretizeGraphics[Normal[gr2 /. (Lighting -> _) :> Lighting -> Automatic]];
polys2 = MeshPrimitives[bdg2, 2];
polygr2 = Graphics3D[Transpose[{RandomColor[MeshCellCount[bdg2, 2]], polys2}]];

Row[{gr2, polygr2}]