February 2012

February 29, 2012

Another really interesting, developer-oriented feature in AutoCAD 2013 is something we’ve been calling “Dynamic .NET”. I don’t know whether that’s official branding, or not – I suspect not – but it’s the moniker we’ve been using to describe this capability internally and to ADN members.

The capability is based on an addition to .NET in version 4.0: to complement (or perhaps just as part of) the integration of the Dynamic Language Runtime into the core .NET Framework, various interfaces – including IDynamicMetaObjectProvider – were provided to developers to let their objects participate in “dynamic” operations.

Providers of APIs therefore now have the option to make certain classes dynamic – for instance, having them implement the IDynamicMetaObjectProvider interface – which means they can effectively be scripted, without binding properties and methods at compile-time. This implements a concept known as duck typing, and is at the core of many scripting language implementations. You may also hear it referred to as late binding: take a look at this excellent post by Eric Lippert if you’d like to get a better understanding of what that term means.

So what have we actually done with AutoCAD 2013? We have made Autodesk.AutoCAD.DatabaseServices.ObjectId support dynamic operations: you can declare an ObjectId using the dynamic keyword in C#, and you can then choose to access any of the properties or methods exposed by the object possessing that ID directly from the ObjectId itself. The ObjectId then takes care of the open & close operations implicitly, so there’s no need to mess around with transactions, etc.

If you’re interested in more background to this, I suggest taking a look at Albert Szilvasy’s class from AU 2009, where he gave us a hint of what was to come (although it was far from being a certainty, at that time – we have since found it to be a popular option via our annual API wishlist surveys).

OK, so what does this mean, in practice? Here’s some statically-typed C# code that iterates through the layer table and prefixes all layers (other than layer “0”) with a string:

publicstaticvoid ChangeLayerNames()

{

Database db = HostApplicationServices.WorkingDatabase;

using (Transaction tr = db.TransactionManager.StartTransaction())

{

LayerTable lt =

(LayerTable)tr.GetObject(

db.LayerTableId, OpenMode.ForRead

);

foreach (ObjectId ltrId in lt)

{

LayerTableRecord ltr =

(LayerTableRecord)tr.GetObject(ltrId, OpenMode.ForRead);

if (ltr.Name != "0")

{

ltr.UpgradeOpen();

ltr.Name = "First Floor " + ltr.Name;

}

}

tr.Commit();

}

}

Here’s how the code could look if using dynamic types. I’ve done my best not to exaggerate the differences in numbers of lines by using different conventions in the two implementations, but I did have to add more line breaks in the above code to have it fit the width of this blog.

publicstaticvoid ChangeLayerNamesDynamically()

{

dynamic layers =

HostApplicationServices.WorkingDatabase.LayerTableId;

foreach (dynamic l in layers)

if (l.Name != "0")

l.Name = "First Floor " + l.Name;

}

Something that needs to be stressed: as the properties and methods are effectively being determined and then accessed/called at run-time, there is no IntelliSense available when coding dynamically (something Python and Ruby developers are used to, but we statically-typed .NET developers are probably less so). It would not be impossible to implement, but it would be very difficult: aside from types being inferred at the time code is being written (I’ll use the term develop-time, although I haven’t seen it used before), there would need to be some work done to effectively populate the list of properties & methods, which would typically mean opening the object and reflecting over its protocol… not at all straightforward to implement, as far as I can tell, especially when you don’t actually have AutoCAD running to help with accessing the objects.

A lack of IntelliSense does mean this style of coding it likely to be less accessible to novice coders, unless a more integrated REPL environment becomes available (as this is something that tends to balance out the lack of other develop-time tooling when scripting – just as people have found with the ability to execute fragments of AutoLISP at AutoCAD’s command-line, over the years).

But that’s all for the future, and is very much my own speculation, at this stage.

One question that has come up in relation to “Dynamic .NET” in AutoCAD is around performance: as we’re effectively opening and closing an object every time we execute a method (and we may actually open and close an object twice – once for read, once for write – when we read one of its properties, modify it and store it back), isn’t there a performance impact? I admit I haven’t performed extensive benchmarking of this, and there is likely to be some overhead if dealing with lots of properties and methods. But in terms of getting something working quickly – which you can then tune for performance later, as needed – this is likely to be a valuable tool for developers, irrespective of the performance impact (which I’d hope to be modest).

And it really comes into its own when used in combination with LINQ. Here’s some code that Stephen Preston put together to test out some of the possibilities with LINQ (I’ve included the CLN and CLN2 commands – shown above – for completeness).

using Autodesk.AutoCAD.Runtime;

using Autodesk.AutoCAD.ApplicationServices;

using Autodesk.AutoCAD.DatabaseServices;

using Autodesk.AutoCAD.Geometry;

using Autodesk.AutoCAD.EditorInput;

using System.Collections.Generic;

using System.Linq;

namespace DynamicTyping

