Through the Interface

AU 2012 class recordings now available online

This week has so far had an AU theme to it, just as last week we talked exclusively about Leap Motion and AutoCAD. Perhaps I’m sub-consciously shifting this blog to a weekly-themed format? Hmm.

Like many of you, I’m sure, I received an email over the weekend to let me know that the recorded sessions from Autodesk University 2012 are now available online (for anyone with a valid AU online account).

I wasn’t sure which of my sessions had made it up there from this last year’s event (thank goodness I rarely have to write cheques anymore), and so was reasonably happy to find out that “Moving Code to the Cloud: It's Easier Than You Think” and “Secure Your AutoCAD® Environment” (unsurprisingly, really, as it’s a virtual session) have both been posted. My best-rated session from last year’s AU (“Developing a Simple Metro-Style Application for Microsoft® Windows® 8”) wasn’t recorded, but that’s how it goes. Maybe that’s why it received better feedback. :-)

There is loads of great content up there on the AU Online site. I strongly recommend checking out the content from the other members of my team who presented at AU 2012, for instance, even if these particular sessions were unfortunately also not recorded:

• Albert Szilvasy’s “The Confluence of AutoCAD® Mobile, Desktop, and the Cloud”
  • Loads of great insights on the future contained in these slides
• Christer Janson’s “GrumpyBlocks: A Tale of Writing a Game in AutoCAD® Using ObjectARX®”
  • Check out the full ObjectARX source code provided (as mentioned back on April 1st, although this time it’s serious :-). Woah!
• Ravi Krishnaswamy’s “PaaSt the Desktop: Implementing Cloud-Based Productivity Solutions with the AutoCAD® ObjectARX® API”
  • Well-documented approach for creating AutoCAD apps that work with Google App Engine

Please feel free to post comments with links to other sessions that people might find of interest (posting links to your own is fine, too, as long as they’re at least vaguely development-focused). I know the ADN team presented many well-received and valuable sessions at AU 2012, for instance.

photo credit: SimonDeanMedia via photopin cc

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

As suggested in the last post, today we’re going to take the results of running the code from that post and use them to generate a hollowed-out sphere. A big thanks to Francesco Tonioni, from our Product Support team in Neuchatel, who spent some time throwing ideas around on a lazy (but very cold) Sunday afternoon, contributing significantly to this post.

A few minor changes to the code were needed: rather than creating the spheres at exactly the size at which they were generated by the F# code, I adjusted the C# code to multiply the radius by 0.98 and use that instead (simply to avoid spheres touching and therefore making it harder for the Autodesk Shape Manager component to effectively perform the Boolean subtract). I also adjusted the code to only make the inversion spheres transparent – not the outer sphere (which also got reduced in radius by 2%). I then manually scaled the outer sphere by 1.05 (which is probably more than was strictly needed, but anyway).

We used AutoCAD’s SUBTRACT command to subtract the results of the packing from the outer sphere. This took some time to complete, but we ended up with a partially hollowed sphere that could then be sectioned to expose its insides:

I don’t have access to a 3D printer (and frankly doubt that current technology would allow this to be printed, given the need for support material when building arches), but I did look a little into the mass properties of this – and the original sphere – using the MASSPROP command. The original sphere has a mass/volume of 4.5639, while our hollowed version has a mass/volume of 1.1752 (i.e. it’s about a quarter of the mass).

Not very surprising – and I’m sure this could be tweaked downwards with a deeper packing and perhaps a smaller percentage change in our spheres’ radii – but nonetheless interesting.

To take this further, I’ve been working on a different packing strategy that provides better overall strength, as well as considering non-spherical volumes. When initially looking into these directions, I was very grateful for pointers from Alex Fielder (who sent me a link to this interesting technique, which generates spheres that are optimally sized for filling an arbitrary volume rather than for strength) and Stephen Preston (who pointed me towards papers on randomly-tessellated foams structures and wasp-inspired construction algorithms, as well as getting me thinking about artificial life-inspired algorithms).

I’ve ended up thinking along different lines – which is actually more in line with the original hyperbolic tessellation approach – but we’ll see more in the next post.

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:

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:

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.

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

So far in this series, we’ve looked at Apollonian circle packing using C# and also F#. The next few posts will look at solving this problem in 3D: performing Apollonian sphere packing. I’ve decided to stay in F# for the algorithmic side of things: it just feels a much cleaner environment for dealing with this kind of problem, and, besides, I’ve been having too much fun with it. :-)

One thing I should mention, as this series is nominally about 3D printing fill algorithms: printing hollow spheres isn’t at all straightforward with today’s 3D printing technology, as support material is needed to keep the internal arches from collapsing. This material can be water soluble, but then the material around the spherical hollows needs to be porous (or at least have some holes) to allow water to pass through it. So this series is still largely theoretical, in many ways, until such a time as hollow spheres can effectively be printed. But anyway – onwards and upwards. I’ve never been one to let feasibility to get in the way of pursuing something interesting.

Before we get onto the technical details, there’s something else I’d like to reiterate: I have been utterly dependent on the goodwill of others for the underlying mathematics/algorithms allowing creation of this geometry. I’ve begged and borrowed (I don’t like to use the term “stolen” ;-) code from a few different of places, and have even reached out to academia for help. More on that later.

When looking into the 3D side of this problem, I found relatively few websites that were ultimately of help.

Paul Bourke has a nice site, although the posted code dealt with 2D rather than 3D. This site provided some useful information – and a simple algorithm – although ultimately didn’t lead to me being able to create working code.

There were also a couple of academic papers that provided some helpful insights.

My first “breakthrough” in terms of getting actual code working was the discovery of The Mandelbot Dazibao, which included some BASIC code creating 3D Apollonian fills. It’s this code that we’ll be looking at, today.

A few words on what’s coming in the coming posts…

Before getting today’s post working, I had also reached out to Thomas Bonner, the author of ApolFrac. ApolFrac is a very cool tool to visualize – but unfortunately not export – Apollonian packings. Thomas was very helpful, pointing me at the algorithm he had used for ApolFrac, which I’d actually already seen (but obviously hadn’t understood) in one of the aforementioned academic papers.

My second breakthrough came from reaching out to Professor Ronald Peikert from ETH Zurich, who very kindly sent some C++ code that essentially implements a version of the algorithm Thomas had pointed me to. We’ll look at that in the next post.

So yes, back to breakthrough #1. Here’s the F# code I generated from The Mandelbrot Dazibao, which I believe uses the Soddy-Gosset theorem as an alternative to (or really an extension of) Descartes’ theorem.

module SpherePackingFs

 

let maxSpheres = 8000

let spheres = Array.create(maxSpheres)(0.,0.,0.,0.,0,0,0,0,0)

let mutable idc = 0

 

let createSphere (x,y,z,r) id1 id2 id3 id4 rank =

  spheres.[idc] <- x,y,z,r,id1,id2,id3,id4,rank

  idc <- idc + 1

  idc - 1

 

let getSpheres() =

  [for i in 0..(idc-1) do

    let (x,y,z,r,id1,id2,id3,id4,rank) = spheres.[i]

    yield ((x,y,z,r),rank)]

 

