[comp.graphics] Ray Traced Bounding Spheres

ritter@versatc.UUCP (Jack Ritter) (05/11/89)

Given a cluster of points in 3 space, is there
a good method for finding the minumum radius
sphere which encloses all the points?  If not
minumum, at least "small"?  Certainly it should
be tighter than the sphere which encloses the
minimum bounding box.

I have a feeling the solution is iterative. If
so, I could provide a good initial guess for the
center & radius.

-- 
	      ->  S C R E W    E X X O N    <-
   Jack Ritter, S/W Eng. Versatec, 2710 Walsh Av, Santa Clara, CA 95051
   Mail Stop 1-7.  (408)982-4332, or (408)988-2800 X 5743
   UUCP: {pyramid,mips,vsi1,arisia}!versatc!ritter

ritter@versatc.UUCP (Jack Ritter) (05/11/89)

In article <17241@versatc.UUCP>, ritter@versatc.UUCP (Jack Ritter) writes:
> Given a cluster of points in 3 space, is there
> a good method for finding the minumum radius
> sphere which encloses all the points? 

Several people have responded to this,
so I'll post this. 
Wolfgang (sohrt@cs.utah.edu) suggested an
algorithm, which I dont think was quite right,
but it was good enough to trigger me to come
up with this algorithm:

Take 1st 2 pts, center the initial sphere at mid pt M of the
2 pts, with radius = dist to either.
For each new pt P, if it's outside current sphere,
the new minimum sphere's center will be the
the point midway between P and the intersection of the
line P->M with the old sphere (on side OPPOSITE P). 
This will make an ossculating
sphere that just kisses new pt and old sphere
on opposite sides.
Thanks to Wolfie.
-- 
	      ->  S C R E W    E X X O N    <-
   Jack Ritter, S/W Eng. Versatec, 2710 Walsh Av, Santa Clara, CA 95051
   Mail Stop 1-7.  (408)982-4332, or (408)988-2800 X 5743
   UUCP:  {ames,apple,sun,pyramid}!versatc!ritter

flynn@pixel.cps.msu.edu (Patrick J. Flynn) (05/11/89)

In article <17241@versatc.UUCP> ritter@versatc.UUCP (Jack Ritter) writes:
>Given a cluster of points in 3 space, is there
>a good method for finding the minumum radius
>sphere which encloses all the points?  If not
>minumum, at least "small"?  Certainly it should
>be tighter than the sphere which encloses the
>minimum bounding box.

