Now it’s time to shed some light on the reason for the code in the last post. I wrote it to help address a question that came in from Elson Brown via our Plugin of the Month feedback alias:
I have a request for an app that will draw the smallest circle around a polyline object. Pick the polyline and the app draws the smallest possible circle around the shape.
It turns out this is a well-known problem (and thanks to Stephen Preston for pointing me in this direction), commonly known as the minimal enclosing circle or smallest circle problem.
This struck me as a fun problem to solve inside AutoCAD (and probably some of our other products, too), so I started looking at the (typically C++) code projects linked to from the Wikipedia page. Having taken a deep breath, I started looking at Bernd Gärtner's smallest enclosing ball of points code, to see whether I could make sense of it. Sadly it was a little too complex for me to bring across to C#, but it started me in a good direction.
I considered – briefly – the idea of just using an open source library, whether CGAL (for computational geometry) or OpenCV (for computer vision), both of which provide functions that solve this problem. It just seemed overkill: I don’t currently need any of the other functions in either of these very impressive libraries, and using them would most probably add complexity when supporting 32- vs. 64-bit platforms, etc.
So back to implementing something myself (or rather porting something from someone else :-). Searching on “miniball”, I came across this C++ implementation from Frank Nielsen, Professor of Computer Science at the LIX (le Laboratoire d’Informatique de l'École Polytechnique) in France and Senior Researcher at Sony CSL FRL (the Fundamental Research Laboratory of Sony Computer Science Laboratories). I managed to implement this in C# – having contacted Frank to get his permission to do so – and it worked well: the only problem was with the fact it’s a recursive algorithm, and calling it with lots of points (extracted from the selected AutoCAD geometry using the technique shown in the last post) leads to a stack overflow, a common phenomenon with deeply recursive code. I did some work to reduce the set of points extracted from the geometry – using LINQ to filter out points near the centre of the circle, for instance – but the approach could still have caused problems when a certain threshold was exceeded. It should be noted that our use of recursion in the last post – where we recurse for each level of containment, whether for blocks or “complex” entities such as polylines – results in very few recursive calls (I’d be surprised if there were 10 additional items in the call stack because of the recursion), whereas the miniball algorithm was recursing for each point (as far as I could tell).
Thankfully Frank came to my rescue, and suggested another, iterative algorithm, based on work by Bădoiu and Clarkson, to approximate a solution to this problem. Thankfully there was also a C++ implementation of the algorithm available, which helped me a lot (working from the mathematical description of an algorithm makes my head hurt – I’m just not used to thinking like that).
Both of these implementations come from the material for one of Frank’s books entitled Visual Computing: Geometry, Graphics, and Vision. Given what I’ve seen so far from this book, I’m certainly going to pick up a copy of it (and hopefully some of you will, too, to repay Frank for his kind contributions to this post :-).
The iterative algorithm worked a treat. If you iterate 10,000 times (which happens quickly enough) there is a 1% margin of error, which is altogether acceptable (especially as the code allows you to add a buffer percentage to the calculated radius, so that we’re not quite touching the enclosed geometry).
It’s worth noting that the algorithm can work on n-dimensional points, so I’ll probably be using it to create a sphere, at some point (I don’t think I’ll personally be needing to solve the problem for any more than three dimensions, though :-).
Here’s the C# code that builds on the code shown in the last post:
using Autodesk.AutoCAD.ApplicationServices;
using Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.EditorInput;
using Autodesk.AutoCAD.Runtime;
using Autodesk.AutoCAD.Geometry;
using System.Collections.Generic;
namespace MinimumEnclosingCircle
{
public class Commands
{
[CommandMethod("MEC", CommandFlags.UsePickSet)]
public void MinimumEnclosingCircle()
{
Document doc =
Application.DocumentManager.MdiActiveDocument;
Database db = doc.Database;
Editor ed = doc.Editor;
// Ask user to select entities
PromptSelectionOptions pso =
new PromptSelectionOptions();
pso.MessageForAdding = "\nSelect objects to enclose: ";
pso.AllowDuplicates = false;
pso.AllowSubSelections = true;
pso.RejectObjectsFromNonCurrentSpace = true;
pso.RejectObjectsOnLockedLayers = false;
PromptSelectionResult psr = ed.GetSelection(pso);
if (psr.Status != PromptStatus.OK)
return;
bool oneCircPerEnt = false;
if (psr.Value.Count > 1)
{
PromptKeywordOptions pko =
new PromptKeywordOptions(
"\nMultiple objects selected: create " +
"individual circles around each one?"
);
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;
oneCircPerEnt = (pkr.StringResult == "Yes");
}
// There may be a SysVar defining the buffer
// to add to our radius
double buffer = 0.0;
try
{
object bufvar =
Application.GetSystemVariable(
"ENCLOSINGCIRCLEBUFFER"
);
if (bufvar != null)
{
short bufval = (short)bufvar;
buffer = bufval / 100.0;
}
}
catch
{
object bufvar =
Application.GetSystemVariable("USERI1");
if (bufvar != null)
{
short bufval = (short)bufvar;
buffer = bufval / 100.0;
}
}
// Get the current UCS
CoordinateSystem3d ucs =
ed.CurrentUserCoordinateSystem.CoordinateSystem3d;
// Collect points on the component entities
Point3dCollection pts = new Point3dCollection();
Transaction tr =
db.TransactionManager.StartTransaction();
using (tr)
{
BlockTableRecord btr =
(BlockTableRecord)tr.GetObject(
db.CurrentSpaceId,
OpenMode.ForWrite
);
for (int i = 0; i < psr.Value.Count; i++)
{
Entity ent =
(Entity)tr.GetObject(
psr.Value[i].ObjectId,
OpenMode.ForRead
);
// Collect the points for each selected entity
Point3dCollection entPts = CollectPoints(tr, ent);
foreach (Point3d pt in entPts)
{
/*
* Create a DBPoint, for testing purposes
*
DBPoint dbp = new DBPoint(pt);
btr.AppendEntity(dbp);
tr.AddNewlyCreatedDBObject(dbp, true);
*/
pts.Add(pt);
}
// Create a circle for each entity (if so chosen) or
// just once after collecting all the points
if (oneCircPerEnt || i == psr.Value.Count - 1)
{
try
{
Circle cir =
CircleFromPoints(pts, ucs, buffer);
btr.AppendEntity(cir);
tr.AddNewlyCreatedDBObject(cir, true);
}
catch
{
ed.WriteMessage(
"\nUnable to calculate enclosing circle."
);
}
pts.Clear();
}
}
tr.Commit();
}
}
private Point3dCollection CollectPoints(
Transaction tr, Entity ent
)
{
// The collection of points to populate and return
Point3dCollection pts = new Point3dCollection();
// We'll start by checking a block reference for
// attributes, getting their bounds and adding
// them to the point list. We'll still explode
// the BlockReference later, to gather points
// from other geometry, it's just that approach
// doesn't work for attributes (we only get the
// AttributeDefinitions, which don't have bounds)
BlockReference br = ent as BlockReference;
if (br != null)
{
foreach (ObjectId arId in br.AttributeCollection)
{
DBObject obj = tr.GetObject(arId, OpenMode.ForRead);
if (obj is AttributeReference)
{
AttributeReference ar = (AttributeReference)obj;
ExtractBounds(ar, pts);
}
}
}
// If we have a curve - other than a polyline, which
// we will want to explode - we'll get points along
// its length
Curve cur = ent as Curve;
if (cur != null &&
!(cur is Polyline ||
cur is Polyline2d ||
cur is Polyline3d))
{
// Two points are enough for a line, we'll go with
// a higher number for other curves
int segs = (ent is Line ? 2 : 20);
double param = cur.EndParam - cur.StartParam;
for (int i = 0; i < segs; i++)
{
try
{
Point3d pt =
cur.GetPointAtParameter(
cur.StartParam + (i * param / (segs - 1))
);
pts.Add(pt);
}
catch { }
}
}
else if (ent is DBPoint)
{
// Points are easy
pts.Add(((DBPoint)ent).Position);
}
else if (ent is DBText)
{
// For DBText we use the same approach as
// for AttributeReferences
ExtractBounds((DBText)ent, pts);
}
else if (ent is MText)
{
// MText is also easy - you get all four corners
// returned by a function. That said, the points
// are of the MText's box, so may well be different
// from the bounds of the actual contents
MText txt = (MText)ent;
Point3dCollection pts2 = txt.GetBoundingPoints();
foreach (Point3d pt in pts2)
{
pts.Add(pt);
}
}
else if (ent is Face)
{
Face f = (Face)ent;
try
{
for (short i = 0; i < 4; i++)
{
pts.Add(f.GetVertexAt(i));
}
}
catch { }
}
else if (ent is Solid)
{
Solid sol = (Solid)ent;
try
{
for (short i = 0; i < 4; i++)
{
pts.Add(sol.GetPointAt(i));
}
}
catch { }
}
else
{
// Here's where we attempt to explode other types
// of object
DBObjectCollection oc = new DBObjectCollection();
try
{
ent.Explode(oc);
if (oc.Count > 0)
{
foreach (DBObject obj in oc)
{
Entity ent2 = obj as Entity;
if (ent2 != null && ent2.Visible)
{
foreach (Point3d pt in CollectPoints(tr, ent2))
{
pts.Add(pt);
}
}
obj.Dispose();
}
}
}
catch { }
}
return pts;
}
private void ExtractBounds(
DBText txt, Point3dCollection pts
)
{
// We have a special approach for DBText and
// AttributeReference objects, as we want to get
// all four corners of the bounding box, even
// when the text or the containing block reference
// is rotated
if (txt.Bounds.HasValue && txt.Visible)
{
// Create a straight version of the text object
// and copy across all the relevant properties
// (stopped copying AlignmentPoint, as it would
// sometimes cause an eNotApplicable error)
// We'll create the text at the WCS origin
// with no rotation, so it's easier to use its
// extents
DBText txt2 = new DBText();
txt2.Normal = Vector3d.ZAxis;
txt2.Position = Point3d.Origin;
// Other properties are copied from the original
txt2.TextString = txt.TextString;
txt2.TextStyleId = txt.TextStyleId;
txt2.LineWeight = txt.LineWeight;
txt2.Thickness = txt2.Thickness;
txt2.HorizontalMode = txt.HorizontalMode;
txt2.VerticalMode = txt.VerticalMode;
txt2.WidthFactor = txt.WidthFactor;
txt2.Height = txt.Height;
txt2.IsMirroredInX = txt2.IsMirroredInX;
txt2.IsMirroredInY = txt2.IsMirroredInY;
txt2.Oblique = txt.Oblique;
// Get its bounds if it has them defined
// (which it should, as the original did)
if (txt2.Bounds.HasValue)
{
Point3d maxPt = txt2.Bounds.Value.MaxPoint;
// Place all four corners of the bounding box
// in an array
Point2d[] bounds =
new Point2d[] {
Point2d.Origin,
new Point2d(0.0, maxPt.Y),
new Point2d(maxPt.X, maxPt.Y),
new Point2d(maxPt.X, 0.0)
};
// We're going to get each point's WCS coordinates
// using the plane the text is on
Plane pl = new Plane(txt.Position, txt.Normal);
// Rotate each point and add its WCS location to the
// collection
foreach (Point2d pt in bounds)
{
pts.Add(
pl.EvaluatePoint(
pt.RotateBy(txt.Rotation, Point2d.Origin)
)
);
}
}
}
}
private Circle CircleFromPoints(
Point3dCollection pts, CoordinateSystem3d ucs, double buffer
)
{
// Get the plane of the UCS
Plane pl = new Plane(ucs.Origin, ucs.Zaxis);
// We will project these (possibly 3D) points onto
// the plane of the current UCS, as that's where
// we will create our circle
// Project the points onto it
List<Point2d> pts2d = new List<Point2d>(pts.Count);
for (int i = 0; i < pts.Count; i++)
{
pts2d.Add(pl.ParameterOf(pts[i]));
}
// Assuming we have some points in our list...
if (pts.Count > 0)
{
// We need the center and radius of our circle
Point2d center;
double radius = 0;
// Use our fast approximation algorithm to
// calculate the center and radius of our
// circle to within 1% (calling the function
// with 100 iterations gives 10%, calling it
// with 10K gives 1%)
BadoiuClarksonIteration(
pts2d, 10000, out center, out radius
);
// Get our center point in WCS (on the plane
// of our UCS)
Point3d cen3d = pl.EvaluatePoint(center);
// Create the circle and add it to the drawing
return new Circle(
cen3d, ucs.Zaxis, radius * (1.0 + buffer)
);
}
return null;
}
// Algorithm courtesy (and copyright of) Frank Nielsen
// http://blog.informationgeometry.org/article.php?id=164
public void BadoiuClarksonIteration(
List<Point2d> set, int iter,
out Point2d cen, out double rad
)
{
// Choose any point of the set as the initial
// circumcenter
cen = set[0];
rad = 0;
for (int i = 0; i < iter; i++)
{
int winner = 0;
double distmax = (cen - set[0]).Length;
// Maximum distance point
for (int j = 1; j < set.Count; j++)
{
double dist = (cen - set[j]).Length;
if (dist > distmax)
{
winner = j;
distmax = dist;
}
}
rad = distmax;
// Update
cen =
new Point2d(
cen.X + (1.0 / (i + 1.0)) * (set[winner].X - cen.X),
cen.Y + (1.0 / (i + 1.0)) * (set[winner].Y - cen.Y)
);
}
}
}
}
One improvement suggested by Elson, after testing my initial implementation, was to provide the choice of creating multiple circles – one per selected entity – in the case where the user has selected multiple entities: a big time-saver if having to draw circles around a number of single polylines or block references.
As mentioned above, I realised quite early on that circles would not always want to be drawn to be touching (or even intersecting, given the margin of error) the enclosed geometry, so the above code supports the idea of a buffer. This is an amount – a percentage – by which to increase the radius when creating the circle. The code picks this value up via a new system variable (which can be added via the Registry from AutoCAD 2011) called ENCLOSINGCIRCLEBUFFER. Here’s the Registry file that can be used to add this system variable – an integer between 0 and 100 with a default value of 3:
Windows Registry Editor Version 5.00
[HKEY_LOCAL_MACHINE\SOFTWARE\Autodesk\AutoCAD\R18.1\ACAD-9001:409\Variables\ENCLOSINGCIRCLEBUFFER]
"StorageType"=dword:00000002
"LowerBound"=dword:00000000
"UpperBound"=dword:00000064
"PrimaryType"=dword:0000138b
@="3"
As not everyone can write to HKEY_LOCAL_MACHINE on their system – or is working with AutoCAD 2011 or higher – the code also will also use the value of the USERI1 system variable for this should it fail to find ENCLOSINGCIRCLEBUFFER.
Let’s take a look at the code working on the contents of Sample/Dynamic Blocks/Architectural – Metric.dwg. We’ll start by changing the buffer to 3% (after it having previously been set to zero in this session).
Command: ENCLOSINGCIRCLEBUFFER
Enter new value for ENCLOSINGCIRCLEBUFFER <0>: 3
Command: MEC
Select objects to enclose: ALL 9 found
Select objects to enclose:
Multiple objects selected: create individual circles around each one? [Yes/No] <No>: Y
Here are the changes to the geometry:
Overall the code should be reasonably straightforward – with the possible exception of the BadoiuClarksonIteration() function, which is certainly a black box to me.
It’s worth noting that while the points we gather from the geometry could be in 3D space – even if we’ve only really written the code to deal with 2D entities – we project them down to the plane of the UCS in order to create our circle there. As mentioned earlier, I do want to develop a 3D version of this – which shouldn’t be at all hard – to create a minimal enclosing sphere. The BadoiuClarksonIteration() function will need very little work – it would take a list of Point3d objects and also handle the Z coordinate of the center point – and we would need to develop the point collection code to also handle various types of solids and surfaces.
Once again, many thanks to Elson Brown for suggesting this topic – and helping with testing – and to Frank Nielsen for allowing me to adapt his code to address this problem.
Update
Added a missing Dispose() call.