let generateMatrix

  (x1:double, y1:double, z1:double, r1:double)

  (x2:double, y2:double, z2:double, r2:double)

  (x3:double, y3:double, z3:double, r3:double)

  (x4:double, y4:double, z4:double, r4:double) =

 

  let a = Matrix.create 5 5 0.0

  a.Item(1,1) <- 2. * (x1 - x2)

  a.Item(2,1) <- 2. * (x2 - x3)

  a.Item(3,1) <- 2. * (x3 - x4)

  a.Item(4,1) <- 2. * (x4 - x1)

 

  a.Item(1,2) <- 2. * (y1 - y2)

  a.Item(2,2) <- 2. * (y2 - y3)

  a.Item(3,2) <- 2. * (y3 - y4)

  a.Item(4,2) <- 2. * (y4 - y1)

 

  a.Item(1,3) <- 2. * (z1 - z2)

  a.Item(2,3) <- 2. * (z2 - z3)

  a.Item(3,3) <- 2. * (z3 - z4)

  a.Item(4,3) <- 2. * (z4 - z1)

 

  a.Item(1,4) <- 2. * (r1 - r2)

  a.Item(2,4) <- 2. * (r2 - r3)

  a.Item(3,4) <- 2. * (r3 - r4)

  a.Item(4,4) <- 2. * (r4 - r1)

 

  a.Item(1,0) <-

    -(x1**2.-x2**2.+y1**2.-y2**2.+z1**2.-z2**2.+r2**2.-r1**2.)

  a.Item(2,0) <-

    -(x2**2.-x3**2.+y2**2.-y3**2.+z2**2.-z3**2.+r3**2.-r2**2.)

  a.Item(3,0) <-

    -(x3**2.-x4**2.+y3**2.-y4**2.+z3**2.-z4**2.+r4**2.-r3**2.)

  a.Item(4,0) <-

    -(x4**2.-x1**2.+y4**2.-y1**2.+z4**2.-z1**2.+r1**2.-r4**2.)

  a

 

let apollonianSphere id1 id2 id3 id4 offset =

  let (x1,y1,z1,r1,id11,id12,id13,id14,rank1) = spheres.[id1]

  let (x2,y2,z2,r2,id21,id22,id23,id24,rank2) = spheres.[id2]

  let (x3,y3,z3,r3,id31,id32,id33,id34,rank3) = spheres.[id3]

  let (x4,y4,z4,r4,id41,id42,id43,id44,rank4) = spheres.[id4]

 

  let rank = List.max [rank1;rank2;rank3;rank4]

 

  if r1 = r2 && r1 = r3 && r1= r4 then

    createSphere

      (0.,0.,(-r1 * sqrt(1./6.) + offset),(r1 * (sqrt(3./2.)-1.)))

      id1 id2 id3 id4 (rank+1)

      |> ignore

  else

 

    (*

     The four equations to be solved are:

     (x-x1)^2+(y-y1)^2+(z-z1)^2 = (r+r1)^2

     (x-x2)^2+(y-y2)^2+(z-z2)^2 = (r+r2)^2

     (x-x3)^2+(y-y3)^2+(z-z3)^2 = (r+r3)^2

     (x-x4)^2+(y-y4)^2+(z-z4)^2 = (r+r4)^2

 

     Subtract them 2 by 2:

     x*2*(x1-x2) + y*2*(y1-y2) + z*2*(z1-z2) + r*2*(r1-r2) =

      x1^2-x2^2 + y1^2-y2^2 + z1^2-z2^2 + r2^2-r1^2

     x*2*(x2-x3) + y*2*(y2-y3) + z*2*(z2-z3) + r*2*(r2-r3) =

      x2^2-x3^2 + y2^2-y3^2 + z2^2-z3^2 + r3^2-r2^2

     x*2*(x3-x4) + y*2*(y3-y4) + z*2*(z3-z4) + r*2*(r3-r4) =

      x3^2-x4^2 + y3^2-y4^2 + z3^2-z4^2 + r4^2-r3^2

     x*2*(x4-x1) + y*2*(y4-y1) + z*2*(z4-z1) + r*2*(r4-r1) =

      x4^2-x1^2 + y4^2-y1^2 + z4^2-z1^2 + r1^2-r4^2

 

     Matrix writing:

     a(1,1)*x+a(1,2)*y+a(1,3)*z+a(1,4)*r+a(1,0) = 0

     a(2,1)*x+a(2,2)*y+a(2,3)*z+a(2,4)*r+a(2,0) = 0

     a(3,1)*x+a(3,2)*y+a(3,3)*z+a(3,4)*r+a(3,0) = 0

     a(4,1)*x+a(4,2)*y+a(4,3)*z+a(4,4)*r+a(4,0) = 0

    *)

 

    let mutable a =

      generateMatrix

        (x1,y1,z1,r1) (x2,y2,z2,r2) (x3,y3,z3,r3) (x4,y4,z4,r4)

 

    if a.Item(3,3) = 0. then

      a <-

        generateMatrix

          (x4,y4,z4,r4) (x2,y2,z2,r2) (x3,y3,z3,r3) (x1,y1,z1,r1)

   (*

     Get x, y and z as functions of r

     a(1,1)*x+a(1,2)*y+a(1,3)*z=-a(1,4)*r-a(1,0)

     a(2,1)*x+a(2,2)*y+a(2,3)*z=-a(2,4)*r-a(2,0)

     a(3,1)*x+a(3,2)*y+a(3,3)*z=-a(3,4)*r-a(3,0)

    *)

 

    let det (a : matrix) m n p q =

      a.Item(m,p) * a.Item(n,q) - a.Item(n,p) * a.Item(m,q)

 

    let D1313 = det a 1 3 1 3

    let D2323 = det a 2 3 2 3

    let D1323 = det a 1 3 2 3

    let D2313 = det a 2 3 1 3

    let DetMajor = D2313 * D1323 - D1313 * D2323

 

    let RX =

      (a.Item(1,4) * a.Item(3,3) * D2323 - a.Item(2,4) *

       a.Item(3,3) * D1323 - a.Item(3,4) * a.Item(1,3) *

       D2323 + a.Item(3,4) * a.Item(2,3) * D1323) / DetMajor

 

    let RY =

      (-a.Item(1,4) * a.Item(3,3) * D2313 + a.Item(2,4) *

       a.Item(3,3) * D1313 - a.Item(3,4) * a.Item(2,3) *

       D1313 + a.Item(3,4) * a.Item(1,3) * D2313) / DetMajor

 

    let RX0 =

      (a.Item(1,0) * a.Item(3,3) * D2323 - a.Item(2,0) *

       a.Item(3,3) * D1323 - a.Item(3,0) * a.Item(1,3) *

       D2323 + a.Item(3,0) * a.Item(2,3) * D1323) / DetMajor

 

    let RY0 =

      (-a.Item(1,0) * a.Item(3,3) * D2313 + a.Item(2,0) *

       a.Item(3,3) * D1313 - a.Item(3,0) * a.Item(2,3) *

       D1313 + a.Item(3,0) * a.Item(1,3) * D2313) / DetMajor

 

    let RZ =

      -a.Item(3,4) / a.Item(3,3) - RX * a.Item(3,1) / a.Item(3,3) -

      a.Item(3,2) / a.Item(3,3) * RY

 

    let RZ0 =

      -a.Item(3,0) / a.Item(3,3) - RX0 * a.Item(3,1) / a.Item(3,3) -

      RY0 * a.Item(3,2) / a.Item(3,3)

 

    let Pp = RX**2. + RY**2. + RZ**2. - 1.

    let Pq =

      2. * (RX * (RX0-x4) + RY * (RY0-y4) + RZ * (RZ0-z4) - r4)

    let Pr =

      (RX0-x4)**2. + (RY0-y4)**2. + (RZ0-z4)**2. - r4**2.

 

    let Delta = Pq**2. - 4. * Pp * Pr

    let Radius = (-Pq - sqrt(Delta)) / 2. / Pp

 

    let x = RX * Radius + RX0

    let y = RY * Radius + RY0

    let z = RZ * Radius + RZ0

 

    createSphere (x,y,z,Radius) id1 id2 id3 id4 (rank+1)

      |> ignore

 

