Hello
after some attempts, here is a not too bad script to compute the Minimum Bounding Sphere. I put 2 approaches for computing this.
One is incremental, it take 2 points to make a first sphere, then try to add another point …Random points is a way to fasten calculation, as calculations are long just one one point is outside bounding sphere. This algorithm also outputs points used to make the sphere.
It take on my PC 0.48 s for 35 000 points (no random) 0.3 s if random points.
The other approach uses Welzl method, with a part of the script from there
BEWARE it crashes rhino if more than 300 points (recursions ?) are chosen.
This site seems more clean
One of the flaws of most algorithm description is that they used circumscribed circle (3 points) or sphere (4 points). On many cases it is not the minimum Bounding Sphere. You’ll see in my script that I used these methods but inside method that compute Minimum Bounding Sphere, for 1, 2, 3, 4 or 5 points.
bounding_sphere_LEGACY.gh (1.4 MB)
private void RunScript(List<Point3d> lst_points, bool random, ref object A, ref object B, ref object C, ref object D)
{
//Laurent Delrieu
//18/12/2020
//https://discourse.mcneel.com/t/c-welzls-minimum-bounding-sphere/109246/7
// Use this one if you want to use Welzl method but beware there is a bug somewhere
//Don't put more than 300 points
// A = MinimumBoundingSphereWelzl.GetMinimumBoundingSphere(lst_points);
// B = lst_points;
List<Point3d> boundary_points;
C = MinimumBoundingSphereIncrement.GetMinimumBoundingSphere(lst_points, out boundary_points, 1000000, random);
D = boundary_points;
Component.Message = "Minimum Bounding Sphere";
}
// <Custom additional code>
/// <summary>
/// Class to calculate the minimum bounding sphere
/// it takes the first and second points to make a sphere then add one point at each increment
///Laurent Delrieu, 18/12/2020
/// </summary>
public static class MinimumBoundingSphereIncrement
{
//Tolerance to be inside the sphere
static double epsilon = 1e-12;
/// <summary>
///
/// </summary>
/// <param name="points">List of points</param>
/// <param name="boundaryPoints">Points that are used to make the sphere</param>
/// <param name="nMax">Max number of iterations, must be suppressed in future release</param>
/// <returns>Minimum sphere</returns>
public static Sphere GetMinimumBoundingSphere(List<Point3d> points, out List<Point3d> boundaryPoints, int nMax, bool random)
{
Random rnd = new Random();
boundaryPoints = new List<Point3d>();
List<Point3d> allPoints = new List<Point3d>(points);
Sphere sphere = Sphere.Unset;
if (allPoints.Count == 1)
{
sphere = new Sphere(allPoints[0], 0);
boundaryPoints.Add(allPoints[0]);
allPoints.RemoveAt(0);
}
else if (points.Count > 1)
{
sphere = Sphere2Points(allPoints[0], allPoints[1]);
boundaryPoints.Add(allPoints[0]);
boundaryPoints.Add(allPoints[1]);
allPoints.RemoveAt(1);
allPoints.RemoveAt(0);
}
int count = 0;
if (random)
{
while (allPoints.Count > 0 && count < nMax)
{
int next = rnd.Next(allPoints.Count);
sphere = AddPoint(sphere, ref boundaryPoints, allPoints[next]);
allPoints.RemoveAt(next);
count++;
}}
else
{
while (allPoints.Count > 0 && count < nMax)
{
sphere = AddPoint(sphere, ref boundaryPoints, allPoints[0]);
allPoints.RemoveAt(0);
count++;
}
}
return sphere;
}
/// <summary>
/// Add a point to the minimum sphere and update boundary points
/// </summary>
/// <param name="initialSphere">Initial bounding sphere</param>
/// <param name="boundaryPoints">List of boundary points as input and output</param>
/// <param name="newPoint">New point to add</param>
/// <returns>Updated or not minimum bounding sphere</returns>
public static Sphere AddPoint(Sphere initialSphere, ref List<Point3d> boundaryPoints, Point3d newPoint)
{
//If new point is inside sphere
Sphere output = initialSphere;
//New Point is not on the sphere
if (newPoint.DistanceTo(output.Center) > output.Radius)
{
// Print(" Ajout d'un point" + newPoint);
boundaryPoints.Add(newPoint);
output = Sphere5Points(boundaryPoints);
for (int i = boundaryPoints.Count - 2; i >= 0; i--)
{
if (ContainsStrict(output, boundaryPoints[i], epsilon))
{
// Print(" Suppression : " + boundaryPoints[i]);
boundaryPoints.RemoveAt(i);
}
}
//Print(" Boundary cout = " + boundaryPoints.Count);
if (boundaryPoints.Count >= 5)
{
boundaryPoints.RemoveAt(0);
}
}
return output;
}
}
/// <summary>
/// Minimum bounding sphere of a maximum of 5 points
/// </summary>
/// <param name="lst_points">List of points, no more than 5 points</param>
/// <returns>Bounding sphere</returns>
public static Sphere Sphere5Points(List<Point3d> lst_points)
{
//Default output
Sphere sphere = Sphere.Unset;
//If there is 0, 1, 2 3 or 4 points
if (lst_points.Count <= 4)
{
switch (lst_points.Count)
{
case 0: return Sphere.Unset;
case 1: return new Sphere(lst_points[0], 0);
case 2: return Sphere2Points(lst_points[0], lst_points[1]);
case 3: return Sphere3Points(lst_points[0], lst_points[1], lst_points[2]);
case 4: return Sphere4Points(lst_points[0], lst_points[1], lst_points[2], lst_points[3]);
}
}
//All poissible spheres around the 4 among the 5 points
Sphere sphere0 = Sphere4Points(lst_points[1], lst_points[2], lst_points[3], lst_points[4]);
Sphere sphere1 = Sphere4Points(lst_points[0], lst_points[2], lst_points[3], lst_points[4]);
Sphere sphere2 = Sphere4Points(lst_points[0], lst_points[1], lst_points[3], lst_points[4]);
Sphere sphere3 = Sphere4Points(lst_points[0], lst_points[1], lst_points[2], lst_points[4]);
Sphere sphere4 = Sphere4Points(lst_points[0], lst_points[1], lst_points[2], lst_points[3]);
//Spheres that contains all 4 points
List<Sphere> lst_spheres = new List<Sphere>();
if (Contains(sphere0, lst_points[0])) lst_spheres.Add(sphere0);
if (Contains(sphere1, lst_points[1])) lst_spheres.Add(sphere1);
if (Contains(sphere2, lst_points[2])) lst_spheres.Add(sphere2);
if (Contains(sphere3, lst_points[3])) lst_spheres.Add(sphere3);
if (Contains(sphere4, lst_points[4])) lst_spheres.Add(sphere4);
//We keep the minimum sphere
if (lst_spheres.Count > 0)
{
sphere = lst_spheres[0];
for (int i = 1; i < lst_spheres.Count; i++)
{
if (lst_spheres[i].Radius < sphere.Radius)
{
sphere = lst_spheres[i];
}
}
}
return sphere;
}
/// <summary>
/// Class to compute the minimum bounding sphere with Welzl method
/// BEWARE if more than 300 points, there is a bug somewhere
/// https://gamedev.stackexchange.com/questions/162731/welzl-algorithm-to-find-the-smallest-bounding-sphere
/// </summary>
public static class MinimumBoundingSphereWelzl
{
public static List<Point3d> _boundaryPoints = null;
static Random random = new Random();
/// <summary>
/// Main method to compute the Minimum Bounding Sphere using Welzl method
/// </summary>
/// <param name="points">List of points</param>
/// <returns>Minimum Bounding Sphere</returns>
public static Sphere GetMinimumBoundingSphere(List<Point3d> points)
{
if (_boundaryPoints == null)
{
_boundaryPoints = new List<Point3d>(4);
}
else
{
_boundaryPoints.Clear();
}
return BallWithBounds(points);
}
/// <summary>
/// Recursive method to calculate a minimum Bounding Sphere
/// </summary>
/// <param name="contained"></param>
/// <returns></returns>
private static Sphere BallWithBounds(List<Point3d> contained)
{
if (contained.Count == 0 || _boundaryPoints.Count == 4)
{
switch (_boundaryPoints.Count)
{
case 0: return Sphere.Unset;
case 1: return new Sphere(_boundaryPoints[0], 0);
case 2: return Sphere2Points(_boundaryPoints[0], _boundaryPoints[1]);
case 3: return Sphere3Points(_boundaryPoints[0], _boundaryPoints[1], _boundaryPoints[2]);
case 4: return Sphere4Points(_boundaryPoints[0], _boundaryPoints[1], _boundaryPoints[2], _boundaryPoints[3]);
}
}
int last = contained.Count - 1;
int removeAt = random.Next(0, contained.Count);
Point3d removed = contained[removeAt];
contained[removeAt] = contained[last];
contained.RemoveAt(last);
var ball = BallWithBounds(contained);
if(!Contains(ball, removed))
{
_boundaryPoints.Add(removed);
ball = BallWithBounds(contained);
_boundaryPoints.RemoveAt(_boundaryPoints.Count - 1);
}
contained.Add(removed);
return ball;
}
}
/// <summary>
/// Test if a point is strictly inside a sphere and not on the sphere
/// </summary>
/// <param name="sphere">Sphere</param>
/// <param name="point">Point</param>
/// <param name="epsilon">Precision or tolerance</param>
/// <returns>Return true if point is inside</returns>
public static bool ContainsStrict(Sphere sphere, Point3d point, double epsilon)
{
bool output = false;
if ((sphere.Radius - sphere.Center.DistanceTo(point)) >= epsilon)
{
output = true;
}
return output;
}
/// <summary>
/// Return true if point is on or inside the sphere
/// </summary>
/// <param name="sphere">Sphere</param>
/// <param name="point">Point</param>
/// <returns>Return true if point is inside or on the sphere</returns>
public static bool Contains(Sphere sphere, Point3d point)
{
bool output = false;
if (sphere.Center.DistanceTo(point) <= sphere.Radius)
{
output = true;
}
return output;
}
/// <summary>
/// Minimum bounding sphere having 2 points
/// </summary>
/// <param name="a">First point</param>
/// <param name="b">Second Point</param>
/// <returns>Minimum bounding sphere</returns>
public static Sphere Sphere2Points(Point3d a, Point3d b)
{
Sphere sphere = Sphere.Unset;
double radius = a.DistanceTo(b) / 2;
if (radius > 1e-12)
{
sphere = new Sphere((a + b) / 2, radius);
}
return sphere;
}
/// <summary>
/// Minimum bounding sphere having 3 points
/// </summary>
/// <param name="a"></param>
/// <param name="b"></param>
/// <param name="c"></param>
/// <returns>Minimum bounding sphere</returns>
public static Sphere Sphere3Points(Point3d a, Point3d b, Point3d c)
{
Sphere sphere3 = GetCircumSphere(a, b, c);
Sphere sphere_ab = Sphere2Points(a, b);
Sphere sphere_ac = Sphere2Points(a, c);
Sphere sphere_bc = Sphere2Points(b, c);
Sphere sphereOK = sphere3;
if ((sphere_ab.IsValid) && sphere_ab.Center.DistanceTo(c) < sphereOK.Radius)
{
if ( sphere_ab.Radius < sphereOK.Radius)
{
sphereOK = sphere_ab;
}
}
if ((sphere_ac.IsValid) && sphere_ac.Center.DistanceTo(b) < sphereOK.Radius)
{
if ( sphere_ac.Radius < sphereOK.Radius)
{
sphereOK = sphere_ac;
}
}
if ((sphere_bc.IsValid) && sphere_bc.Center.DistanceTo(a) < sphereOK.Radius)
{
if ( sphere_bc.Radius < sphereOK.Radius)
{
sphereOK = sphere_bc;
}
}
return sphereOK;
}
/// <summary>
/// Minimum bounding sphere having 4 points
/// </summary>
/// <param name="a"></param>
/// <param name="b"></param>
/// <param name="c"></param>
/// <param name="d"></param>
/// <returns>Minimum bounding sphere</returns>
public static Sphere Sphere4Points(Point3d a, Point3d b, Point3d c, Point3d d)
{
Sphere sphere4 = GetCircumSphere(a, b, c, d);
Sphere sphere_abc = Sphere3Points(a, b, c);
Sphere sphere_abd = Sphere3Points(a, b, d);
Sphere sphere_acd = Sphere3Points(a, c, d);
Sphere sphere_bcd = Sphere3Points(b, c, d);
Sphere sphereOK = sphere4;
if (sphere_abc.IsValid && sphere_abc.Center.DistanceTo(d) < sphereOK.Radius)
{
if ( sphere_abc.Radius < sphereOK.Radius)
{
sphereOK = sphere_abc;
}
}
if (sphere_abd.IsValid && sphere_abd.Center.DistanceTo(c) < sphereOK.Radius)
{
if ( sphere_abd.Radius < sphereOK.Radius)
{
sphereOK = sphere_abd;
}
}
if (sphere_acd.IsValid && sphere_acd.Center.DistanceTo(b) < sphereOK.Radius)
{
if ( sphere_acd.Radius < sphereOK.Radius)
{
sphereOK = sphere_acd;
}
}
if (sphere_bcd.IsValid && sphere_bcd.Center.DistanceTo(a) < sphereOK.Radius)
{
if ( sphere_bcd.Radius < sphereOK.Radius)
{
sphereOK = sphere_bcd;
}
}
return sphereOK;
}
/// <summary>
///Calculate the circumscribed sphere of 3 points, it is not the minimal bounding sphere
/// https://gamedev.stackexchange.com/questions/60630/how-do-i-find-the-circumcenter-of-a-triangle-in-3d
/// </summary>
/// <param name="a">First point</param>
/// <param name="b">Second point</param>
/// <param name="c">Third point</param>
/// <returns></returns>
public static Sphere GetCircumSphere(Point3d a, Point3d b, Point3d c)
{
Vector3d ac = c - a;
Vector3d ab = b - a;
Vector3d abXac = Vector3d.CrossProduct(ab, ac);
Sphere sphere = Sphere.Unset;
if (abXac.Length > 1e-12)
{
// this is the vector from a TO the circumsphere center
Vector3d toCircumsphereCenter = (Vector3d.CrossProduct(abXac, ab) * ac.SquareLength + Vector3d.CrossProduct(ac, abXac) * ab.SquareLength) / (2.0 * abXac.SquareLength);
double circumsphereRadius = toCircumsphereCenter.Length;
sphere = new Sphere(a + toCircumsphereCenter, circumsphereRadius);
}
else
{
Sphere sphere_ab = new Sphere((a + b) / 2, a.DistanceTo(b) / 2);
Sphere sphere_ac = new Sphere((a + c) / 2, a.DistanceTo(c) / 2);
Sphere sphere_bc = new Sphere((b + c) / 2, b.DistanceTo(c) / 2);
sphere = sphere_ab;
if (sphere_ac.Radius > sphere.Radius) sphere = sphere_ac;
if (sphere_bc.Radius > sphere.Radius) sphere = sphere_bc;
}
return sphere;
}
/// <summary>
///Calculate the circumscribed sphere of 4 points, it is not the minimal bounding sphere
///
/// https://mathworld.wolfram.com/Circumsphere.html
/// </summary>
/// <param name="p1">First point</param>
/// <param name="p2">Second Point</param>
/// <param name="p3">Third Point</param>
/// <param name="p4">Fourh Point</param>
/// <returns>Minimum bounding sphere</returns>
public static Sphere GetCircumSphere(Point3d p1, Point3d p2, Point3d p3, Point3d p4)
{
//TODO of point aligned or in a plane =>
Sphere sphere = Sphere.Unset;
Point3d center = Point3d.Origin;
double radius = 1;
//Matrix a = new Matrix(4, 4);
Transform a = new Transform(1);
//row colum,)
a[0, 0] = p1.X;
a[0, 1] = p1.Y;
a[0, 2] = p1.Z;
a[0, 3] = 1;
a[1, 0] = p2.X;
a[1, 1] = p2.Y;
a[1, 2] = p2.Z;
a[1, 3] = 1;
a[2, 0] = p3.X;
a[2, 1] = p3.Y;
a[2, 2] = p3.Z;
a[2, 3] = 1;
a[3, 0] = p4.X;
a[3, 1] = p4.Y;
a[3, 2] = p4.Z;
a[3, 3] = 1;
//Matrix dx = new Matrix(4, 4);
Transform dx = new Transform(1);
dx[0, 0] = p1.X * p1.X + p1.Y * p1.Y + p1.Z * p1.Z;
dx[0, 1] = p1.Y;
dx[0, 2] = p1.Z;
dx[0, 3] = 1;
dx[1, 0] = p2.X * p2.X + p2.Y * p2.Y + p2.Z * p2.Z;
dx[1, 1] = p2.Y;
dx[1, 2] = p2.Z;
dx[1, 3] = 1;
dx[2, 0] = p3.X * p3.X + p3.Y * p3.Y + p3.Z * p3.Z;
dx[2, 1] = p3.Y;
dx[2, 2] = p3.Z;
dx[2, 3] = 1;
dx[3, 0] = p4.X * p4.X + p4.Y * p4.Y + p4.Z * p4.Z;
dx[3, 1] = p4.Y;
dx[3, 2] = p4.Z;
dx[3, 3] = 1;
Transform dy = new Transform(1);
dy[0, 0] = p1.X * p1.X + p1.Y * p1.Y + p1.Z * p1.Z;
dy[0, 1] = p1.X;
dy[0, 2] = p1.Z;
dy[0, 3] = 1;
dy[1, 0] = p2.X * p2.X + p2.Y * p2.Y + p2.Z * p2.Z;
dy[1, 1] = p2.X;
dy[1, 2] = p2.Z;
dy[1, 3] = 1;
dy[2, 0] = p3.X * p3.X + p3.Y * p3.Y + p3.Z * p3.Z;
dy[2, 1] = p3.X;
dy[2, 2] = p3.Z;
dy[2, 3] = 1;
dy[3, 0] = p4.X * p4.X + p4.Y * p4.Y + p4.Z * p4.Z;
dy[3, 1] = p4.X;
dy[3, 2] = p4.Z;
dy[3, 3] = 1;
Transform dz = new Transform(1);
dz[0, 0] = p1.X * p1.X + p1.Y * p1.Y + p1.Z * p1.Z;
dz[0, 1] = p1.X;
dz[0, 2] = p1.Y;
dz[0, 3] = 1;
dz[1, 0] = p2.X * p2.X + p2.Y * p2.Y + p2.Z * p2.Z;
dz[1, 1] = p2.X;
dz[1, 2] = p2.Y;
dz[1, 3] = 1;
dz[2, 0] = p3.X * p3.X + p3.Y * p3.Y + p3.Z * p3.Z;
dz[2, 1] = p3.X;
dz[2, 2] = p3.Y;
dz[2, 3] = 1;
dz[3, 0] = p4.X * p4.X + p4.Y * p4.Y + p4.Z * p4.Z;
dz[3, 1] = p4.X;
dz[3, 2] = p4.Y;
dz[3, 3] = 1;
Transform c = new Transform(1);
c[0, 0] = p1.X * p1.X + p1.Y * p1.Y + p1.Z * p1.Z;
c[0, 1] = p1.X;
c[0, 2] = p1.Y;
c[0, 3] = p1.Z;
c[1, 0] = p2.X * p2.X + p2.Y * p2.Y + p2.Z * p2.Z;
c[1, 1] = p2.X;
c[1, 2] = p2.Y;
c[1, 3] = p2.Z;
c[2, 0] = p3.X * p3.X + p3.Y * p3.Y + p3.Z * p3.Z;
c[2, 1] = p3.X;
c[2, 2] = p3.Y;
c[2, 3] = p3.Z;
c[3, 0] = p4.X * p4.X + p4.Y * p4.Y + p4.Z * p4.Z;
c[3, 1] = p4.X;
c[3, 2] = p4.Y;
c[3, 3] = p4.Z;
center = new Point3d(dx.Determinant / (2 * a.Determinant), -dy.Determinant / (2 * a.Determinant), dz.Determinant / (2 * a.Determinant));
double dx2 = dx.Determinant * dx.Determinant;
double dy2 = dy.Determinant * dy.Determinant;
double dz2 = dz.Determinant * dz.Determinant;
double ac = a.Determinant * c.Determinant;
radius = Math.Sqrt(dx2 + dy2 + dz2 - 4 * ac) / (2 * Math.Abs(a.Determinant));
return new Sphere(center, radius);
}
