// ------------------------------------------------------------------------
  // This program is complementary material for the book:
  //
  // Frank Nielsen
  //
  // Visual Computing: Geometry, Graphics, and Vision
  //
  // ISBN: 1-58450-427-7
  //
  // Charles River Media, Inc.
  //
  //
  // All programs are available at www.charlesriver.com/visualcomputing/
  //
  // You may use this program for ACADEMIC and PERSONAL purposes ONLY. 
  //
  //
  // The use of this program in a commercial product requires EXPLICITLY
  // written permission from the author. The author is NOT responsible or 
  // liable for damage or loss that may be caused by the use of this program. 
  //
  // Copyright (c) 2005. Frank Nielsen. All rights reserved.
  // ------------------------------------------------------------------------
   
  // ------------------------------------------------------------------------
  // File: SmallEnclosingBall.cpp
  // 
  // Description: Approximate the smallest enclosing ball of a high-dimensional
  // point set using Badoui and Clarkson's algorithm.
  // ------------------------------------------------------------------------
  
  include <stdafx.h>
  include <math.h>
  
  // dimension
  const int d=100;
  
  using namespace std;
  
  inline double drand()
  {
  return (double)rand()/(double)RAND_MAX;
  }
  
  // d-dimensional point set
  class PointdD{
  public:
          double coord[d];
  
          double Distance(PointdD p)
          {
                  double result=0.0;
                  for(int i=0;i<d;i++)
                                  result+= ((coord[i]-p.coord[i])*(coord[i]-p.coord[i]));
                  return sqrt(result);
          }
  
          void InitializeRandom()
          {
                  for(int i=0;i<d;i++) 
                                  coord[i]=drand();
          }
  
          // initialize vertex of a cube
          void InitializeCube()
          {
          int setd=rand()\%d;
          for(int i=0;i<d;i++) 
                                  coord[i]=0;
          coord[setd]=1.0;
          }
  
          void PrintPoint()
          {
          cout<<"[";
          for(int i=0;i<d;i++)  cout<<coord[i]<<" ";
          cout<<"]"<<endl;
          }
  };
  
  // Paper is available here:
  // cm.bell-labs.com/who/clarkson/coresets1.html
  void  BadoiuClarsonIteration(PointdD *setcore,  int core, PointdD &center, double &radius, int nbiter)
  {
  int i,iter,winner;
  double distmax,dist;
  
  // choose any point of the set as the initial circumcenter
  center=setcore[0];
  
  for(iter=0;iter<nbiter;iter++)
  {
          winner=0;
          distmax=center.Distance(setcore[0]);
  
          // Maximum distance point
          for(i=1;i<core;i++)
                  {
                  dist=center.Distance(setcore[i]);
                  if (dist>distmax) {winner=i;distmax=dist;}
          }
          radius=distmax;
          // Update
          for(i=0;i<d;i++)
          {
          center.coord[i]+=(1.0/(iter+1.0))*(setcore[winner].coord[i]-center.coord[i]);
          }
  
  } // for
  }
  
  int _tmain(int argc, _TCHAR* argv[])
  {
  PointdD *set, center;
  int i,n,nb;
  double epsilon=0.01,radius;
  
  cout<<"Visual Computing: Geometry, Graphics, and Vision (ISBN:1-58450-427-7)"<<endl;
  cout<<"Demo program\n\n"<<endl;
  
  srand(2005);
     
          n=2*d;
  
          set=new PointdD[n];
  
          for(i=0;i<n;i++)
          {
          set[i].InitializeRandom();
          //set[i].InitializeCube();
          }
  
          cout<<"Implementation of Badoiu-Clarkson algorithm for computing"<<endl;
          cout<<"a small enclosing ball in large dimensions."<<endl;
          cout<<endl;
          cout<<"Create a set of "<<n<<" points in dimension "<<d<<endl;
  
          nb=(int)ceil(1.0/(epsilon*epsilon));
          cout<<"Epsilon="<<epsilon<<" Nb of iterations:"<<nb<<endl;
  
          cout<<"Be patient..."<<endl;
          BadoiuClarsonIteration(set,n,center,radius,nb);
  
          cout<<"The radius r of the ball is:"<<radius<<endl;
          cout<<"r is guaranteed to be less than "<<1+epsilon<<"r*, where r* is the optimal radius"<<endl;
          cout<<"Circumcenter:";
          center.PrintPoint();
  
          delete [] set;
  
          
  cout<<"Press Return key"<<endl;
  char line[256];
  gets(line);
  
          return 0;
  }