Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

EDIT: (my conclusion and thank you note) I want to thank you all guys for this unexpected intellectual and artistic journey. Hope you had fun and enjoyed it the same as I did.


I would like to generate a circle pack that mimics this: (don't pay attention on numbers, and colors at all, at the moment I am interested in circle positions and radii only)

enter image description here

or this:

enter image description here

I am new in Mathematica, could you give me some guidance? Thnx.

EDIT: The question was strictly for planar case (and remains so), however I see @Jacob Akkerboom in his answer added a solution for 3D generalization (thanks!), and, speaking of that, I just want here to bring to your attention this picture:

enter image description here

EDIT 2: There are some applications of circle packing in irregular shapes, like this: (author Jerome St Claire, handpainted)

enter image description here

... and a font called Dotted: (author Maggie Janssen)

enter image description here

... and some logos:

enter image description here

... and these hypnotic images: (from percolatorapp)

enter image description here

enter image description here

enter image description here

share|improve this question
4  
You can take a look here and here as a start. –  Öskå Jan 13 at 13:00
7  
Somewhat related is a word cloud. Functions used there could be easily adapted to disks. –  István Zachar Jan 13 at 13:27
1  
randomly place the largest circles, then place successively smaller ones where they can fit. Its quite straigntforward, and really this feels like a "give me some code" question.. –  george2079 Jan 13 at 13:49
4  
Thanks for links above, I'll try to find some ideas there. @george2079, no, its not that straightforward. I have a version in another language that is my test bed, it produces a solid output, but I want a little intellectual stimulation from other people here, in order to create a working and beautiful solution at the end. –  VividD Jan 13 at 14:09
    
I'm not sure the irregular shapes and fonts are not hand-made... –  István Zachar Jan 14 at 8:41

6 Answers 6

up vote 42 down vote accepted

replacing RandomReal function in István's code with

u = RandomVariate[UniformDistribution[{0,1 - ((1 - 2 min)/(max - min) (r - min) + 2 min)}]]

leads to non-uniform distribution enter image description here

Randomization for the angle can also be non-uniform:

randomPoint = 
  Compile[{{r, _Real}}, 
   Module[{u = 
      RandomVariate[
       UniformDistribution[{0, 
         1 - (-((1 - 2 min)/(max - min)) (r - min) + 1)}]], 
     a = RandomVariate[
       UniformDistribution[{π/(max - min)^(1/10) (r - min)^(1/10),
          2 π - π/(max - min)^(1/10) (r - min)^(
           1/10)}]]}, {Sqrt@u*Cos[a 2 Pi], Sqrt@u*Sin[a 2 Pi]}], 
   Parallelization -> True, CompilationTarget -> "C", 
   RuntimeOptions -> "Speed"];

enter image description here

The same applies to color. I think that after playing for long enough with these distributions you may even get some beautiful shapes. The real challenge would be to build new compilable distributions based on some graphics, like the figures in your example, or even some edge-detected pictures.

Edit (thanks to Simon Woods). The same idea may be implemented much easier using Simon's approach. We just have to make the radius choice dependent on the distance to the border. Inside the main loop replace the definition of r:

r = Min[max, d, m Exp[-(d/m)^0.2]]

This way the code respects fine details of the shape. You can see the the elephant's tail is drawn in small circles, which is common sense.

enter image description here

And it takes about 40 seconds to render all the zigs and zags of Norway's shoreline (set imagesizeto 500, max=10, min=0.5, pad=0.2).

enter image description here

Further, changing Simon's definition of m by adding a background value we can create distinguishable shapes in a pool of small circles:

distance = 
  Compile[{{pt, _Real, 1}, {centers, _Real, 2}, {radii, _Real, 1}}, 
    (Sqrt[Abs[First@# - First@pt]^2 + Abs[Last@# - Last@pt]^2] & /@ centers) - radii,
      Parallelization -> True, CompilationTarget -> "C", RuntimeOptions -> "Speed"];

max = 20;(*largest disk radius*)
min = 0.5;(*smallest disk radius *)
pad = 1;(*padding between disks*)
color = ColorData["DeepSeaColors"];
timeconstraint = 10;
shape = Binarize@ColorNegate@ImageCrop@Rasterize@Style["A", FontSize -> 1000];
centers = radii = {};
Module[{dim, dt, pt, m, d, r},
 dim = ImageDimensions[shape];
 dt = DistanceTransform[shape];
 TimeConstrained[While[True,
   While[
    While[
     pt = RandomReal[{1, #}] & /@ dim;
     (m = 3 + ImageValue[dt, pt]) < min];
    (d = Min[distance[pt, centers, radii]] - pad) < min];
   r = Min[max, d, m];
   centers = Join[centers, {pt}];
   radii = Join[radii, {r}]
   ], timeconstraint]]

Graphics[{color@RandomReal[], #} & /@ MapThread[Disk, {centers, radii}]]

enter image description here

And after that we can finally get to coloring (again, this is a modification of Simon's code):

distance = Compile[{{pt, _Real, 1}, {centers, _Real, 2}, {radii, _Real,1}}, (Sqrt[Abs[First@# - First@pt]^2 + Abs[Last@# - Last@pt]^2] & /@centers) - radii, Parallelization -> True, CompilationTarget -> "C", RuntimeOptions -> "Speed"];
max = 8;(*largest disk radius*)
min = 2;(*smallest disk radius*)
pad = 1.5;(*padding between disks*)
color1 = ColorData["Aquamarine"];
color2 = ColorData["SunsetColors"];
timeconstraint = 10;
background = 7;

shape = Binarize@ColorNegate@Rasterize@Style["74", Bold, Italic, FontFamily -> "PT Serif", FontSize -> 250];
centers = radii = colors = {};
Module[{dim, dt, pt, m, d, r}, dim = ImageDimensions[shape];
 dt = DistanceTransform[shape];
 TimeConstrained[
 While[True, While[While[pt = RandomReal[{1, #}] & /@ (2 dim);
    (m = If[Norm[pt - dim] < 200, background, 0] + If[pt[[1]] < dim[[1]] 3/2 && pt[[1]] > dim[[1]]/2 && pt[[2]] < dim[[2]] 3/2 && pt[[2]] > dim[[2]]/2, ImageValue[dt, pt - dim/2], 0]) < min];
    (d = Min[distance[pt, centers, radii]] - pad) < min];
    r = Min[max, d, m ];
    centers = Join[centers, {pt}];
    radii = Join[radii, {r}];
    colors =Join[colors, {Blend[{color2@RandomReal[{0.4, 0.7}], color1@RandomReal[{0.4, 0.7}]}, Piecewise[{{1/max*(m - background), m < background + max/2}, {1, m >= background + max/2}}]]}]];, timeconstraint]]

Graphics[MapThread[{#1, Disk[#2, #3]} &, {colors, centers, radii}]]

enter image description here

share|improve this answer
    
I like the first one, it reminds me on a scene from movie "The Lawnmower Man". –  VividD Jan 14 at 17:43
    
Wow, for A in a sea of small circles! –  VividD Jan 15 at 11:21
1  
Now make the background of the 74 a Disk and you win :) –  Öskå Jan 15 at 14:07
2  
@Öskå you got it, Sir! –  level1807 Jan 15 at 14:50

This methods relies on generating random circles and then removing circles that overlap with circles that were found earlier.

I suppose one should really divide the surface into bins and only check for overlaps between subsets of the circles. Especially if there is an upper bound to the size of a circle (the bound in my code is 1, which is not practical). This would be an improvement.

Here is the code

nn = 1*^5;
randRsPrefilter = RandomReal[1, nn];

randRadiiPrefilter = RandomVariate[BetaDistribution[8, 200], nn];

filter = Compile[{{radii1, _Real, 1}, {radii2, _Real, 1}}, 
   Block[{len, remaining}, len = Length@radii1;
    remaining = ConstantArray[1, len];
    Do[If[radii1[[i]] + radii2[[i]] > 1., remaining[[i]] = 0], {i, 
      len}];
    remaining]];

filt = filter[randRsPrefilter, randRadiiPrefilter];

randRs = Pick[randRsPrefilter, filt, 1];
randRadii = Pick[randRadiiPrefilter, filt, 1];

randAngles = RandomReal[2 Pi, Length@randRs];

toCoords = 
  Compile[{{l1, _Real, 1}, {l2, _Real, 1}}, 
   Table[l1[[i]] {Cos[l2[[i]]], Sin[l2[[i]]]}, {i, Length@l1}]];

coords = toCoords[randRs, randAngles];

overlapFilter = 
  Compile[{{coords, _Real, 2}, {radii, _Real, 1}, {start, _Integer}}, 
   Block[{res, curr, len, remaining, test, j}, len = Length@coords;
    remaining = ConstantArray[1, len];
    j = 1;
    test = True;
    Do[If[remaining[[i]] == 1, test = True;
      j = Max[start, i + 1];
      While[test, 
       If[Sqrt[Total[(coords[[i]] - coords[[j]])^2]] < 
         radii[[i]] + radii[[j]] + 0.01, remaining[[j]] = 0];
       If[j == len, test = False, j++;]]], {i, 1, len - 1}];
    remaining], CompilationTarget -> "C"];

overlapFilt = overlapFilter[coords, randRadii, 1];

goodCoords = Pick[coords, overlapFilt, 1];
goodRadii = Pick[randRadii, overlapFilt, 1];

goodPairs = Transpose[{goodCoords, goodRadii}];

Graphics[Disk @@@ goodPairs]

Output

enter image description here

3D version output

This only requires a slight a slight modification of the code, one only has to convert spherical coordinates to euclidean, rather than polar to euclidean. The distance function used in the overlapFilter function is sufficiently abstracted to deal with this.

enter image description here

share|improve this answer
    
I like this method, because it is somehow fair. I had some distribution in mind and I didn't cheat by putting new circles between old ones in an artificial way. But I guess that has it's cost, Count[overlapFilt[[-3000 ;;]], 1] gives 2, which means that in the end almost all circles get rejected. –  Jacob Akkerboom Jan 13 at 16:26
    
I like your method too! This will help me a lot, since I am new in Mathematica. The algorithm is short, simple, easy to modify, perfect for me. Thanks! –  VividD Jan 13 at 16:31
    
@VividD Please do not edit answers to ask a question or leave a comment. Writing it in this space (as a comment) is the right thing to do. –  rm -rf Jan 13 at 20:03
    
@rm -rf OK, no problem... :) –  VividD Jan 13 at 20:03

A simple algorithm that measures the distance of existing disks from a new, candidate disk, while decreasing radius size.

The following two functions generate a random point in the unit disk and measures the distance to all existing disks.

randomPoint = Compile[{{r, _Real}}, Module[
   {u = RandomReal@{0, 1 - 2 r}, a = RandomReal@{0, 2 Pi}},
   {Sqrt@u*Cos[a 2 Pi], Sqrt@u*Sin[a 2 Pi]}],
   Parallelization -> True, CompilationTarget -> "C", RuntimeOptions -> "Speed"];
distance = Compile[{{pt, _Real, 1}, {centers, _Real, 2}, {radii, _Real, 1}},
   (Sqrt[Abs[First@# - First@pt]^2 + Abs[Last@# - Last@pt]^2] & /@ centers) - radii,
   Parallelization -> True, CompilationTarget -> "C", RuntimeOptions -> "Speed"];


max = .08; (* largest disk radius *)
min = 0.005; (* smallest disk radius and step size *)
pad = 0.005; (* padding between disks *)
tolerance = 1000; (* wait this many rejections before decreasing radius *)
color = ColorData["BlueGreenYellow"];

centers = radii = {};
Do[failed = 0;
 While[failed < tolerance,
  pt = randomPoint@r;
  dist = distance[pt, centers, radii];
  If[Min@dist > r + pad,
   centers = Join[centers, {pt}];
   radii = Join[radii, {r}];,
   failed++;
   ]];, {r, max, min, -min}]

Graphics[{color@RandomReal[], #} & /@ MapThread[Disk, {centers, radii}],
  AspectRatio -> 1, Frame -> False, PlotRange -> {{-1, 1}, {-1, 1}},
  PlotRangePadding -> Scaled@.02, Axes -> False, ImageSize -> 400]   

Mathematica graphics

By increasing tolerance one can achieve more dense packings, using of course more time. With various min/max radius values and paddings, I got the following packings:

Mathematica graphics


Concerning other, possibly irregular shapes

Since OP requested other shapes, here is my solution for any, possibly irregular polygon. While the distance function remains intact, this approach requires a new randomPoint function that draws random points from the $(x, y)$ range of the polygon coordinates from inside the shape (thanks to rm -rf).

The function randomPoint expects a single polygon (with no holes) or a list of polygons (where the first is the outer boundary shape and the rest are the holes):

randomPoint[r_, {in_Polygon, ex___Polygon}] := Module[{p, range},
   range = {Min@#, Max@#} & /@ Transpose@First@in;
   While[(p = RandomReal /@ range; Not@And[
       Graphics`Mesh`InPolygonQ[in, p], 
       And @@ (Not@Graphics`Mesh`InPolygonQ[#, p] & /@ {ex})]
     ),]; p];
randomPoint[r_, poly_Polygon] := randomPoint[r, {poly}];

max = 1; (* largest disk radius *)
d = 0.01; (* smallest disk radius and step size *)
pad = 0; (* padding between disks *)
tolerance = 300; (* wait for this many rejections before decreasing radius *)
color = ColorData@"BlueGreenYellow";
shape = Polygon@N@(First@First@CountryData["Australia", "SchematicPolygon"]);

centers = radii = {};
Do[failed = 0;
 While[failed < tolerance,
  pt = randomPoint[r, shape];
  dist = distance[pt, centers, radii];
  If[Min@dist > r + pad,
   centers = Join[centers, {pt}];
   radii = Join[radii, {r}];,
   failed++;
   ]];, {r, max, d, -d}];

{
 Graphics[{EdgeForm@Gray, FaceForm@White, shape}, ImageSize -> 300],
 Graphics[{color@RandomReal[], #} & /@ MapThread[Disk, {centers, radii}],
    ImageSize -> 300]}

Mathematica graphics

For shapes with holes, I used Szabolcs's conversion to polygons:

shape = Block[{fun, g, xmin, xmax, ymin, ymax},
   fun = ListInterpolation@
     Rasterize[Style[Rotate["β", -Pi/2], FontSize -> 24, 
        FontFamily -> "Times"], "Data", ImageSize -> 300][[All, All, 1]];
   {{xmin, xmax}, {ymin, ymax}} = fun@"Domain";
   g = RegionPlot[fun[x, y] < 128, {x, xmin, xmax}, {y, ymin, ymax}, 
     PlotPoints -> 50, AspectRatio -> Automatic];
   Cases[Normal@g, Line[x___] :> Polygon@x, Infinity]
   ];

Result with {max = 3, d = .05} is:

Mathematica graphics

share|improve this answer
    
Thanks! It looks great, it gives natural, even, feel. However, I believe it would look even better if there is always at least a small padding (2 or 3 pixels) between two circles. –  VividD Jan 13 at 17:27
1  
@VividD You can easily have that by increasing pad. –  István Zachar Jan 13 at 17:42
    
Beautiful! Will check it out! –  VividD Jan 13 at 17:48
    
Can you also implement a region function so that you can pack the disks into any shape? –  rm -rf Jan 13 at 18:06
1  
Why not use a built-in inPolyQ from here: How to check if a 2D point is in a polygon? ? :) –  rm -rf Jan 14 at 21:16

Here's a pretty general way of filling an arbitrary shape with circle packs that starts with an arbitrary black and white image.

Start with a grid of points and perturb them (the size of the perturbation will dictate the variability in the sizes of the circles). Then find how far it is to the nearest point -- each point will then be grown into a disk with diameter equal to this distance. First for a square region:

n = 10;
tab = Flatten[Table[{i, j}, {i, -n, n}, {j, -n, n}], 1];
pts = tab + RandomReal[{-0.3, 0.3}, {Length[tab], 2}];
nf[x_] := Nearest[pts, x, 2];
radii = EuclideanDistance[nf[pts[[#]]][[1]], nf[pts[[#]]][[2]]] & /@ 
    Range[Length[pts]]/2;
Graphics[Table[{RGBColor[RandomReal[], RandomReal[], RandomReal[]], 
   Disk[pts[[i]], radii[[i]]]}, {i, 1, Length[pts]}]]

enter image description here

To control the shape, begin with an image that is white wherever we wish the circles to be. For example, consider the horse

img = Import["http://i.stack.imgur.com/H6OUl.png"]

enter image description here

It is easy to apply the above method to a simplified (downsampled) version of the horse

imgSimp = Downsample[ImageData[img], 10];
imgDim = Dimensions[imgSimp];
tab = Flatten[Table[{i, j}, {i, 1, imgDim[[1]]}, {j, 1, imgDim[[2]]}],
    1];
mask = Flatten[Partition[tab, imgDim[[2]]] img, 1];
pts = mask + RandomReal[{-0.3, 0.3}, {Length[mask], 2}];
nf[x_] := Nearest[pts, x, 2];
radii = EuclideanDistance[nf[pts[[#]]][[1]], nf[pts[[#]]][[2]]] & /@ 
    Range[Length[pts]]/2;
Rotate[Graphics[
  Table[{RGBColor[g = RandomReal[], g, g], 
    Disk[pts[[i]], radii[[i]]]}, {i, 1, Length[pts]}]], -Pi/2]

enter image description here

share|improve this answer

The idea is to

  1. Choose a random white spot on a binarized image.
  2. Try to fit as large circle as possible (of the givens ones) in this position.
  3. If no circle fits at this position, choose a new random position. If there are n positions in a row for which you cannot fit a circle then terminate the process.

Instead of working with graphics primitives and the position of circles and their radii I work with the image matrix directly.

Code:

padding = 1;
circles = {Position[DiskMatrix[# + padding], 1] - # - padding, 
     Position[DiskMatrix[#], 1] - #} & /@ {15, 12, 10, 6, 3, 1};

shape = DiskMatrix[250];
space = Position[shape, 1];
i = 0;

While[i < 1000,
 pt = RandomChoice[space]; placed = False;
 Do[
   occupied = pt + # & /@ c[[1]];
   If[Length@occupied == Length@Intersection[space, occupied],
    space = Complement[space, occupied];
    shape = 
     ReplacePart[shape, 
      pt + # & /@ c[[2]] -> ColorData["BlueGreenYellow"]@RandomReal[] /. 
       RGBColor[r_, g_, b_] :> {r, g, b}];
    placed = True; i = 0; Break[]
    ], {c, circles}]
  If[! placed, i++];
 ]

shape = ReplacePart[shape, Position[shape, 1, {2}] -> {1, 1, 1}];
shape = ReplacePart[shape, Position[shape, 0, {2}] -> {1, 1, 1}];

shape // Image

Adjustable parameters are the paddings, the size of the circles and for how long it continues to try to pack the shape.

With the color scheme BlueGreenYellow and shape = DiskMatrix[250] we get

demo1

With the color scheme DeepSeaColors and

shape = ImageData@ColorNegate@ImageCrop@Binarize@Rasterize@Style["A", FontSize -> 1000];

we get

demo2

Finally, this is a circle packed map of Sweden using the color scheme DarkTerrain:

shape = ImageData@ColorNegate@Binarize[Rasterize@Show[CountryData["Sweden", "Shape"], ImageSize -> 200],0.99];

demo3

The careful observer will note that the smallest objects are not actually round. Don't worry about this, it's because the smallest object is just one pixel and you can't make a circle out of that. It saved me time to generate the graphics like this, it can easily be fixed by making a larger image and setting the smallest circle to a radius of say three or five, then shrinking the image to whatever size one wants.

share|improve this answer
    
Thanks, amazing!! What I would also like to see is Norway, detailed map, I am interested in how your method handle fjords. –  VividD Jan 14 at 21:23
    
@VividD The code is available, you can try it out. Just change "Sweden" to "Norway." The larger the image the more detail you can capture, afterwards you can scale it down. –  Pickett Jan 14 at 21:30
    
@VividD see the edit in my post. Seems to do Norway just fine! –  level1807 Jan 15 at 9:51

Here's another shape-packing one, with a binary image used to define the shape to be filled. I use a DistanceTransform on the image, which provides a convenient way to measure the distance from any point to the boundary of the shape.

I've used Istvan's distance function, but instead of choosing a spot size and then locating somewhere to put it, I choose a location and then determine the spot size subject to the constraints.

The packing continues for a fixed time, using TimeConstrained - the longer you allow the code to run the more densely packed the shape will be.

distance = 
  Compile[{{pt, _Real, 1}, {centers, _Real, 2}, {radii, _Real, 1}}, 
    (Sqrt[Abs[First@# - First@pt]^2 + Abs[Last@# - Last@pt]^2] & /@ centers) - radii,
      Parallelization -> True, CompilationTarget -> "C", RuntimeOptions -> "Speed"];

max = 5;(*largest disk radius*)
min = 1;(*smallest disk radius *)
pad = 1;(*padding between disks*)
color = ColorData["CandyColors"];
timeconstraint = 10;

shape = Binarize@ColorNegate@Import["http://i.stack.imgur.com/wtJoA.png"]

enter image description here

centers = radii = {};
Module[{dim, dt, pt, m, d, r},
 dim = ImageDimensions[shape];
 dt = DistanceTransform[shape];
 TimeConstrained[While[True,
   While[
    While[
     pt = RandomReal[{1, #}] & /@ dim;
     (m = ImageValue[dt, pt]) < min];
    (d = Min[distance[pt, centers, radii]] - pad) < min];
   r = Min[max, d, m];
   centers = Join[centers, {pt}];
   radii = Join[radii, {r}]
   ], timeconstraint]]

Graphics[{color@RandomReal[], #} & /@ MapThread[Disk, {centers, radii}]]

enter image description here

share|improve this answer
6  
This one is beautiful. Contours of the original are clearly visible. –  VividD Jan 14 at 23:11

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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