{

publicclassCommands

{

[CommandMethod("CLN")]

publicstaticvoid ChangeLayerNames()

{

Database db = HostApplicationServices.WorkingDatabase;

using (

Transaction tr = db.TransactionManager.StartTransaction()

)

{

LayerTable lt =

(LayerTable)tr.GetObject(

db.LayerTableId, OpenMode.ForRead

);

foreach (ObjectId ltrId in lt)

{

LayerTableRecord ltr =

(LayerTableRecord)tr.GetObject(ltrId, OpenMode.ForRead);

if (ltr.Name != "0")

{

ltr.UpgradeOpen();

ltr.Name = "First Floor " + ltr.Name;

}

}

tr.Commit();

}

}

[CommandMethod("CLN2")]

publicstaticvoid ChangeLayerNamesDynamically()

{

dynamic layers =

HostApplicationServices.WorkingDatabase.LayerTableId;

foreach (dynamic l in layers)

if (l.Name != "0")

l.Name = "First Floor " + l.Name;

}

// Adds a custom dictionary in the extension dictionary of

// selected objects. Uses dynamic capabilities, but not LINQ.

[CommandMethod("AddMyDict")]

publicvoid AddMyDict()

{

Document doc = Application.DocumentManager.MdiActiveDocument;

Database db = doc.Database;

Editor ed = doc.Editor;

PromptSelectionResult res = ed.GetSelection();

if (res.Status != PromptStatus.OK)

return;

foreach (dynamic ent in res.Value.GetObjectIds())

{

dynamic dictId = ent.ExtensionDictionary;

if (dictId == ObjectId.Null)

{

ent.CreateExtensionDictionary();

dictId = ent.ExtensionDictionary();

}

if (!dictId.Contains("MyDict"))

dictId.SetAt("MyDict", newDBDictionary());

}

}

// Adds a custom dictionary in the extension dictionary of

// selected objects, but this time using dynamic capabilities

// with LINQ. Conceptually simpler, but with some performance

// overhead, as we're using two queries: one to get entities

// without extension dictionaries (and then add them) and the

// other to get entities with extension dictionaries.

[CommandMethod("AddMyDict2")]

publicvoid AddMyDict2()

{

PromptSelectionResult res =

Application.DocumentManager.MdiActiveDocument.Editor.

GetSelection();

if (res.Status != PromptStatus.OK)

return;

// Query for ents in selset without ExtensionDictionaries

var noExtDicts =

from ent in res.Value.GetObjectIds().Cast<dynamic>()

where ent.ExtensionDictionary == ObjectId.Null

select ent;

// Add extension dictionaries

foreach (dynamic ent in noExtDicts)

ent.CreateExtensionDictionary();

// Now we've added the ext dicts, we add our dict to each

var noMyDicts =

from ent in res.Value.GetObjectIds().Cast<dynamic>()

where !ent.ExtensionDictionary.Contains("MyDict")

select ent.ExtensionDictionary;

foreach (dynamic dict in noMyDicts)

dict.SetAt("MyDict", newDBDictionary());

}

// Access various bits of information using LINQ

[CommandMethod("IUL")]

publicvoid InfoUsingLINQ()

{

Document doc = Application.DocumentManager.MdiActiveDocument;

Database db = doc.Database;

Editor ed = doc.Editor;

dynamic bt = db.BlockTableId;

// Dynamic .NET loop iteration

ed.WriteMessage("\n*** BlockTableRecords in this DWG ***");

foreach (dynamic btr in bt)

ed.WriteMessage("\n" + btr.Name);

// LINQ query - returns startpoints of all lines

ed.WriteMessage(

"\n\n*** StartPoints of Lines in ModelSpace ***"

);

dynamic ms = SymbolUtilityServices.GetBlockModelSpaceId(db);

var lineStartPoints =

from ent in (IEnumerable<dynamic>)ms

where ent.IsKindOf(typeof(Line))

select ent.StartPoint;

foreach (Point3d start in lineStartPoints)

ed.WriteMessage("\n" + start.ToString());

// LINQ query - all entities on layer '0'

ed.WriteMessage("\n\n*** Entities on Layer 0 ***");

var entsOnLayer0 =

from ent in (IEnumerable<dynamic>)ms

where ent.Layer == "0"

select ent;

foreach (dynamic e in entsOnLayer0)

ed.WriteMessage(

"\nHandle=" + e.Handle.ToString() + ", ObjectId=" +

((ObjectId)e).ToString() + ", Class=" + e.ToString()

);

ed.WriteMessage("\n\n");

// Using LINQ with selection sets

PromptSelectionResult res = ed.GetSelection();

if (res.Status != PromptStatus.OK)

return;

// Select all entities in selection set that have an object

// called "MyDict" in their extension dictionary

var extDicts =

from ent in res.Value.GetObjectIds().Cast<dynamic>()

where ent.ExtensionDictionary != ObjectId.Null &&

ent.ExtensionDictionary.Contains("MyDict")

select ent.ExtensionDictionary.Item("MyDict");

// Erase our dictionary

foreach (dynamic myDict in extDicts)

myDict.Erase();

}

}

}

You’ll see that certain operations become very straightforward when combining Dynamic .NET and LINQ. This could well be the “sweet spot” of the dynamic capability – it certainly makes LINQ a more natural choice when accessing various types of data in AutoCAD drawings.

That’s it for today’s post. Next time we’ll summarize the remaining API enhancements coming in AutoCAD 2013, before looking at the migration steps needed for .NET applications.

February 27, 2012

In the first of this week’s posts on the new, developer-oriented features in AutoCAD 2013, we’re going to take a look at the AutoCAD Core Console. I predict that this one feature alone is going to be worth its weight in gold to many customers and developers.

Picture this: a stripped down version of AutoCAD – but only from a UI perspective – that can be executed from a standard Command Prompt. We’re talking about a version with almost all the capabilities of full AutoCAD – with the ability to load certain applications and execute scripts – but launches in a couple of seconds (if you have a slow machine :-) and is absolutely perfect for batch operations.

[Note: you can only run this executable on systems that have AutoCAD 2013 (or one of its verticals) installed – you shouldn’t expect to be able to run it elsewhere.]

And as well as being good at performing a set of sequential operations, the Core Console is a great way to get process-level parallelism (that harnesses all those spare cores many of us now have): you can launch multiple instances of the executable to run concurrently (or – for that matter – simply launch one to implement a “background” version of a command that would otherwise hog the editor during a time-consuming operation).

