August 2014

Sun Mon Tue Wed Thu Fri Sat
          1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30
31            








« Gathering points defining 2D AutoCAD geometry using .NET | Main | Gathering points defining 3D AutoCAD geometry using .NET »

February 09, 2011

Creating the smallest possible circle around 2D AutoCAD geometry using .NET

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.

Smallest circles around polylines

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

Visual Computing - Geometry, Graphics, and VisionBoth 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:

Circles enclosing various block references

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.

blog comments powered by Disqus

Feed/Share

10 Random Posts