This is a 3D generalization of the `smallest enclosing circle' problem, which
is covered in Preparata and Shamos' book on computational geometry (see pp.
248ff.). In 2D, you can find a smallest bounding circle for n points in
n log n time.  I don't know how that scales to 3D.

>I have a feeling the solution is iterative.

Hmm, it shouldn't be.  It's a discrete problem.  See P&S for a discussion
of the (mis-)formulation of the bounding-circle problem as a continuous problem.
------------------------------------------------------------------------------
Pat Flynn, CS, Mich. State U | "In your heart, you know it's biscuit shaped."
flynn@cps.msu.edu            |                                    -Me
(517) 353-4638               |

jep@oink.UUCP (James E. Prior) (05/12/89)

In article <17241@versatc.UUCP> ritter@versatc.UUCP (Jack Ritter) writes:
>Given a cluster of points in 3 space, is there
>a good method for finding the minumum radius
>sphere which encloses all the points?

I'll try my hand at it:
1) Figure out which pair of points are furthest away from each other.
	(This need not be unique, i.e. this algorithm will work even
	if 2 or 3 of the points are at the same location.)

2) Define Sphere A as having its center halfway between the furthest points,
	and its radius as half the distance between the furthest points.

If the other point lies inside or on Sphere A, then Sphere A is the solution.

Otherwise, the other point lies outside Sphere A, 
and the solution is Sphere B defined as follows:

Sphere B has the same center and radius as the circle defined by the
three points.

It feels right, but I don't guarantee it.  This is the result of my
own mental hacking, and not a reference book.  
-- 
Jim Prior    jep@oink    osu-cis!n8emr!oink!jep    N8KSM

jep@oink.UUCP (James E. Prior) (05/12/89)

In article <
7241@versatc.UUCP> ritter@versatc.UUCP (Jack Ritter) writes:
>Given a cluster of points in 3 space, is there
>a good method for finding the minumum radius
>sphere which encloses all the points?

I'll try my hand at it:
1) Figure out which pair of points are furthest away from each other.
	(This need not be unique, i.e. this algorithm will work even
	if 2 or 3 of the points are at the same location.)

2) Define Sphere A as having its center halfway between the furthest points,
	and its radius as half the distance between the furthest points.

If the other point lies inside or on Sphere A, then Sphere A is the solution.

Otherwise, the other point lies outside Sphere A, 
and the solution is Sphere B defined as follows:

Sphere B has the same center and radius as the circle defined by the
three points.

It feels right, but I don't guarantee it.  This is the result of my
own mental hacking, and not a reference book.  
-- 
Jim Prior    jep@oink    osu-cis!n8emr!oink!jep    N8KS

ECULHAM@UALTAVM.BITNET (Earl Culham) (05/12/89)

In article <17241@versatc.UUCP>, ritter@versatc.UUCP (Jack Ritter) writes:
 
>Given a cluster of points in 3 space, is there
>a good method for finding the minumum radius
>sphere which encloses all the points?  If not
>minumum, at least "small"?  Certainly it should
>be tighter than the sphere which encloses the
>minimum bounding box.
>
>I have a feeling the solution is iterative. If
>so, I could provide a good initial guess for the
>center & radius.
>
 
The solution is not iterative. A simple solution is available in
a two step process, and is characterised by whether three, or only two
of the points are on the surface of the optimal sphere.
 
Given the points, A, B, and C.
Searching for the center of the optimal sphere, P, with the smallest
radius, R.
 
Checking the midpoints.
 
set P = point halfway between A and B
set R = 1/2 of length from A to B
If length from C to P is less than or equal to R ---> done
Repeat with line A-C, and B-C if needed.
 
Extending the midpoints at right angles.
 
build the line which intersects A-B at a right angle through the midpoint
  of A-B, call the line D
build the line which intersects A-C at a right angle through the midpoint
  of A-C, call it E.
P is the point where D and E intersect. (Solve the simultaneous equations
R is the length A-P.
 

mplevine@phoenix.Princeton.EDU (Marshall Polimer Levine) (05/14/89)

In article <17241@versatc.UUCP>, ritter@versatc.UUCP (Jack Ritter) writes:
> Given a cluster of points in 3 space, is there
> a good method for finding the minumum radius
> sphere which encloses all the points? 


Well, Jack, I just so happened to write exactly the algorithm you are looking
for one week ago.  It is fully tested and FULLY LINEAR.  I have included the
PASCAL source.  If anyone has any questions, feel free to ask me.
I have included all testing and timing results after the source.
Here we go:


<================================ CUT HERE ==============================>

{ Written by Marshall Levine     5/2/89

mplevine@phoenix.princeton.edu

This program searches for a cluster of 10 stars with the smallest possible
radius (centered around one star) in a space full of points.

First, I implemented what I call a pri list.  This is a linked
list of real numbers (representing distances).  The list is
kept in order.  However, when a number that would go farther into
the list than the number of points per sphere (NUMINSPHERE) tries to go into
the list, the insert procedure stops it.  This is done because the distance
is useless.  For example, if NUMINSPHERE is 5 and a number would be inserted
into the 7th slot in the list, it is not inserted.  The minimum radius of a
sphere with 5 points would be determined by the 5th element of the list
(not including the header), so any number inserted after the 5th element
is useless and is therefore not inserted.  If there are not NUMINSPHERE
elements in the pri, then there were not enough points to fill the sphere.

The brute-force algorithm loops through every point in space.  For each
point, the algorithm finds the distance between that point and every other
point and puts that distance in the pri.  When all points have been compared
against this point, the NUMINSPHERE'th element is taken to be the minimum
radius of a sphere centered at this point containing NUMINSPHERE points.
However, points are not compared against themselves, so the exact number of
comparisons is N^2-N, making this an N^2 algorithm.

My efficient algorithm divides the space up into a 3-dimensional grid.  If
the grid were divided into NUMPOINTS/NUMINSPHERE sectors, then at least one
of the sectors would have at least NUMINSPHERE points.  Now, make spheres
with the same volume as the sectors.  At least one sphere surrounding one
point will have at least NUMINSPHERE points.  It turns out that the tradeoff
between fewer computations and more overhead is minimized by choosing the
grid to have enough sectors such that each sector is r/2 on each side (where
r is the radius of the aforementioned sphere).  My algorithm starts with a
sphere radius equal to the distance from one corner of the unit cube to
another (3^.5).  Given the first point in the list, we compare that point
against every other point in sectors touching the sphere (In this case,
every other point in space!)  By storing the distances and
then taking the NUMINSPHERE'th number from the pri list, as in
the brute algorithm, we probably reduce the size of the sphere.  Then, we
check the next point with the new, smaller sphere size and continue in this
way until we have tested every point.  As we go along, the minimum sphere
size keeps shrinking until for any given point, we only check a few
neighboring sectors, if any.  In practice, this radius shrinks so quickly
that my algorithm displays LINEAR BEHAVIOR!!!  See the second program and
its output for clarification.
}

program clusters(input,output);

const
  MAXNUMPOINTS = 501;    { The maximum # of points we can handle  }
  NUMINSPHERE = 10;      { # stars to find inside sphere          }
  INFINITY = 999999.9;   { Larger than largest distance possible  }
  MAXUSESPACE = 20;      { Maximum length per edge of space-grid  }
  PI = 3.1415926535;

type
  datatype = real;
  point = record         { The type of a point }
            x : real;
            y : real;
            z : real;
            data : datatype;
          end;
  ptr = ^node;
  node = record          { Linked list for a distances list called "pri" }
           data : real;
           next : ptr;
         end;
  sptr = ^spacenode;     { Linked list for each sector in the space-grid }
  spacenode = record
                index : integer; { Stores index of that point in points[] }
                next : sptr;
              end;


var
  rndnm : integer;       { Needed for the random number generator }
  points : array [1..MAXNUMPOINTS] of point;   { All points in space }
  listhead : ptr;        { List head for distances list called "pri" }
  space : array[0..MAXUSESPACE, 0..MAXUSESPACE, 0..MAXUSESPACE] of sptr;
                         { The space-grid (hereafter called 'grid') }
  spacesize, usespace : integer;  { Size per edge of grid }
  NUMPOINTS : integer;   { The number of points we have in space }


{ **************** Support routines for random generators ************** }

procedure seed;        { Seed the random number generator }
begin
  writeln('seed:');
  readln(rndnm);
end;

function rndom(scale : integer) : real; { Make random real from 0 to scale }
begin
  rndnm := abs(abs((rndnm*921+1)) mod 32749);
  rndom := (rndnm*scale/32749)
end;

procedure randompoint(var pt : point);  { Generate a random point within }
begin                                   {   a unit cube.                 }
  pt.x := rndom(1);
  pt.y := rndom(1);
  pt.z := rndom(1)
end;

procedure generatepoints;           { Generate NUMPOINTS points in space }
var x : integer;
begin
  for x := 1 to NUMPOINTS do
    randompoint(points[x])
end;


{ *************** Support routines for the "pri" list ******************** }

procedure initprilist;    { Initialize the pri list }
begin
  new(listhead);
  listhead^.data := 0.0;
  new(listhead^.next);
  listhead^.next^.data := INFINITY;
  listhead^.next^.next := nil
end;

procedure clearprilist;   { Clear the pri list }
var p,oldp : ptr;
begin
  p := listhead;
  while p <> nil do
  begin
    oldp := p;
    p := p^.next;
    dispose(oldp)
  end;
  new(listhead);
  listhead^.data := 0.0;
  new(listhead^.next);
  listhead^.next^.data := INFINITY;
  listhead^.next^.next := nil
end;


procedure insertprilist(r : real);  { Insert a distance into pri list    }
var p,oldp,temp : ptr;       { "pri" is just a linked list of distances  }
    x : integer;             { kept in low -> high order. The catch is   }
begin                        { that if a number should be inserted after }
  x := 1;                    { the NUMINSPHERE'th node, we don't bother  }
  p := listhead^.next;       { inserting it, because it isn't in the     }
  oldp := listhead;          { smallest sphere with NUMINSPHERE points.  }
  while (r > p^.data) and (x <= NUMINSPHERE) do
  begin
    oldp := p;
    p := p^.next;
    x := x + 1
  end;
  if x <= NUMINSPHERE then
  begin
    new(temp);
    temp^.data := r;
    temp^.next := p;
    oldp^.next := temp
  end
end;

function getbiggestinsphere : real;  { Returns value of the NUMINSPHERE'th }
var x : integer;                     { element in pri list, or INFINITY    }
    p : ptr;                         { if the list isn't that long.        }
begin
  x := 1;
  p := listhead^.next;
  while (x < NUMINSPHERE) and (p <> nil) do
  begin
    x := x + 1;
    p := p^.next
  end;
  if (x < NUMINSPHERE) or (p = nil) then getbiggestinsphere := INFINITY
  else getbiggestinsphere := p^.data
end;

procedure printprilist;              { Print the pri list, for debugging }
var p : ptr;
begin
  p := listhead;  { DO print the head }
  while p <> nil do
  begin
    writeln(p^.data:15:10);
    p := p^.next
  end;
  writeln('nil')
end;


{ ******************* Miscellaneous support routines ******************** }

procedure printpoint(pt : point);   { Print out a point }
begin
  writeln('(',pt.x:8:5,',',pt.y:8:5,',',pt.z:8:5,')')
end;


function cube(x : real) : real;     { Return cube root of a number }
begin
  cube := exp((1/3)*ln(x))
end;


{ *********************** Brute Force algorithm ************************* }

procedure brutecluster;    { Find minimum sphere containing NUMINSPHERE }
                           {   points by testing the distance between   }
                           {   every point.                             }
var distx,disty,distz,dist : real;      { Find distance between two points }
    bestsphere,trysphere : real;        { Find minimum sphere              }
    numcomps : integer;                 { # comparisons                    }
    thispoint,againstpoint : integer;   { Counters                         }
begin
  clearprilist;                           { Kill the priority list          }
  bestsphere := INFINITY;
  numcomps := 0;
  for thispoint := 1 to NUMPOINTS do      { Test every point...             }
  begin
    clearprilist;
    for againstpoint := 1 to NUMPOINTS do { ...against every other point    }
      if thispoint <> againstpoint then   { Don't compare point against self}
      begin
        distx := points[thispoint].x - points[againstpoint].x;
        disty := points[thispoint].y - points[againstpoint].y;
        distz := points[thispoint].z - points[againstpoint].z;
        dist := sqrt(distx*distx + disty*disty + distz*distz);
        numcomps := numcomps + 1;
        if dist < bestsphere then       { If dist less than smallest sphere,}
          insertprilist(dist)           {   insert distance into pri list   }
      end;
    trysphere := getbiggestinsphere;   { Get 'NUMINSPHERE'th item from list }
    if trysphere < bestsphere then     { If this radius is the smallest yet,}
      bestsphere := trysphere;         {   then remember it.                }
  end;
  writeln('Brute-force:');
  writeln('   Best sphere: ',bestsphere:15:10);
  writeln('   Number of comparisons: ',numcomps:8)
end;


{ **************************** My algorithm *********************** }

procedure makespace;        { Build the space-grid.  See documentation at }
var x,y,z : integer;        { beginning of program for details.           }
    temp : sptr;
    thispoint : integer;
begin
  spacesize := trunc(cube(8*PI*NUMPOINTS/NUMINSPHERE));
  usespace := spacesize-1;
  if usespace > MAXUSESPACE then writeln('****** NOT ENOUGH MEMORY FOR GRID');
  for x := 0 to usespace do
    for y := 0 to usespace do
      for z := 0 to usespace do
        space[x,y,z] := nil;     { Clear the grid }
  for thispoint := 1 to NUMPOINTS do     { Go through every point... }
  begin
    new(temp);
    temp^.index := thispoint;
    x := trunc(points[thispoint].x * spacesize);
    y := trunc(points[thispoint].y * spacesize);
    z := trunc(points[thispoint].z * spacesize);
    temp^.next := space[x,y,z];          { Put this point into proper }
    space[x,y,z] := temp;                {   sector in grid.          }
  end
end;


procedure elegantcluster;    { Find smallest sphere containing NUMINSPHERE }
                             {   points by looping through every point,    }
                             {   checking ROUGHLY only the points within   }
                             {   a radius less than or equal to the        }
                             {   minimum radius found so far.              }
var bestsphere,trysphere : real;
    xmin,xmax,ymin,ymax,zmin,zmax : integer; { Dimensions of box to check }
    thispoint : integer;              { The current point to test against }
    x,y,z : integer;                  { The current grid we are testing   }
    distx,disty,distz,dist : real;    { For computing distances           }
    numcomps : integer;               { # comparisons                     }
    cpindex : sptr;          { Pointer into point list for a grid sector  }
begin
  makespace;
  bestsphere := 1.732050808;    { Start with radius of distance from one }
  numcomps := 0;                {   corner of unit cube to other: 3^.5   }
  for thispoint := 1 to NUMPOINTS do    { Loop for every point }
  begin
    clearprilist;
    xmin := trunc((points[thispoint].x - bestsphere) * spacesize);
    xmax := trunc((points[thispoint].x + bestsphere) * spacesize);
    ymin := trunc((points[thispoint].y - bestsphere) * spacesize);
    ymax := trunc((points[thispoint].y + bestsphere) * spacesize);
    zmin := trunc((points[thispoint].z - bestsphere) * spacesize);
    zmax := trunc((points[thispoint].z + bestsphere) * spacesize);
    if xmin < 0 then xmin := 0;
    if ymin < 0 then ymin := 0;               { Get dimensions of box      }
    if zmin < 0 then zmin := 0;               { containing every sector in }
    if xmax > usespace then xmax := usespace; { grid that we want to check }
    if ymax > usespace then ymax := usespace; { against the current point  }
    if zmax > usespace then zmax := usespace;
    for x := xmin to xmax do
      for y := ymin to ymax do
        for z := ymin to ymax do   { Loop through every sector in this box }
        begin
          cpindex := space[x,y,z];
          while cpindex <> nil do  { Test against every point in this sector}
          begin
            if thispoint <> cpindex^.index then  { Don't test point against }
            begin                                {   itself.                }
              distx := points[thispoint].x - points[cpindex^.index].x;
              disty := points[thispoint].y - points[cpindex^.index].y;
              distz := points[thispoint].z - points[cpindex^.index].z;
              dist := sqrt(distx*distx + disty*disty + distz*distz);
              numcomps := numcomps + 1;
              if dist < bestsphere then  { If dist less than smallest sphere}
                insertprilist(dist)      {   insert distance into pri list  }
            end;
            cpindex := cpindex^.next     { Get next point in this sector }
          end
        end;
    trysphere := getbiggestinsphere;
    if trysphere < bestsphere then
      bestsphere := trysphere
  end;
  writeln('Efficient Algorithm:');
  writeln('   Best sphere: ',bestsphere:15:10);
  writeln('   Number of comparisons: ',numcomps:8)
end;


begin
  seed;
  writeln('How many points?');
  readln(NUMPOINTS);
  if NUMPOINTS < NUMINSPHERE then
    writeln('***** Must have at least ',NUMINSPHERE:1,' points.')
  else
  begin
    writeln('Testing ',NUMPOINTS:1,' points.');
    initprilist;
    generatepoints;
    brutecluster;
    elegantcluster
  end
end.

<================================== CUT HERE ==============================>
<-------------------------  Speed and accuracy test results --------------->

seed:
How many points?
Testing 50 points.
Brute-force:
   Best sphere:    0.3067962678
   Number of comparisons:     2450
Efficient Algorithm:
   Best sphere:    0.3067962678
   Number of comparisons:      946

Profile:-----------------------------------------
 %time  cumsecs  #call  ms/call  name
  31.6     0.10      1   100.01  _brutecluster        <-----------
  15.8     0.15                  libm$dsqrt_r5
  10.5     0.18      1    33.34  _elegantcluster      <-----------
   5.3     0.20    730     0.02  _DISPOSE
   5.3     0.22                  _FCALL
   5.3     0.23    784     0.02  _NEW
   5.3     0.25    101     0.17  _clearprilist
   5.3     0.27    581     0.03  _insertprilist



seed:
How many points?
Testing 300 points.
Brute-force:
   Best sphere:    0.1569231423
   Number of comparisons:    89700
Efficient Algorithm:
   Best sphere:    0.1569231423
   Number of comparisons:     5617

Profile:--------------------------------------
 %time  cumsecs  #call  ms/call  name
  44.2     3.27      1  3267.00  _brutecluster    <------------------
  31.7     5.61                  libm$dsqrt_r5
   7.0     6.13                  mcount
   6.4     6.60  95317     0.00  _sqrt
   2.9     6.82      1   216.69  _elegantcluster  <------------------
   1.4     6.92   3343     0.03  _malloc
   1.1     7.00   2358     0.04  _insertprilist




                  Bruteforce:                    My algorithm:
                                               (Average of 3 runs)
Points:   #Comparisons:  comps/points:     #Comparisons:  comps/points:
     50            2450         49.000               958         19.160
     75            5550         74.000              1241         16.547
    100            9900         99.000              2111         21.110
    150           22350        149.000              2785         18.567
    200                                             3689         18.445
    250                                             5120         20.480
    300                                             6010         20.033
    350                                             7149         20.426
    400                                             7658         19.145
    600                                            11404         19.007
    800                                            16781         20.976
   1000                                            20438         20.438

                                                                   ^
                                                                   |
                                                             Look, LINEAR!!!



Good luck.

-- Marshall Levine
mplevine@phoenix.princeton.edu

sarrel@galley.cis.ohio-state.edu (Marc Sarrel) (05/15/89)

Well, I don't argue that the posted algorithm works (although I'm not
sure about its linearity), I don't think that it answers the question
exactly.

Here is my solution, which I _know_ to be linear in the number of
points.

Problem:
Given a set of points in three space, find the two that are farthest
apart and use them to define a minimum bounding sphere.

Solution:
1) Pick a point, any point, call it p1
2) Find the farthest point from p1, call it p2.
3) Find the farthest point from p2, call it p3.
4) Use p2, p3 and the distance between them to define your minimum
   bounding sphere.

That's it.  It is guaranteed to be linear, too.

The reason that this algorithm works is as follows.  We know that
there are two points such that the distance between them is greater
than the distance between any other two points, call them p2 and p3
(as above).  Call the line between p2 and p3 d1 (for diameter).  Call
the point midway along d1 c1 (for center).  Put a plane perpendicular
to d1 through the point c1.  This divides the points into two sets.
For all the points on one side of the plane, point p2 is the farthest
away and for all the points on the other side, p3 is farthest away.
Now, we can see that step two above finds one of the points that
define the bounding sphere and that step three must find the other.
QED (for an informal argument).

If you don't beleive that it works, try a few examples from the 2D
case to convince yourself.  I don't know if this has been published
anywhere (I'm almost sure that it has), but I can't remember reading
about it.  It was inspired by a "Computer Games" (or whatever it was
called at the time) column in an old issue of Scientific American.

The main problem with implementing this algorithm is that the square
and square root function needed to calculate distances are expensive.
In practice, the square root can be ignored except to calculate the
final distance for the diameter.  However, there may be more efficient
algorithms that don't bear any relation to this one.
-=-
"Master, why is the letter 'i' the symbol for current?"  "Because there is
no letter 'i' in the word 'current'."  "Master, why do we use the letter
'j' for sqrt(-1)?"  "Because we use the letter 'i' for current."  Whereupon
the Master struck the Disciple, and the Disciple became enlightened.

turk@media-lab.media.mit.edu (Matthew Turk) (05/15/89)

In article <SARREL.89May14223947@galley.cis.ohio-state.edu> sarrel@galley.cis.ohio-state.edu (Marc Sarrel) writes:
>
>  Here is my solution, which I _know_ to be linear in the number of
>  points.
>
>  Problem:
>  Given a set of points in three space, find the two that are farthest
>  apart and use them to define a minimum bounding sphere.
>
>  Solution:
>  1) Pick a point, any point, call it p1
>  2) Find the farthest point from p1, call it p2.
>  3) Find the farthest point from p2, call it p3.
>  4) Use p2, p3 and the distance between them to define your minimum
>     bounding sphere.
>
>  That's it.  It is guaranteed to be linear, too.
>
Yes, but unfortunately it's *not* guaranteed to work!

>  The reason that this algorithm works is as follows.  We know that
>  there are two points such that the distance between them is greater
>  than the distance between any other two points, call them p2 and p3
>  (as above).  Call the line between p2 and p3 d1 (for diameter).  Call
>  the point midway along d1 c1 (for center).  Put a plane perpendicular
>  to d1 through the point c1.  This divides the points into two sets.
>  For all the points on one side of the plane, point p2 is the farthest
>  away and for all the points on the other side, p3 is farthest away.
>  Now, we can see that step two above finds one of the points that
>  define the bounding sphere and that step three must find the other.
>  QED (for an informal argument).
>
This reasoning is incorrect.  Your p2 and p3 are the points farthest
away on each side of a plane -- they are not related to the bounding
sphere.  As a counter example to your algorithm, note this figure:

			   |
                           3        point 1 - (0,0)
                           |	       	  2 - (10,0)
                           |	       	  3 - (0,9)
                           |	       	  4 - (0,-8)
                           |
              -------------1------------2--
                           |
                           |
                           |
                           4
			   |

Step (1) - I'll pick p1, at (0,0).  Step (2) - the farthest point is
p2, at (10,0).  Step (3) - the farthest from p2 is p3, at (0,9).  Step (4)
- the minimum sphere defined by p2 and p3 is centered at (5,4.5), and
has a radius of 6.727.  But since the distance between the sphere's
center and p4(0,-8) is 13.463, p4 is not in the minimum sphere.

Better luck next time!


	Matthew Turk
	Vision Science Group, MIT Media Lab
	turk@media-lab.media.mit.edu

ken@uf.msc.umn.edu (Kenneth Chin-Purcell) (05/15/89)

In article <SARREL.89May14223947@galley.cis.ohio-state.edu> 
      sarrel@galley.cis.ohio-state.edu (Marc Sarrel) writes:
>Problem:
>Given a set of points in three space, find the two that are farthest
>apart and use them to define a minimum bounding sphere.

You are assuming that a line between the two farthest points 
defines the axis of the minimum bounding sphere.  I don't think
this is true.  Take these three points:

                           P3
                        xxxxxxx
                      xx       xx
                     x           x
           ---------P1---axis----P2---------
                     x           x
                      xx       xx
                        xxxxxxx

P1 and P2 are the farthest apart.  Use line P1-P2 to define the axis
of sphere x.  There is a region, near the equator of the sphere, 
where a point (P3) could be outside the sphere, yet less than P1-P2
distance from either point.

To visualize this, think of two larger spheres, each with *radius*
P1-P2, one centered at P1, the other at P2.  The intersection of
these spheres is the region where P3 can be while keeping P1 and P2
the farthest points.  Sphere x is just a subset of this region.

Since in general it takes 4 points to define a sphere, my
intuition leads me to think that the true general solution will
involve 4 bounding points.

-- 
Ken Chin-Purcell                                 (aka ken@msc.umn.edu)
Minnesota Supercomputer Center                   "Seymore doesn't
1200 Washinton Ave., Minneapolis, MN 55415        believe in paging"

guilford@rpics (Jim Guilford) (05/15/89)

In article <17241@versatc.UUCP>, ritter@versatc.UUCP (Jack Ritter) writes:
> Given a cluster of points in 3 space, is there
> a good method for finding the minumum radius
> sphere which encloses all the points?

In article <SARREL.89May14223947@galley.cis.ohio-state.edu> sarrel@galley.cis.ohio-state.edu (Marc Sarrel) writes:
>Here is my solution, which I _know_ to be linear in the number of
>points.
>
>Problem:
>Given a set of points in three space, find the two that are farthest
>apart and use them to define a minimum bounding sphere.
>
>Solution:
>1) Pick a point, any point, call it p1
>2) Find the farthest point from p1, call it p2.
>3) Find the farthest point from p2, call it p3.
>4) Use p2, p3 and the distance between them to define your minimum
>   bounding sphere.
>
>That's it.  It is guaranteed to be linear, too.

I'm not sure if your problem statement is actually the problem
statement, or a partial solution statement as well, but if you are
trying to solve the earlier problem (bounding sphere of n points),
then I'm afraid that doesn't work.

Consider three points in a plane arranged as an equilateral triangle: 

         (1,sqrt(3))
            A
           / \
          /   \
         /     \
        /       \
       B---------C
     (0,0)      (2,0)

If you (for simplicity) pick B and C as the farthest points, and then
center a sphere at (1,0) with a radius of 1 (so that it just encloses
B and C), then A will be outside the sphere. In general, the center
will not lie on the line between any two points.

Consider also the case of four points arranged in a tetrahedron.

--JimG
guilford@turing.cs.rpi.edu

tada@athena.mit.edu (Michael Zehr) (05/15/89)

In article <17241@versatc.UUCP>, ritter@versatc.UUCP (Jack Ritter) writes:
> Given a cluster of points in 3 space, is there
> a good method for finding the minumum radius
> sphere which encloses all the points?

In article <SARREL.89May14223947@galley.cis.ohio-state.edu> sarrel@galley.cis.ohio-state.edu (Marc Sarrel) writes:
>Here is my solution, which I _know_ to be linear in the number of
>points.
>
>Problem:
>Given a set of points in three space, find the two that are farthest
>apart and use them to define a minimum bounding sphere.
>
>Solution: [...]

People seem to be having a problem understanding what the problem is.
The problem is not limited to 3 points (like one poster's solution), nor
is it to find which pair of points are farthest apart in a set of
points.

The real problem is to speed up a ray-tracing program by providing some
sort of bounding sphere for an object.  This is really 2 problems:

1. given a set of points, find the minimum-bounding sphere
2. given a set of points, quickly find a "good" bounding sphere

if the data is static, but you're going to view it from many directions
(or lighting models, etc), then you want 1.  it doesn't really matter
how long it takes, because it's done once and used many times.  but you
want it to be minimal, to save as much time as possible.

if the data isn't static (the points move with respect to each other
between each run of the program) then you want to find the bounding
sphere quickly.  and it's better to quickly find a sphere that's
slightly larger than what's needed if it takes a less time to compute
it.

In conclusion...  be sure you understand the problem before giving a
solution.  try and be more certain of your "knowing" the solution.  but
do go ahead and post a possible solution.  (i'm not trying to discourage
people from posting, i'm trying to prevent people from using a solution
that isn't.)

-michael j zehr

ECULHAM@UALTAVM.BITNET (Earl Culham) (05/16/89)

>Here is my solution, which I _know_ to be linear in the number of
>points.
>
>Problem:
>Given a set of points in three space, find the two that are farthest
>apart and use them to define a minimum bounding sphere.
>
>Solution:
>1) Pick a point, any point, call it p1
>2) Find the farthest point from p1, call it p2.
>3) Find the farthest point from p2, call it p3.
>4) Use p2, p3 and the distance between them to define your minimum
>   bounding sphere.
>
>That's it.  It is guaranteed to be linear, too.
>
 
This algorithm does not necessarily find the two points furthest apart.
 
Even if it did, the two points furthest apart do not need to be on
opposite sides of the minimum bounding sphere.
 
Counter Example (in only 2D)
A = (0,0)
B = (0,10)
C = (6,5)
D = (-6,5)
 
Choose P1 to be A.
P2 becomes B.
P3 becomes A.
 
This did not find the two furthest points C and D.
 
If we forget about the point D, then this algorithm does find the two
furthest. But the circle delimited by these two furthest points, does
not contiain the point C.

jep@oink.UUCP (James E. Prior) (05/16/89)

In article <12851@umn-cs.CS.UMN.EDU> ken@uf.msc.umn.edu (Kenneth Chin-Purcell) writes:
...
>Since in general it takes 4 points to define a sphere, my
>intuition leads me to think that the true general solution will
>involve 4 bounding points.

You are correct in that a general solution requires 4 points, but we
aren't faced with a general problem.  We are faced with a special one,
one that needs the smallest sphere.  The requirement of the smallest 
sphere is information itself.  Enough so that one doesn't need 4 points
to define the sphere.  
-- 
Jim Prior    jep@oink    osu-cis!n8emr!oink!jep    N8KSM

jep@oink.UUCP (James E. Prior) (05/16/89)

In article <11444@bloom-beacon.MIT.EDU> tada@athena.mit.edu (Michael Zehr) writes:
...
<The real problem is to speed up a ray-tracing program by providing some 
<sort of bounding sphere for an object.  This is really 2 problems: 
< 
<1. given a set of points, find the minimum-bounding sphere 
<2. given a set of points, quickly find a "good" bounding sphere 
...
<if the data isn't static (the points move with respect to each other 
<between each run of the program) then you want to find the bounding 
<sphere quickly.  and it's better to quickly find a sphere that's 
<slightly larger than what's needed if it takes a less time to compute 
<it. 

I would always want the closest fit, not merely a close fit, regardless of
how many frames I was calculating.  Ray traced images are so expensive
computationally, that spending the time at the beginning of each frame 
for the closest fit will speed it up.  If it's a good optimization for
a single frame, you should expect the improvement for doing many frames.
-- 
Jim Prior    jep@oink    osu-cis!n8emr!oink!jep    N8KSM

mplevine@phoenix.Princeton.EDU (Marshall P. Levine) (05/18/89)

In article <SARREL.89May14223947@galley.cis.ohio-state.edu> sarrel@galley.cis.ohio-state.edu (Marc Sarrel) writes:
>Well, I don't argue that the posted algorithm works (although I'm not
>sure about its linearity), I don't think that it answers the question
>exactly.


Actually, Marc, I carefully analyzed my program.  It finds the minimum spanning
sphere on the average at a linear time.  It is possible that its time will not
be spectacular (given a really bad distribution of points) but it should never
be n^2.  I don't think most people who are posting replies to my program 
understand what I did.  My algorithm, given a space of randomly distributed 
points, will find the sphere of minimum radius, centered on ANY POINT IN THE
SPACE, that contains NUMINSPHERE other points.  The key here is that the sphere
can be centered on ANY POINT IN THE SPACE.  Most of the ideas that I have seen
posted involve comparing every point against every other point.  That is 
guaranteed to be n^2.  In fact, my BRUTEFORCE algorithm included in my posting
will do precisely that.  If you are skeptical about the accuracy of my algo,
compile the source youself.  I wrote it in standard pascal, and you should have
no problems compiling it.  As for speed tweeking, obviously the square-root
function is unnecessary.  I am not stupid.  While my algos are sorting the
distances that they consider, they do not need the square roots.  The only 
square root that needs to be calculated is the square root of the final number
produced by the algos, which will be the answer to the problem.  I included
the square-root functions because I figured that it would be easier for people
to understand what I did.  Once again, if there are any questions, feel free
to contact me.  Good luck.


Best wishes,
Marshall Levine
mplevine@phoenix.princeton.edu