# Generating visually pleasing circle packs

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)

or this:

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:

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

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

... and some logos:

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

-
You can take a look here and here as a start. –  Öskå Jan 13 '14 at 13:00
Somewhat related is a word cloud. Functions used there could be easily adapted to disks. –  István Zachar Jan 13 '14 at 13:27
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 '14 at 13:49
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 '14 at 14:09
I'm not sure the irregular shapes and fonts are not hand-made... –  István Zachar Jan 14 '14 at 8:41

replacing RandomReal function in István's code with

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


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"];


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.

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

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"];

min = 0.5;(*smallest disk radius *)
color = ColorData["DeepSeaColors"];
timeconstraint = 10;
shape = Binarize@ColorNegate@ImageCrop@Rasterize@Style["A", FontSize -> 1000];
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];
r = Min[max, d, m];
centers = Join[centers, {pt}];
], timeconstraint]]



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"];
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];
r = Min[max, d, m ];
centers = Join[centers, {pt}];
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]]



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

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;
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

With the color scheme DeepSeaColors and

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


we get

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];


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.

-
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 '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 '14 at 21:30
@VividD see the edit in my post. Seems to do Norway just fine! –  level1807 Jan 15 '14 at 9:51

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 *)
tolerance = 1000; (* wait this many rejections before decreasing radius *)
color = ColorData["BlueGreenYellow"];

Do[failed = 0;
While[failed < tolerance,
pt = randomPoint@r;
centers = Join[centers, {pt}];
failed++;
]];, {r, max, min, -min}]

AspectRatio -> 1, Frame -> False, PlotRange -> {{-1, 1}, {-1, 1}},
PlotRangePadding -> Scaled@.02, Axes -> False, ImageSize -> 400]


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:

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[
GraphicsMeshInPolygonQ[in, p],
And @@ (Not@GraphicsMeshInPolygonQ[#, p] & /@ {ex})]
),]; p];
randomPoint[r_, poly_Polygon] := randomPoint[r, {poly}];

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

Do[failed = 0;
While[failed < tolerance,
pt = randomPoint[r, shape];
centers = Join[centers, {pt}];
failed++;
]];, {r, max, d, -d}];

{
Graphics[{EdgeForm@Gray, FaceForm@White, shape}, ImageSize -> 300],
ImageSize -> 300]}


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:

-
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 '14 at 17:27
@VividD You can easily have that by increasing pad. –  István Zachar Jan 13 '14 at 17:42
Beautiful! Will check it out! –  VividD Jan 13 '14 at 17:48
Can you also implement a region function so that you can pack the disks into any shape? –  rm -rf Jan 13 '14 at 18:06
Why not use a built-in inPolyQ from here: How to check if a 2D point is in a polygon? ? :) –  rm -rf Jan 14 '14 at 21:16

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"];

min = 1;(*smallest disk radius *)
color = ColorData["CandyColors"];
timeconstraint = 10;

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


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];
r = Min[max, d, m];
centers = Join[centers, {pt}];
], timeconstraint]]



-
This one is beautiful. Contours of the original are clearly visible. –  VividD Jan 14 '14 at 23:11

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[]],


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"]


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];
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]


-

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];

remaining = ConstantArray[1, len];
len}];
remaining]];

randRs = Pick[randRsPrefilter, 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]] <
If[j == len, test = False, j++;]]], {i, 1, len - 1}];
remaining], CompilationTarget -> "C"];

goodCoords = Pick[coords, overlapFilt, 1];

Graphics[Disk @@@ goodPairs]


Output

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.

-
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 '14 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 '14 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 '14 at 20:03
@rm -rf OK, no problem... :) –  VividD Jan 13 '14 at 20:03