October 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  










« Sphere packing in AutoCAD: creating an Apollonian packing using F# – Part 1 | Main | Sphere packing in AutoCAD: creating an Apollonian packing using F# – Part 3 »

February 20, 2012

Sphere packing in AutoCAD: creating an Apollonian packing using F# – Part 2

This post continues the series on fill algorithms for 3D printing by looking specifically at an Apollonian sphere packing. In the last post we got most of the way there, but today we’re going to introduce a more elegant algorithm for solving the problem (with pretty impressive results :-).

Many thanks to Professor Ronald Peikert from ETH Zurich for kindly providing the C++ code used to generate the F# code in today’s post.

The original algorithm was outlined in this paper co-authored by Professor Peikert, under the section The “Inversion Algorithm” and – as you might divine from the name – uses inversion to simplify the challenge of identifying mutually tangent spheres.

Professor Peikert has also provided notes describing the steps to derive the special coordinates he has used to simplify working with large sets of spheres, taking us through the transformation between Cartesian, inversive and then "special” coordinates.

Here is the F# code:

module SpherePackingInversionFs

 

let mutable minRad = 0.01

let mutable maxGen = 6

let mutable clans = false  // assign color to clans, not generations

 

let draw generation clan A B C D E =

 

  // returns false if the radius is less than minRad

 

  // Convert special coords to inversive coords

 

  let rt2o2 = sqrt(2.)/2.

  let a = (-B + C + D - E) * rt2o2

  let b = ( B - C + D - E) * rt2o2

  let c = ( B + C - D - E) * rt2o2

  let d = A - B - C - D - E

  let e = ( B + C + D + E) * sqrt(6.)/2.

 

  // Convert inversive coords to Cartesian coords

 

  // Positive radius required for drawing

 

  let r = if e-d > 0. then 1./(e-d) else 1./(d-e)

  let x = a*r

  let y = b*r

  let z = c*r

 

  if r < minRad then

    []

  else

    let inversionCircle = generation < 0

    let outerCircle = (generation = 0 && clan = 0)

    let transparent = inversionCircle || outerCircle

 

    let colorIndex = if clans then clan else (generation + 1)

 

    if transparent then

      []

    else

      // Check for z < 0 to draw lower half of the packing only

      [((x, y, z, r), colorIndex)]

 

let generateItem item itemVal i A B C D E =

  match i with

  | 0 -> if item = i then -A else A + itemVal

  | 1 -> if item = i then -B else B + itemVal

  | 2 -> if item = i then -C else C + itemVal

  | 3 -> if item = i then -D else D + itemVal

  | _ -> if item = i then -E else E + itemVal

 

let rec generate generation clan i A B C D E =

  let A1 = generateItem 0 A i A B C D E

  let B1 = generateItem 1 B i A B C D E

  let C1 = generateItem 2 C i A B C D E

  let D1 = generateItem 3 D i A B C D E

  let E1 = generateItem 4 E i A B C D E

 

  // Avoid duplicates

 

  if (i=0 && (        B1< A1||C1< A1||D1< A1||E1< A1) ||

      i=1 && (A1<=B1||        C1< B1||D1< B1||E1< B1) ||

      i=2 && (A1<=C1||B1<=C1||        D1< C1||E1< C1) ||

      i=3 && (A1<=D1||B1<=D1||C1<=D1||        E1< D1) ||

      i=4 && (A1<=E1||B1<=E1||C1<=E1||D1<=E1       )) then

    []

  else

    let res = draw generation clan A1 B1 C1 D1 E1

    if res = [] then

      []

    else

      res @

      if generation < maxGen then

        List.concat

          (List.map

            (fun j ->

              if j <> i then

                generate (generation+1) clan j A1 B1 C1 D1 E1

              else

                [])

            [0..4])

      else

        []

 

type Packer() =

  static member ApollonianGasket steps minR cl =

 

    maxGen <- steps

    minRad <- minR

    clans <- cl

 

    // Draw initial circles and generate all circles of their clans

 

    let init =

      draw 0 0 1. 0. 0. 0. 0. @ generate 1 0 0 1. 0. 0. 0. 0. @

      draw 0 1 0. 1. 0. 0. 0. @ generate 1 1 1 0. 1. 0. 0. 0. @

      draw 0 2 0. 0. 1. 0. 0. @ generate 1 2 2 0. 0. 1. 0. 0. @

      draw 0 3 0. 0. 0. 1. 0. @ generate 1 3 3 0. 0. 0. 1. 0. @

      draw 0 4 0. 0. 0. 0. 1. @ generate 1 4 4 0. 0. 0. 0. 1.

 

    // Draw inversion spheres

 

    let Q = sqrt(3.)/6.

 

    init @

    draw -1 -1 (-2.*Q) Q Q Q Q @

    draw -1 -1 Q (-2.*Q) Q Q Q @

    draw -1 -1 Q Q (-2.*Q) Q Q @

    draw -1 -1 Q Q Q (-2.*Q) Q @

    draw -1 -1 Q Q Q Q (-2.*Q)

 

