Path: utzoo!attcan!utgpu!jarvis.csri.toronto.edu!mailrus!tut.cis.ohio-state.edu!rutgers!njin!princeton!phoenix!mplevine
From: mplevine@phoenix.Princeton.EDU (Marshall Polimer Levine)
Newsgroups: comp.graphics
Subject: Re: Ray Traced Bounding Spheres
Summary: Algorithm for finding minimum sphere enclosing N points - linear!
Keywords: Ray trace, bounding volumes, intersection, sphere, linear, enclosing, algorithm
Message-ID: <8411@phoenix.Princeton.EDU>
Date: 14 May 89 16:08:32 GMT
References: <17241@versatc.UUCP> <17249@versatc.UUCP>
Reply-To: mplevine@phoenix.Princeton.EDU (Marshall Polimer Levine)
Organization: Princeton University, NJ
Lines: 441

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