A method to draw a surface around a set of regularly spaced points

This is mostly an artistic endeavour.

I have a 3D array of points, 1 representing a point within the surface, -1 without. The shape is not convex. I would like to produce a surface which encloses the points. One simple way is ListContourPlot:

The surface is not very pleasing though; it is similar to the result of simply building the shape out of cubes without using any interpolation. Perhaps as a result of this, it is also very large, taking up 2gb of RAM.

What is a better way of doing this?

• Oh, I literally posted the same question just a few minutes ago. I am stuck with ListContourPlot3D right now, but there might be some options I do not know of for different "intepolation" methods. Oct 21, 2015 at 17:09
• What if you use ListSurfacePlot3D[] instead? Oct 21, 2015 at 17:10
• Actually I believe this paper is more relevant: Lempitsky, "Surface Extraction from Binary Volumes with Higher-Order Smoothness" (2010). That said, I got a nice-looking plot from your sample data by doing ListContourPlot3D[Downsample[GaussianFilter[data, 5], 2], Contours -> {0}].
– user484
Oct 21, 2015 at 19:09
• I sure hope not. I'm trying to whip up some fancy graphics for a grant. Oct 21, 2015 at 19:28
• Is the binary data all you have? If it was, for example, created by binarizing a continuous scalar field, you would get better results by running ListContourPlot3D on the pre-binarized data.
– user484
Oct 22, 2015 at 3:46

So the algorithm in the paper I linked to in a comment, "Surface Extraction from Binary Volumes with Higher-Order Smoothness" by Lempitsky (2010), turned out to be pretty easy to implement (though for speed I changed eq. (10a) to a difference of Gaussians). And it works much better than my attempt, so I'm replacing that with this.

Build a signed distance field (SDF):

dOut = ImageClip[
ImageSubtract[DistanceTransform[Image3D[-data]], 0.5], {0, 1*^6}];
dIn = ImageClip[
ImageSubtract[DistanceTransform[Image3D[data]], 0.5], {0, 1*^6}];
sdf = ImageSubtract[dOut, dIn];

Define lower and upper bounds for the smoothed SDF:

l = ImageApply[Which[# >= 0, Max[# - 1, 0], True, -1*^6] &, sdf];
u = ImageApply[Which[# <= 0, Min[# + 1, 0], True, 1*^6] &, sdf];

Define the filtering operation:

filter[r_][sdf_] :=
ImageApply[
Clip[#1, {#2, #3}] &, {ImageSubtract[
ImageMultiply[GaussianFilter[sdf, r], 4/3],
ImageMultiply[GaussianFilter[sdf, 2 r], 1/3]], l, u}]

And that's it!

If you don't have much time, use a large radius and a handful of iterations. Otherwise, use a small radius and a large number of iterations for higher-quality results.

draw[sdf_] :=
ListContourPlot3D[ImageData[sdf], Contours -> {0},
ContourStyle -> White, Mesh -> None]
draw[sdf]
Print[draw[filter[4][sdf]]]; // Timing
Print[draw[Nest[filter[2], sdf, 10]]]; // Timing
Print[draw[Nest[filter[1.2], sdf, 100]]]; // Timing

(* {6.74007, Null} *)

(* {39.5372, Null} *)

(* {365.001, Null} *)
• @Rahul Your previous solution using a GaussianFilter has produced quite excellent results, although it is a touch slow. I'd like to accept it as an answer but it's just a comment. Oct 29, 2015 at 14:38
• @user5751: The reason I didn't recommend GaussianFilter in my answer is that it can completely lose narrow features: i.sstatic.net/S9M8m.png becomes i.sstatic.net/zCmLW.png
– user484
Oct 29, 2015 at 23:28
• @user5751: See my update (though it's even slower than a plain old GaussianFilter).
– user484
Oct 30, 2015 at 2:52
• I dunno, sometimes one just has to spend more time in exchange for higher quality. It's too bad I can't upvote again. Oct 30, 2015 at 3:24