Creating the smallest possible sphere around 3D AutoCAD geometry using .NET

mccad

Given the last few posts – where we gathered 2D points, created a minimal circle around them, and then gathered 3D points – the topic of this post shouldn’t come as much of a surprise. :-)

As suggested, the changes needed to support 3D in the algorithm we used to solve the
smallest circle problem were trivial: we had to adjust the function protocol to pass in a list of 3D points and the implementation to deal with the Z coordinates provided. Looking at the code again, I did consider adjusting it to pass in a Point3dCollection rather than a List<Point3d> (I has chosen to use a generic list when implementing the original recursive algorithm and had just kept it with the newer iterative one), but then decided to leave that particular aspect alone.

Here’s the updated C# code:

using Autodesk.AutoCAD.ApplicationServices;

using Autodesk.AutoCAD.DatabaseServices;

using
DbSurface =

  Autodesk.AutoCAD.DatabaseServices.Surface;

using
DbNurbSurface =

  Autodesk.AutoCAD.DatabaseServices.NurbSurface;

using Autodesk.AutoCAD.EditorInput;

using Autodesk.AutoCAD.Runtime;

using Autodesk.AutoCAD.Geometry;

using System.Collections.Generic;

using System;

namespace MinimumEnclosingSphere

{

public
class
Commands

  {

    [CommandMethod("MES", CommandFlags.UsePickSet)]

public
void MinimumEnclosingSphere()

    {

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 oneSpherePerEnt = false;

if (psr.Value.Count > 1)

      {

PromptKeywordOptions pko =

new
PromptKeywordOptions(

"\nMultiple objects selected: create " +

"individual spheres 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;

        oneSpherePerEnt = (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;

        }

      }

// 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.ObjectId,

OpenMode.ForRead

            );

// Collect the points for each selected entity

          CollectPoints(tr, ent, pts);

// Create a sphere for each entity (if so chosen) or

// just once after collecting all the points

if (oneSpherePerEnt || i == psr.Value.Count - 1)

          {

try

            {

Solid3d sol = SphereFromPoints(pts, buffer);

              btr.AppendEntity(sol);

              tr.AddNewlyCreatedDBObject(sol, true);

            }

catch

            {

              ed.WriteMessage(

"\nUnable to calculate enclosing sphere."

              );

            }

            pts.Clear();

          }

        }

        tr.Commit();

      }

    }

private
Solid3d SphereFromPoints(

Point3dCollection pts, double buffer

    )

    {

// Copy the points across

List<Point3d> pts3d = new
List<Point3d>(pts.Count);

for (int i = 0; i < pts.Count; i++)

      {

        pts3d.Add(pts);

      }

// Assuming we have some points in our list...

if (pts.Count > 0)

      {

// We need the center and radius of our sphere

Point3d center;

double radius = 0;

// Use our fast approximation algorithm to

// calculate the center and radius of our

// sphere to within 1% (calling the function

// with 100 iterations gives 10%, calling it

// with 10K gives 1%)

        BadoiuClarksonIteration(

          pts3d, 10000, out center, out radius

        );

// Create the sphere and return it

Solid3d sol = new
Solid3d();

        sol.CreateSphere(radius * (1.0 + buffer));

        sol.TransformBy(

Matrix3d.Displacement(center - Point3d.Origin)

        );

return sol;

      }

return
null;

    }

private
void CollectPoints(

Transaction tr, Entity ent, Point3dCollection pts

    )

    {

// 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);

          }

        }

      }

// For surfaces we'll do something similar: we'll

// collect points across its surface (by first getting

// NURBS surfaces for the surface) and will still

// explode the surface later to get points from the

// curves

DbSurface sur = ent as
DbSurface;

if (sur != null)

      {

DbNurbSurface[] nurbs = sur.ConvertToNurbSurface();

foreach (DbNurbSurface nurb in nurbs)

        {

// Calculate the parameter increments in u and v

double ustart = nurb.UKnots.StartParameter,

                uend = nurb.UKnots.EndParameter,

                uinc = (uend - ustart) / nurb.UKnots.Count,

                vstart = nurb.VKnots.StartParameter,

                vend = nurb.VKnots.EndParameter,

                vinc = (vend - vstart) / nurb.VKnots.Count;

// Pick up points across the surface

for (double u = ustart; u <= uend; u += uinc)

          {

for (double v = vstart; v <= vend; v += vinc)

            {

              pts.Add(nurb.Evaluate(u, v));

            }

          }

        }

      }

