Circle packing theorem states:

For every connected simple planar graph G there is a circle packing in the plane whose intersection graph is (isomorphic to) G.

How to use Mathematica to draw examples of such pairs (planar graph, circle packing)?

More specifically, the requirements of my question could be formulated as follows:

Generate a random connected simple planar graph with given number of nodes $N$ (this can be done lets say using adjacency matrix; in a way another parameter affecting resulting graph could be "density" of links). Draw its usual graphical representation containing nodes and links. Then draw its circle packing representation.

In other words, a solution could contain a slider for $N$, and a slider for the "density" (those two variables are parameters). Based on values of these sliders, get regular drawing of a random graph with given $N$ and density, and also drawing of its circle pack.


To my knowledge, this is not a trivial problem, and there is no built-in Mathematica solution. Though, there are several questions and answers here. One of best is Create triangular mesh from random list of points :

pts = RandomReal[{0, 12}, {100, 2}];
dt = DelaunayTriangulation[pts];
dt // Column
toPairs[{m_, ns_List}] := Map[{m, #} &, ns];
edges = Flatten[Map[toPairs, dt], 1];
Graphics[GraphicsComplex[pts, {Line[edges], 
  Red, PointSize[Large], Point[pts]}]]


See: Collins, Stephenson: A circle packing algorithm.

Visual examples:

enter image description here

enter image description here

enter image description here

This question is not an assignment. It is an intellectual experiment, related to aestetics and structure. I am not requesting code. I am seeking ideas, approaches, and solutions.

  • 2
    $\begingroup$ This question is too broad. Please be more specific. $\endgroup$
    – Yves Klett
    May 12 '14 at 7:09
  • 10
    $\begingroup$ And please, this is not a feature/code-request site. It would be much more motivating for others if you showed what you already tried. Plus, if it is an assignment, add the appropriate tags. $\endgroup$
    – Yves Klett
    May 12 '14 at 7:17
  • 8
    $\begingroup$ It was not intended to be accusatory, but personally I feel it is polite to show previous efforts. Just to clarify: is it an assignment or not? Just to get the tags right. If it is not, please accept my apologies, but verbatim pasted paragraphs often carry that notion. $\endgroup$
    – Yves Klett
    May 12 '14 at 7:34
  • 3
    $\begingroup$ @VividD I'd like to see you answer some questions so that I know you are not here only to take. Though I speak only for myself I imagine others feel the same way. (Nevertheless it is an interesting question, with or without effort shown, and it has my vote.) $\endgroup$
    – Mr.Wizard
    May 12 '14 at 16:03
  • 4
    $\begingroup$ @rm-rf: Actually, looking at VividD's contributions to the site so far, I see zero evidence that they are capable of writing Mathematica code or have any interest in learning to do so. $\endgroup$
    – user484
    May 13 '14 at 9:41

I bought 'Introduction to Circle Packing' by Kenneth Stephenson a few years ago, and have been trying to understand it since. I found theFindMaximumanswer by @SimonWoods very good, but larger graphs slow that algorithm down, as he said.

@VividD said "I am not requesting code" and "I am seeking ideas, approaches, and solutions". I finally have the code, but there is a lot of it. Therefore, here is a description of one approach based on my interpretation of Stephenson.

The key is to break the problem into two halves.

First: Find radii without placing any circles.

Find the radii of all circles such that the combinatorics of the graph are satisfied. That is, find radii such that the distance between adjacent nodes i and j is radii[[i]]+radii[[j]], for all i and j, i $\ne$ j.

Stephenson (page 244) describes the uniform neighbour algorithm to iteratively find these radii. For any node $V$ in the graph, adjust its current radius by multiplying it by $b (1-d) /(d (1-b))$, where $b=\sin(t/(2k))$, $d=\sin(\pi/k)$, $t$ is the angle sum, and $k$ is the number of neighbouring nodes surrounding $V$. The angle sum (pages 15, 57) is the total angle subtended at the central node $V$ by adjacent pairs of surrounding circles labelled with their current, as yet approximate, radii. An angle sum greater than $2\pi$ means the central circle radius is too small, causing the surrounding circles to overlap. An angle sum less than $2\pi$ means the central circle radius is too large, producing gaps between the surrounding circles. An angle sum of $2\pi$ is just right, the $k$ surrounding circles are pair-wise tangent.

The uniform neighbour algorithm of Stephenson cycles through each interior node, adjusting the radius corresponding to that node. Multiple cycles through the interior nodes enables convergence to a stable set of radii.

Second: Find the circle centres corresponding to their radii.

Once the radii are found, place any two adjacent circles, a and b, such that their centres are separated by the sum of their radii, that is, they are tangent. Then place a (positively oriented) triangular face by finding the position of the centre of the third circle c which makes a mutually tangent triple with these two placed circles. This process generates two new pairs of placed circles, (c,b) and (a,c), leading to placements of new third circles, triangular face by face.

Following this outline, plus some compiling and memoization, I was able to construct circle packings of 1500 circles in about 125 s. Finding the radii is by far the most time-consuming part.

Rendering is the easy part:


  • $\begingroup$ I wonder if this is the same as the algorithm described in the Wikipedia article which is the first link in the question? $\endgroup$
    – user484
    Jun 18 '14 at 23:32
  • $\begingroup$ It is indeed. Have I wasted bandwidth? I hope there is some value added with the explicit expression for the modified value of the central radius r... $\endgroup$ Jun 19 '14 at 0:00
  • $\begingroup$ No no, I was just curious if it was the same or whether there was some subtlety that I missed. $\endgroup$
    – user484
    Jun 19 '14 at 0:03
  • $\begingroup$ This is great! This answer required a lot of effort and expertise! There is a big difference between reading algorithm description (in Wikipedia or elsewhere) and knowing that it actually works, how it works, where its strong and weak points are, etc! I wonder if you could post some more results, perhaps for graphs mentioned here. Maybe you could export results and place them in a dropbox shared folder? $\endgroup$
    – VividD
    Jun 20 '14 at 8:24
  • $\begingroup$ @KennyColnago, I have some suggestions: It would be nice if you take all planar graphs from Mathematica curated data (using GraphData[]; there are several thousand such graphs, I believe), and export both original graph drawing, and the circle pack obtained by your code. $\endgroup$
    – VividD
    Jun 21 '14 at 5:26

Two versions of a solution to this question will be presented. Both are applications of built-in Mathematica optimization tools. The key point of both version is the choice of optimization goal and constraints (or, in another words, model formulation).

Interested reader may find article Configuration Analysis and Design by Using Optimization Tools in Mathematica by FRANK J. KAMPAS and JÁNOS D. PINTÉR useful. It is not about the problem from this question, but deals with solving geometric problems using Mathematica's optimization tools. There are also several related questions on thus site.

First version

This will surely be hopeless for large graphs, but here's a crude approach using NMinimize. The aim is to minimize the difference between the inter-vertex distance and the sum of their radii, and also the offset distance between the vertex positions in the original and final graph.

Let's say heptahedral-5 graph is chosen as the original graph:

g = GraphData[{"Heptahedral", 5}]

enter image description here

This code will produce its circle packing, using NMinimize in appropriate setting:

n = VertexCount[g];    
pts = Table[PropertyValue[{g, i}, VertexCoordinates], {i, n}];    
rads = Array[r, n];
xs = Array[x, n];
ys = Array[y, n];

f1 = Total[(r[#1] + r[#2] - EuclideanDistance[{x[#1], y[#1]}, {x[#2], y[#2]}])^2 & @@@ EdgeList[g]];
f2 = Total[MapThread[({#1, #2} - #3)^2 &, {xs, ys, pts}], -1];

rules = Last@NMinimize[f1 + 0.1 f2, Join[rads, xs, ys]];

circles = MapThread[Circle[{#1, #2}, #3] &, {xs, ys, rads} /. rules];
g2 = Fold[SetProperty[{##}, VertexCoordinates -> {x[#2], y[#2]} /. rules] &, g, Range[n]];

Show[Graphics[circles], g2]

enter image description here

Improved version

This is an improved version which seems to cope with larger graphs. This is still based on black box optimization but with some changes:

  • Use FindMaximum instead of NMaximize to find a local maximum.

  • Circles which should be touching (i.e. graph edges) are put in as constraints, i.e. sum of radii == distance between centres.

  • All other pairs of vertices (non-edges) are placed as far apart as possible - the sum of these distances is what we maximize.

Also, in this example, I've used a Delaunay triangulation (as suggested in the question) to create a simple random planar graph (however, using an undocumented DelaunayMesh function.

Other steps are also marked and described by comments in the code:

dist[{a_, b_}] := Sqrt[(a - b).(a - b)];
n = 30;

(* get a random mesh *)
pts = RandomReal[1, {n, 2}];
m = DelaunayMesh[pts]["MeshObject"];

(* get edges and non-edges *)
edges = Sort /@ m["Edges"];
nonedges = Complement[Subsets[Range[n], {2}], edges];

(* variable lists *)
rads = Array[r, n];
pos = Array[{x[#], y[#]} &, n];

(* set up constraints for edges *)
cons = (Total@rads[[#]] == dist@pos[[#]]) & /@ edges;
cons = Join[cons, Thread[rads > 0]];

(* define function to maximize *)
f = Total[dist@pos[[#]] & /@ nonedges];

(* variable initializations *)
vars = Join[rads, Flatten[{pos, pts}, {2, 3}]];

(* do the maximization *)
{posv, radsv} = {pos, rads} /. Last@FindMaximum[{f, cons}, vars];

(* result *)
Graphics[{MapThread[Circle, {posv, radsv}],
  Opacity[0.5, Green], Line[posv[[#]] & /@ edges]}]

enter image description here

  • $\begingroup$ Great application of Mathematica built-in optimisations! Thanks a lot! I am also curious why do you think that in FindMaximum vs NMaximize, FindMaximum wins for larger graphs? (Or maybe this is not the case?) $\endgroup$
    – VividD
    May 14 '14 at 9:40
  • $\begingroup$ for your first method, I tried g = GraphData[{"BananaTree", {7, 4}}] , with all other code the same, the resulting circle pac is not entirely correct for some reason... $\endgroup$
    – VividD
    May 14 '14 at 10:21
  • $\begingroup$ The original algorithm is not very good - it doesn't constrain the circles on connected vertices to touch, it just attempts to minimize the gaps. At the same time it is trying to minimize the movement of the vertices from their initial positions (which is not really what we want to do). So it reaches a point where it will sacrifice the circle-packing in order to avoid rearranging the vertex positions. $\endgroup$ May 14 '14 at 10:57
  • $\begingroup$ The second version improves things by using constraints to insist that the circles on connected vertices touch each other. The original vertex positions go in as a starting point to FindMaximum, which seeks a local maximum (whereas NMaxmimize looks for a global optimum). It can still give unsatisfactory results, with overlapping circles for non-edges. I think this means that FindMaximum has found a local maximum which is not a global optimum. I guess this could be remedied with more constraints. $\endgroup$ May 14 '14 at 10:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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