let fleshOut maxRank offset =

  for curRank in 1..maxRank do

    for id in 1..(idc-1) do

      let (x,y,z,r,id1,id2,id3,id4,rank) = spheres.[id]

      if rank = curRank then

        apollonianSphere id id1 id2 id3 offset

        apollonianSphere id id1 id2 id4 offset

        apollonianSphere id id1 id3 id4 offset

        apollonianSphere id id2 id3 id4 offset

 

let exteriorGasket c0 c1 c2 c3 c4 offset steps =

 

  idc <- 0

 

  let c0 = createSphere c0 0 0 0 0 0

  let c1 = createSphere c1 0 0 0 0 0

  let c2 = createSphere c2 0 0 0 0 0

  let c3 = createSphere c3 0 0 0 0 0

  let c4 = createSphere c4 0 0 0 0 0

 

  apollonianSphere c0 c1 c2 c3 offset

  apollonianSphere c0 c1 c2 c4 offset

  apollonianSphere c0 c1 c3 c4 offset

  apollonianSphere c0 c2 c3 c4 offset

  apollonianSphere c1 c2 c3 c4 offset

 

  fleshOut steps offset

 

  getSpheres()

 

let interiorGasket c0 c1 c2 c3 offset steps =

 

  idc <- 0

 

  let c0 = createSphere c0 0 0 0 0 0

  let c1 = createSphere c1 0 0 0 0 0

  let c2 = createSphere c2 0 0 0 0 0

  let c3 = createSphere c3 0 0 0 0 0

 

  apollonianSphere c0 c1 c2 c3 offset

 

  fleshOut steps offset

 

  getSpheres()

 

type Packer() =

 

  static member ExteriorApollonianGasket steps =

    let size = 500.

    let rt1o6 = sqrt(1./6.)

    let rt3 = sqrt(3.)

    let rt2o3 = sqrt(2./3.)

    let offset = size * rt1o6

    let outerRad = size * (2. + (1. / (2. * rt1o6)))

 

    exteriorGasket

      (0., 0., (offset - (size * rt1o6)), outerRad)

      (((2. * size) / rt3), 0., offset, size)

      (0., 0., (offset - (2. * size * rt2o3)), size)

      ((-size / rt3), size, offset, size)

      ((-size / rt3), -size, offset, size)

      offset steps

 

  static member InteriorApollonianGasket steps =

    let size = 500.

    let rt3 = sqrt(3.)

    let rt2o3 = sqrt(2./3.)

    let offset = size * 0.5

 

    interiorGasket

      ((2. * size / rt3), 0., offset, size)

      (0., 0., (offset - (2. * size * rt2o3)), size)

      ((-size / rt3), size, offset, size)

      ((-size / rt3), -size, offset, size)

      offset steps

As the original code was in BASIC, I’ve found it simplest to keep certain approaches in place: rather than building – and returning – a list of spheres to create, the code builds an fixed-size array and populates it with the results. Which means there are fundamental limits to the number of spheres that can be generated (and if you increase the number of levels to generate, you also need to be aware you may need to increase the size of the array). This isn’t ideal, by any means, but I’ve left the implementation fairly similar to the original, especially as it’s ultimately being superseded by the code in the next post.

To make use of the F# code, I’ve added this C# loader to the project (which I’ll post next time):

using Autodesk.AutoCAD.ApplicationServices;

using Autodesk.AutoCAD.DatabaseServices;

using Autodesk.AutoCAD.Geometry;

using Autodesk.AutoCAD.Runtime;

 

namespace SpherePackingLoader

{

  public class Commands

  {

    [CommandMethod("AGFSE")]

    public static void ExteriorApollonianGasket()

    {

      ApollonianGasket(true);

    }

 

    [CommandMethod("AGFSI")]

    public static void InteriorApollonianGasket()

    {

      ApollonianGasket(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 exterior)

    {

      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 =

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

      }

    }

  }

}

It defines two commands, AGFSI and AGFSE, creating interior and exterior gaskets, respectively. To aid with visualizing the various levels, the C# code places the resultant spheres on corresponding layers.

Running the AGFSI command – and turning off layer 0, after changing the current later, of course – generates a nice interior Apollonian gasket:

Running the AGFSE command creates an external packing, although I’m somewhat loathe to call it Apollonian. I suspect through some transcription error on my part it’s not a perfect packing – you can see some intersecting spheres, for instance. I also had to manually delete some outer spheres to get to this geometry.

Anyway – rather than spend time perfecting the irrelevant, I’m leaving this implementation as it stands. After the weekend I’ll post a different – and more elegant – implementation that does a much better job of generating an Apollonian packing in 3D.

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:

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

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#.

Circle packing in AutoCAD: creating an Apollonian gasket using .NET

To follow on from the recent series on using hyperbolic tessellation to generate patterns that might be used for 3D printing, I decided to research a slightly different approach.

While I found hyperbolic tessellation reasonably straightforward for generating 2D patterns, it was much harder to adapt to 3D, mainly because we’d need to create irregular polyhedra rather than just irregular 2D polygons. I had enough trouble getting my head around the 2D side of things that I decided not to go down that path for now, at least. The good news from Alex Fielder is that he’s been making good progress getting the hyperbolic tessellation code working inside Inventor (very cool), focusing – as far as I know – on using the 2D pattern to generate ribs rather than polyhedra (a very pragmatic approach).

I do think the type of pattern created through hyperbolic tessellation is probably more interesting from a structural perspective (having smaller versions at the edges with larger ones at the middle just makes a lot of sense), especially if working with polyhedra rather than ribs.

I expect to come back to something similar, in due course, but for now have been looking into a different approach: thinking about more regular (and therefore simpler) geometric forms by looking at circle packing before moving on to the very much related sphere packing.

There are a number of interesting strategies for circle packing: the one I like best is called the Apollonian gasket (which has been covered on a number of different Wikipedia pages, as well as on Wolfram MathWorld), which is also known as Leibniz packing. A related, specific case is that of Ford circles. The idea is that after creating three initial circles (each tangent to each other and to the outer circle), you draw the largest possible circle in the remaining space, continuing to sub-divide fractally as far as you’d like.

[As an aside… while researching this topic, I came across this interesting site by Takaya Iwamoto, a former colleague at Autodesk, who clearly shares my interest in fractals.]

And this of course also applies to 3D, although it may require a slightly different approach to get there.

I really like this pattern: while it may not create an internal structure optimised for strength, it should still be pretty interesting from a material savings perspective (short of having a completely hollow shell, of course). And I fully expect to looking into more structural filling techniques, further down the line.

The main problem to solve when generating the above pattern is to find a fourth mutually tangent circle (in 2D) or a fifth mutually tangent sphere (in 3D). It’s possible to use Descartes’ theorem in 2D or the Soddy-Gosset theorem in 3D (and beyond, for the truly fearless).

My mathematics being as they are (i.e rusty at best), I did my best to find someone else’s code to adapt for the various posts in this series. I started with a C# implementation (it also has a VB.NET sibling), and reworked it to generate AutoCAD geometry.

using Autodesk.AutoCAD.ApplicationServices;

using Autodesk.AutoCAD.DatabaseServices;

using Autodesk.AutoCAD.EditorInput;

using Autodesk.AutoCAD.Geometry;

using Autodesk.AutoCAD.Runtime;

using System.Collections.Generic;

using System;

 

namespace CirclePacking

{

  public class Commands

  {

    static double _minRad = 0.01;

    static short _color = 5;

 

    public class MyCircle

    {

      public Point2d Center;

      public double Radius;

      public short ColorIndex;

    }

 

    [CommandMethod("AGC")]

    public static void ApollonianGasket()

