Circle packing in AutoCAD: creating an Apollonian gasket using F# – Part 1

To carry on from the last post in this series, today’s post is looking at a simple, initial attempt to pack circles into a space using F#.

Rather than starting from the C# code in the previous post, I decided to look for a solution that makes better use of F#’s mathematical capabilities. I came across this simple Common LISP implementation, which creates a subset of a full Apollonian gasket.

[Aside from the links in the previous post, this page may also provide additional insights into programmatic approaches for solving this problem.]

Here’s my equivalent F# code:

module CirclePackingFs

 

open System.Numerics;

 

// Use Descartes theorem to calculate the radius/position

// of the 4th circle

// k4 = k1 + k2 +/- sqrt(k1k2 + k2k3 + k3k1)

 

let solve a b c =

  let D = (a*b + b*c + a*c) in

  [(+);(-)] |> List.map (fun f -> (f (a + b + c) (2.0 * (sqrt D))))

 

let multComplex a (c: Complex) = Complex(a * c.Real, a * c.Imaginary)

 

let solveComplex (a : Complex) (b : Complex) (c : Complex) =

  let D =

    Complex.Multiply(a,b) +

    Complex.Multiply(b,c) +

    Complex.Multiply(a,c) in

  [(+);(-)] |>

    List.map

      (fun f -> (f (a + b + c) (multComplex 2.0 (Complex.Sqrt(D)))))

 

let fourthCircle (x1, y1, k1) (x2, y2, k2) (x3, y3, k3) =

  let pos1 = Complex(x1,y1)

  let pos2 = Complex(x2,y2)

  let pos3 = Complex(x3,y3)

  let r1 = 1. / k1

  let r2 = 1. / k2

  let r3 = 1. / k3

  let kz1 = multComplex k1 pos1

  let kz2 = multComplex k2 pos2

  let kz3 = multComplex k3 pos3

  let k4 = List.head(solve k1 k2 k3)

  let r4 = 1.0 / k4

  let pos4 = multComplex r4 (List.head(solveComplex kz1 kz2 kz3))

  pos4.Real, pos4.Imaginary, k4

 

let rec apollonianGasketInner c1 c2 c3 steps =

  match steps with

  | 0 -> []

  | _ ->

    let c4 = fourthCircle c1 c2 c3

    [(c4, steps)] @

      (apollonianGasketInner c1 c2 c4 (steps - 1) @

        (apollonianGasketInner c2 c3 c4 (steps - 1) @

          (apollonianGasketInner c3 c1 c4 (steps - 1))))

 

type Packer() =

  static member ApollonianGasket outerRad steps =

 

    let rt3 = sqrt(3.)

    let size = outerRad * 2.

    let a = 1. + 2. / rt3

    let innerRad = outerRad / a

    let innerCur = 1. / innerRad

    let h = innerRad * rt3

 

    let c1 = outerRad, size - innerRad, innerCur

    let c2 = outerRad + innerRad, size - (h + innerRad), innerCur

    let c3 = outerRad - innerRad, size - (h + innerRad), innerCur

 

    [(c1,steps);(c2,steps);(c3,steps)] @

    apollonianGasketInner c1 c2 c3 steps

 

I like the succinctness of this: in fact, if it wasn’t for the fact that F# is pretty strict about types (even if they’re inferred at runtime), the code would be even smaller, as we could do away with duplicated functions for handling complex numbers (which are used to perform Descartes’ theorem on our x,y coordinates – a very elegant way of doing things).

The implementation is admittedly a bit different from the original LISP in a couple of areas: rather than creating the geometry in the code, it returns a list of circles for our calling function to process. This allows us to effectively separate our algorithm from AutoCAD. The second difference is in the representation of our circles: rather than returning a point with a radius, we’re returning the point with the curvature of the circle (which is simply one over the radius: curvature = 1 / radius). This isn’t really needed for this post, but it will help us when we come to the next post, which has a full Apollonian gasket implementation.

