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.