Before we get into what it means to be able to load “certain” applications, let’s talk a little about the background to this tool. It was made possible by a multi-year architectural project known internally as The Big Split: to make AutoCAD truly – and sustainably – cross-platform, we needed a core component that could be compiled for multiple (for now that means two) platforms. The majority of AutoCAD’s capabilities – including the drawing canvas and graphics sub-system – would be contained “inside” (if you think about it in logical – rather than physical – terms) this component.

This component is AcCore.dll on the Windows platform. If you take a look at acad.exe for AutoCAD 2013, for instance, you’ll see it is now about 5.7 MB, with AcCore.dll weighing in at 13.4 MB. Compare this with AutoCAD 2012, where we had a monolithic acad.exe of 17.5 MB.

Subsetting out this core functionality means we can share a common codebase for this OS-independent set of functionality, and expose an appropriate UI to it for whichever platform we’re targetting. For AutoCAD for Windows, this means using some MFC, some WPF, etc., while for AutoCAD for Mac this means using Cocoa.

And for the simplest possible, command-line only UI, we have the 24 KB AcCoreConsole.exe.

So what kinds of application can be loaded into the Core Console? It can certainly load DBX modules – which are coded against ObjectDBX/RealDWG, which is a subset of the “core” capabilities – but it can also load CRX modules and .NET DLLs that have been coded against the new AcCoreMgd.dll (rather than AcMgd.dll).

You can think of CRX modules as ARX modules developed to run against the core functionality: many of your commands – especially if using the command-line and jigs for their user-input – can be ported to CRX modules (or their .NET equivalent). Ideally you’d then have the GUI integration alone in your ARX modules (or .NET DLLs using AcMgd.dll) which then call into the “core” implementations.

This would allow them to be loaded into the Core Console and even in other environments, moving forwards (watch this space for more on that :-).

You can also load AutoLISP files into the Core Console – just as you can execute scripts – you just need to be aware that certain (mostly GUI-oriented) capabilities will not be available to you.

Let’s now take a look at a very concrete use for this tool, along with some code to drive it.

A very common requirement is for some kind of server-resident process publishing DWFs or PDFs from DWGs. Balaji Ramamoorthy, from DevTech India, wrote an elegant batch script that can sit on a server and use the Core Console to generate PDFs from DWGs (actually it’s a couple of batch files and an AutoCAD script, but anyway).

The main batch file sits and watches a particular “in” folder (which would ideally be network-accessible to be of much use). If it’s empty, the batch file waits for 10 seconds before checking again (to avoid unnecessary thrashing). If there’s something in the folder, the main batch file calls a secondary batch file (GenPDF.bat) to process the file(s) and create the results in the “out” folder.

Until now, the focus on the series has largely been around relatively pure mathematical approaches to fill simple (in our case circular or spherical) areas/volumes: we’ve so far avoided having to adapt the filling/packing algorithm depending on the form selected.

Today’s post changes that: it uses a very different approach to fill a selected Solid3d with spheres. A bit like the “hyperbolic tessellation” patterns we saw, early on, we want to create smaller spheres at the surface of the solid, with the radius gradually increasing, layer upon layer, as we get towards the centre.

Something like this sketch, created by my colleague Francesco Tonioni (as mentioned in the previous post) as we brainstormed optimal filling strategies (from a strength – rather than material – perspective):

To do this for an arbitrary solid, we’re using this basic approach:

Section a solid along an appropriate plane (one that’s defined by two user-selected points and the view direction) and create a ring of (initially very small) spheres along the resultant curve.

Move the points along in one direction (normal to the plane) by the diameter of the sphere. Repeat step 1 (and then 2) until you run out of solid to section.

Starting from the initially-created section, move the points in the other direction repeating the process until you hit the other end of the solid.

Now that you’ve got a complete layer of spheres at the surface of the solid, offset the solid by the diameter of the spheres and start again with a slightly larger sphere diameter.

There are a host of devils in the detail, of course: to reduce the chances of any of the outer layer of spheres intersecting the surface of the solid, we need to know the curvature of the surface at a particular point (and while we get that right for relatively uniform surfaces, I’ll bet it’ll fail for highly irregular ones), which allows us to adjust the depth of the sphere at that point to make sure it fits.

In fact, here are the back-of-the-envelope sketches I drew to work out some of the trigonometry needed to find out how adjacent rings of spheres should be positioned:

These calculations are far from foolproof, so we iterate adjusting the values to see if we can get a better fit. Once again, this probably works well with solids with a uniform shape along the direction of our slices (the normal to the section vector) – so rotations, etc., can work very well, but I’d expect less uniform shapes to cause problems.

Here’s the C# code, also contained in the project with the various code files for this series. I ended up using C# for this post, as it had a lot to do with accessing and managing AutoCAD geometry: the times I’ve used F# of late have been when I’ve had more “pure” calculations to perform. There’s no strict need for this division, but I do like having a split between core algorithms and AutoCAD-specific code. And while it might seem a bit artificial to use a language barrier to enforce it, it does feel like it plays to the strengths of the respective languages.

using Autodesk.AutoCAD.ApplicationServices;

using Autodesk.AutoCAD.DatabaseServices;

using Autodesk.AutoCAD.EditorInput;

using Autodesk.AutoCAD.Geometry;

using Autodesk.AutoCAD.Runtime;

using System;

namespace SpherePacking