    {

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

        }

 

        List<MyCircle> circles =

          FindApollonianPacking(cir.Radius * 2);

 

        // Offset each of the returned circles to fit our outer one

 

        Vector3d offset =

          cir.Center - Point3d.Origin -

          new Vector3d(cir.Radius, cir.Radius, 0);

 

        // Add each of the returned circles to the current space

 

        foreach (MyCircle myc in circles)

        {

          Point3d cen =

            new Point3d(myc.Center.X, myc.Center.Y, 0) + offset;

          Circle c = new Circle(cen, Vector3d.ZAxis, myc.Radius);

          c.ColorIndex = myc.ColorIndex;

 

          btr.AppendEntity(c);

          tr.AddNewlyCreatedDBObject(c, true);

        }

 

        tr.Commit();

      }

    }

 

    // Find the Apollonian packing

 

    private static List<MyCircle> FindApollonianPacking(

      double width

    )

    {

      List<MyCircle> circles = new List<MyCircle>();

 

      // Create the three central tangent circles

 

      double radius = width * 0.232f,

             x = width / 2,

             gaskHeight = 2 * (radius + 2 * radius / Math.Sqrt(3)),

             y = (width - gaskHeight) / 2 + radius;

 

      MyCircle cir0 = CreateCircle(x, y, radius);

      circles.Add(cir0);

 

      x -= radius;

      y += radius * Math.Sqrt(3);

 

      MyCircle cir1 = CreateCircle(x, y, radius);

      circles.Add(cir1);

 

      x += 2 * radius;

 

      MyCircle cir2 = CreateCircle(x, y, radius);

      circles.Add(cir2);

 

      // Find the circle that contains them all

 

      MyCircle outer =

        FindApollonianCircle(cir0, cir1, cir2, -1, -1, -1);

      circles.Add(outer);

 

      // Set level to smaller values such as 3 to see partially

      // drawn gaskets

 

      int level = 10000;

 

      // Find the central circle

 

      FindCircleOutsideAll(circles, level, cir0, cir1, cir2);

 

      // Find circles tangent to the big circle

 

      FindCircleOutsideTwo(circles, level, cir0, cir1, outer);

      FindCircleOutsideTwo(circles, level, cir1, cir2, outer);

      FindCircleOutsideTwo(circles, level, cir2, cir0, outer);

 

      return circles;

    }

 

    // Draw a circle tangent to these three circles and that is

    // outside all three

 

    private static void FindCircleOutsideAll(

      List<MyCircle> circles,

      int level, MyCircle cir0, MyCircle cir1, MyCircle cir2

    )

    {

      MyCircle newCir =

        FindApollonianCircle(cir0, cir1, cir2, 1, 1, 1);

      if (newCir == null || newCir.Radius < _minRad)

        return;

      circles.Add(newCir);

 

      if (--level > 0)

      {

        FindCircleOutsideAll(circles, level, cir0, cir1, newCir);

        FindCircleOutsideAll(circles, level, cir0, cir2, newCir);

        FindCircleOutsideAll(circles, level, cir1, cir2, newCir);

      }

    }

 

    // Draw a circle tangent to these three circles and that is

    // outside two of them

 

    private static void FindCircleOutsideTwo(

      List<MyCircle> circles, int level,

      MyCircle cir0, MyCircle cir1, MyCircle contains

    )

    {

      MyCircle newCir =

        FindApollonianCircle(cir0, cir1, contains, 1, 1, -1);

      if (newCir == null || newCir.Radius < _minRad)

        return;

      circles.Add(newCir);

 

      if (--level > 0)

      {

        FindCircleOutsideTwo(circles, level, newCir, cir0, contains);

        FindCircleOutsideTwo(circles, level, newCir, cir1, contains);

        FindCircleOutsideAll(circles, level, cir0, cir1, newCir);

      }

    }

 

    private static MyCircle FindApollonianCircle(

      MyCircle c1, MyCircle c2, MyCircle c3, int s1, int s2, int s3

    )

    {

      // Make sure c2 doesn't have the same X or Y coordinate as

      // the others

 

      double tiny = Tolerance.Global.EqualPoint;

 

      if (

        Math.Abs(c2.Center.X - c1.Center.X) < tiny ||

        Math.Abs(c2.Center.Y - c1.Center.Y) < tiny

      )

      {

        MyCircle tempCir = c2;

        c2 = c3;

        c3 = tempCir;

        int temp_s = s2;

        s2 = s3;

        s3 = temp_s;

      }

 

      if (

        Math.Abs(c2.Center.X - c3.Center.X) < tiny ||

        Math.Abs(c2.Center.Y - c3.Center.Y) < tiny

      )

      {

        MyCircle tempCir = c2;

        c2 = c1;

        c1 = tempCir;

        int temp_s = s2;

        s2 = s1;

        s1 = temp_s;

      }

 

      double x1 = c1.Center.X,

             y1 = c1.Center.Y,

             r1 = c1.Radius,

             x2 = c2.Center.X,

             y2 = c2.Center.Y,

             r2 = c2.Radius,

             x3 = c3.Center.X,

             y3 = c3.Center.Y,

             r3 = c3.Radius;

 

      double v11 = 2 * x2 - 2 * x1,

             v12 = 2 * y2 - 2 * y1,

             v13 =

               x1 * x1 - x2 * x2 + y1 * y1 -

               y2 * y2 - r1 * r1 + r2 * r2,

             v14 = 2 * s2 * r2 - 2 * s1 * r1;

 

      double v21 = 2 * x3 - 2 * x2,

             v22 = 2 * y3 - 2 * y2,

             v23 =

               x2 * x2 - x3 * x3 + y2 * y2 - y3 * y3 -

               r2 * r2 + r3 * r3,

             v24 = 2 * s3 * r3 - 2 * s2 * r2;

 

      double w12 = v12 / v11,

             w13 = v13 / v11,

             w14 = v14 / v11;

 

      double w22 = v22 / v21 - w12,

             w23 = v23 / v21 - w13,

             w24 = v24 / v21 - w14;

 

      double P = -w23 / w22,

             Q = w24 / w22,

             M = -w12 * P - w13,

             N = w14 - w12 * Q;

 

      double a = N * N + Q * Q - 1,

             b =

               2 * M * N - 2 * N * x1 + 2 * P * Q -

               2 * Q * y1 + 2 * s1 * r1,

             c =

               x1 * x1 + M * M - 2 * M * x1 + P * P +

               y1 * y1 - 2 * P * y1 - r1 * r1;

 

      // Find roots of a quadratic equation

 

      double[] solutions = QuadraticSolutions(a, b, c);

 

      if (solutions.Length < 1)

        return null;

 

      double rs = solutions[0],

             xs = M + N * rs,

             ys = P + Q * rs;

 

      if (

        Math.Abs(xs) < tiny ||

        Math.Abs(ys) < tiny ||

        Math.Abs(rs) < tiny

      )

        return null;

 

      return CreateCircle(xs, ys, Math.Abs(rs));

    }

 

    // Return solutions to a quadratic equation

 

    private static double[] QuadraticSolutions(

      double a, double b, double c

    )

    {

      double discriminant = b * b - 4 * a * c;

 

      // See if there are no real solutions

 

      if (discriminant < 0)

        return new double[] { };

 

      // See if there is one solution

 

      if (discriminant < Tolerance.Global.EqualPoint)

        return new double[] { -b / (2 * a) };

 

      // There are two solutions

 

      return new double[]

      {

        (-b + Math.Sqrt(discriminant)) / (2 * a),

        (-b - Math.Sqrt(discriminant)) / (2 * a),

      };

    }

 

    private static MyCircle CreateCircle(

      double x, double y, double radius

    )

    {

      MyCircle cir = new MyCircle();

      cir.Center = new Point2d(x, y);

      cir.Radius = radius;

      cir.ColorIndex = _color;

      return cir;

    }

  }

}

