# Showing the empty space in a Graphics3D plot

I would like to be able to visualize the empty space in a Graphics3D plot. For example, if I plot hexagonally close packed spheres using something like:

Show[Graphics3D[Sphere[#, Sqrt[2]/2] & /@ Select[Tuples[{-1, 0, 1}, 3],
Mod[Total[#], 2] == 0 &], Boxed -> False]]


Then how would I show the space between the spheres? That is, I would like the space occupied by the spheres to be empty and the space between them to be filled in.

• The empty space between them. I just found out it is often called an inverse opal structure. See, for example, https://duckduckgo.com/?q=inverse+opal+structure Commented May 25, 2014 at 15:48
• I appreciate how the answers might do what the OP wants, but the question seems to make no sense. The spheres as closed surfaces divide the space into regions. Using even/odd crossing one can identify two regions, one of which is the space outside all the spheres. Its visualization would look like a solid block with holes inside it, which might be shown with opacity < 1. What does "between" the spheres mean -- inside the convex hull? Commented May 26, 2014 at 3:45

You may try this

eq = And @@ (Total[({x, y, z} - #)^2] > 1/2 & /@
Select[Tuples[{-1, 0, 1}, 3], Mod[Total[#], 2] == 0 &])
RegionPlot3D[eq, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, Mesh -> None,
PlotPoints -> 150]


Notice there are small holes at points of contact between the spheres.

You can also "bound" by a sphere instead of a cuboid, with

eq = And @@
Prepend[Total[({x, y, z} - #)^2] > 1/2 & /@
Select[Tuples[{-1, 0, 1}, 3], Mod[Total[#], 2] == 0 &],
x^2 + y^2 + z^2 < 1]


• This is very helpful. Thank you! I am going to play around with it a bit, but it does essentially what I need. Commented May 25, 2014 at 16:59

Using some tricks (for speed and show-off :)

b = Normal@LatticeData["FaceCenteredCubic", "Basis"]
l1 = Flatten[Table[{i, j, k}.b, {i, -1, 1}, {j, -1, 1}, {k, -1, 1}], 2];
f = Nearest[l1];

RegionPlot3D[
Length@f[{x, y, z}, {1, 1/N@Sqrt@2}] < 1,
{x, -1, 1}, {y, -1, 1}, {z, -1, 1},
Mesh -> None, PlotStyle -> Directive[Yellow, Opacity[0.5]],
PlotPoints -> 100]


I thought a cute idea would be to take the rectangle formed by four points of contact between the spheres and blow it out like a bubble:

It seemed a simple idea, but the code to put together the surface of the space bounded by the spheres and the planes through the points of contact was a little longer than I expected.

Graphics3D @ space[]


The edges of the "triangles" and "squares" with arcs for sides lie on planes. We can solve for these and sort the points (some of the edges will be reversed, and it's easier just to sort them).

### Code dump

Because I wanted to do the animation, I began saving bits of data. The function  spacesetup sets up the all the data the first time space[] is called. There is an option PlotPoints that controls how finely the spherical segment is subdivided into polygons. The data for the various properties may accessed with calls of the form space[prop, PlotPoints -> plotPoints].

ClearAll[space, spacesetup];
Options[space] = {PlotPoints -> 15};
space["Properties"] = {"fullsurface", "allpoints", "spherepoints",
"rectangle", "rectpoints", "allspheres"};

space[property_String: "fullsurface", opts : OptionsPattern[]] /;
MemberQ[space["Properties"], property] &&
spacesetup[OptionValue[space, opts, PlotPoints]] :=
space[property, PlotPoints -> OptionValue[PlotPoints]];

spacesetup[Automatic] := spacesetup[PlotPoints /. Options[space]];
spacesetup[plotPoints_] :=
Module[{get, vertices, pairs, edges, rectpoints, p, q, r, s, spacexf},
get[prop_] := space[prop, PlotPoints -> plotPoints];
vertices = DeleteCases[
Select[Tuples[{-1, 0, 1}, 3], Mod[Total[#], 2] == 0 &],
{0, 0, 0}];
pairs = Subsets[vertices, {2}];
edges = Pick[pairs, #, Min[#]] &@ N[EuclideanDistance @@@ pairs];
{p, q, r, s} = # - First@vertices & /@
Mean /@ Select[edges, MemberQ[#, First@vertices] &][[{1, 2, 4, 3}]];
rectpoints =
N @ Flatten[
Table[(1 - u) p + u q + (1 - v) p + v s - p,
{u, 0, 1, 1/(plotPoints - 1)}, {v, 0, 1, 1/(plotPoints - 1)}],
1];
space["spherepoints", PlotPoints -> plotPoints] =
First@vertices + # & /@ (Normalize /@ rectpoints/Sqrt[2]);
space["rectpoints", PlotPoints -> plotPoints] =
First@vertices + # & /@ rectpoints;

spacexf = Flatten[Table[
RotationMatrix[r2, {0, 1, 0}].RotationMatrix[r1, {1, 0, 0}],
{r1, 0, Pi, Pi/2}, {r2, 0, 3 Pi/2, Pi/2}],
1];
space["allpoints", PlotPoints -> plotPoints] = Flatten[
Transpose[
spacexf.Transpose@get["spherepoints"],
{1, 3, 2}], 1];
space["allnormals", PlotPoints -> plotPoints] = Flatten[
Transpose[
spacexf.Transpose[-(First@vertices - # & /@ (get["spherepoints"]))],
{1, 3, 2}], 1];

space["allspheres", PlotPoints -> plotPoints] = Polygon[Flatten[
Table[plotPoints^2 k + {
{i + 1 + plotPoints*(j - 1),  i + plotPoints*(j - 1),  i + plotPoints*j},
{i + 1 + plotPoints*j,    i + 1 + plotPoints*(j - 1),  i + plotPoints*j}},
{k, 0, 11}, {i, plotPoints - 1}, {j, plotPoints - 1}],
3]];
space["alltriangles", PlotPoints -> plotPoints] = Polygon[
Map[
Function[pat,
With[{proj = Orthogonalize[Prepend[IdentityMatrix[3], pat]][[2 ;; 3]]},
SortBy[Flatten @ Position[get["allpoints"], _?(#.pat >= 2 &)],
ArcTan @@ (proj.get["allpoints"][[#]]) &]]
],
Tuples[{-1, 1}, 3]
],
VertexNormals -> None];
space["allsquares", PlotPoints -> plotPoints] = Polygon[
Map[
Function[pat,
SortBy[Flatten @ Position[get["allpoints"], pat],
ArcTan @@
get["allpoints"][[#,
Flatten @ Position[pat, Except[1. | -1.], Heads -> False]]] &]],
(Join[#, -#] &@N@IdentityMatrix[3] /. 0. -> _)
],
VertexNormals -> None];
space["fullsurface", PlotPoints -> plotPoints] =
GraphicsComplex[get["allpoints"],
{EdgeForm[],
get["allspheres"],
get["alltriangles"],
get["allsquares"]
},
VertexNormals -> get["allnormals"]
];

spacesetup[plotPoints] = True  (* Unset to reset *)
];


For the animation:

vertices = DeleteCases[
Select[Tuples[{-1, 0, 1}, 3], Mod[Total[#], 2] == 0 &],
{0, 0, 0}];
With[{plotPoints = 15},
anim = Table[
Graphics3D[{EdgeForm[],
{GraphicsComplex[(1 - t) space["rectpoints"] + t space["spherepoints"],
space["allspheres"][[All, 1 ;; 2*(plotPoints-1)^2]]
],
Red, Point[vertices]}
}, ViewPoint -> {-3, 1, 0}/2],
{t, 0, 1, 1/15}
]
];

Export["Example.gif", anim]
`