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?
Here is a subset of my data in .MAT format.
Answer
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} *)




Comments
Post a Comment