The AGC command allows you to select an existing circle to fill or it creates one at the origin and fills that. One thing to bear in mind: with this implementation the size of the circle does matter: as the fractal stops at circles of radius 0.01, an outer circle of 500 radius is going to lead to more geometry that one with a radius of 10 (for instance).

The implementation is reasonably straightforward – the main function, FindApollonianCircle(), is probably the hardest to understand – but I tried not to change anything very much with the algorithm (which means, just like last time, there may well be more elegant ways to solve the problem).

Throughout this series of investigations, I’ve been itching to switch across to F# to see what’s possible, there. In the next post or two in the series, we’ll take a look at implementing an Apollonian gasket in F# before moving on to the ultimately more relevant world of 3D.

Scripting the generation of hyperbolic tessellations inside AutoCAD

To follow on from this recent topic, today’s post looks at a simple script to generate various hyperbolic tessellations, laying them out in an order that makes some sense of the progressive nature of the patterns that can be generated using the HT command.

Here’s an AutoCAD script (which can be saved as an .scr and executed using the SCRIPT command) to generate all the valid {n k} patterns where n <= 11 and k <= 7 (n being the number of sides in each polygon, k is the number of polygons that meet at each vertex). The application module implementing the HT command (using the code in the previous post, must be NETLOADed before this script will function properly, of course).

_zoom _w 0,0 50,25

; k=3 series

_text _j _m 25,5 0.2 0 {7 3}

_circle 25,5 2 _ht l 7 3 _s

_text _j _m 30,5 0.2 0 {8 3}

_circle 30,5 2 _ht l 8 3 _s

_text _j _m 35,5 0.2 0 {9 3}

_circle 35,5 2 _ht l 9 3 _s

_text _j _m 40,5 0.2 0 {10 3}

_circle 40,5 2 _ht l 10 3 _s

_text _j _m 45,5 0.2 0 {11 3}

_circle 45,5 2 _ht l 11 3 _s

; k=4 series

_text _j _m 15,10 0.2 0 {5 4}

_circle 15,10 2 _ht l 5 4 _s

_text _j _m 20,10 0.2 0 {6 4}

_circle 20,10 2 _ht l 6 4 _s

_text _j _m 25,10 0.2 0 {7 4}

_circle 25,10 2 _ht l 7 4 _s

_text _j _m 30,10 0.2 0 {8 4}

_circle 30,10 2 _ht l 8 4 _s

_text _j _m 35,10 0.2 0 {9 4}

_circle 35,10 2 _ht l 9 4 _s

_text _j _m 40,10 0.2 0 {10 4}

_circle 40,10 2 _ht l 10 4 _s

_text _j _m 45,10 0.2 0 {11 4}

_circle 45,10 2 _ht l 11 4 _s

; k=5 series

_text _j _m 10,15 0.2 0 {4 5}

_circle 10,15 2 _ht l 4 5 _s

_text _j _m 15,15 0.2 0 {5 5}

_circle 15,15 2 _ht l 5 5 _s

_text _j _m 20,15 0.2 0 {6 5}

_circle 20,15 2 _ht l 6 5 _s

_text _j _m 25,15 0.2 0 {7 5}

_circle 25,15 2 _ht l 7 5 _s

_text _j _m 30,15 0.2 0 {8 5}

_circle 30,15 2 _ht l 8 5 _s

_text _j _m 35,15 0.2 0 {9 5}

_circle 35,15 2 _ht l 9 5 _s

_text _j _m 40,15 0.2 0 {10 5}

_circle 40,15 2 _ht l 10 5 _s

_text _j _m 45,15 0.2 0 {11 5}

_circle 45,15 2 _ht l 11 5 _s

; k=6 series

_text _j _m 10,20 0.2 0 {4 6}

_circle 10,20 2 _ht l 4 6 _s

_text _j _m 15,20 0.2 0 {5 6}

_circle 15,20 2 _ht l 5 6 _s

_text _j _m 20,20 0.2 0 {6 6}

_circle 20,20 2 _ht l 6 6 _s

_text _j _m 25,20 0.2 0 {7 6}

_circle 25,20 2 _ht l 7 6 _s

_text _j _m 30,20 0.2 0 {8 6}

_circle 30,20 2 _ht l 8 6 _s

_text _j _m 35,20 0.2 0 {9 6}

_circle 35,20 2 _ht l 9 6 _s

_text _j _m 40,20 0.2 0 {10 6}

_circle 40,20 2 _ht l 10 6 _s

_text _j _m 45,20 0.2 0 {11 6}

_circle 45,20 2 _ht l 11 6 _s

; k=7 series

_text _j _m 5,25 0.2 0 {3 7}

_circle 5,25 2 _ht l 3 7 _s

_text _j _m 10,25 0.2 0 {4 7}

_circle 10,25 2 _ht l 4 7 _s

_text _j _m 15,25 0.2 0 {5 7}

_circle 15,25 2 _ht l 5 7 _s

_text _j _m 20,25 0.2 0 {6 7}

_circle 20,25 2 _ht l 6 7 _s

_text _j _m 25,25 0.2 0 {7 7}

_circle 25,25 2 _ht l 7 7 _s

_text _j _m 30,25 0.2 0 {8 7}

_circle 30,25 2 _ht l 8 7 _s

_text _j _m 35,25 0.2 0 {9 7}

_circle 35,25 2 _ht l 9 7 _s

_text _j _m 40,25 0.2 0 {10 7}

_circle 40,25 2 _ht l 10 7 _s

_text _j _m 45,25 0.2 0 {11 7}

_circle 45,25 2 _ht l 11 7 _s

_zoom _e

Just looking at the output results (extracted manually from the command-line output), we can see the level used and the number of polygons created for each tessellation:

{7 3} (level 5) => 617 polygons.

{8 3} (level 4) => 609 polygons.

{9 3} (level 3) => 271 polygons.

{10 3} (level 3) => 421 polygons.

{11 3} (level 3) => 617 polygons.

 

{5 4} (level 4) => 761 polygons.

{6 4} (level 3) => 505 polygons.

{7 4} (level 2) => 127 polygons.

{8 4} (level 2) => 177 polygons.

{9 4} (level 2) => 235 polygons.

{10 4} (level 2) => 301 polygons.

{11 4} (level 2) => 375 polygons.

 

{4 5} (level 4) => 913 polygons.

{5 5} (level 3) => 841 polygons.

{6 5} (level 2) => 199 polygons.

{7 5} (level 2) => 295 polygons.

{8 5} (level 2) => 409 polygons.

{9 5} (level 2) => 541 polygons.

{10 5} (level 2) => 691 polygons.

{11 5} (level 2) => 859 polygons.

 

{4 6} (level 3) => 673 polygons.

{5 6} (level 2) => 221 polygons.

{6 6} (level 2) => 361 polygons.

{7 6} (level 2) => 533 polygons.

{8 6} (level 2) => 737 polygons.

{9 6} (level 2) => 973 polygons.

{10 6} (level 1) => 41 polygons.

{11 6} (level 1) => 45 polygons.

 

{3 7} (level 4) => 496 polygons.

{4 7} (level 2) => 181 polygons.

{5 7} (level 2) => 351 polygons.

{6 7} (level 2) => 571 polygons.

{7 7} (level 2) => 841 polygons.

{8 7} (level 1) => 41 polygons.

{9 7} (level 1) => 46 polygons.

{10 7} (level 1) => 51 polygons.

{11 7} (level 1) => 56 polygons.