{

publicclassPackingCommands

{

[CommandMethod("PACKSPHERES")]

staticpublicvoid PackSolidWithSpheres()

{

Document doc =

Application.DocumentManager.MdiActiveDocument;

Database db = doc.Database;

Editor ed = doc.Editor;

// We need a solid

PromptEntityOptions peo =

newPromptEntityOptions("\nSelect solid to fill");

peo.SetRejectMessage("\nMust be a Solid3d.");

peo.AddAllowedClass(typeof(Solid3d), false);

PromptEntityResult per = ed.GetEntity(peo);

if (per.Status != PromptStatus.OK)

return;

ObjectId solId = per.ObjectId;

// And an initial section plane

// (the view direction is used to define the third point)

PromptPointOptions ppo =

newPromptPointOptions(

"\nSelect first point for initial section plane"

);

PromptPointResult ppr = ed.GetPoint(ppo);

if (ppr.Status != PromptStatus.OK)

return;

Point3d first = ppr.Value;

ppo.Message =

"\nSelect second point for initial section plane";

ppo.BasePoint = first;

ppo.UseBasePoint = true;

ppo.UseDashedLine = true;

ppr = ed.GetPoint(ppo);

if (ppr.Status != PromptStatus.OK)

return;

Point3d second = ppr.Value;

// We either generate spheres or subtract them from the

// selected solid

PromptKeywordOptions pko =

newPromptKeywordOptions(

"\nSubtract spheres from original?"

);

pko.AllowNone = true;

pko.Keywords.Add("Yes");

pko.Keywords.Add("No");

pko.Keywords.Default = "No";

PromptResult pkr = ed.GetKeywords(pko);

if (pkr.Status != PromptStatus.OK)

return;

bool subtract = (pkr.StringResult == "Yes");

// We need a transaction

Transaction tr =

db.TransactionManager.StartTransaction();

using (tr)

{

// And we also need the modelspace open for write

BlockTable bt =

(BlockTable)tr.GetObject(

db.BlockTableId,

OpenMode.ForRead,

false

);

BlockTableRecord btr =

(BlockTableRecord)tr.GetObject(

bt[BlockTableRecord.ModelSpace],

OpenMode.ForWrite,

false

);

// Get the view direction

ViewportTableRecord vtr =

(ViewportTableRecord)tr.GetObject(

ed.ActiveViewportId, OpenMode.ForRead

);

Vector3d viewDir = vtr.ViewDirection;

// Open the selected solid for write

Solid3d outer =

tr.GetObject(solId, OpenMode.ForWrite) asSolid3d;

bool cancelled = false;

if (outer != null)

{

try

{

// Create our sphere packer and use it to fill the

// selected solid with spheres

SpherePacker sp =

newSpherePacker(

db, tr, btr, outer, first, second, viewDir, subtract

);

sp.FillWithSpheres();

}

catch

{

cancelled = true;

}

}

if (!cancelled)

tr.Commit();

}

}

}

publicclassSpherePacker

{

Database _db = null;

Transaction _tr = null;

BlockTableRecord _btr = null;

Solid3d _parent = null, _sol = null;

Point3d _first, _second;

Vector3d _viewDir;

bool _sub = false;

ProgressMeter _pm = null;

// Constructor

public SpherePacker(

Database db, Transaction tr, BlockTableRecord btr,

Solid3d parent, Point3d first, Point3d second,

Vector3d viewDir, bool subtract

)

{

// Set the member data from the arguments

_db = db;

_tr = tr;

_btr = btr;

_parent = parent;

_sub = subtract;

_first = first;

_second = second;

_viewDir = viewDir;

// Make a copy of our outer solid, which we will

// repeatedly shrink for each layer of spheres

_sol = newSolid3d();

_sol.CopyFrom(_parent);

// Add it to the drawing

AddToDrawing(_sol);

// Make sure we have target layers in place for each layer

// of spheres

CreateLayers();

}

// Get the name for a layer based on its integer level

// (a very simple function, but nice to have it centralised

// in case we want to do things differently)

privatestring GetLayer(int layer)

{

return layer.ToString();

}

// Create a set of layers for our spheres

internalvoid CreateLayers()

{

// Open the layer table for write

LayerTable lt =

(LayerTable)_tr.GetObject(

_db.LayerTableId, OpenMode.ForWrite

);

// Loop through, creating 10 layers

for (short i = 1; i <= 10; i++)

{

string name = GetLayer(i);

if (!lt.Has(name))

{

LayerTableRecord ltr = newLayerTableRecord();

ltr.Color =

Autodesk.AutoCAD.Colors.Color.FromColorIndex(

Autodesk.AutoCAD.Colors.ColorMethod.ByAci, i

);

ltr.Name = name;

lt.Add(ltr);

_tr.AddNewlyCreatedDBObject(ltr, true);

}

}

}

// Our main method

internalvoid FillWithSpheres()

{

// Create and start a progress meter (the limit is slightly

// arbitrary - it works well for spheres, but is likely to

// work less well for other forms)

_pm = newProgressMeter();

_pm.SetLimit(180);

_pm.Start(

_sub ?

"Subtracting spheres from selected solid" :

"Filling solid with spheres"

);

try

{

// Start with layer 1 and a radius of 0.2 (which should

// get overwritten immediately of the solid has valid

// extents)

int layer = 1;

double radius = 0.2;

if (_sol.Bounds.HasValue)

{

// Set the radius relative to the solid's extents

Extents3d ext = _sol.Bounds.Value;

Vector3d vec = ext.MaxPoint - ext.MinPoint;

double mag = vec.Length;

radius = mag / 150;

}

// Shrink the body by half the radius of the first layer,

// just to give us a little gap between the spheres

// and the outer solid

_sol.OffsetBody(-0.5 * radius);

// Create layers of spheres, starting with the initial

// radius value and multiplying by a factor for each

// new layer

for (double r = radius; r < radius * 20; r *= 1.3)

{

GenerateSectionsAcrossSolid(r, layer++);

// Shrink the body by the diameter of the layer's spheres

try

{

_sol.OffsetBody(-2 * r);

}

catch

{

break;

}

// We won't attempt to fill volumes that aren't large

// enough (we check the radius against the volume)

if (_sol.MassProperties.Volume < 350 * r)

break;

}

// Set our shrinking solid to be the next level down

// (via its layer)

_sol.Layer = GetLayer(layer);

// If subtracting, do so for this remaining solid

if (_sub)

{

_parent.BooleanOperation(

BooleanOperationType.BoolSubtract, _sol

);

}

}

catch

{

// Rethrow any exception, letting is get caught in our

// command and dealt with there

throw;

}

finally

{

// Make sure we stop our progress meter

_pm.Stop();

}

}

// Generate a layer of spheres across the complete solid

internalvoid GenerateSectionsAcrossSolid(

double radius, int layer

)

{

// Use a coordinate system object to work out the direction

// vector for us to move the section plane

// (should be perpendicular to the section plane, basically)

CoordinateSystem3d cs =

newCoordinateSystem3d(

Point3d.Origin, _viewDir, _second - _first

);

Vector3d unitVec = cs.Zaxis / cs.Zaxis.Length;

// Local variables that get adjusted by each call to our

// GenerateSpheresAlongSection method

Point3d pt1 = _first,

pt2 = _second,

prevCen = Point3d.Origin;

Curve prevCur = null;

double prevDep = 0.0;

// Create a ring of spheres around our initial section curve

GenerateSpheresAlongSection(

ref pt1, ref pt2, ref prevCen, ref prevCur, ref prevDep,

unitVec, radius, 0, layer

);

if (prevCur != null)

prevCur.Dispose();

prevCur = null;

bool done = false;

Point3d initCen = prevCen;

double initDep = prevDep;

// Now go all the way across in one direction

do

{

done =

!GenerateSpheresAlongSection(

ref pt1, ref pt2, ref prevCen, ref prevCur, ref prevDep,

unitVec, radius, 1, layer

);

}

while (!done);

if (prevCur != null)

prevCur.Dispose();

prevCur = null;

// And then go all the way in the other direction

pt1 = _first;

pt2 = _second;

prevCen = initCen;

prevDep = initDep;

do

{

done =

!GenerateSpheresAlongSection(

ref pt1, ref pt2, ref prevCen, ref prevCur, ref prevDep,

unitVec, radius, -1, layer

);

}

while (!done);

if (prevCur != null)

prevCur.Dispose();

}

// Generate a ring of spheres along a single section

privatebool GenerateSpheresAlongSection(

refPoint3d pt1, refPoint3d pt2, refPoint3d prevCen,

refCurve prevCur, refdouble prevDep,

Vector3d unitVec, double radius,

int direction, int layer

)

{

// Tick our progress meter and check for user break

_pm.MeterProgress();

System.Windows.Forms.Application.DoEvents();

if (HostApplicationServices.Current.UserBreak())

{

if (prevCur != null)

prevCur.Dispose();

thrownew System.Exception("User break");

}

// Take local copies of the points defining the section

Point3d first = pt1, second = pt2;

// Approximate the distance from the outer curve for the first

// pass. Set some limits for our target distance from the last

// curve

double bestDist = radius,

bestDepth = radius,

upperDistLimit = 2.1 * radius,

lowerDistLimit = 2 * radius;

// We'll count the number of iterations

int tries = 0;

bool res = true;

Curve cur = null;

do

{

// If we're creating ring in a direction along the surface,

// we need to transform the points defining the section

if (direction != 0)

{

Matrix3d disp =

Matrix3d.Displacement(direction * unitVec * bestDist);

pt1 = pt1.TransformBy(disp);

pt2 = pt2.TransformBy(disp);

}

// Get the curve from this section

cur = CurveFromSection(pt1, pt2);

if (cur == null)

{

FillRemainingCurve(prevCur, radius, prevDep, layer);

res = false;

break;

}

// Get the distance between the proposed position of the

// first sphere and the first sphere of the previous ring

Point3d cen = GetFirstCenter(cur, bestDepth);

double center2center = (cen - prevCen).Length;

// If this is either the first time run (i.e. no previous

// sphere) or the spheres are at an acceptable distance,

// skip the trigonometric calculations

if (

prevCen.DistanceTo(Point3d.Origin) >

Tolerance.Global.EqualPoint &&

(center2center < lowerDistLimit ||

center2center > upperDistLimit)

)

{

// Only do the calculations if this is the first run,

// otherwise we will just adjust the values

if (tries == 0)

{

// Calculate the distance to the next section

// Get slope at the current point on the surface

// (ideal would be to get the tangent at the mid-

// point of the spheres, but we have to approximate)

Vector3d slope =

GetSlopeAtPoint(pt1, pt2, unitVec, radius);

if (slope.Length < Tolerance.Global.EqualPoint)

{

res = false;

break;

}

// Get the "down" vector based on the curve direction

Vector3d down =

cur.GetPointAtParameter(cur.EndParam / 2) -

cur.GetPointAtParameter(cur.StartParam);

// Get the viewing plane

using (Plane p = newPlane(Point3d.Origin, _viewDir))

{

// Project the various vectors onto the plane

Vector2d unit2d = unitVec.Convert2d(p),

slope2d = slope.Convert2d(p),

down2d = down.Convert2d(p);

// Get the angle of the slope to the offset direction

double theta = unit2d.GetAngleTo(slope2d);

// Get the angle of the slope to the down direction

double theta2 = down2d.GetAngleTo(slope2d);

// Calculate the distance and the depth

bestDist = Math.Abs(2 * radius * Math.Cos(theta));

bestDepth = Math.Abs(radius / Math.Sin(theta2));

}

}

else

{

// Incrementally adjust the distance upwards or

// downwards, depending on which side of ideal we are

if (center2center < lowerDistLimit)

bestDist *= 1.05;

elseif (center2center > upperDistLimit)

bestDist *= 0.95;

}

// Reset the points, so that we start from scratch

// but with (hopefully) a better distance value

pt1 = first;

pt2 = second;

tries++;

// Max out at 20 iterations (gives much better coverage

// than 10)

if (tries > 20)

{

double curLen =

cur.GetDistanceAtParameter(cur.EndParam);

double prevLen =

prevCur.GetDistanceAtParameter(prevCur.EndParam);

// Only fill previous curves when the current one

// is much shorter than the previous one

if (curLen < prevLen / 4)

FillRemainingCurve(prevCur, radius, prevDep, layer);

res = false;

break;

}

}

else

{

// We have spheres at an acceptable distance from

// each other (or this is the first one)

prevCen = cen;

GenerateSpheresAlongCurve(

cur, radius, bestDepth, layer

);

res = true;

break;

}

}

while (true);

if (prevCur != null)

prevCur.Dispose();

prevCur = cur;

return res;

}

privatevoid FillRemainingCurve(

Curve cur, double radius, double initDepth, int layer

)

{

Curve inner = null;

try

{

inner =

(initDepth == 0.0 ? cur : GetInnerCurve(cur, initDepth));

if (inner == null)

return;

// Fill the remaining space as best we can

double inc = 2 * radius,

depth = inc;

bool done = false;

while (!done)

{

done =

!GenerateSpheresAlongCurve(

inner, radius, depth, layer

);

depth += inc;

}

}

catch { }

finally

{

if (inner != null && initDepth > 0.0)

inner.Dispose();

}

}

// Return the slope of the solid's surface at a particular

// section location

privateVector3d GetSlopeAtPoint(

Point3d pt1, Point3d pt2, Vector3d unitVec, double radius

)

{

// To get the slope of the surface at the point we care

// about, we'll take a section either side of the

// point - using 5% of the radius as the distance in

// each direction

constdouble factor = 0.05;

// Transform the points defining our section line in one

// direction

Matrix3d disp =

Matrix3d.Displacement(radius * factor * unitVec);

Point3d pt11 = pt1.TransformBy(disp);

Point3d pt12 = pt2.TransformBy(disp);

// And then the other direction

disp = Matrix3d.Displacement(-radius * factor * unitVec);

Point3d pt21 = pt1.TransformBy(disp);

Point3d pt22 = pt2.TransformBy(disp);

// Get the defining curves of these two sections

Curve c1 = CurveFromSection(pt11, pt12);

Curve c2 = CurveFromSection(pt21, pt22);

// Assume their start points are at the same relative point

if (c1 != null && c2 != null)

{

// Return the angle between the two

return c1.StartPoint - c2.StartPoint;

}

returnnewVector3d(0, 0, 0);

}

privateCurve CurveFromSection(Point3d pt1, Point3d pt2)

{

// Make a collection from our defining points

Point3dCollection pts =

newPoint3dCollection(newPoint3d[2] { pt1, pt2 });

// Declare the arrays for our results

Array fillEnts, bgEnts, fgEnts, fTangEnts, cTangEnts;

// Create the section and generate the geometry

using (Section sec = newSection(pts, _viewDir))

{

sec.GenerateSectionGeometry(

_sol, out fillEnts, out bgEnts,

out fgEnts, out fTangEnts, out cTangEnts

);

// We only care about the fillEnts - dispose of the rest

foreach (Entity e in bgEnts)

e.Dispose();

foreach (Entity e in fTangEnts)

e.Dispose();

foreach (Entity e in cTangEnts)

e.Dispose();

if (fillEnts.Length == 0)

returnnull;

// Return the first curve, dispose of any remaining

Curve c = fillEnts.GetValue(0) asCurve;

for (int i = 1; i < fillEnts.Length; i++)

{

DBObject o = fillEnts.GetValue(i) asDBObject;

if (o != null)

o.Dispose();

}

return c;

}

}

// Get an inner curve offset at a certain depth

privatestaticCurve GetInnerCurve(Curve c, double depth)

{

try

{

// We only return a single curve, if there is one.

// If there are more, we return null

DBObjectCollection inners = c.GetOffsetCurves(-depth);

if (inners.Count == 1)

return inners[0] asCurve;

foreach (DBObject o in inners)

o.Dispose();

}

catch { }

returnnull;

}

// Return the center of the 1st sphere to be placed on a curve

privatestaticPoint3d GetFirstCenter(Curve c, double radius)

{

// Only handle closed curves

if (c.Closed)

{

// Get the inner curve at the depth of the radius

using (Curve inner = GetInnerCurve(c, radius))

{

// Return the start point of the curve (which we assume

// will be at a consistent point across the various

// section profiles)

if (inner != null)

return inner.StartPoint;

}

}

returnPoint3d.Origin;

}

// Generate adjacent spheres along a curve

privatebool GenerateSpheresAlongCurve(

Curve curve, double radius, double depth, int layer

)

{

// Only deal with closed curves

if (!curve.Closed)

returnfalse;

try

{

// Get the inner curve at a certain depth

using (Curve inner = GetInnerCurve(curve, depth))

{

// Return null if we didn't get a curve or if it's too

// short to realistically fill with spheres

if (

inner == null ||

inner.GetDistanceAtParameter(inner.EndParam) <

8.85 * radius

)

{

if (inner != null)

{

// For shorter curves we'll create one last

// sphere at the center of the curve, assuming

// it's a circle or an ellipse (spheres

// generally generate ellipse section curves,

// it seems)

if (inner isEllipse || inner isCircle)

{

Point3d lastCen =

(inner isEllipse ?

((Ellipse)inner).Center :

((Circle)inner).Center);

CreateSphereAtPoint(lastCen, radius, layer);

}

}

returnfalse;

}

// Get the curve's plane

Plane p = inner.GetPlane();

// Get the center of the first sphere (the curve's start

// point) and make a local copy of it

Point3d cen = inner.StartPoint,

initial = cen;

// We'll step along the curve parametrically

double curParam = inner.StartParam;

// A flag to indication completion

bool done = false;

// We'll loop until we can't create more spheres

while (!done)

{

// Create a circle at the first location

CreateSphereAtPoint(cen, radius, layer);

// Get the center of the next sphere: cen and curParam

// will get updated, and lastRadius will only get set

// in the case where there's only space to create a

// smaller sphere

double lastRadius;

done =

!GetNextSphereCenter(

inner, p, radius, initial, ref cen, ref curParam,

out lastRadius

);

// If we're at the end of the ring, create the last

// sphere

if (done && lastRadius > 0.0)

CreateSphereAtPoint(cen, lastRadius, layer);

}

returntrue;

}

}

catch

{

returnfalse;

}

}

// Get the next sphere's center along the ring

privatebool GetNextSphereCenter(

Curve curve, Plane p, double radius, Point3d initial,

refPoint3d cen, refdouble curParam, outdouble lastRadius

)

{

// We may not return a lastRadius, so default to 0

lastRadius = 0.0;

// Make a local copy of the previous sphere's center

Point3d prevCen = cen;

// We're going to use a circle on the plane of the curve

// to find the center of our next sphere

using (Circle c = newCircle(cen, p.Normal, 2 * radius))

{

// The first - at the previous sphere's center - gets

// intersected with the curve

Point3dCollection pts = newPoint3dCollection();

curve.IntersectWith(

c, Intersect.OnBothOperands, pts,

IntPtr.Zero, IntPtr.Zero

);

// For each of the intersection points, we get the parameter

foreach (Point3d pt in pts)

{

double param = curve.GetParameterAtPoint(pt);

// If the point's parameter is greater than the current

// one, we'll use it

if (param > curParam)

{

double distFromInitial = (pt - initial).Length;

if (distFromInitial >= (2 * radius * 0.98))

{

cen = pt;

curParam = param;

returntrue;

}

else

{

// Fill the gap at the end of a ring by reducing

// the radius and adjusting the center to be the

// mid-point on the curve between the first and

// the previous spheres' centers

curParam =

curParam + ((curve.EndParam - curParam) / 2);

cen = curve.GetPointAtParameter(curParam);

double distFromInitial2 = (cen - initial).Length;

lastRadius = distFromInitial2 - radius;

returnfalse;

}

}

break;

}

}

returnfalse;

}

// Create a sphere at the given point

privatevoid CreateSphereAtPoint(

Point3d cen, double radius, int layer

)

{

// Just make sure the center isn't at the origin (which

// we've naughtily been using to denote an error) and

// that the radius is greater than zero

if (

cen.DistanceTo(Point3d.Origin) > Tolerance.Global.EqualPoint

&& radius > Tolerance.Global.EqualPoint

)

{

// Create and position our sphere

Solid3d sol = newSolid3d();

// If we're subtracting, make the sphere 2% smaller, so

// they don't touch each other

sol.CreateSphere(_sub ? radius * 0.98 : radius);

// Move the sphere to the desired location

Matrix3d disp = Matrix3d.Displacement(cen - Point3d.Origin);

sol.TransformBy(disp);

// Set its layer

sol.Layer = GetLayer(layer);

// Add it to the modelspace and the transaction

AddToDrawing(sol);

// If we're subtracting from the outer solid, do so

if (_sub)

{

_parent.BooleanOperation(

BooleanOperationType.BoolSubtract, sol

);

}

}

}

// Add an entity to the open block table record and the

// current transaction

privatevoid AddToDrawing(Entity ent)

{

_btr.AppendEntity(ent);

_tr.AddNewlyCreatedDBObject(ent, true);

}

}

}