The C# loader code shown last time has been updated to add a new command (AGFS) calling this new F# implementation.

using Autodesk.AutoCAD.ApplicationServices;

using Autodesk.AutoCAD.DatabaseServices;

using Autodesk.AutoCAD.Geometry;

using Autodesk.AutoCAD.Runtime;

 

namespace SpherePackingLoader

{

  public class Commands

  {

    [CommandMethod("AGFS")]

    public static void ApollonianGasket()

    {

      ApollonianGasket(true);

    }

 

    [CommandMethod("AGFSE")]

    public static void ExteriorApollonianGasket()

    {

      ApollonianGasket(false, true);

    }

 

    [CommandMethod("AGFSI")]

    public static void InteriorApollonianGasket()

    {

      ApollonianGasket(false, false);

    }

 

    public static void CreateLayers(Database db, Transaction tr)

    {

      LayerTable lt =

        (LayerTable)tr.GetObject(db.LayerTableId, OpenMode.ForWrite);

      for (short i = 1; i <= 20; i++)

      {

        string name = i.ToString();

        if (!lt.Has(name))

        {

          LayerTableRecord ltr = new LayerTableRecord();

          ltr.Color =

            Autodesk.AutoCAD.Colors.Color.FromColorIndex(

              Autodesk.AutoCAD.Colors.ColorMethod.ByAci, i

            );

          ltr.Name = name;

          lt.Add(ltr);

          tr.AddNewlyCreatedDBObject(ltr, true);

        }

      }

    }

 

    public static void ApollonianGasket(

      bool inverse, bool exterior = false

    )

    {

      Document doc =

        Application.DocumentManager.MdiActiveDocument;

      Database db = doc.Database;

 

      Transaction tr =

        doc.TransactionManager.StartTransaction();

      using (tr)

      {

        CreateLayers(db, tr);

 

        BlockTable bt =

          (BlockTable)tr.GetObject(

            db.BlockTableId, OpenMode.ForRead

          );

        BlockTableRecord btr =

          (BlockTableRecord)tr.GetObject(

            bt[BlockTableRecord.ModelSpace], OpenMode.ForWrite

          );

 

        var res =

          inverse ?

            SpherePackingInversionFs.Packer.ApollonianGasket(

              7, 0.001, false

            )

            :

              exterior ?

                SpherePackingFs.Packer.ExteriorApollonianGasket(5)

                :

                SpherePackingFs.Packer.InteriorApollonianGasket(5);

 

        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 radius = System.Math.Abs(tup.Item1.Item4);

            if (radius > 0.0)

            {

              // x, y & z are items 1, 2 & 3 (in the nested tuple)

              // respectively

 

              Solid3d s = new Solid3d();

              s.CreateSphere(radius);

              Point3d center =

                new Point3d(

                  tup.Item1.Item1, tup.Item1.Item2, tup.Item1.Item3

                );

              Vector3d disp = center - Point3d.Origin;

              s.TransformBy(Matrix3d.Displacement(disp));

 

              // The Layer (and therefore the colour) will be based

              // on the "level" of each sphere

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

 

              s.Layer = tup.Item2.ToString();

 

              btr.AppendEntity(s);

              tr.AddNewlyCreatedDBObject(s, true);

            }

          }

        }

        tr.Commit();

      }

    }

  }

}

When we run the AGFS command, we see a beautiful Apollonian packing gets generated in 3D:

Our full 3D Apollonian packing

And if we turn off layer “1”, we effectively hide the four primary internal spheres, allowing a deeper view into the guts of our packing:

And with the layer with the four primary internal spheres turned off

If you’re interested in looking at the results without running the code, here they are as an AutoCAD 2010 drawing without the enclosing sphere (6.16 MB).

In the next post, we’ll take the results of this packing algorithm and subtract them from an enclosing sphere, to see what kind of results we get.

There’s one last implementation that I’d like to present before publishing the complete source project (this will hopefully come at the end of the week). The project will also contain a version of today’s post that runs in parallel (as it’s a highly parallelizable algorithm :-). While it’s probably a good topic for a future post, I left the code in there for anyone who’s interested, in the meantime.

blog comments powered by Disqus

Feed/Share

10 Random Posts