Circle packing – theory and practice

I am spending a few months in Göttingen as a Courant Distinguished Visiting Professor, and talking a bit to Laurent Bartholdi about rational functions — i.e. holomorphic maps from the Riemann sphere \widehat{\mathbb C} to itself. A rational function is determined (up to multiplication by a constant) by its zeroes and poles, and can therefore generically be put in the form f:z \to P(z)/Q(z) where P and Q are polynomials of degree d. If d=1 then f is invertible, and is called a fractional linear transformation (or, sometimes, a Mobius transformation). The critical points are the zeroes of P'Q-Q'P; note that this is a polynomial of degree \le 2d-2 (not 2d-1) and the images of these points under f are the critical values. Again, generically, there will be 2d-2 critical values; let’s call them V. Precomposing f with a fractional linear transformation will not change the set of critical values.

The map f cannot usually be recovered from V (even up to precomposition with a fractional linear transformation); one needs to specify some extra global topological information. If we let \overline{C} denote the preimage of V under f, and let C denote the subset consisting of critical points, then the restriction f:\widehat{\mathbb C} - \overline{C} \to \widehat{\mathbb C} - V is a covering map of degree d, and to specify the rational map we must specify both V and the topological data of this covering. Let’s assume for convenience that 0 is not a critical value. To specify the rational map is to give both V and a representation \rho:\pi_1(\widehat{\mathbb C} - V,0) \to S_d (here S_d denotes the group of permutations of the set \lbrace 1,2,\cdots,d\rbrace) which describes how the branches of f^{-1} are permuted by monodromy about V. Such a representation is not arbitrary, of course; first of all it must be irreducible (i.e. not conjugate into S_e \times S_{d-e} for any 1\le e \le d-1) so that the cover is connected. Second of all, the cover must be topologically a sphere. Let’s call the (branched) cover \Sigma for the moment, before we know what it is. The Riemann-Hurwitz formula lets one compute the Euler characteristic of \Sigma from the representation \rho. A nice presentation for \pi_1(\widehat{\mathbb C}-V,0) has generators e_i represented by small loops around the points v_i \in V, and the relation \prod_{i=1}^{|V|} e_i = 1. For each e_i define o_i to be the number of orbits of \rho(e_i) on the set \lbrace 1,2,\cdots,d\rbrace. Then

\chi(\Sigma) = d\chi(S^2) - \sum_i (d-o_i)

If each \rho(e_i) is a transposition (i.e. in the generic case), then o_i=d-1 and we recover the fact that |V|=2d-2.

This raises the following natural question:

Basic Question: Given a set of points V in the Riemann sphere, and an irreducible representation \rho:\pi_1(\widehat{\mathbb C} - V,0) \to S_d satisfying \sum_i (d-o_i) = 2d-2, what are the coefficients of the rational function z \to P(z)/Q(z) that they determine (up to precomposition by a fractional linear transformation)?

Note that we would like to recover the coefficients numerically (i.e. as numbers with decimal points). And we are really interested in finding a practical algorithm, and then implementing it on computer. One obvious (and bad) idea is to just solve for the coefficients of P and Q subject to the constraint that the set of critical values is V (after normalizing so that three of the critical points are 0, 1 and infinity to remove the ambiguity of the precomposition). The problem is that the number of such solutions is exponential in the degree, and although Newton’s method will quickly find some solution, it is very, very unlikely to be the solution with the correct combinatorial data.

Another idea — and one that leads to the point of this blog post — is to try to build the branched cover \Sigma directly as a Riemann surface together with a holomorphic map with the correct critical values and combinatorics, and then uniformize it as \widehat{\mathbb C} to determine the numerical location of the zeros and poles. This sounds more promising, since there is an obvious way to build \Sigma piecewise from copies of regions in \widehat{\mathbb C}- V glued together by very explicit maps. The problem is that (numerical) uniformization itself is very slow, at least if one wants any kind of accuracy. On the other hand, we do not need to know the values of the uniformization map everywhere, only the locations of the zeros and poles. So we can try to ask for a fast and approximate method of uniformization which gives sufficiently accurate values of these numbers, that they can then be adjusted quickly to very accurate values by Newton’s method.

One potential idea is to use the method of circle packing. A circle packing is a (rigid) configuration of round circles with disjoint interiors and prescribed combinatorial pattern of tangencies. Abstractly, the circle packing determines a triangulation, with one vertex for each circle, one edge for each tangency, and one triangle for each triple of mutually tangent circles. Implicitly, by using the term “round circle”, the domain in which the circles are packed should be a Riemann surface together with a complex projective structure; for example, the Riemann sphere, or the Euclidean or hyperbolic planes. Given a projective surface and a triangulation satisfying mild topological conditions, such a circle packing exists and is unique; this is known as the Circle Packing Theorem (aka the Koebe-Andreev-Thurston theorem). One can also solve for configurations of circles intersecting at prescribed angles; for instance, one can look for a configuration of round circles with prescribed combinatorics, and meeting always at right angles. Such a configuration in the Riemann sphere can be interpreted as the boundaries of a collection of geodesic planes in hyperbolic 3-space meeting in right angles, and cutting out a compact right-angled polyhedron. The existence of such a polyhedron is the base step in Thurston’s inductive proof of geometrization for Haken 3-manifolds.

