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

Following on from the previous post in this series, today’s post completes the implementation to create a full Apollonian gasket in AutoCAD using F#.

As a comment on the original Common LISP implementation, someone had contributed a more complete version which allowed me to complete today’s F# version.

Here’s the additional F# file for the project (which I’ll be providing in full at the end of the series):

module CirclePackingFullFs
 
open System.Numerics;
 
// Use Descartes' theorem to calculate the radius/position
// of the 4th circle
// k4 = k1 + k2 + k3 +/- sqrt(k1k2 + k2k3 + k3k1)
 
let solve a b c =
  [(+);(-)] |>
    List.map
      (fun f -> (f (a + b + c) (2. * (sqrt (a*b + b*c + a*c)))))
 
let multComplex a (c : Complex) =
  Complex(a * c.Real, a * c.Imaginary)
 
let solveComplex (a : Complex) (b : Complex) (c : Complex) =
  [(+);(-)] |>
    List.map
      (fun f ->
        (f (a + b + c)
           (multComplex 2.
             (Complex.Sqrt
               (Complex.Multiply(a,b) +
                Complex.Multiply(b,c) +
                Complex.Multiply(a,c))))))
 
// Return a list of all solution circles (positive and negative)
 
let apollonianCircle (x1,y1,k1) (x2,y2,k2) (x3,y3,k3) =
  let p1 = Complex(x1,y1)
  let p2 = Complex(x2,y2)
  let p3 = Complex(x3,y3)
  let kz1 = multComplex k1 p1
  let kz2 = multComplex k2 p2
  let kz3 = multComplex k3 p3
  let curSols = solve k1 k2 k3
  let posSols = solveComplex kz1 kz2 kz3
  List.map2
    (fun cur pos ->
      let r4 = 1. / cur
      let pos4 = multComplex r4 pos
      pos4.Real, pos4.Imaginary, cur)
    curSols posSols
 
// Only return certain circles - for use with "hard" gaps
 
let apollonianCircle2 (x1,y1,k1) (x2,y2,k2) (x3,y3,k3) =
  let p1 = Complex(x1,y1)
  let p2 = Complex(x2,y2)
  let p3 = Complex(x3,y3)
  let kz1 = multComplex k1 p1
  let kz2 = multComplex k2 p2
  let kz3 = multComplex k3 p3
  let cur = List.head (solve k1 k2 k3)
  let pos = List.nth (solveComplex kz1 kz2 kz3) 1
  let r4 = 1. / cur
  let pos4 = multComplex r4 pos
  pos4.Real, pos4.Imaginary, cur
 
let rec easyGap c1 c2 c3 steps =
  match steps with
  | 0 -> []
  | _ ->
    let c4 = List.head(apollonianCircle c1 c2 c3)        
    [(c4,steps-1)] @
    easyGap c1 c2 c4 (steps-1) @
    easyGap c2 c3 c4 (steps-1) @
    easyGap c3 c1 c4 (steps-1)
 
let hardGap c1 c2 c3 steps =
  match steps with
  | 0 -> []
  | _ ->
    let g1 = apollonianCircle2 c1 c2 c3       
    [(g1,steps-1)] @
    easyGap c2 c3 g1 (steps-1) @
    easyGap c2 c1 g1 (steps-1) @
    if steps = 1 then
      []
    else
      let g2 = apollonianCircle2 c1 c3 g1
      [(g2,steps-2)] @
      easyGap c3 g1 g2 (steps-2) @
      easyGap c1 g1 g2 (steps-2) @
      if steps = 2 then
        []
      else
        let g3 = apollonianCircle2 c1 c3 g2
        [(g3,steps-3)] @
        easyGap c1 g2 g3 (steps-3) @
        easyGap c1 c3 g3 (steps-3) @
        if steps = 3 then
          []
        else
          let g4 = apollonianCircle2 c3 g2 g3
          [(g4,steps-4)] @
          easyGap c3 g2 g4 (steps-4) @
          easyGap g2 g3 g4 (steps-4)
 
 
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 c0 = outerRad, outerRad, -1. / outerRad
    let c1 = outerRad, size - innerRad, innerCur
    let c2 = outerRad + innerRad, size - (h + innerRad), innerCur
    let c3 = outerRad - innerRad, size - (h + innerRad), innerCur
 
    let c4 = List.head(apollonianCircle c1 c2 c3)
 
    let nextStep = steps-1
 
    // Main circles
    // (Assumption is that the outer circle already exists)
 
    [(c1,steps);(c2,steps);(c3,steps);(c4,nextStep)] @
 
    // Inner gaps
 
    easyGap c1 c2 c4 nextStep @
    easyGap c2 c3 c4 nextStep @
    easyGap c3 c1 c4 nextStep @
 
    // Peripheral gaps
 
    easyGap c1 c2 c0 nextStep @
    easyGap c2 c3 c0 nextStep @
    hardGap c3 c1 c0 nextStep
 

The hardGap function does a fair amount of manual filling of the top-left gap… I suspect this could also have been addressed by mirroring or rotating the results from one of the easyGap calls, but there you go. Ultimately the code performs well enough.

The corresponding C# loader file clearly also needs updating to define an additional command calling the F# code:

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("AGFC")]
    public static void ApollonianGasket()
    {
      PackApollonianCircles(true);
    }
 
    [CommandMethod("AGFCI")]
    public static void InnerApollonianGasket()
    {
      PackApollonianCircles(false);
    }
 
    public static void PackApollonianCircles(bool full)
    {
      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 =
            full ?
            CirclePackingFullFs.Packer.ApollonianGasket(
              cir.Radius, numSteps
            )
            :
            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();
      }
    }
  }
}