Topic: Creating an Apollonian gasket using F#

[kaefer]

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}

[vegbruiser]

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]

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]

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]

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]

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