In order to use our core algorithm’s F# implementation in AutoCAD, I created a simple C# loader that defines a command and calls into the F# library:

using Autodesk.AutoCAD.ApplicationServices;

using Autodesk.AutoCAD.DatabaseServices;

using Autodesk.AutoCAD.EditorInput;

using Autodesk.AutoCAD.Geometry;

using Autodesk.AutoCAD.Runtime;

 

namespace CirclePackingLoader

{

  public class Commands

  {

    [CommandMethod("AGFCI")]

    public static void InnerApollonianGasket()

    {

      Document doc =

        Application.DocumentManager.MdiActiveDocument;

      Database db = doc.Database;

      Editor ed = doc.Editor;

 

      PromptEntityOptions peo =

        new PromptEntityOptions("\nSelect circle: ");

      peo.AllowNone = true;

      peo.SetRejectMessage("\nMust be a circle.");

      peo.AddAllowedClass(typeof(Circle), false);

 

      PromptEntityResult per = ed.GetEntity(peo);

 

      if (per.Status != PromptStatus.OK &&

          per.Status != PromptStatus.None

        )

        return;

 

      Transaction tr =

        doc.TransactionManager.StartTransaction();

      using (tr)

      {

        BlockTableRecord btr =

          (BlockTableRecord)tr.GetObject(

            db.CurrentSpaceId, OpenMode.ForWrite

          );

 

        // If the user selected a circle, use it

 

        Circle cir;

        if (per.Status != PromptStatus.None && per.ObjectId != null)

        {

          cir =

            tr.GetObject(per.ObjectId, OpenMode.ForRead) as Circle;

        }

        else

        {

          // Otherwise create a new one at the default location

 

          cir = new Circle(Point3d.Origin, Vector3d.ZAxis, 500);

          btr.AppendEntity(cir);

          tr.AddNewlyCreatedDBObject(cir, true);

        }

 

        int numSteps = 11;

 

        var res =

          CirclePackingFs.Packer.ApollonianGasket(

            cir.Radius, numSteps

          );

 

        Vector3d offset =

          cir.Center - new Point3d(cir.Radius, cir.Radius, 0.0);

 

        if (res.Length > 0)

        {

          // We get an F# list of tuples containing our circle

          // definitions (no need for a special class)

 

          foreach (var tup in res)

          {

            // Our circles are defined in terms of position (x,y)

            // and curvature (the 3rd item in the nested tuple)

 

            double curvature = System.Math.Abs(tup.Item1.Item3);

            if (1.0 / curvature > 0.0)

            {

              // x and y are Items 1 and 2 (in the nested tuple)

              // respectively

 

              Circle c =

                new Circle(

                  new Point3d(

                    tup.Item1.Item1, tup.Item1.Item2, 0.0

                  ) + offset,

                  Vector3d.ZAxis,

                  1.0 / curvature

                );

 

              // Color index will be based on the "level" of each

              // circle (item 2 in the top-level tuple)

 

              c.ColorIndex = (numSteps - tup.Item2) + 1;

 

              btr.AppendEntity(c);

              tr.AddNewlyCreatedDBObject(c, true);

            }

          }

        }

        tr.Commit();

      }

    }

  }

}

In order to use this, our project needs to reference the F# class library containing our other code as well as having references to FSharp.Core.dll, which lets us process F# lists and tuples without having to fool around with some shared class or structure (something I love about F#).

When we run the AGFCI command, we see a portion of the Apollonian gasket has been created:

Here’s a closer look at the inner geometry:

You’ll see the geometry is no longer in blue: we’re assigning a colour to each circle, to allow us to inspect the levels that are being created (the lowest level chosen has a pink colour, which ultimately gets picked up as the dominant colour of the pattern when zoomed out).

In the next post, we’ll extend this implementation to generate a full Apollonian gasket using F#.