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[i].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[i]);
}
// 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.
