   WelzlSphere algorithm.  How do I create a sphere given 3 points?
  I'm using the following algorithm to create a minimum bounding sphere for a set of points.  It is written in C# and based off a method found in Christer Ericson's awesome orange book:   private static CD_BoundingSphere WelzlSphere(Vector3[] points, int numPts, Vector3[] sos, int numSos)   {  // If no input points, the recursion has bottomed out.   // Now compute an exact sphere based on points in Set Of Support (zero through four points)  if (numPts == 0)  { switch (numSos) {case 0: return new CD_BoundingSphere();   //    Radius = 0.0fcase 1: return new CD_BoundingSphere(sos[0]);  // Centre = p   Radius = 0.0fcase 2: return new CD_BoundingSphere(sos[0], sos[1]);case 3: return new CD_BoundingSphere(sos[0], sos[1], sos[2]);case 4: return new CD_BoundingSphere(sos[0], sos[1], sos[2], sos[3]); }  }    // Pick a point at "random" (here just the last point of the input set)  int index = numPts - 1;    // Recursively compute the smallest bounding sphere of the remaining points  CD_BoundingSphere smallestSphere = WelzlSphere(points, numPts - 1, sos, numSos);    // If the selected point lies inside this sphere, it is indeed the smallest  if (PointInsideSphere(points[index], smallestSphere))  { return smallestSphere;  }  // Otherwise, update set of support to additionally contain the new point  sos[numSos] = points[index];  // Recursively compute the smallest sphere of remaining points with new s.o.s.  return WelzlSphere(points, numPts - 1, sos, numSos + 1);   }   public static bool PointInsideSphere(Vector3 point, CD_BoundingSphere sphere)   {  // Relative position of point from sphere centre  Vector3 d = point - sphere.Centroid;  return Vector3.Dot(d, d) <= sphere.Radius * sphere.Radius;   }The problem is that I can create a sphere using 1, 2 and 4 points but how can I create a sphere from 3 points?   Do I need to create 3 planes and find out where they intersect to find the sphere's centre?Also, is the following code: (2 * a1 * (c2 * d3 + d2 * b3 + b2 * c3 - d2 * c3 - b2 * d3 - c2 * b3) + b1 * (d2 * c3 + c2 * a3 + a2 * d3 - c2 * d3 - a2 * c3 - d2 * a3) + c1 * (a2 * b3 + b2 * d3 + d2 * a3 - b2 * a3 - d2 * b3 - a2 * d3) + d1 * (b2 * a3 + a2 * c3 + c2 * b3 - a2 * b3 - c2 * a3 - b2 * c3));the same as: (2(a1)(c2d3 + d2b3 + b2c3 - d2c3 - b2d3 - c2b3) +   (b1)(d2c3 + c2a3 + a2d3 - c2d3 - a2c3 - d2a3) +   (c1)(a2b3 + b2d3 + d2a3 - b2a3 - d2b3 - a2d3) +   (d1)(b2a3 + a2c3 + c2b3 - a2b3 - c2a3 - b2c3))The brackets confused me a little, so I wasn't completely sure.   
  The minimum sphere passing through 3 given points (presumably what you are asking) can be computed by looking at the problem on the plane that contains the three points. Simply compute the center of the circle that contains those three points, and that's the center of the sphere you are looking for.You can also look at it as the intersection of the plane that contains the three points with the three planes of equidistance to each pair of points.;  Ok so I've calculated the plane on which the 3 points lie.  How can I find the further planes in order find the point of intersection?I'm not sure how I can remove a dimension otherwise to use Vector3's to calculate Vector2 lines for a 2D circle.   ;Quote:Original post by Spa8nkyOk so I've calculated the plane on which the 3 points lie.  How can I find the further planes in order find the point of intersection?Plane plane_of_equidistant_points(Point A, Point B) {  Point P = (A+B)/2;  Vector normal = A-B;  return plane_from_point_and_normal(P, normal);};  Thanks alvaro, the sphere from 3 points constructor now looks like the following:   public CD_BoundingSphere(Vector3 p0, Vector3 p1, Vector3 p2)   {  // Calculate plane on which all 3 points lie  CD_Plane plane0 = new CD_Plane(p0, p1, p2);  // Calculate two planes from the 3 points  CD_Plane plane1 = new CD_Plane(p0 - p1, (p0 + p1) * 0.5f);  CD_Plane plane2 = new CD_Plane(p1 - p2, (p1 + p2) * 0.5f);  // Sphere centre is the point of intersection for all 3 planes  centre = CD_Plane.IntersectionOf3Planes(plane0, plane1, plane2);  position = centre;  // Calculate squared radius  Vector3 d = centre - p0;  radius = Vector3.Dot(d, d);   }Is there a way of working out the maximum number for the set of support points given a set of points from which the minimum bounding sphere is calculated?  Is it just the same as the number of points?   public static CD_BoundingSphere WelzlSphere(Vector3[] points)   {  int numPoints = points.Length;  Vector3[] supportPoints = new Vector3[numPoints];  // Number of support points is intialised at 0  return WelzlSphere(points, numPoints, supportPoints, 0);   };Quote:Original post by Spa8nkyIs there a way of working out the maximum number for the set of support points given a set of points from which the minimum bounding sphere is calculated?Did you mean "is to be calculated"?Quote:Is it just the same as the number of points?No, you cannot know beforehand how many points the set of support points will contain. You just know that it will contain at most 4 points.As far as I remember, you start by computing a random permutation of the set of points. You initialize your set of support points with the first point of the permutation and you compute the minimum sphere containing it.Then,Real minimum_radius = infinite;Sphere min_spheres[number_of_spheres];int index_of_the_smallest_sphere = -1;int count_of_sphere = -1;for each point P:    if P is not contained in the current_minimum sphere   for each support point SP  sphere = minimum_sphere( SP, P )  increment count_of_sphere  min_spheres[count_of_sphere] = sphere  if sphere contains all the remaining support points if sphere.radius < minimum_radiusminimum_radius = sphere.radiusindex_of_the_smallest_sphere = count_of_sphere   for each combination of support points ( SP1, SP2 )  sphere = compute_minimum_sphere( SP1, SP2, P )  increment count_of_sphere  min_spheres[count_of_sphere] = sphere  if sphere contains all the remaining support points if sphere.radius < minimum_radiusminimum_radius = sphere.radius;index_of_the_smallest_sphere = count_of_sphere   for each combination of support points ( SP1, SP2, SP3 )  sphere = compute_minimum_sphere( SP1, SP2, SP3, P )  increment count_of_sphere  min_spheres[count_of_sphere] = sphere  if sphere contains all the remaining support points if sphere.radius < minimum_radiusminimum_radius = sphere.radiusindex_of_the_smallest_sphere = count_of_sphereswitch( index_of_the_smallest_sphere ){    case 0:    .    .    case i:   you know how min_spheres[index_of_the_smallest_sphere] has been   computed if you hard-coded all the possible combinations(good job!!!).   You can because you know depending on the number of support points   (n_support_points) how to combine n_support_points - 1   (think of "n chose k")   so you know what point(s) to remove from the support point   remove them and add point P    .    .}You will only compute minimum spheres of combinations of three supporting points plus the point not contained if the set of supporting points contains four points.You will only compute minimum spheres of combination of two supporting points plus the point not contained if the set of supporting points contains at least three points.So you'll have updating routines depending on the number of supporting points when a point is not contained in the minimum sphere built with the supporting points.P.S: Dave Eberly has an implementation in his site geometrictools.comIt will save you time.[Edited by - johnstanp on May 25, 2010 1:12:27 PM];Quote:P.S: Dave Eberly has an implementation in his site geometrictools.comIt will save you time.Thank you.  I think I've cracked it though:   public static CD_BoundingSphere WelzlSphere(Vector3[] points)   {  // Maximum of 4 support points  Vector3[] supportPoints = new Vector3[4];  // Number of support points is intialised at 0  CD_BoundingSphere s = WelzlSphere(points, points.Length, supportPoints, 0);    // Now that we have a sphere, calculate the real radius of the sphere using an expensive square root  s.Radius = (float)Math.Sqrt(s.Radius);  return s;   }   private static CD_BoundingSphere WelzlSphere(Vector3[] points, int numPts, Vector3[] sos, int numSos)   {  // If no input points, the recursion has bottomed out.   // Now compute an exact sphere based on points in Set Of Support (zero through four points)  if (numPts == 0)  { switch (numSos) {case 0: return new CD_BoundingSphere();   //    Radius = 0.0fcase 1: return new CD_BoundingSphere(sos[0]);  // Centre = p   Radius = 0.0fcase 2: return new CD_BoundingSphere(sos[0], sos[1]);case 3: return new CD_BoundingSphere(sos[0], sos[1], sos[2]);case 4: return new CD_BoundingSphere(sos[0], sos[1], sos[2], sos[3]); }  }  // Pick a point at "random" (here just the last point of the input set)  int index = numPts - 1;  // Recursively compute the smallest bounding sphere of the remaining points  CD_BoundingSphere smallestSphere = WelzlSphere(points, numPts - 1, sos, numSos);  // If the selected point lies inside this sphere, it is indeed the smallest  // Sphere radius is already squared so no need to square it here, just square the distance  if (Vector3.DistanceSquared(points[index], smallestSphere.Centroid) < smallestSphere.Radius)  { return smallestSphere;  }  // Otherwise, update set of support to additionally contain the new point  sos[numSos] = points[index];  // Recursively compute the smallest sphere of remaining points with new s.o.s.  return WelzlSphere(points, numPts - 1, sos, numSos + 1);   }I haven't posted the sphere constructors but everything else is there.  It appears to work absolutely fine.  If you can see any potential pitfalls or gremlins then please shout.EDIT: It seems that numSos can be greater than 4!  Why is this?   ;Quote:Original post by Spa8nkyQuote:P.S: Dave Eberly has an implementation in his site geometrictools.comIt will save you time.Thank you.  I think I've cracked it though:*** Source Snippet Removed ***I haven't posted the sphere constructors but everything else is there.  It appears to work absolutely fine.  If you can see any potential pitfalls or gremlins then please shout.EDIT: It seems that numSos can be greater than 4!  Why is this?Well, with 4 points, you can find the equation of a sphere: adding points located on the surface of that sphere is redundant. You could afford such redundance (to a certain precision) by including in your support set, points that are very close to each other, or very close to the sphere.You could have in your set, Math::abs( squared_length( substract( points[5], points[10] ) ) ) <= fepsMath::abs( squared_length( substract( points[5], points[7] ) ) ) <= fepsMath::abs( squared_length( substract( points[11], sphere.get_center() ) ) - Math::pow( sphere.get_radius(), 2 ) ) <= fepspoints[5], points[10] and points[7] are so close( you define how much ), that you can safely assume that they are equivalent.points[11] is so close to the sphere defined by the other support points,that you can assume it lies on the sphere.P.S: I'll take a close look to your algorithm later. It seems very interesting.I didn't do it and I am leaving.[Edited by - johnstanp on May 26, 2010 12:32:31 AM];  Thanks John.  I'm pm'ing you the sphere constructors as it is a lot of code to post here.The methods for creating the spheres from 0 to 4 points all work perfectly with random points.  The problem lies with the algorithm finding more than 4 support points and causing a stack overflow error in the process.   ;  Another method I, which I am inclined to stick with as it is less error prone, is as follows:- Find the minimum AABB of all the model vertices- Create a minimum bounding sphere using the AABB position and SQRT3 * the maximum AABB extent.I'm finding the minimum AABB using this method:   public static CD_AABB CreateFromPoints(Vector3[] points)   {  int numPoints = points.Length;  if (numPoints == 0)  { return null;  }  Vector3 min = new Vector3(float.MaxValue);  Vector3 max = new Vector3(float.MinValue);  // No need to project points onto an axis as the values are axis aligned  for (int i = 0; i < numPoints; ++i)  { Vector3 point = points; Vector3.Min(ref min, ref point, out min); Vector3.Max(ref max, ref point, out max);  }  Vector3 centre = (min + max) * 0.5f;  Vector3 dim = max - min;  return new CD_AABB(centre, dim.X, dim.Y, dim.Z);   }I presume the dimension calculation is "max - min"?  If you can see any problems with the algorithm or improvements to be made, please shout.   