// For 3D solids we'll fire a number of rays from the

// centroid in random directions in order to get a

// sampling of points on the outside

Solid3d sol = ent as
Solid3d;

if (sol != null)

      {

for (int i = 0; i < 500; i++)

        {

Solid3dMassProperties mp = sol.MassProperties;

using (Plane pl = new
Plane())

          {

            pl.Set(mp.Centroid, randomVector3d());

using (Region reg = sol.GetSection(pl))

            {

using (Ray ray = new
Ray())

              {

                ray.BasePoint = mp.Centroid;

                ray.UnitDir = randomVectorOnPlane(pl);

                reg.IntersectWith(

                  ray, Intersect.OnBothOperands, pts,

IntPtr.Zero, IntPtr.Zero

                );

              }

            }

          }

        }

      }

// Now we start the terminal cases - for basic objects -

// before we recurse for more complex objects (including

// the ones we've already collected points for above).

// 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 s = (Solid)ent;

try

        {

for (short i = 0; i < 4; i++)

          {

            pts.Add(s.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)

              {

                CollectPoints(tr, ent2, pts);

              }

              obj.Dispose();

            }

          }

        }

catch { }

      }

    }

// Return a random 3D vector on a plane

private
Vector3d randomVectorOnPlane(Plane pl)

    {

// Create our random number generator

Random ran = new
Random();

// First we get the absolute value

// of our x, y and z coordinates

double absx = ran.NextDouble();

double absy = ran.NextDouble();

// Then we negate them, half of the time

double x = (ran.NextDouble() < 0.5 ? -absx : absx);

double y = (ran.NextDouble() < 0.5 ? -absy : absy);

Vector2d v2 = new
Vector2d(x,y);

return
new
Vector3d(pl, v2);

    }

// Return a random 3D vector

private
Vector3d randomVector3d()

    {

// Create our random number generator

Random ran = new
Random();

// First we get the absolute value

// of our x, y and z coordinates

double absx = ran.NextDouble();

double absy = ran.NextDouble();

double absz = ran.NextDouble();

// Then we negate them, half of the time

double x = (ran.NextDouble() < 0.5 ? -absx : absx);

double y = (ran.NextDouble() < 0.5 ? -absy : absy);

double z = (ran.NextDouble() < 0.5 ? -absz : absz);

return
new
Vector3d(x, y, z);

    }

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)

              )

            );

          }

        }

      }

    }

// Algorithm courtesy (and copyright of) Frank Nielsen

// http://blog.informationgeometry.org/article.php?id=164

public
void BadoiuClarksonIteration(

List<Point3d> set, int iter,

out
Point3d 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
Point3d(

            cen.X + (1.0 / (i + 1.0)) * (set[winner].X - cen.X),

            cen.Y + (1.0 / (i + 1.0)) * (set[winner].Y - cen.Y),

            cen.Z + (1.0 / (i + 1.0)) * (set[winner].Z - cen.Z)

          );

      }

    }

  }

}

We’ll run the new MES (Minimal Enclosing Sphere) command against the 3D geometry we saw in the last post…

… after setting ENCLOSINGCIRCLEBUFFER (we haven’t bothered implementing a separate ENCLOSINGSPHEREBUFFER, but we very easily could, if needed) to zero.

Command: ENCLOSINGCIRCLEBUFFER

Enter new value for ENCLOSINGCIRCLEBUFFER <3>: 0

Command: MES

Select objects to enclose: ALL

9 found

Select objects to enclose:

Multiple objects selected: create individual spheres around each one? [Yes/No] <No>: Y

Here are the results:

Which are perhaps more visible with the x-ray visual style applied:

Such operations certainly have the potential to be time-consuming… the code should probably be adjusted to show a progress meter or – at the very least – check HostApplicationServices.Current.UserBreak() from time to time. Something I’d certainly do myself, if preparing this code for production use, but for now left as an exercise for the reader.