Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Problem: Smallest sphere including all 3D points

• To: mathgroup at smc.vnet.net
• Subject: [mg16621] Re: Problem: Smallest sphere including all 3D points
• From: Martin Kraus <mkraus at theorie3.physik.uni-erlangen.de>
• Date: Fri, 19 Mar 1999 12:53:48 -0500
• Organization: Regionales Rechenzentrum Erlangen, Germany
• References: <7cd869$o91@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

Barthelet, Luc wrote:
> 
> 
>  I was just asked by a friend how to find the smallest sphere that
>  would include all points from a set of 3D points.
>  It feels like finding a 3D convex hull and then finding the best
>  sphere (??)
>  I would appreciate any complete solution or best set of pointers...
> 
>         Thank you,
> 
>         Luc Barthelet
>         General manager, Maxis
>         http://www.simcity.com <http://www.simcity.com>  (we are #1!)
> 
> 

For any practical purpose your suggestions (realized with a
numerical minimization) is probably the best.
However, we found an algorithm for an "exact" solution.
The problem is in fact almost the same in two dimensions;
therefore, I will describe the algorithm in two dimensions first:

1) Take all pairs of points and calculate the minimal circle
covering each pair. (The radius of such a circle is just half
the distance between the pair of points, the center is the
midpoint between them.)

2) Take all triples of points and calculate the minimal circle
covering each triple. (That is just the circum circle.)

3) Throw away all circles which do not cover all points.

4) From the remaining circles take the one with the minimal
radius.

This algorithm guarantees to give the correct result. 
(Problems might arise from numerics, however.)

Here is a simple implementation in Mathematica:


<<DiscreteMath`Combinatorica`;

includesAllPoints[{radius_,center_},points_,epsilon_:10^-8]:=
  Module[{square=radius^2},
    Apply[And,Map[((#-center).(#-center)<=square+epsilon)&,points]]]

circumCircle[{p1_,p2_}]:=Module[{},{Sqrt[(p1-p2).(p1-p2)]/2.,(p1+p2)/2.}]

circumCircle[{p1_,p2_,p3_}]:=
 
Module[{n1=p1-p2,b1=N[(p1+p2)/2],n2=p3-p2,b2=N[(p3+p2)/2],nn1,nn2,x,y},
    nn1=N[n1/Sqrt[n1.n1]];
    nn2=N[n2/Sqrt[n2.n2]];{x,
       
y}={x,y}/.First[Solve[{({x,y}-b1).n1==0,({x,y}-b2).n2==0},{x,y}]];{
      Sqrt[(p2-{x,y}).(p2-{x,y})],{x,y}}]

minimalCircumCircle[points_]:=
  Module[{pairs,triangles,circles,coveringCircles},
    pairs=KSubsets[N[points],2];
			triangles=KSubsets[N[points],3];
    circles=Map[circumCircle,Join[triangles,pairs]];
   
coveringCircles=DeleteCases[circles,c_/;(!includesAllPoints[c,points]),1];
    If[Length[coveringCircles]>0,Return[First[Sort[coveringCircles]]],
      Print["Found no solution. (This is an error.)"];Abort[]]]



minimalCircumCircle[<list of points>] will return
the wanted circle in the form {<radius>, <center>}.
We can check the result with a simple function:

maxDistance[point_,points_]:=Max[Map[Sqrt[(point-#).(point-#)]&,points]]

maxDistance[<center>, <list of points>] should give the 
returned radius.

Warnings arise from identical points or three points on one line.
However, these warnings can be ignored.

The algorithm for three dimensions needs a small extension
but the basic idea is the same:

1) Take all pairs of points and calculate the minimal sphere
covering each pair. (The radius of such a sphere is just half
the distance between the pair of points, the center is the
midpoint between them.)

2) Take all triples of points and calculate the minimal sphere
covering each triple. (That is just the three dimensional 
continuation of the circum circle of the flat polygon.)

3) Take all quadruples of points and calculate the minimal sphere
covering each quadruple. (That is just the circum sphere.)

4) Throw away all spheres which do not cover all points.

5) From the remaining spheres take the one with the minimal
radius.


The implementation is:

<<DiscreteMath`Combinatorica`;

includesAllPoints[{radius_,center_},points_,epsilon_:10^-8]:=
  Module[{square=radius^2},
    Apply[And,Map[((#-center).(#-center)<=square+epsilon)&,points]]]

circumSphere[{p1_,p2_}]:=Module[{},{Sqrt[(p1-p2).(p1-p2)]/2.,(p1+p2)/2.}]

circumSphere[{p1_,p2_,p3_}]:=
 
Module[{n1=p1-p2,b1=N[(p1+p2)/2],n2=p3-p2,b2=N[(p3+p2)/2],n3,b3,nn1,nn2,nn3,
      x,y,z},n3=Cross[n1,n2];b3=N[p2];nn1=N[n1/Sqrt[n1.n1]];
    nn2=N[n2/Sqrt[n2.n2]];nn3=N[n3/Sqrt[n3.n3]];x;y;
    z;{x,y,z}={x,y,z}/.First[
         
Solve[{({x,y,z}-b1).n1==0,({x,y,z}-b2).n2==0,({x,y,z}-b3).n3==0},{x,
              y,z}]];{Sqrt[(p2-{x,y,z}).(p2-{x,y,z})],{x,y,z}}]

circumSphere[{p1_,p2_,p3_,p4_}]:=
  Module[{x,y,
      z},{x,y,z}={x,y,z}/.First[
          Solve[{({x,y,z}-p1).({x,y,z}-p1)==({x,y,z}-p2).({x,y,z}-
                     
p2),({x,y,z}-p1).({x,y,z}-p1)==({x,y,z}-p3).({x,y,z}-
                     
p3),({x,y,z}-p1).({x,y,z}-p1)==({x,y,z}-p4).({x,y,z}-
                      p4)}]];{Sqrt[({x,y,z}-p1).({x,y,z}-p1)],{x,y,z}}]

minimalCircumSphere[points_]:=
  Module[{pairs,triangles,tetrahedra,spheres,coveringSpheres},
    pairs=KSubsets[N[points],2];
			triangles=KSubsets[N[points],3];tetrahedra=KSubsets[N[points],4];
    spheres=Map[circumSphere,Join[tetrahedra,triangles,pairs]];
   
coveringSpheres=DeleteCases[spheres,s_/;(!includesAllPoints[s,points]),1];
    If[Length[coveringSpheres]>0,Return[First[Sort[coveringSpheres]]],
      Print["Found no solution. (This is an error.)"];Abort[]]]


Checks of the result can be done as above. Note that both 
algorithms need extremely long for more than just a few (about 9)
points.
Therefore, this algorithm has probably no real world applications.

(This solution was created yesterday evening by Oliver Jahn, 
Martin Kraus, Verena Schoen and Christian Winkler; Mathematica 
implementation by Martin Kraus. No optimization was done; better 
algorithms and implementations are surely available.)

Greetings

Martin Kraus