10
$\begingroup$

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$.

$\endgroup$
  • $\begingroup$ Interesting problem. Could you explain the context of the problem? Just curiosity? $\endgroup$ – Dr. belisarius Oct 12 '12 at 1:37
  • $\begingroup$ How many spheres? $\endgroup$ – Mark McClure Oct 12 '12 at 1:38
  • $\begingroup$ I'm trying to solve a packing density problem with mathematica. Eventually, I want to fill the cube to a volume fraction of 1/3 of the spheres. Right now I'm struggling to plot the cube and hard spheres. $\endgroup$ – Jen Oct 12 '12 at 1:43
  • $\begingroup$ So you are just trying to plot right now? You can do cube = {Opacity[.3], Cuboid[{0, 0, 0}, {20, 20, 20}]}; spheres = Table[{Black, Sphere[RandomReal[{.35, 20 - .35}, 3], .7]}, {j, 1, 10}]; Graphics3D[{cube, spheres}] $\endgroup$ – David Slater Oct 12 '12 at 1:54
  • $\begingroup$ It's easier to pack them than to place them randomly, especially if you want to assume that they cannot sit still, suspended by the air. $\endgroup$ – DavidC Oct 12 '12 at 5:44
9
$\begingroup$

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]

1

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.

$\endgroup$
  • 1
    $\begingroup$ It's a good way to get the OP started. An extension would be to allow for wiggling the locations. One could also treat it as like an electromagnetics problem, for instance — think of the spheres as all having unit charge and reposition the spheres just as they would when another charged sphere is added to the existing bunch, subject to the constraint that none leave the box. $\endgroup$ – rm -rf Oct 12 '12 at 3:36
  • $\begingroup$ The alternate type approach is to randomly place spheres so that each exactly touches three others. Much tougher to implement but it will get you ultimately to higher density (50+%). $\endgroup$ – george2079 Oct 12 '12 at 4:01
  • 1
    $\begingroup$ A relatively simple scheme could be to select a random sphere in the list, then creating a new sphere that touches it in a random direction. This got me to 1856 spheres in a few minutes. $\endgroup$ – wxffles Oct 12 '12 at 4:54
  • $\begingroup$ Why is the standing assumption that the metric space is Euclidean (L2-norm)? $\endgroup$ – alancalvitti Oct 12 '12 at 6:22
6
$\begingroup$

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.

nonoverlapping spheres in a cube

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.

nonoverlapping spheres in a cube

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

nonoverlapping spheres in a cube


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.

nonoverlapping spheres in a cube

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

nonoverlapping spheres in a cube

$\endgroup$
  • $\begingroup$ Why do you need to check if the center AND the hull are inside the cube? $\endgroup$ – Dr. belisarius Oct 12 '12 at 12:14
  • 1
    $\begingroup$ @bel: I'm exploiting the fact that And[] does short-circuit evaluation; if it fails the easy test, then the other tests are not touched. On the other hand, I suppose I can merge those first two tests; wait... $\endgroup$ – J. M. is away Oct 12 '12 at 12:19
5
$\begingroup$

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.

densely packed spheres

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

densely packed spheres

$\endgroup$
  • $\begingroup$ being nit picky, if i understand this scheme correctly it seems you are introducing a certain non-randomness in that you have a layered structure with constant z planes. That may be important depending on what you want to do with the result. Perhaps throw in a small random z within each layer. Might need to check two layers below if you do that though. $\endgroup$ – george2079 Oct 15 '12 at 12:25
  • $\begingroup$ No @george, the heights at each layer (except the first) are not constant; that is, unless z = RandomReal[{j, j + 1} r] is doing something peculiar that you know and I don't. Yes, I check two layers below as well; that's what the If[j > 2,(* stuff *)] sections of the seeding and layer-building parts are for. $\endgroup$ – J. M. is away Oct 15 '12 at 12:30
  • $\begingroup$ sorry..in my head that was outside the loop. Now i'm not sure that seeding the layer is necessary, you could find the mean z for each layer and move up a fixed amount. In the end i expect that seeding is a small part of the time though. $\endgroup$ – george2079 Oct 15 '12 at 13:08
1
$\begingroup$

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:

packed cube

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.

$\endgroup$

protected by Community Dec 31 '18 at 20:35

Thank you for your interest in this question. Because it has attracted low-quality or spam answers that had to be removed, posting an answer now requires 10 reputation on this site (the association bonus does not count).

Would you like to answer one of these unanswered questions instead?

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