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.