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.
Subscribe via RSS