RE: Problem: Smallest sphere including all 3D points
- To: mathgroup at smc.vnet.net
- Subject: [mg16638] RE: Problem: Smallest sphere including all 3D points
- From: "Barthelet, Luc" <lucb at ea.com>
- Date: Fri, 19 Mar 1999 12:53:56 -0500
- Sender: owner-wri-mathgroup at wolfram.com
It works. Thanks.
I sped it up slightly by testing for pairs, then triangles then tetra...
minimalCircumSphere[points_] :=
Module[{pairs, triangles, tetrahedra, spheres, coveringSpheres},
pairs = KSubsets[N[points], 2];
spheres = Map[circumSphere, pairs];
coveringSpheres =
DeleteCases[spheres, s_ /; (! includesAllPoints[s, points]), 1];
If[Length[coveringSpheres] > 0, Return[First[Sort[coveringSpheres]]],
Print["Found no solution for pairs. "];];
triangles = KSubsets[N[points], 3];
spheres = Map[circumSphere, triangles];
coveringSpheres =
DeleteCases[spheres, s_ /; (! includesAllPoints[s, points]), 1];
If[Length[coveringSpheres] > 0, Return[First[Sort[coveringSpheres]]],
Print["Found no solution for triangles. "];];
tetrahedra = KSubsets[N[points], 4];
spheres = Map[circumSphere, tetrahedra];
coveringSpheres =
DeleteCases[spheres, s_ /; (! includesAllPoints[s, points]), 1];
If[Length[coveringSpheres] > 0, Return[First[Sort[coveringSpheres]]],
Print["Found no solution for tetrahedra. This is a bug "];]
]
it could be sped up further. you only need to test one sphere in each group
for covering, the largest.
with our data set, so far all solutions have been found in triangles.
This along with the convexHull3D from wouter fixed the problem. The
solutions we tested were 3 percent better than the approximation we were
using, we will probably switch to that method.
thanks to all of you for your help. (I receive about 20 emails on this from
people suggesting solutions, cool).
Luc Barthelet
General Manager, Maxis
http://www.simcity.com (SimCity 3000 is #1 Game on the PC in february
worldwide!)
-----Original Message-----
From: Martin Kraus [mailto:mkraus at theorie3.physik.uni-erlangen.de]
To: mathgroup at smc.vnet.net
Subject: [mg16638] Re: Problem: Smallest sphere including all 3D points
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