Randomly packing spheres of fixed radius within a cube

Question

How can I have Mathematica randomly place spheres in a cube so they won't overlap? The cube is $20 \times 20 \times 20$, and the spheres have a radius of $0.7$.

Answer 1

Here's a brute force and soporifically slow method.

Start with an empty list, and something to monitor progress:

spheres = {};
Dynamic[Length[spheres]]

Then create a new random sphere, and if it's not too close to any existing sphere add it to the list:

While[Length[spheres] < 1856, s = RandomReal[{0.7, 20 - 0.7}, 3]; 
  If[And @@ (Norm[# - s] > 1.4 & /@ spheres), AppendTo[spheres, s]]];

According to my calculations, you need 1856 spheres of radius 0.7 to fill 1/3 of the volume of a 20 x 20 x 20 cube.

To display:

cube = {Opacity[0.3], Cuboid[{0, 0, 0}, {20, 20, 20}]};
Graphics3D[{cube, Sphere[#, 0.7] & /@ spheres}, Boxed -> False]

As I stated at the beginning, this is very slow and inefficient. The above image is from when I gave up after 1500 spheres. You can Alt. and resume at your leisure. It's entirely possible, what with the randomness, that you will never fit the 1856 spheres in.

Answer 2

Here's a brute-force method (based on a bounded random walk) that attempts to generate n nonintersecting spheres of radius r within a cube of side length s, as well as reporting the actual number of spheres it managed to generate after lim unsuccessful iterations:

sphrPak[s_, r_, n_Integer?Positive, lim_Integer?Positive] := 
  Module[{sh = s/2, k, p0, pc, sphList},

   (* generate initial sphere within box *)
   While[
         p0 = RandomReal[{-sh, sh}, 3];
         And @@ Thread[Map[Min[sh - #, sh + #] &, p0] < r]];
   sphList = {p0};

   k = 0;
   While[Length[sphList] < n && k < lim,
         (* center for new sphere chosen in random direction *)
         pc = p0 + 2 r Normalize[RandomVariate[NormalDistribution[], 3]];
         If[(* is sphere within the cube? *)
            PolyhedronData["Cube", "RegionFunction"] @@ 
            ScalingTransform[ConstantArray[1/(s - 2 r), 3]][pc] &&
            (* does the sphere not overlap with other spheres? *)
            (And @@ Thread[Map[EuclideanDistance[pc, #] &, sphList] >= 2 r]),
            k = 0; AppendTo[sphList, pc]; p0 = pc,
            (* else *) ++k];
         ];

   If[k == lim, Print[StringForm["Only `1` spheres were generated.", Length[sphList]]]];

   Graphics3D[Sphere[sphList, r], Axes -> Automatic,
              PlotRange -> {{-sh, sh}, {-sh, sh}, {-sh, sh}}]];

sphrPak[s_, r_, n_Integer] := sphrPak[s, r, n, Quotient[2 n, 3]]

Try it out:

BlockRandom[SeedRandom[1024, Method -> "MersenneTwister"]; (* for reproducibility *)
            sphrPak[20, 0.7, 1500]]
   Only 725 spheres were generated.

Of course, a different attempt might yield a result with more spheres within the cube.

---

More attempts:

BlockRandom[SeedRandom[4092, Method -> "MersenneTwister"];
            sphrPak[20, 0.7, 1500]]
   Only 805 spheres were generated.

BlockRandom[SeedRandom[2012, Method -> "MersenneTwister"];
            sphrPak[20, 0.7, 1500]]
   Only 932 spheres were generated.

---

Here's a variation that attempts to do more than one bounded random walk, where the starting point of the $n$th random walk is a randomly chosen point from the $(n-1)$-th random walk:

sphrPak2[s_, r_, n_Integer?Positive, walks_Integer?Positive, its_Integer?Positive] :=
Module[{sh = s/2, j, k, p0, pc, sphList},
       While[
             p0 = RandomReal[{-sh, sh}, 3];
             And @@ Thread[Map[Min[sh - #, sh + #] &, p0] < r]];
       sphList = {p0}; j = 0;
       While[Length[sphList] < n && j < walks,
             k = 0;
             While[Length[sphList] < n && k < its,
                   pc = p0 + 2 r Normalize[RandomVariate[NormalDistribution[], 3]];
                   If[(PolyhedronData["Cube", "RegionFunction"] @@
                      ScalingTransform[ConstantArray[1/(s - 2 r), 3]][pc]) &&
                      (And @@ Thread[Map[EuclideanDistance[pc, #] &, sphList] >= 2 r]),
                      k = 0; AppendTo[sphList, pc]; p0 = pc,
                      ++k];
                   ];
             If[Length[sphList] < n, 
                p0 = RandomChoice[Map[EuclideanDistance[#, sphList[[1]]] &,
                                      sphList] -> sphList]];
             ++j;
             ];
       If[j == walks, Print[StringForm["Only `1` spheres were generated.",
                                       Length[sphList]]]];
       Graphics3D[Sphere[sphList, r], Axes -> Automatic, 
                  PlotRange -> {{-sh, sh}, {-sh, sh}, {-sh, sh}}]]

Try it out:

BlockRandom[SeedRandom[85, Method -> "Rule50025CA"]; 
            sphrPak2[20, 0.7, 1500, 1, 1000]]
   Only 662 spheres were generated.

BlockRandom[SeedRandom[85, Method -> "Rule50025CA"]; 
            sphrPak2[20, 0.7, 1500, 20, 1000]]
   Only 1335 spheres were generated.

Answer 3

I will present another brute-force method for densely packing spheres that works by randomly depositing the spheres in successive "layers". The first layer is built by solving the two-dimensional analog of the problem (that is, packing disks in a square). Having constructed the first layer, successive layers are built up by dropping new spheres, checking that the incoming spheres do not intersect with spheres in previous layers (as well as with other spheres in the same layer). For the second layer, the new spheres are compared against spheres in the first layer, while for the third and succeeding layers, the new sphere is tested for intersections with the two previous layers. All this is done up until we hit the top of the cube.

Here is a Mathematica implementation of the idea presented above:

sphrPak[s_, r_, its_Integer: 2000] := 
 Module[{j, k, lap, lap2, layer, layers, p0, pc, sphrs, z},

  (* pack disks into square *)

  (* pick first disk *)
  While[
        p0 = RandomReal[{r, s - r}, 2];
        And @@ Thread[Map[Min[#, s - #] &, p0] < r]];

  (* successively pack disks into square *)
  layer = {p0}; k = 0;
  While[k < its,
        pc = RandomReal[{r, s - r}, 2]; (* candidate disk *)
        If[And @@ Thread[Map[EuclideanDistance[pc, #] &, layer] >= 2 r],
           k = 0; AppendTo[layer, pc];,
           ++k];
        ];

  (* build and store first layer *)
  layers = {PadRight[layer, {Length[layer], 3}, r]};

  (* build successive layers *)
  j = 2;
  While[True,

        (* previous layers *)
        lap = layers[[j - 1]]; If[j > 2, lap2 = layers[[j - 2]]];

        (* find "seed" sphere for layer *)
        p0 = {}; k = 0;
        While[k < its,
              z = RandomReal[{j, j + 1} r];
              pc = Append[RandomReal[{r, s - r}, 2], z]; ++k; (* candidate seed *)
              If[z <= s - r &&
                 (And @@ Thread[Map[EuclideanDistance[pc, #] &, lap] >= 2 r]),
                 If[j > 2,
                    If[And @@ Thread[Map[EuclideanDistance[pc, #] &, lap2] >= 2 r], 
                       p0 = pc; Break[]],
                    p0 = pc; Break[]]]];

        If[p0 === {}, Break[]]; (* no seed found *)

        (* start building new layer *)
        layer = {p0}; k = 0;
        While[k < its,
              z = RandomReal[{j, j + 1} r];
              pc = Append[RandomReal[{r, s - r}, 2], z]; (* candidate point *)
              If[z <= s - r &&
                 (And @@ Thread[Map[EuclideanDistance[pc, #] &, layer] >= 2 r]) &&
                 (And @@ Thread[Map[EuclideanDistance[pc, #] &, lap] >= 2 r]),
                 If[j > 2,
                    If[And @@ Thread[Map[EuclideanDistance[pc, #] &, lap2] >= 2 r],
                       k = 0; AppendTo[layer, pc], ++k],
                    k = 0; AppendTo[layer, pc]], ++k];
              ];
        AppendTo[layers, layer];
        ++j];
  Print[StringForm["`1` spheres generated.", Length[sphrs = Join @@ layers]]];
  Graphics3D[Sphere[sphrs, r], Axes -> Automatic, PlotRange -> {{0, s}, {0, s}, {0, s}}]]

Try it out:

BlockRandom[SeedRandom[2012, Method -> "MersenneTwister"]; (* for reproducibility *)
            sphrPak[20, 0.7]]
   1996 spheres generated.

BlockRandom[SeedRandom[1024, Method -> "MersenneTwister"];
            sphrPak[20, 0.7]]
   2041 spheres generated.

Answer 4

So using another random walk-like method (well I tried two which I realized are equivalent) I was able to do this pretty quickly and compactly:

generateRandomAnchoring[] :=
  
  With[{p = RandomReal[{r, m - r}, 3], c = RandomChoice@points, 
    d = (2 r)},
   c + d*Normalize@(c - p)
   ];
generateRandomWalking[] :=
  
  With[{p = RandomChoice[points], \[Theta] = 
     RandomReal[2 \[Pi]], \[Phi] = RandomReal[2 \[Pi]], d = (2 r)},
   p + {d*Cos[\[Phi]] Sin[\[Theta]], d*Sin[\[Theta]] Sin[\[Phi]], 
     d*Cos[\[Theta]]}
   ];
tryPoint[s_] := 
  If[AllTrue[s, r < # < m - r &] && 
    AllTrue[points, (Norm[# - s] >= 2 r) &], AppendTo[points, s]];

Both of these take a random pregenerated point and attach a point to it at a random position, just in one case it generates a spherical polar coordinate and in the other it just attaches a randomly generated point. They're entirely equivalent as best I can reason out.

Then we'll just initialize the state and update the point generating code from wxffles answer so that it also spits out the time elapsed:

r = .7;
m = 20;
points = {};
tryPoint@RandomReal[{r, m - r}, 3];

With[{now = Now},
 Monitor[While[Length[points] < 1825,
   tryPoint@generateRandomAnchoring[]
   ], Column@{Now - now, Length@points}]
 ]

It runs for about a minute before hitting the 1825 mark. And we'll use the same viewer code as wxffles too:

cube = {Opacity[0.3], Cuboid[{0, 0, 0}, {20, 20, 20}]};
Graphics3D[{cube, Sphere[#, 0.7] & /@ points}, Boxed -> False]

And this is what we see:

Then just adjust the number of spheres and timing to your liking and you have a your packed cube.

Edit

Note that one could also add a weight-tracking scheme for the random choice of points (in the second of the methods). I.e., points that have been chosen and for which the addition failed decrease in weight by some percentage. The likelihood that this will be a significant improvement is low, but it could eke out some performance gains.