August 2014

Sun Mon Tue Wed Thu Fri Sat
          1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30
31            








« Circle packing in AutoCAD: creating an Apollonian gasket using F# – Part 1 | Main | DevCamps 2012 »

February 13, 2012

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();

      }

    }

  }

}

When we run our new AGFC command, we now see a full Apollonian gasket generated:

A full F#-generated 2D Apollonian gasket

Over the next few posts we’ll make the transition across to 3D, looking at sphere packing.

blog comments powered by Disqus

Feed/Share

10 Random Posts