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

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

typedef double Real;

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

class point2d {public: Real x,y;};

// 
// Solve decision problems by dichotomy (w/o bootstrapping)
void seb2d(point2d *P, int n, Real e,  point2d &center, Real& radius)
{
Real r,rs,l,x,ym,yM,common,xmin,xmax,xm,xM,a,b,distsqr,d,maxd; 
int i,d1,d2; 
bool ebpierceable,qdisjoint;
	 
	// Initialization
	//
	xmin=P[0].x;xmax=xmin;maxd=0.0;
	
	for(i=1;i<n;i++)
	{xmin=min2(xmin,P[i].x);xmax=max2(xmax,P[i].x);
	 d=(P[0].x-P[i].x)*(P[0].x-P[i].x)+(P[0].y-P[i].y)*(P[0].y-P[i].y);
	 maxd=max2(maxd,d);}
	maxd=realsqrt(maxd);

	b=maxd;a=b/2.0;  // 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); 
	radius=1.0e8;center.x=0;

	while(b-a>e)
		{
		r=(a+b)/2.0;rs=r*r;
		// candidate strip
		//
		xm=max2(xmax-r,center.x-radius*e);
		xM=min2(xmin+r,center.x+radius*e); // 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(rs-(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(rs-(l-P[i].x)*(l-P[i].x));
				//
				// those tests are of degree 4 on Integers
				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) {
				center.x=l;center.y=ym;
				ebpierceable=true; 
				radius=r;
			}
			else{
				distsqr=(P[d1].x-P[d2].x)*(P[d1].x-P[d2].x)+(P[d1].y-P[d2].y)*(P[d1].y-P[d2].y);
				distsqr-=4.0*(r-e)*(r-e);
				if (distsqr>0) {qdisjoint=true;}
				else 
					{// choose on which side to recurse	
					x=(P[d1].x+P[d2].x)/2.0; 
					if (x>l) xm=l; else xM=l;
				}
				}
			}// while decision problem
		if (ebpierceable) {b=r;
		//FilterBall(P,n,center,r/(b-a)-(b-a)*r);
		} else {a=r;}
		}// while r
}

int N=100000;
Real epsilon=1.0e-3; 

int main(int argc, char* argv[])
{
	point2d *S,center;
	Real radius,radiusopt,distsq,app,theta;
	double ct,avgt,mint,maxt;
	LARGE_INTEGER startcounter,stopcounter,startt,stopt,freq;
	int i,trial,maxtrial=100,distribution;
	
	printf("Smallest Enclosing Disk Approximation.\n\n");
	S=new point2d[N]; if (S==NULL) {printf("Not enough memory!\n");return -1;}
	QueryPerformanceFrequency(&freq);
	SetThreadPriority(GetCurrentThread(), THREAD_PRIORITY_TIME_CRITICAL);
	printf("#Trial=%d, N=%d, Epsilon=%lf\n",maxtrial,N,epsilon);
	QueryPerformanceCounter(&startt);

	for(distribution=1;distribution<=2;distribution++)
	{mint=1.0e8;maxt=0.0;avgt=0.0;

	for(trial=0;trial<maxtrial;trial++)
	{
	srand(2003+trial); 
	// Draw a uniform sample
	//
	if (distribution==1)
	{for(i=0;i<N;i++){S[i].x=RandReal(-0.5,0.5);S[i].y=RandReal(-0.5,0.5);}
	i=rand()%N; S[i].x=S[i].y=-0.5;i=(i+rand())%N; S[i].x=S[i].y=0.5;
	radiusopt=0.5*realsqrt(2.0);
	}
	// Draw a ring distribution
	//
	if (distribution==2)
	{for(i=0;i<N;i++){radius=RandReal(0.99,1.0);theta=RandReal(-M_PI,M_PI);S[i].x=radius*cos(theta);S[i].y=radius*sin(theta);}
	i=rand()%N; S[i].x=0.0;S[i].y=1.0;i=(i+rand())%N; S[i].x=0.0;S[i].y=-1.0;
	radiusopt=1.0;
	}

	app=(1.0+epsilon)*radiusopt;app*=app;

	QueryPerformanceCounter(&startcounter);
		if (distribution==1) seb2d(S,N,epsilon,center,radius);
		if (distribution==2) seb2d(S,N,epsilon,center,radius);
	QueryPerformanceCounter(&stopcounter);

	ct=((double)stopcounter.LowPart-(double)startcounter.LowPart)/(double)freq.LowPart;
	mint=min2(mint,ct);maxt=max2(maxt,ct);avgt+=ct;
	} // number of trials
		printf("Distribution %d Average time %lf max %lf min %lf\n",distribution,avgt/(double)(trial+1),maxt,mint);
}// distribution

	delete [] S;
	QueryPerformanceCounter(&stopt);
	ct=((double)stopt.LowPart-(double)startt.LowPart)/(double)freq.LowPart;
	printf("All tests performed in %lf seconds.\n",ct);

	return 0;
}