// ------------------------------------------------------------------------
  // 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: MiniBall.cpp
  // 
  // Description: Randomized Recursive Incremental Construction of the Smallest
  // Enclosing Ball of a Set of 2D Points.
  // ------------------------------------------------------------------------
  
  // Miniball: this is a pedagogical implementation. 
  // For efficient implementation, rather update the enclosing balls and consider square radii
  
  include <stdafx.h>
  include <math.h>
  
  using namespace std;
  
  typedef  long double number;
  
  class point {
  public: 
  number x,y;
  
  inline number Distance(point q)
  {return sqrt((q.x-x)*(q.x-x) + (q.y-y)*(q.y-y));}
  };
  
  class disc{
  public: 
  point circumcenter;
  number radius;
  };
  
  inline double drand(){return (double) rand () / RAND_MAX ;}
  
  inline SwapPoint(point& a, point& b)
  {point tmp; tmp=a; a=b; b=tmp;}
  
  // The unique circle passing through exactly three non-collinear points  p1, p2 ,p3
  void SolveDisc3(point p1, point p2, point p3, point & center, number &rad)
  {
  number a = p2.x - p1.x;
  number b = p2.y - p1.y;
  number c = p3.x - p1.x;
  number d = p3.y - p1.y;
  number e = a*(p2.x + p1.x)*0.5 + b*(p2.y + p1.y)*0.5;
  number f = (c*(p3.x + p1.x)*0.5) + (d*(p3.y + p1.y)*0.5);
  number det = a*d - b*c;    
  
  center.x = (d*e - b*f)/det;   
  center.y = (-c*e + a*f)/det;
  
  rad =p1.Distance(center);
  }
  
  // The smallest circle passing through two points
  void SolveDisc2(point p1, point p2, point & center, number &rad)
  {
  center.x = 0.5*(p1.x+p2.x);   
  center.y = 0.5*(p1.y+p2.y);
  
  rad =p1.Distance(center);
  }
  
  void MiniDisc(point* set, int n, point* basis, int b, point& center, number& rad)
  {
  point cl; number rl;
  int k;
  
  // Terminal cases
  if (b==3)  SolveDisc3(basis[0], basis[1], basis[2], cl, rl);
  
  if ((n==1)&&(b==0)) {rl=0.0; cl=set[0];}
  if ((n==2)&&(b==0))  SolveDisc2(set[0],set[1], cl, rl); 
  if ((n==0)&&(b==2))  SolveDisc2(basis[0],basis[1], cl, rl); 
  if ((n==1)&&(b==1))  SolveDisc2(basis[0], set[0], cl, rl); 
  
  // General case
  if ((b<3)&&(n+b>2))
          {
          // Randomization: choosing a pivot
          k=rand()\%n; // between 0 and n-1
          if (k!=0) SwapPoint(set[0],set[k]);
  
          MiniDisc(&set[1],n-1,basis,b,cl,rl);
  
          if (cl.Distance(set[0])>rl)
                  {
                  // Then set[0] necessarily belongs to the basis.
                          basis[b++]=set[0];
                          MiniDisc(&set[1],n-1,basis,b,cl,rl);
                  }        
          }
          center=cl;rad=rl;
  }
  
  // Example
  int _tmain(int argc, _TCHAR* argv[])
  {
  point *set, * basis, center;
  int i, b=0;
  number rad;
  int N=10000;
  
  cout<<"Visual Computing: Geometry, Graphics, and Vision (ISBN:1-58450-427-7)"<<endl;
  cout<<"Demo program\n\n"<<endl;
  
  cout<<"Compute the smallest enclosing ball of a 2D point set."<<endl;
  
  center.x=center.y=0;rad=0;
  set=new point[N];
  basis=new point[3];
  
  for(i=0;i<N;i++) {set[i].x=drand();set[i].y=drand(); 
                                  }
  
  cout <<"Randomized optimization: MiniDisc"<<endl;
  MiniDisc(set,N,basis,b, center, rad);
  
  cout <<"Center ("<<center.x<<","<<center.y<<") rad="<<rad<<endl;
  
  delete [] set;
  
  char line[256];
  cout<<"Press Return key"<<endl;
  gets(line);
  
          return 0;
  }