Let’s see what we get when we run the PACKSPHERES command selecting a sphere and a section line right through the middle.

At the top-left of the below image is the wireframe view of the sphere prior to packing, top-middle is the packed sphere, also in wireframe (and prior to graphics regeneration, which would lead to the inner layers not showing as well), and the subsequent shots are of the various layers of the solid in realistic view (one more layer of the onion is peeled off in each shot).

I didn’t bother hiding the last layer, as we can – in any case – see the inner solid and I prefer the symmetry of the below image as it stands. The inner solid is the one we have been using to generate each layer, and has been created by iterative offsets from the original sphere: by design we choose to stop the process - leaving the inner body – when the radius of the spheres we would be creating would be relatively too large to create a useful layer.

If you look at the model from a different angle, you’ll see some quirks: the poles have some gaps, for instance, and as spheres of a particular radius are not guaranteed to fill a ring completely, we fill the remainder with smaller spheres, where needed.

There are alternative approaches which would avoid this “fault line” of small spheres: we could adjust the radius of the spheres very slightly to make sure they are even over the length of the curve, or we could rotate the whole curve by a random amount. Each approach ultimately has advantages and disadvantages, so I’ve left it as is for this initial implementation.

One benefit of an approach that takes sections across the length of a solid is that it can pack different forms full of spheres, too. For instance, here’s a bottle:

And here’s a tree-like form (the algorithm needs some work to deal with sudden changes in size of the section curve, such as when we go from the lower part of the “leafy part” to the trunk):

The PACKSPHERES command does have an option we haven’t looked at in this post, which is to automatically subtract the inner spheres from the outer solid. This can take a long time to work, but generates results that are quite like those in the last post (where we performed more-or-less the same process manually).

Here’s a section view of the resultant hollowed bottle when that option is selected (the size of the individual holes clearly depends on where the spheres happen to be relative to the section plane):

AutoCAD’s STLOUT command generated a 257 MB STL file for the above solid. Thanks to Alex Fielder for suggesting netfabb as a tool for viewing STL files (with the caveat that he hadn’t tested it on files quite that large): while netfabb Studio Basic managed to load the file, it unfortunately had difficulty slicing it.

Next week we’re going to have a change of pace: as the news aboutAutoCAD 2013 is starting to hit the web, I’ll spend the week talking about the coming developer-oriented features in the 2013 release of AutoCAD.

February 22, 2012

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

February 20, 2012

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 :-).

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

letmutable minRad = 0.01

letmutable maxGen = 6

letmutable 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

letrec 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() =

staticmember 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

{

publicclassCommands

{

[CommandMethod("AGFS")]

publicstaticvoid ApollonianGasket()

{

ApollonianGasket(true);

}

[CommandMethod("AGFSE")]

publicstaticvoid ExteriorApollonianGasket()

{

ApollonianGasket(false, true);

}

[CommandMethod("AGFSI")]

publicstaticvoid InteriorApollonianGasket()

{

ApollonianGasket(false, false);

}

publicstaticvoid 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 = newLayerTableRecord();

ltr.Color =

Autodesk.AutoCAD.Colors.Color.FromColorIndex(

Autodesk.AutoCAD.Colors.ColorMethod.ByAci, i

);

ltr.Name = name;

lt.Add(ltr);

tr.AddNewlyCreatedDBObject(ltr, true);

}

}

}