And here are the patterns, themselves:

Just to see the difference, here are the same tessellations with curved geometry (as per the implementation in the second post in the series). Created using a simple search & replace in the above script, from “_s” to “_c”.

Thinking about future directions for this investigation… (some of which have been suggested by Alex Fielder, after previewing drafts of these posts…)

• Clearly it’d be great to get to 3D and get the results integrated with Inventor and Autodesk’s Simulation (FEA) tools.
• Remaining in 2D, it’d be interesting to see about supporting non-circular boundaries.
  • One option for this would be to use the minimal enclosing circle around the geometry as the boundary, and then trim back the results or filter them based on whether they’re inside or outside/intersecting the enclosing geometry.
  • Ideally the fill pattern would be generated in a way that followed the shape being filled, rather than relying on a circle at all. That seems more of a stretch, though, given the nature of the geometry.

If you have your own ideas on uses and directions for this implementation, please post a comment!

Generating hyperbolic tessellations inside AutoCAD using .NET

Following on from the post introducing this series, and the last post focused on generating simple hyperbolic geometry, this post looks at generating hyperbolic tessellations inside AutoCAD.

Having “borrowed” some C++ code, last time, today we’re going to borrow some Java. That’s one of the great things about the C family of languages: the relatively small amount of work that’s often needed to move code between environments. :-)

This article, written by David E. Joyce, Professor of Mathematics and Computer Science at Clark University, has provided the best information I could find with regards to a hyperbolic tiling algorithm (especially thanks to its posted source code).

The algorithm works by creating a central polygon of the specified number of sides, and then – based on the other primary parameter into the algorithm, the number of polygons that meet at each vertex – it mirrors this polygon again and again in the hyperbolic plane (which means it gets smaller and a bit squished as it gets mirrored out towards the edge of the circle) in order to make up our pattern. The number of mirroring operations – and the number of resultant polygons – will depend on the number of levels (or layers): 0 means we have just the central polygon, 1 means we have that plus the layer around it, etc.

The number of levels needed to effectively fill a circle depends on the resolution needed but also greatly on the pattern. Interestingly the more even patterns (such as 7-sided polygons meeting 3 per vertex) create fewer polygons per level, but need more levels to fill a circle adequately. The patterns with larger central polygons create many more polygons per level (which is probably not that surprising, thinking about it), and so need fewer levels to fill a circle.

Rather than ask the user for the number of levels, this implementation has a coded limit on the number of polygons (which, of course, could also be asked of the user, if needed). The code starts with a default level of 6 (which is pretty high for most patterns), but quickly reduces the level by one if it finds it would result in more than 1000 polygons being created. Once it has found a level number that creates a small enough of set of polygons, it goes ahead and uses that.

In addition to the HyperbolicEngine we saw in the last post, we’re now introducing a PolygonEngine to generate polygons for a particular pattern. Here’s the C# file:

using Autodesk.AutoCAD.Geometry;

using System;

 

namespace HyperbolicGeometry

{

  public class Polygon

  {

    public Point2dCollection Vertices = new Point2dCollection();

  }

 

  public class PolygonEngine

  {

    // The list of polygons

 

    Polygon[] _polys;

 

    // Previously created neighbors for the polygons

 

    int[] _rules;

 

    // The total number of polygons in all the layers

 

    int _totalPolygons;

 

    // The number through one less layer

 

    int _innerPolygons;

 

    public Polygon[] GetPolygons(

      int n, int k, ref int layNum, int max

    )

    {

      do

      {

        CountPolygons(layNum, n, k, max);

        if (_totalPolygons > max)

          layNum--;

      }

      while (_totalPolygons > max);

 

      DeterminePolygons(n, k);

 

      return _polys;

    }

 

    private static Polygon ConstructCenterPolygon(int n, int k)

    {

      // Initialize P as the center polygon in an n-k regular tiling

      // Let ABC be a triangle in a regular (n,k0-tiling, where

      //    A is the center of an n-gon (also center of the disk),

      //    B is a vertex of the n-gon, and

      //    C is the midpoint of a side of the n-gon adjacent to B

 

      double angleA = Math.PI / n;

      double angleB = Math.PI / k;

      double angleC = Math.PI / 2.0;

 

      // For a regular tiling, we need to compute the distance s

      // from A to B

 

      double sinA = Math.Sin(angleA);

      double sinB = Math.Sin(angleB);

      double s =

        Math.Sin(angleC - angleB - angleA) /

        Math.Sqrt(1.0 - sinB * sinB - sinA * sinA);

 

      // Determine the coordinates of the n vertices of the n-gon.

      // They're all at distance s from center of the Poincare disk

 

      Polygon p = new Polygon();

 

      for (int i = 0; i < n; ++i)

      {

        p.Vertices.Add(

          new Point2d(

            s * Math.Cos((3 + 2 * i) * angleA),

            s * Math.Sin((3 + 2 * i) * angleA)

          )

        );

      }

      return p;

    }

 

    private void CountPolygons(int layer, int n, int k, int max)

    {

      // Determine

      //   totalPolygons:  the number of polygons there are through

      //                   that many layers

      //   innerPolygons:  the number through one less layer

 

      _totalPolygons = 1;  // count the central polygon

      _innerPolygons = 0;

 

      int a = n * (k - 3); // polygons in 1st layer joined by vertex

      int b = n;           // polygons in 1st layer joined by edge

 

      int next_a, next_b;

      for (int l = 1; l <= layer; ++l)

      {

        _innerPolygons = _totalPolygons;

        if (k == 3)

        {

          next_a = a + b;

          next_b = (n - 6) * a + (n - 5) * b;

        }

        else

        { // k > 3

          next_a =

            ((n - 2) * (k - 3) - 1) * a +

            ((n - 3) * (k - 3) - 1) * b;

          next_b = (n - 2) * a + (n - 3) * b;

        }

        _totalPolygons += a + b;

 

        // If we've reached our maximum threshold, just return,

        // as we won't be using this layer, anyway

 

        if (_totalPolygons > max)

          break;

 

        a = next_a;

        b = next_b;

      }

    }

 

    /* Rule codes

     *   0:  initial polygon.  Needs neighbors on all n sides

     *   1:  polygon already has 2 neighbors, but one less

     *       around corner needed

     *   2:  polygon already has 2 neighbors

     *   3:  polygon already has 3 neighbors

     *   4:  polygon already has 4 neighbors

     */

 

    private void DeterminePolygons(int n, int k)

    {

      _polys = new Polygon[_totalPolygons];

      _rules = new int[_totalPolygons];

      _polys[0] = ConstructCenterPolygon(n, k);

      _rules[0] = 0;

      int j = 1; // index of the next polygon to create

      for (int i = 0; i < _innerPolygons; ++i)

      {

        j = ApplyRule(i, j, n, k);

      }

    }

 

    private int ApplyRule(int i, int j, int n, int k)

    {

      int r = _rules[i];

      bool special = (r == 1);

      if (special) r = 2;

      int start = (r == 4) ? 3 : 2;

      int quantity = (k == 3 && r != 0) ? n - r - 1 : n - r;

      for (int s = start; s < start + quantity; ++s)

      {

        // Create a polygon adjacent to P[i]

 

        _polys[j] = CreateNextPolygon(_polys[i], s % n, n, k);

        _rules[j] = (k == 3 && s == start && r != 0) ? 4 : 3;

        j++;

        int m;

        if (special) m = 2;

        else if (s == 2 && r != 0) m = 1;

        else m = 0;

        for (; m < k - 3; ++m)

        {

          // Create a polygon adjacent to P[j-1]

 

          _polys[j] = CreateNextPolygon(_polys[j - 1], 1, n, k);

          _rules[j] = (n == 3 && m == k - 4) ? 1 : 2;

          j++;

        }

      }

      return j;

    }

 

