Author Topic: Creating an Apollonian gasket using F#  (Read 2380 times)

0 Members and 1 Guest are viewing this topic.

kaefer

  • Guest
Creating an Apollonian gasket using F#
« on: February 07, 2012, 06:37:17 PM »
For those interested in the unfolding story at Kean's blog and here, I'd like to humbly submit my take on this in F#, with a few minor tweaks.

The progression is 2*3^(n-1)+2 circles for n steps.

Doesn't need .Net4.0 for Complex, but a reference to FSharp.PowerPack.dll.

Code - F#: [Select]
  1. // F# adaptation of Kean Walmsley's Apollonian gasket code - 2012-02-06
  2. // http://through-the-interface.typepad.com/through_the_interface/2012/02/circle-packing-in-autocad-creating-an-apollonian-gasket-using-net.html
  3. // with Gaston Nunez's comment (#1) taken into account, see also
  4. // http://en.wikibooks.org/wiki/Fractals/Apollonian_fractals#Soddy.27s_algorithm
  5.  
  6. module CirclePacking
  7.  
  8. open Autodesk.AutoCAD.EditorInput
  9. open Autodesk.AutoCAD.DatabaseServices
  10. open Autodesk.AutoCAD.Geometry
  11. open Autodesk.AutoCAD.Runtime
  12.  
  13. type acadApp = Autodesk.AutoCAD.ApplicationServices.Application
  14.  
  15. type MyCircle = {
  16.     Center : Point2d
  17.     Radius : float
  18.     ColorIndex : int } with
  19.     static member Create(x, y, radius) = {
  20.         Center = new Point2d(x, y)
  21.         Radius = radius
  22.         ColorIndex = 5 }
  23.  
  24. type TangentType = OutsideAll | OutsideTwo | Outer
  25.  
  26. let findApollonianCircle c1 c2 c3 tanType =
  27.  
  28.     let (k1, k2, k3) =
  29.         match tanType with
  30.         | OutsideAll ->  1. / c1.Radius,  1. / c2.Radius,  1. / c3.Radius
  31.         | OutsideTwo ->  1. / c1.Radius,  1. / c2.Radius, -1. / c3.Radius
  32.         | Outer ->      -1. / c1.Radius, -1. / c2.Radius, -1. / c3.Radius
  33.     let k4 = k1 + k2 + k3 + 2. * sqrt(k1 * k2 + k1 * k3 + k2 * k3)
  34.  
  35.     let z1 = Complex.Create(c1.Center.X, c1.Center.Y)
  36.     let z2 = Complex.Create(c2.Center.X, c2.Center.Y)
  37.     let z3 = Complex.Create(c3.Center.X, c3.Center.Y)
  38.     let z4 = k1 * z1 + k2 * z2 + k3 * z3
  39.     let z4 = z4 + 2. * sqrt(k1 * k2 * z1 * z2 + k1 * k3 * z1 * z3 + k2 * k3 * z2 * z3)
  40.     let z4 = z4 / Complex.Create(k4, 0.)
  41.  
  42.     MyCircle.Create(z4.r, z4.i, 1. / k4)
  43.  
  44. // Draw a circle tangent to these three circles and that is
  45. // outside all three
  46. let rec findCircleOutsideAll level cir0 cir1 cir2 = seq{
  47.    
  48.     let newCir = findApollonianCircle cir0 cir1 cir2 OutsideAll
  49.     yield newCir
  50.     if level > 1 then
  51.         yield! findCircleOutsideAll (level - 1) cir0 cir1 newCir
  52.         yield! findCircleOutsideAll (level - 1) cir0 cir2 newCir
  53.         yield! findCircleOutsideAll (level - 1) cir1 cir2 newCir }
  54.  
  55. // Draw a circle tangent to these three circles and that is
  56. // outside two of them
  57. let rec findCircleOutsideTwo level cir0 cir1 contains = seq{
  58.    
  59.     let newCir = findApollonianCircle cir0 cir1 contains OutsideTwo
  60.     yield newCir
  61.     if level > 1 then
  62.         yield! findCircleOutsideTwo (level - 1) newCir cir0 contains
  63.         yield! findCircleOutsideTwo (level - 1) newCir cir1 contains
  64.         yield! findCircleOutsideAll (level - 1) cir0 cir1 newCir }
  65.  
  66. // Find the Apollonian packing
  67. let findApollonianPacking level =
  68.    
  69.     // Create the three central tangent circles
  70.     let radius = 1. / (2. + 4. / sqrt 3.)
  71.    
  72.     // The centroid of the central circles lies at point (1, 1)
  73.     let (x0, y0) = 1.0, 0.5 + radius
  74.     let (x1, y1) = x0 - radius, y0 + radius * sqrt 3.
  75.     let (x2, y2) = x0 + radius, y1
  76.  
  77.     let cir0 = MyCircle.Create(x0, y0, radius)
  78.     let cir1 = MyCircle.Create(x1, y1, radius)
  79.     let cir2 = MyCircle.Create(x2, y2, radius)
  80.  
  81.     // Find the circle that contains them all, we know already
  82.     // its center (1, 1) and diameter (1)
  83.     let outer = findApollonianCircle cir0 cir1 cir2 Outer
  84.  
  85.     seq{
  86.         yield outer
  87.         yield cir0
  88.         yield cir1
  89.         yield cir2
  90.         if level > 1 then
  91.             // Find the central circle
  92.             yield! findCircleOutsideAll (level - 1) cir0 cir1 cir2
  93.  
  94.             // Find circles tangent to the big circle
  95.             yield! findCircleOutsideTwo (level - 1) cir0 cir1 outer
  96.             yield! findCircleOutsideTwo (level - 1) cir1 cir2 outer
  97.             yield! findCircleOutsideTwo (level - 1) cir2 cir0 outer }
  98.  
  99. let drawCircles cirId level =
  100.     let db = HostApplicationServices.WorkingDatabase
  101.     use tr = db.TransactionManager.StartTransaction()
  102.     let btr =
  103.         tr.GetObject(
  104.             db.CurrentSpaceId, OpenMode.ForWrite )
  105.                 :?> BlockTableRecord
  106.  
  107.     let cir =
  108.         if cirId = ObjectId.Null then
  109.             // Create a new circlee at the default location
  110.             let cir =
  111.                 new Circle(Point3d.Origin, Vector3d.ZAxis, 500.)
  112.             btr.AppendEntity cir |> ignore
  113.             tr.AddNewlyCreatedDBObject(cir, true)
  114.             cir
  115.         else
  116.             // Otherwise use the selected circle
  117.             tr.GetObject(cirId, OpenMode.ForRead) :?> Circle
  118.  
  119.     let mat =
  120.         Matrix3d.Displacement(cir.Center.GetAsVector()) *
  121.         Matrix3d.Scaling(cir.Diameter, Point3d.Origin) *
  122.         Matrix3d.WorldToPlane(cir.Normal).Inverse() *
  123.         Matrix3d.Displacement(new Vector3d(-1., -1., 0.))
  124.  
  125.     // Add each of the returned circles to the current space
  126.     let cnt =
  127.         findApollonianPacking level
  128.         |> Seq.fold(fun n myC ->
  129.             let c =
  130.                 new Circle(
  131.                     new Point3d(myC.Center.X, myC.Center.Y, 0.),
  132.                     Vector3d.ZAxis, myC.Radius,
  133.                     ColorIndex = myC.ColorIndex )
  134.             c.TransformBy mat
  135.             btr.AppendEntity c |> ignore
  136.             tr.AddNewlyCreatedDBObject(c, true)
  137.             n + 1 ) 0
  138.     tr.Commit()
  139.     cnt
  140.  
  141. [<CommandMethod "myAG">]
  142. let apollonianGasket() =
  143.     let ed = acadApp.DocumentManager.MdiActiveDocument.Editor
  144.     let peo =
  145.         new PromptEntityOptions(
  146.             "\nSelect circle: ",
  147.             AllowNone = true )
  148.     peo.SetRejectMessage "Must be a circle."
  149.     peo.AddAllowedClass(typeof<Circle>, false)
  150.     let per = ed.GetEntity peo
  151.  
  152.     if  per.Status = PromptStatus.OK ||
  153.         per.Status = PromptStatus.None then
  154.  
  155.         let pio =
  156.             new PromptIntegerOptions(
  157.                 "\nEnter maximum level: ",
  158.                 AllowNone = true,
  159.                 LowerLimit = 1,  // 2*3^0+2 = 4 circles
  160.                 UpperLimit = 13, // 2*3^12+2 = 1062884 circles
  161.                 DefaultValue = 8 )
  162.         let pir = ed.GetInteger pio
  163.  
  164.         if  pir.Status = PromptStatus.OK ||
  165.             pir.Status = PromptStatus.None then
  166.  
  167.             drawCircles
  168.                 (   if per.Status = PromptStatus.OK then per.ObjectId
  169.                     else ObjectId.Null )
  170.                 (   if pir.Status = PromptStatus.OK then pir.Value
  171.                     else pio.DefaultValue )
  172.             |> fun cnt -> ed.WriteMessage("\n{0} Circles. ", cnt)