publicstaticvoid 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 = newSolid3d();

s.CreateSphere(radius);

Point3d center =

newPoint3d(

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:

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.

February 16, 2012

So far in this series, we’ve looked at Apollonian circle packing using C# and alsoF#. 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 academicpapers 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)

letmutable 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

*)

letmutable 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() =

staticmember 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

staticmember 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

{

publicclassCommands

{

[CommandMethod("AGFSE")]

publicstaticvoid ExteriorApollonianGasket()

{

ApollonianGasket(true);

}

[CommandMethod("AGFSI")]

publicstaticvoid InteriorApollonianGasket()

{

ApollonianGasket(false);

}

publicstaticvoid 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 = newLayerTableRecord();

ltr.Color =

Autodesk.AutoCAD.Colors.Color.FromColorIndex(

Autodesk.AutoCAD.Colors.ColorMethod.ByAci, i

);

ltr.Name = name;

lt.Add(ltr);

tr.AddNewlyCreatedDBObject(ltr, true);

}

}

}

publicstaticvoid 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 = newSolid3d();

s.CreateSphere(radius);

Point3d center =

newPoint3d(

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.

February 14, 2012

After initially holding DevCamps annually in 2007 and 2008, we’ve now settled down to the more regular rhythm of every other year (the last set was in 2010). Which means, of course, it’s time to do them again. :-)