One can also think of the circle packing as a discrete version of a conformal structure; at a talk at the conference in 1985 celebrating de Branges’ proof of the Bieberbach conjecture, Thurston proposed using circle packings to approximate conformal mappings. One starts with a region U in the complex plane, and packs it nearly tightly with a hexagonal packing of small round circles. Together with the boundary of U one obtains a topological circle packing; the round circle packing with the same combinatorics can be thought of as a packing of the unit disk. One therefore obtains a coarse “map” from U to the disk, taking each round circle in the packing to a round circle in the disk, and one should think of this as a discrete version of a conformal map. As the mesh size goes to zero, Thurston conjectured these maps should converge to the uniformizing map. This conjecture was proved by Rodin-Sullivan.

In the context of our Basic Question we can try to find our desired rational map as follows. Starting with the collection of points V in the Riemann sphere, we build an (almost) round circle packing in such a way that one of the centers is at 0 and infinity and at each point of V. One should probably choose the mesh size quite small near these special points, since the derivative of f^{-1} is going to blow up near V. This determines a topological graph (the 1-skeleton of the triangulation described above). The branching data defines a new graph in an obvious way, and we can find the circle packing associated to that graph. The “preimages” of 0 and infinity are given as the centers of d circles in this new circle packing, and we can take these as our approximate zeros and poles for f.

Anyway, this is all preamble to explaining that I wrote a little code to perform circle packing with prescribed combinatorial data, and in case I don’t do anything else with it (which is quite likely) I thought it might be amusing to post the code and some of the pictures it produces. Note that very sophisticated and highly optimized code for circle packing is already available from many other places; for example, Ken Stephenson has an amazing collection of resources (both theoretical and computational) on his website.

The latest code is, and will continue to be, available at github here:

(link to github repository)

Download the source as a zip file, then unzip and type “make” to make. The program is written in C++, and outputs graphics to the screen using Xlib, and to an .eps file. It can either be run interactively (without an argument) or non-interactively, taking a file name (containing a topological circle packing) as a parameter.

A topological circle packing is written in a (ASCII text) file with the following structure:

number of vertices

valence of vertex 0; list of adjacent vertices in circular order

valence of vertex 1; list of adjacent vertices in circular order (etc.)

valence of vertex n-1; list of adjacent vertices in circular order

initial radius of vertex 0

initial radius of vertex 1 (etc)

initial radius of vertex n-1.

By convention, vertex 0 is the “outer” circle (centered at infinity). The program doesn’t check that the adjacency data that it is given is consistent, or that it gives rise to a topological sphere.

If you try this out, please let me know what you think could be improved or fixed. Feel free to modify or change the code however you like (subject to the usual GPL license conditions). A wishlist would include to add the functionality I sketched above, i.e. to find approximate rational maps with prescribed critical values and branching data; any reader with some time on their hands is warmly invited to do this!


Some screen shots of the X-windows output:

and some sample .eps output (converted to .jpg for wordpress):

The effect of taking a double branched cover over two circles:

This entry was posted in Complex analysis, Visualization and tagged , . Bookmark the permalink.

6 Responses to Circle packing – theory and practice

  1. Chris Jerdonek says:

    You should put the code on GitHub! It’s pretty easy to use. It would be easier for people to see and communicate about any changes that way (and to collaborate on making them, too), and it would give the code more visibility.

  2. Jonathon says:

    I really like what you guys are usually up too.
    This sort of clever work and exposure! Keep up the superb works guys I’ve incorporated you guys to my blogroll.

  3. Torsten says:

    Hi this is nice work,
    I tried your algorithm for some simple examples, but it gave no out output, meaning it didn’t stop. Does it work only only for hexgonal circle configurations?

    For example the following vertices-edge file produces no result:
    4 1 2 3 4
    7 0 3 4 5 6 7 8
    7 0 3 4 17 18 19 20
    7 0 1 2 5 9 13 17
    7 0 1 2 8 12 16 20
    4 1 3 6 9
    4 1 5 7 10
    4 1 6 8 11
    4 1 4 7 12
    4 3 5 10 13
    4 6 9 11 14
    4 7 10 12 15
    4 4 8 11 16
    4 3 9 14 17
    4 10 13 15 18
    4 11 14 16 19
    4 4 12 15 20
    4 2 3 13 18
    4 2 14 17 19
    4 2 15 18 20
    4 2 4 16 19

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s