Edit: Reference to PowerPack
« Last Edit: February 07, 2012, 06:42:28 PM by kaefer »

vegbruiser

  • Guest
Re: Creating an Apollonian gasket using F#
« Reply #1 on: February 08, 2012, 06:39:10 AM »
That's some good work Kaefer! Do you have any idea how much of a speed increase it achieves (if any) over and above the C# version Kean wrote?

kaefer

  • Guest
Re: Creating an Apollonian gasket using F#
« Reply #2 on: February 08, 2012, 09:26:47 AM »
Do you have any idea how much of a speed increase it achieves (if any) over and above the C# version Kean wrote?

Basically, f stands for functional, not for fast. Even if most of the algorithmic complexity is now hidden away in the Complex module, those floating point operations have still to be performed. Compared to mutable collections, using sequence expressions to generate the objects recursively incurs a performance hit too; albeit one you also could have in C# (with yield return).

So I wouldn't expect any increase in speed, only in readability.


kaefer

  • Guest
Re: Creating an Apollonian gasket using F#
« Reply #3 on: February 10, 2012, 03:23:21 AM »
Now that the head of the cat's out of the bag, over at Through The Interface, many thanks to Kean, I came to realize that even with Descartes there's still something to abstract away. Here's his fourth circle problem solver with the formula inlined, which will work on floats and complex numbers regardless.