    // Reflect P thru the point or the side indicated by the side s

    // to produce the resulting polygon Q

 

    private Polygon CreateNextPolygon(

      Polygon P, int s, int n, int k

    )

    {

      Point2d start = P.Vertices[s],

              end = P.Vertices[(s + 1) % n];

 

      Polygon Q = new Polygon();

      for (int i = 0; i < n; i++)

        Q.Vertices.Add(Point2d.Origin);

 

      for (int i = 0; i < n; ++i)

      {

        // Reflect P[i] thru C to get Q[j]}

 

        int j = (n + s - i + 1) % n;

        Q.Vertices[j] = Reflect(start, end, P.Vertices[i]);

      }

      return Q;

    }

 

    // Reflect the point R across the line through A to B

 

    public Point2d Reflect(Point2d A, Point2d B, Point2d R)

    {

      // Find a unit vector D in the direction of the line

 

      double den = A.X * B.Y - B.X * A.Y;

      bool straight = Math.Abs(den) < Tolerance.Global.EqualPoint;

 

      if (straight)

      {

        // The less interesting case - we could do this less

        // manually using AcGe

 

        Point2d P = A;

        den =

          Math.Sqrt(

            (A.X - B.X) * (A.X - B.X) + (A.Y - B.Y) * (A.Y - B.Y)

          );

        Point2d D =

          new Point2d(

            (B.X - A.X) / den,

            (B.Y - A.Y) / den

          );

 

        // Reflect method

 

        double factor =

          2.0 * ((R.X - P.X) * D.X + (R.Y - P.Y) * D.Y);

        return

          new Point2d(

            2.0 * P.X + factor * D.X - R.X,

            2.0 * P.Y + factor * D.Y - R.Y

          );

      }

      else

      {

        // Find the center of the circle through these points

        // (this is what we really need to use this for)

 

        double s1 = (1.0 + A.X * A.X + A.Y * A.Y) / 2.0;

        double s2 = (1.0 + B.X * B.X + B.Y * B.Y) / 2.0;

 

        Point2d C =

          new Point2d(

            (s1 * B.Y - s2 * A.Y) / den,

            (A.X * s2 - B.X * s1) / den

          );

        double r = Math.Sqrt(C.X * C.X + C.Y * C.Y - 1.0);

 

        // Reflect method

 

        double factor =

          r * r / (

            (R.X - C.X) * (R.X - C.X) + (R.Y - C.Y) * (R.Y - C.Y)

          );

        return

          new Point2d(

            C.X + factor * (R.X - C.X),

            C.Y + factor * (R.Y - C.Y)

          );

      }

    }

  }

}

Here’s the updated command class that includes a new command, HT, to generate polygons for a hyperbolic tessellation with a given set of parameters:

using Autodesk.AutoCAD.ApplicationServices;

using Autodesk.AutoCAD.Colors;

using Autodesk.AutoCAD.DatabaseServices;

using Autodesk.AutoCAD.EditorInput;

using Autodesk.AutoCAD.Geometry;

using Autodesk.AutoCAD.Runtime;

using System;

 

namespace HyperbolicGeometry

{

  public class Commands

  {

    [CommandMethod("HT")]

    public static void CreateHyperbolicTiles()

    {

      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;

 

      int n, k;

      bool done;

 

      do

      {

        PromptIntegerOptions pio =

          new PromptIntegerOptions(

            "\nNumber of sides for polygons"

          );

        pio.AllowZero = false;

        pio.AllowNegative = false;

        pio.AllowNone = true;

        pio.DefaultValue = 7;

        pio.UseDefaultValue = true;

        pio.LowerLimit = 3;

        pio.UpperLimit = 20;

 

        PromptIntegerResult pir = ed.GetInteger(pio);

 

        if (pir.Status != PromptStatus.None &&

            pir.Status != PromptStatus.OK)

          return;

 

        n = pio.DefaultValue;

 

        if (pir.Status == PromptStatus.OK)

          n = pir.Value;

 

        pio.Message =

          "\nNumber of polygons meeting at each vertex";

        pio.DefaultValue = 3;

 

        pir = ed.GetInteger(pio);

 

        if (pir.Status != PromptStatus.None &&

            pir.Status != PromptStatus.OK)

          return;

 

        k = pio.DefaultValue;

 

        if (pir.Status == PromptStatus.OK)

          k = pir.Value;

 

        done = ((n - 2) * (k - 2) > 4);

 

        if (!done)

        {

          ed.WriteMessage(

            "\nPlease increase values for a successful tiling."

          );

        }

      }

      while (!done);

 

      PromptKeywordOptions pko =

        new PromptKeywordOptions(

          "\nUse curved or straight segments"

        );

      pko.Keywords.Add("Curved");

      pko.Keywords.Add("Straight");

      pko.Keywords.Default = "Straight";

 

      PromptResult pr = ed.GetKeywords(pko);

      if (pr.Status != PromptStatus.Keyword &&

          pr.Status != PromptStatus.OK)

        return;

 

      bool straight = pr.StringResult == "Straight";

 

      Transaction tr =

        doc.TransactionManager.StartTransaction();

      using (tr)

      {

        LayerTable lt =

          (LayerTable)tr.GetObject(

            db.LayerTableId, OpenMode.ForWrite

          );

 

        // Points and supporting lines will be on layers that

        // are grey and initially frozen

 

        ObjectId ptLay = CreateLayer(tr, lt, "HG_PTS", 9, false);

        ObjectId lnLay = CreateLayer(tr, lt, "HG_LNS", 9, false);

        ObjectId htLay = CreateLayer(tr, lt, "HG_CNT", 6, true);

 

        Circle cir = null;

 

        if (per.Status == PromptStatus.OK)

        {

          cir =

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

        }

 

        // Create our hyperbolic engine

 

        HyperbolicEngine he =

          new HyperbolicEngine(

            tr, db, cir, straight, ptLay, lnLay, htLay

          );

 

        PolygonEngine pe = new PolygonEngine();

 

        if (cir == null)

        {

          he.DrawUnitCircle();

        }

 

        // We'll default to level 6, but this will get reduced

        // if resulting in > 1000 polygons

 

        int level = 6;

        Polygon[] polys = pe.GetPolygons(n, k, ref level, 1000);

 

        if (HostApplicationServices.Current.UserBreak())

          return;

 

        ed.WriteMessage(

          "\n{{{0} {1}}} (level {2}) => {3} polygons.",

          n, k, level, polys.Length

        );

 

        foreach (Polygon p in polys)

        {

          he.DrawPolygon(p.Vertices.ToArray());

          p.Vertices.Clear();

          p.Vertices.Capacity = 0;

        }

 

        tr.Commit();

      }

    }

 

    [CommandMethod("HGT")]

    public static void CreateHyperbolicTriangles()

    {

      Document doc =

        Application.DocumentManager.MdiActiveDocument;

      Database db = doc.Database;

 

      Transaction tr =

        doc.TransactionManager.StartTransaction();     

      using (tr)

      {

        LayerTable lt =

          (LayerTable)tr.GetObject(

            db.LayerTableId, OpenMode.ForWrite

          );

 

        // Points and supporting lines will be on layers that

        // are grey and initially frozen

 

        ObjectId ptLay = CreateLayer(tr, lt, "HG_PTS", 9, false);

        ObjectId lnLay = CreateLayer(tr, lt, "HG_LNS", 9, false);

        ObjectId htLay = CreateLayer(tr, lt, "HG_CNT", 6, true);

 

        // Create our hyperbolic engine

 

        HyperbolicEngine he =

          new HyperbolicEngine(

            tr, db, null, false, ptLay, lnLay, htLay

          );

 

        he.DrawUnitCircle();

 

        he.DrawTriangle(

          new Point2d(0.8, -0.3),

          new Point2d(0.2, -0.6),

          new Point2d(-0.61, -0.6),

          true, true

        );

 

        tr.Commit();

      }

    }

 

