// Frank Nielsen and Richard Nock
// Smallest Enclosing Disk Approximation Algorithm for Disk Sets
//			(November 2003. v1.0)

#include "stdafx.h"
#include <windows.h>
#include <stdlib.h>
#include <math.h>

typedef float Real;

#define min2(a,b) ((a)<(b) ? (a): (b))
#define max2(a,b) ((a)>(b) ? (a): (b))
#define realsqrt(x) sqrtf(x)
#define RandReal(a,b) ((a)+(((b)-(a))*((Real)rand()/(Real)(RAND_MAX))))

class disk2d {public: Real x,y,r;};

// 
// Solve decision problems by dichotomy (w/o bootstrapping)
void sebdisk(disk2d *P, int n, Real e,  disk2d &seb)
{
Real r,rs,l,x,ym,yM,common,xmin,xmax,xm,xM,a,b,dist,adist,d,maxd; 
int i,d1,d2; 
bool ebpierceable,qdisjoint;
	 
	// Initialization
	//

	xmin=P[0].x-P[0].r;xmax=P[0].x+P[0].r;maxd=P[0].r;
	
	for(i=1;i<n;i++)
	{xmin=min2(xmin,P[i].x-P[i].r);xmax=max2(xmax,P[i].x+P[i].r);
	 d=realsqrt((P[0].x-P[i].x)*(P[0].x-P[i].x)+(P[0].y-P[i].y)*(P[0].y-P[i].y))+P[i].r;
	 maxd=max2(maxd,d);}

	b=maxd;a=max2(b/2.0,P[0].r);  // optimal radius lies in range [a,b] always greater than (xmax-xmin)/2
	e*=(b/2.0); // Convert to absolute epsilon precision. Should be e*=(b/2.0); 

	while(b-a>e)
		{
		r=(a+b)/2.0;rs=r*r;
		// candidate strip
		//
		xM=xmin+r;xm=xmax-r;// band intersection all disks	
		ebpierceable=false;qdisjoint=false;
		//
		// dichotomy for the decision problem
		while((xM>xm)&&(xM-xm>e)&&(!ebpierceable)&&(!qdisjoint))
			{	
			l=(xm+xM)/2.0;d1=d2=0;
			common=realsqrt((r-P[0].r)*(r-P[0].r)-(l-P[0].x)*(l-P[0].x));
			ym=P[0].y-common;yM=P[0].y+common;i=1;
				while((i<n)&&(ym<=yM))
				{
				common=realsqrt((r-P[i].r)*(r-P[i].r)-(l-P[i].x)*(l-P[i].x));	
				//
				// those tests are of degree 4
				if (ym<P[i].y-common) {d1=i;ym=P[i].y-common;}
				if (yM>P[i].y+common) {d2=i;yM=P[i].y+common;}
				i++;
				}

			if (yM>=ym) {
				seb.x=l;seb.y=ym;
				ebpierceable=true; 
			}
			else{// choose on which side to recurse	
				dist=realsqrt((P[d2].x-P[d1].x)*(P[d2].x-P[d1].x)+(P[d2].y-P[d1].y)*(P[d2].y-P[d1].y));
				adist=(P[d1].r*P[d1].r-P[d2].r*P[d2].r+dist*dist)/(2.0*dist);
				x=P[d1].x+(adist/dist)*(P[d2].x-P[d1].x);
				if (x>l) xm=l; else xM=l;
				}
			}// while decision problem
		if (ebpierceable) {b=r;seb.r=r;} else a=r;
		}// while r
}



#define N 100000
#define epsilon 0.001

int main(int argc, char* argv[])
{
	disk2d *S,seb;
	int i,seed=23,trial,maxtrial=100;
	Real radiusopt,distsq,distd,app;
	double ct;
	LARGE_INTEGER startcounter,stopcounter,freq;
	
	printf("Smallest Enclosing Disk Approximation.\n");
	S=new disk2d[N]; if (S==NULL) {printf("Not enough memory!\n");return -1;}
	
	QueryPerformanceFrequency(&freq);
	
	for(trial=0;trial<maxtrial;trial++){
	srand(seed++); 
	// Draw a uniform sample
	//
	for(i=0;i<N;i++){S[i].x=RandReal(-1.0,1.0);S[i].y=RandReal(-1.0,1.0);S[i].r=RandReal(0.0,1.0);}
	i=rand()%N; S[i].x=S[i].y=-1.0;S[i].r=1.0;i=(i+rand())%N; S[i].x=S[i].y=1.0;S[i].r=1.0;
	radiusopt=1.0+realsqrt(2.0);app=(1.0+epsilon)*radiusopt;

	SetThreadPriority(GetCurrentThread(), THREAD_PRIORITY_TIME_CRITICAL);
	QueryPerformanceCounter(&startcounter);
		sebdisk(S,N,epsilon,seb);
	QueryPerformanceCounter(&stopcounter);

	ct=((double)stopcounter.LowPart-(double)startcounter.LowPart)/(double)freq.LowPart;
	printf("N=%d\tepsilon=%lf\n[Time %lf sec.] Center=(%lf,%lf), Radius=%lf\n",N,epsilon,ct,seb.x,seb.y,seb.r);
	} // trial
	return 0;
}