Code - F#: [Select]
  1. let inline solve a b c =
  2.     let (z, d) = a + b + c, sqrt(a * b + b * c + a * c)
  3.     z + d + d // Discard the second solution z - d - d
  4.  
  5. // Custom operator to multiply a real
  6. let ( *!) a b = Complex.Create(a, 0.) * b
  7.  
  8. let fourthCircle (x1, y1, k1) (x2, y2, k2) (x3, y3, k3) =
  9.     let pos1 = Complex.Create(x1, y1)
  10.     let pos2 = Complex.Create(x2, y2)
  11.     let pos3 = Complex.Create(x3, y3)
  12.     let (kz1, kz2, kz3) = k1 *! pos1, k2 *! pos2, k3 *! pos3
  13.     let k4 = solve k1 k2 k3
  14.     let pos4 = 1. / k4 *! solve kz1 kz2 kz3
  15.     pos4.r, pos4.i, k4

The image shows the consequences of straying too far from the first quadrant to the origin.

Kean

  • Newt
  • Posts: 48
Re: Creating an Apollonian gasket using F#
« Reply #4 on: February 14, 2012, 11:04:17 AM »
Interesting - thanks, kaefer!

I'm glad you found a way to factor away the two sets of functions - that was really annoying me.

Kean

kaefer

  • Guest
Re: Creating an Apollonian gasket using F#
« Reply #5 on: February 14, 2012, 03:44:19 PM »
While we're at it, we can also eliminate the awkward multiplication of real and complex. Sorry, I didn't see that before.

Code - F#: [Select]
  1. let inline solve a b c =
  2.     let (z, d) = a + b + c, sqrt(a * b + b * c + a * c)
  3.     z + d + d // Discard the second solution z - d - d
  4.  
  5. let fourthCircle (x1, y1, k1) (x2, y2, k2) (x3, y3, k3) =
  6.     let kz1 = Complex.Create(x1 * k1, y1 * k1)
  7.     let kz2 = Complex.Create(x2 * k2, y2 * k2)
  8.     let kz3 = Complex.Create(x3 * k3, y3 * k3)
  9.     let (k4, kz4)  = solve k1 k2 k3, solve kz1 kz2 kz3
  10.     kz4.r / k4, kz4.i / k4, k4