    [CommandMethod("HGS")]

    public static void CreateHyperbolicSpiro()

    {

      Document doc =

        Application.DocumentManager.MdiActiveDocument;

      Database db = doc.Database;

 

      Transaction tr =

        doc.TransactionManager.StartTransaction();

      using (tr)

      {

        LayerTable lt =

          (LayerTable)tr.GetObject(

            db.LayerTableId, OpenMode.ForWrite

          );

 

        // Points and supporting lines will be on layers that

        // are grey and initially frozen

 

        ObjectId ptLay = CreateLayer(tr, lt, "HG_PTS", 9, false);

        ObjectId lnLay = CreateLayer(tr, lt, "HG_LNS", 9, false);

        ObjectId htLay = CreateLayer(tr, lt, "HG_CNT", 6, true);

 

        // Create our hyperbolic engine

 

        HyperbolicEngine he =

          new HyperbolicEngine(

            tr, db, null, false, ptLay, lnLay, htLay

          );

 

        he.DrawUnitCircle();

 

        int max = 13;

        for (int k = 0; k < (max - 1); k++)

        {

          he.DrawTriangle(

            new Point2d(1, 0),

            new Point2d(

              Math.Cos(k * 134 / (double)max),

              Math.Sin(k * 134 / (double)max)

            ),

            new Point2d(

              Math.Cos((k + 1) * 134 / (double)max),

              Math.Sin((k + 1) * 134 / (double)max)

            ),

            true, true

          );

        }

 

        tr.Commit();

      }

    }

 

    private static ObjectId CreateLayer(

      Transaction tr, LayerTable lt,

      string lname, short col, bool on

    )

    {

      if (lt.Has(lname))

        return lt[lname];

 

      LayerTableRecord ltr = new LayerTableRecord();

      ltr.Color = Color.FromColorIndex(ColorMethod.ByAci, col);

      ltr.Name = lname;

      ltr.IsFrozen = !on;

 

      ObjectId layId = lt.Add(ltr);

      tr.AddNewlyCreatedDBObject(ltr, true);

 

      return layId;

    }

  }

}

When we run the HT command we have a number of options. We can either select a circle to fill, or hit Enter to have it create a standard circle at the origin, and then are asked to enter the number of sides of the polygons, the number meeting at each vertex, and the number of levels. We’re also asked whether to straight or curved segments (and it’s in the case of curved segments that the code we showed in the last post gets exercised).

Let’s see the results of running the HT command with the default option values selected:

Command: HT

Select circle:

Number of sides for polygons <7>:

Number of polygons meeting at each vertex <3>:

Use curved or straight segments [Curved/Straight] <Straight>:

{7 3} (level 5) => 617 polygons.

We can see that our {7 3} tessellation has resulted in 617 polygons being created. The various patterns generated by this algorithm have a single, central polygon. Other algorithms can also generate a group of polygons centered on the circle – this shouldn’t be overly difficult to implement, but has been left as an exercise for the engaged reader. :-)

Just to see, we can try the same operation with curved segments instead of straight, which will exercise the code in the last post but not actually make a very noticeable difference to the results (the resultant arcs mostly have very low curvature).

I’ll follow up with one more post in this series, to script the HT command to create a varied range of the patterns that can be created with the current implementation.

Also, if people are interested in thinking about how this might be taken to the next dimension (i.e. 3D), then this document should be interesting: it even looks at the possibilities around printing 3D hyperbolic tessellations, although it tends to do so more from an artistic – rather than structural – perspective.

An interesting challenge: generating variable density fill patterns for 3D printing

A regular follower of this blog and someone I now consider a friend through our online interactions, Alex Fielder, recently laid down the gauntlet for the topic of the coming series of posts.

He started with Twitter…

… and then moved on to The Swamp (which he thankfully also brought to my attention using Twitter).

Alex has become interested in an area that I can see becoming increasingly relevant: given the fact that 3D printing is an additive, rather than subtractive or mold-based, process, shouldn’t we take advantage of this fact to optimise the composition of an object’s internal structure? This optimisation could bring benefits in terms of cost and weight (as we can use less material) or strength (as we can add structure that increases tensile strength). And it seems he’s not alone, in thinking along these lines.

The more I think about it, myself – and I admit this is a long way from being my field of specialization – this feels like it has the potential to be a real game-changer. And when you can hook up simulation and analysis tools to provide information on optimal fill patterns, things become really interesting.

But anyway – first things first. Clearly the ideal would be to tackle this kind of challenge in 3D (I don’t know of any 3D printers that take third angle projection input ;-), but you have to start somewhere.

Alex’s post pointed us, quite specifically, at hyperbolic tessellation as a technique for generating fill patterns (see this article for an idea of the possibilities):

Alex had suggested my Spiro app (which is also posted as a Plugin of the Month on Autodesk Labs), as being a potential starting point for this.

After looking into it, it seemed to make more sense to attack this a bit differently than my work on Spiro (which used an approach of collecting a set of points to generate a single polyline).

Also, getting my head around the underlying mathematics was going to be a challenge, so I knew I was going to need to borrow some code from somewhere for this, anyway. Alex also suggested that my series of posts on turtle graphics might also be a good starting point – and this also occurred to me – although a lot would depend on the ultimate algorithm I ended up finding (and whether this kind of tiling could be generated in the hyperbolic plane using L-systems).

To get started on this, I decided to look into two main areas, each of which will be addressed by posts next week:

1. Creating hyperbolic geometry inside AutoCAD
2. Creating hyperbolic tessellations inside AutoCAD

There seem to be a few approaches for defining hyperbolic geometry: the most popular appears to be the Poincaré disk model – which connects points using circular arcs – and the Klein disk model – which connects points using straight lines (at least that appears to be the difference from my – probably ill-informed – viewpoint… I have no doubt there are more subtle differences that I’m missing).

We don’t care very much about the underlying mathematics, at the end of the day, as this is really about the generated pattern. In the next post, we’ll look at getting a straightforward Poincaré disk (something akin to this one) implemented inside AutoCAD, and then (most likely) disregard it in the following post, where we generate hyperbolic tessellations using straight-line segments to connect our points. At least that’s my expectation. :-)

Incidentally, for those interested in the domain of 3D printing, there’s a Singularity University live webcast coming up on January 24, 2012 that you’re likely to find interesting (it also has Carl Bass, Autodesk’s CEO, as guest speaker):

Advances in fields such as artificial intelligence, robotics, and digital manufacturing are undoubtedly going to revolutionize manufacturing during this decade, enabling us to design and “print” complex products and “manufacture” these in our own homes. Exponentially advancing technologies will provide major new opportunities for entrepreneurs to create world-changing technologies. But they also may threaten industries and jobs across the world.

What are these technologies, where exactly do the opportunities lie, and where will the jobs go? These are some of the topics that Autodesk CEO Carl Bass and Singularity University VP of Academics and Innovation Vivek Wadhwa will discuss in the next Which Way Next?

Which Way Next?

A Singularity University Live Webcast
Guest: Carl Bass, CEO of Autodesk Inc.
Commentator: Vivek Wadhwa, Singularity University
Topic: AI, Robotics and Digital Manufacturing
Date: Tuesday, January 24, 2012
Time: LIVE at 12:30 PM (Pacific)
Log-in to participate & submit questions: www.singularityu.org