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:

... garden design:

... infographics:

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

• You can take a look here and here as a start.
– Öskå
Jan 13, 2014 at 13:00
• Somewhat related is a word cloud. Functions used there could be easily adapted to disks. Jan 13, 2014 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.. Jan 13, 2014 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. Jan 13, 2014 at 14:09
• Off-topic here admittedly, but I'm posting this Gist just in case anyone finds this page when looking for a JS solution to the same problem (as I did). This script loads an SVG shape and generates the circles to fit within it. gist.github.com/gouldingken/8d0b7a05b0b0156da3b8 Jul 3, 2015 at 14:09

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". Jan 14, 2014 at 17:43
• Wow, for A in a sea of small circles! Jan 15, 2014 at 11:21
• Now make the background of the 74 a Disk and you win :)
– Öskå
Jan 15, 2014 at 14:07
• @Öskå you got it, Sir! Jan 15, 2014 at 14:50

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. Jan 13, 2014 at 17:27
• @VividD You can easily have that by increasing pad. Jan 13, 2014 at 17:42
• Beautiful! Will check it out! Jan 13, 2014 at 17:48
• Can you also implement a region function so that you can pack the disks into any shape?
– rm -rf
Jan 13, 2014 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, 2014 at 21:16

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. Jan 13, 2014 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! Jan 13, 2014 at 16:31
– rm -rf
Jan 13, 2014 at 20:03
• @rm -rf OK, no problem... :) Jan 13, 2014 at 20:03
• Jacob, check your gmail :)
– Kuba
Jun 20, 2017 at 8:46

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. Jan 14, 2014 at 23:11

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

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. The mask removes all the circles in the black area leaving only the circles in the white.

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 Flatten[imgSimp], imgDim[[2]]], 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]


• What's your meaning of the mask = Flatten[Partition[tab, imgDim[[2]]] img, 1];And your imgDim[[2]]] img is a typo?
– yode
Sep 28, 2015 at 11:41
• If change the mask tomask=Pick[Partition[tab,imgDim[[2]]],Map[Flatten,imgSimp,{-3}],1.]//Level[#,{2}]&;.The code will work in my computer. And I found a strange thing that if you run{Downsample[ImageData[img],3]//Dimensions,Downsample[ImageData[img],4]//Dimensions},you'll get {{77,80,2},{58,60,1}}
– yode
Sep 28, 2015 at 13:09
• @yode -- I'm not sure what changed, but I have now fixed the mask definition so that it removes all the circles in the black area (leaving only the circles in the white). Thanks for alerting me to the problem! Sep 29, 2015 at 17:52

Share a solution just applicable to 10.4 or later version.The code is very terse,but it seem to have some mysterious bug in it.And I have post it as a discuss.Firstly I make a function name of diskMake

diskMake[region_, n_] :=
Module[{p, rad, dist, temRegion = region}, SeedRandom[1];
Reap[Do[p = RandomPoint[temRegion];
rad = If[(dist = Abs[SignedRegionDistance[temRegion, p]]) < .2,
dist, RandomReal[{.2, Min[{dist, .3}]}]];
temRegion =
n]][[-1, -1]]]


Generate some disk bounded by a certain region

wordRegion =
BoundaryDiscretizeGraphics[
Text[Style["21", FontFamily -> "Arial"]], _Text];
diskLarge =
BoundaryDiscretizeRegion@
RegionDifference[BoundingRegion[wordRegion, "FastBall"],
wordRegion];
wordDisk = diskMake[wordRegion, 200];
largerDisk = diskMake[diskLarge, 600];


Draw it with diffrence color

Graphics[MapThread[
Transpose[{RandomColor[
Hue[#1, NormalDistribution[.6, .2],
NormalDistribution[.6, .07]], #2 // Length], #2}] &, {{1/3,
1/2}, {wordDisk, largerDisk}}]]