Path: utzoo!attcan!utgpu!jarvis.csri.toronto.edu!clyde.concordia.ca!uunet!munnari.oz.au!bruce!dbrmelb!davidp
From: davidp@dbrmelb.dbrhi.oz (David Paterson)
Newsgroups: comp.graphics
Subject: Re: smallest sphere enclosing a set of points
Summary: A FUN solution !
Message-ID: <687@dbrmelb.dbrhi.oz>
Date: 15 Dec 89 07:42:14 GMT
References: <Nov.28.22.18.43.1989.8200@caip.rutgers.edu> <664@dbrmelb.dbrhi.oz>
Organization: CSIRO, Div. Building Constr. and Eng'ing, Melb., Australia
Lines: 87


We've had lots of solutions to this problem posted. I've posted two already
(see above reference). The person who asked is probably sick to death of
them now. I'll test his patience a little bit more by posting this 'fun'
solution.

First a comment on the two solutions that I posted earlier. The first
(brute force) method has convergence order m^(n+1). It can be speeded up
considerably by first finding the convex hull. The second (bottom up)
method has convergence m^2 at best. It can be made watertight in n
dimensions by inserting a bit of the brute force method when it is
needed to find the smallest n-sphere containing n+2 points when this is
required.

The 'fun' solution below is watertight in n dimensions and has convergence
order m. It's not the fastest algorithm of order m but it sure is the
shortest. Here it is:

      PARAMETER(MAXDIMENS=100,MAXPOINTS=100000)
      DIMENSION POINTS(MAXDIMENS,MAXPOINTS),P0(MAXDIMENS),P1(MAXDIMENS)
      READ *, IDIMENSION,IPOINTS,NSTEPS
      READ *, ((POINTS(I,J),I=1,IDIMENSION),J=1,IPOINTS)
      DO 10 I=1,IDIMENSION
 10   P0(I)=0.
      DISTMAX=0.
      DO 30 J=1,IPOINTS
      DISTANCE=0.
      DO 20 I=1,IDIMENSION
 20   DISTANCE=DISTANCE+(P0(I)-POINTS(I,J))**2
 30   IF(DISTANCE.GT.DISTMAX) DISTMAX=DISTANCE
      DO 80 N=1,NSTEPS
      DO 40 I=1,IDIMENSION
 40   P1(I)=P0(I)+RAND()-.5
      XLENMAX=0.
      DO 60 J=1,IPOINTS
      XLENGTH=0.
      DO 50 I=1,IDIMENSION
 50   XLENGTH=XLENGTH+(P1(I)-POINTS(I,J))**2
 60   IF(XLENGTH.GT.XLENMAX) XLENMAX=XLENGTH
      IF(XLENMAX.LT.DISTMAX) THEN
	DISTMAX=XLENMAX
	DO 70 I=1,IDIMENSION
 70     P0(I)=P1(I)
      END IF
 80   CONTINUE
      PRINT *,'RADIUS=',SQRT(DISTMAX)
      PRINT *,'CENTRE=',(P0(I),I=1,IDIMENSION)
      STOP
      END

Here P0 approaches the centre of the n-sphere as N tends to infinity.
NSTEP controls the speed and accuracy of the solution. RAND() denotes
a function that returns a different random number between 0 and 1 each
time it is called. If you don't already have such a function try the
following:

      FUNCTION RAND()
      PARAMETER(K=5701,J=3612,M=566927)
      RM=M
      IRAND=MOD(J*INT(RAND*RM)+K,M)
      RAND=(IRAND+.5)/RM
      RETURN
      END

The algorithm, although it is of order m, is woefully inefficient. The
'fun' part consists of playing around with ways to speed it up. For
example, we could replace:
 10   P0(I)=0.
by
 10   P0(I)=POINT(I,1)
or set P0(I) to be the mean of all points or something else. We could
replace:
 40   P1(I)=P0(I)+RAND()-.5
by
 40   P1(I)=P0(I)+(RAND()-.5)*.1*SQRT(DISTMAX)
or
 40   P1(I)=P0(I)+(RAND()-.5)*SQRT(DISTMAX)/LOG(N+5.)
or move a random distance in a direction towards the point furthest
from P0(I) or something else. We could replace:
      DO 80 N=1,NSTEPS
by a better convergence criterion. The possibilities for playing around
with this method are literally endless.

------------------------------------------------------------------------
David Paterson, CSIRO Division of Building, Construction and Engineering,
Highett, Victoria, AUSTRALIA 3190.
------------------------------------------------------------------------