Convert binary voxel image to geometric region

Question

I have a binary voxel image which described a contiguous region in space.

Example:

arr = Array[
   Function[{x, y, z}, (1 + Sin[5 x] Sin[5 y] Sin[5 z]/2) - (x^2 + y^2 + z^2)],
    {100, 100, 100}, {{-1.25, 1.25}, {-1.25, 1.25}, {-1.25, 1.25}}];

img = Image3D[UnitStep[arr], "Bit"]

I need to convert this into a geometric region (MeshRegion?) that describes the same solid volume.

Can this be done in Mathematica 10.2?

I already found out that if I were doing it from scratch, I could use the marching cubes algorithm. If this is not built in (but maybe it is?), it would have to be done in C I assume. Though it's not clear if this would get me there all the way, or whether it only creates a surface mesh without information about the solid volume.

My actual data is larger (256 by 256 by 128), messier, and even less convex.

Answer 1

You could use ListContourPlot3D and DiscretizeGraphics:

Quiet @ DiscretizeGraphics @ Normal @ ListContourPlot3D[arr, Contours -> {0}, Mesh -> None]

Answer 2

First, here is how I understand your question:

How do I take a region described in a voxel array (binary 3D array), and convert it to a polygonal mesh (obtain a the boundary of the region as a list of polygons) in Mathematica?

Here is my code that takes a binary voxel array and converts it to a list of polygons / Graphics3D object.

Note that this allows you to export voxel graphics to STL (and other 3D formats) in Mathematica, which Raster3D won't do.

g3d[L_]:=Graphics3D[Polygon[L]];
Move[L_,v_]:=Map[#+v&,L];
Dirs={{1,0,0},{0,1,0},{0,0,1},{-1,0,0},{0,-1,0},{0,0,-1}};
CubeFaces = {{{1,1,1},{1,-1,1},{1,-1,-1},{1,1,-1}},{{1,1,1},{1,1,-1},{-1,1,-1},{-1,1,1}},{{1,1,1},{-1,1,1},{-1,-1,1},{1,-1,1}},{{-1,1,1},{-1,1,-1},{-1,-1,-1},{-1,-1,1}},{{1,-1,1},{-1,-1,1},{-1,-1,-1},{1,-1,-1}},{{1,1,-1},{1,-1,-1},{-1,-1,-1},{-1,1,-1}}};
InBounds[v_,dim_]:=Fold[#1&&#2&,True,MapThread[#1>0 && #1<=#2&,{v,dim}]];
NotEmpty[A_,pos_]:=Extract[A,pos]>0;
HasNeighbor[A_,pos_,di_]:= InBounds[pos+Dirs[[di]],Dimensions[A]]&&NotEmpty[A,pos+Dirs[[di]]];
GetFace[A_,pos_,di_]:=If[NotEmpty[A,pos] &&!HasNeighbor[A,pos,di],{Move[CubeFaces[[di]],2pos]},{}]
Voxelize[A_]:=Flatten[Table[ GetFace[A,{i,j,k},l],{i,Dimensions[A][[1]]},{j,Dimensions[A][[2]]},{k,Dimensions[A][[3]]},{l,1,6}],4];

Here is how to use it on your example:

data=UnitStep[Array[Function[{x,y,z},(1+Sin[5 x] Sin[5 y] Sin[5 z]/2)-(x^2+y^2+z^2)],{30,30,30},{{-1.25,1.25},{-1.25,1.25},{-1.25,1.25}}]];
G=g3d[Voxelize[data]]afÞfÞ C@`AÐH–@ÒäP#™Pg›+4@
Export["data.STL",G];

Here is the output (the 3D binary array displayed as voxels):

As a useful side effect, this view is way snappier in Mathematica than Image3D /Raster3D.

The last line exports the 3D polygonal data to STL. You can then use the STL file, say, for 3D-printing your voxel data from Mathematica.

I already found out that if I were doing it from scratch, I could use the marching cubes algorithm.

You can, but this might be an overkill in many cases. Also, it kills the voxels. Not useful if your voxel data is exact (e.g. a rendering of the Menger sponge).

Hope this helps!

Answer 3

As of version 11.0.0, the built-in function ImageMesh may be used