As already reported by Geoff, Heidi, Jeremy and Shaan, here is some information regarding this year’s DevCamps:

Cost: early bird discount fee of US $500 valid to April 30, 2012, thereafter US $600

You probably don’t need me to extoll the virtues of these events – hopefully those are evident from the other posts – but I will just say that they are really great opportunity both to acquire technical knowledge and to network, whether with other people interested in software development, with members of the ADN team, with the product managers and software developers defining and building our products, or even with Autodesk’s divisional heads.

A particular area of focus, this year, will be on web and mobile technology, which reflects the shift our industry is going through. That doesn’t mean we’ll have less focus on our desktop software, though: expect excellent information on developing for those, too, aimed both at beginners and experts.

To stay abreast of the event details as they get nailed down (such as the specific agenda, registration information, etc.) be sure to email Melrose Ross, stating your full name and the company name.

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

letrec 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() =

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

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

letrec 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() =

staticmember 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

{

publicclassCommands

{

[CommandMethod("AGFCI")]

publicstaticvoid InnerApollonianGasket()

{

Document doc =

Application.DocumentManager.MdiActiveDocument;

Database db = doc.Database;

Editor ed = doc.Editor;

PromptEntityOptions peo =

newPromptEntityOptions("\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) asCircle;

}

else

{

// Otherwise create a new one at the default location

cir = newCircle(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 - newPoint3d(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 =

newCircle(

newPoint3d(

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

February 08, 2012

Hot on the heels of the Revit and Inventor installments in this series of guides (OK, OK – perhaps not exactly hot, but then “warm on the heels” doesn’t have quite the same ring to it ;-), I’m happy to announce that the previouslymentioned “My First AutoCAD Plug-in” guide is now available.

Stephen Preston, the author of this excellent guide, has done a great job of creating a compelling sample – which keeps your block attributes horizontal, irrespective of the block rotation – and a series of straightforward explanations to help power-users and programming “newbies” take the plunge into the world of AutoCAD programming.