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 contour between -1's and 1's in my data

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?

Here is a subset of my data in .MAT format.

  • $\begingroup$ 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. $\endgroup$ – Wizard Oct 21 '15 at 17:09
  • $\begingroup$ What if you use ListSurfacePlot3D[] instead? $\endgroup$ – J. M.'s torpor Oct 21 '15 at 17:10
  • 5
    $\begingroup$ 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}]. $\endgroup$ – user484 Oct 21 '15 at 19:09
  • 1
    $\begingroup$ I sure hope not. I'm trying to whip up some fancy graphics for a grant. $\endgroup$ – Crêpo Oct 21 '15 at 19:28
  • 1
    $\begingroup$ 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. $\endgroup$ – user484 Oct 22 '15 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_] := 
  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]
Print[draw[filter[4][sdf]]]; // Timing
Print[draw[Nest[filter[2], sdf, 10]]]; // Timing
Print[draw[Nest[filter[1.2], sdf, 100]]]; // Timing

enter image description here

enter image description here

(* {6.74007, Null} *)

enter image description here

(* {39.5372, Null} *)

enter image description here

(* {365.001, Null} *)
  • $\begingroup$ @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. $\endgroup$ – Crêpo Oct 29 '15 at 14:38
  • 3
    $\begingroup$ @user5751: The reason I didn't recommend GaussianFilter in my answer is that it can completely lose narrow features: i.stack.imgur.com/S9M8m.png becomes i.stack.imgur.com/zCmLW.png $\endgroup$ – user484 Oct 29 '15 at 23:28
  • $\begingroup$ @user5751: See my update (though it's even slower than a plain old GaussianFilter). $\endgroup$ – user484 Oct 30 '15 at 2:52
  • $\begingroup$ I dunno, sometimes one just has to spend more time in exchange for higher quality. It's too bad I can't upvote again. $\endgroup$ – J. M.'s torpor Oct 30 '15 at 3:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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