#import "VoronoiView.h"
#import "Thinker.h"
#import <dpsclient/wraps.h>
#import <math.h>
#import <appkit/color.h>
#import <appkit/nextstd.h>
#import <libc.h>
#import <time.h>

#include "voronoi.h"   /* for the voronoi code */
#include "headers.h" /* more voronoi code */

/*
 * how many pixels wide should the sites and polygon edges be?
 */
#define WIDTH 2.0

/*
 * what is the maximum number of sites to be generated on any iteration?
 */
#define BS_MAXSITES 100

/* 
 * length of idle period in seconds
 */
#define IDLESECONDS 30

@implementation VoronoiView

/*
 * allocate a random number of sites (2D points) around which the
 * Voronoi diagram will be constructed.
 */

- initFrame:(NXRect *)frameRect
{
  self=[super initFrame:frameRect];
  nTotal=nLeft=0;
  sites=(NXPoint *)0;
  srandom(time(0));
  return self;
}

- setUpVoronoi
{
  NXRect r;
  int i;
  [self getBounds:&r];
  if (sites) NX_FREE(sites);
  nTotal=10+random()%(BS_MAXSITES-10);
  nLeft=nTotal-1;
  NX_MALLOC(sites,NXPoint,nTotal);
  for(i=0;i<nTotal;i++) {
     sites[i].x = 20 + r.origin.x + random()%((int)(r.size.width-40));
     sites[i].y = 20 + r.origin.y + random()%((int)(r.size.height-40));
     };
  /* set up voronoi diagram */
  if (load_vsites(nTotal,sites,r.origin.x+5,r.origin.y+5,
                  r.origin.x+r.size.width-5,r.origin.y+r.size.height-5)) {
    nLeft=-1;
    };
  return self;
}

- drawPolygon:(int)index
{
  NXPoint *temp;
  int nv,i;
  NX_MALLOC(temp,NXPoint,nTotal);
  nv=find_vregion(index,temp);
  if (nv < 1) {
    fprintf(stderr,"find_vregion %d failed\n",nLeft);
    nLeft=-1;
    return self;
    };
  NXSetColor(NXConvertHSBToColor(1.0-(1.0*index)/(nTotal-1.0),1.0,1.0));
  /* move the vertices of the voronoi polygon toward the site by 2%.
     This allows edges drawn multiple times to show up in multiple
     colors. */
  for(i=0;i<nv;i++) {
    temp[i].x += 0.02 * (sites[index].x - temp[i].x);
    temp[i].y += 0.02 * (sites[index].y - temp[i].y);
    };
  PSsetlinewidth(WIDTH);
  PSnewpath();
  PSmoveto(temp[0].x,temp[0].y);
  for(i=1;i<nv;i++) PSlineto(temp[i].x,temp[i].y);
  PSclosepath();
  PSstroke(); /* or even PSfill(); */
  NX_FREE(temp);
  return self;
}

- drawSites
{
  int i;
  /* draw the sites */
  for(i=0;i<nTotal;i++) {
    PSnewpath();
    NXSetColor(NXConvertHSBToColor((1.0*i)/(nTotal-1.0),1.0,1.0));
    PSarc(sites[i].x,sites[i].y,WIDTH,0.0,360.0);
    PSfill();/* or even PSstroke();*/
    };
  return self;
}


- drawSelf:(const NXRect *)rects :(int)rectCount
{
  int i;
  if (!rects || !rectCount) return self;
  PSsetgray(0);
  NXRectFill(rects);

  if (!sites) [self setUpVoronoi];

  [self drawSites];

  /* if some voronoi polygons have already been drawn, redraw them.
     this can happen if a window is moved around during or after the
     `active' phase of the oneStep method. */
  if (nLeft<nTotal-1) 
    for(i=nTotal-1;((i>nLeft)&&(i>=0));--i) [self drawPolygon:i];

  return self;
}

/* each time through oneStep, we draw one Voronoi polygon, unless they've
   already all been drawn, in which case we simply sleep for 1/10 second.
   note: nLeft does double duty.  When non-negative, it's the index of the
   polygon to draw (umm, it's actually the index of the SITE whose polygon
   is to be drawn, but you get what I mean.)  If negative, it counts the number
   of idle trips.  When it gets sufficiently negative, we restart the whole
   works (i.e. new sites, new Voronoi polygons, etc.)
*/
- oneStep
{
  if (nLeft>= 0)
    [self drawPolygon:nLeft--];
  else {
    if (nLeft-- <= -(IDLESECONDS*10)) {
      [[self setUpVoronoi] display];
      return self;
      };
    usleep(100000); /* sleep 0.1 second */
    };
  return self;
}

- (BOOL)useBufferedWindow {return NO;}

/**********************************
 **********************************
 *** voronoi code begins here
 *** This code was written by Seth Teller (at SGI?) and is
 *** available from many archive sites.  I squished it into
 *** one source file because I'm lazy.
 **********************************
 **********************************/

PQinsert(he, v, offset)
struct Halfedge *he;
struct Site *v;
float 	offset;
{
    struct Halfedge *last, *next;

    he->vertex = v;
    ref(v);
    he->ystar = v->coord.y + offset;
    last = &PQhash[PQbucket(he)];
    while ((next = last->PQnext) != (struct Halfedge *) NULL &&
      (he->ystar  > next->ystar  ||
	(he->ystar == next->ystar && v->coord.x > next->vertex->coord.x))) {
	last = next;
    }
    he->PQnext = last->PQnext; 
    last->PQnext = he;
    PQcount += 1;
}

PQdelete(he)
struct Halfedge *he;
{
    struct Halfedge *last;

    if(he-> vertex != (struct Site *) NULL) {	
	last = &PQhash[PQbucket(he)];
	while (last->PQnext != he) 
	    last = last->PQnext;
	last->PQnext = he->PQnext;
	PQcount -= 1;
	deref(he->vertex);
	he->vertex = (struct Site *) NULL;
    }
}

int PQbucket(he)
struct Halfedge *he;
{
    int bucket;

    bucket = (he->ystar - ymin)/deltay * PQhashsize;
    if (bucket<0) 
	bucket = 0;
    if (bucket>=PQhashsize) 
	bucket = PQhashsize-1 ;
    if (bucket < PQmin) 
	PQmin = bucket;
    return(bucket);
}

int PQempty()
{
    return(PQcount==0);
}

struct Point PQ_min()
{
    struct Point answer;

    while(PQhash[PQmin].PQnext == (struct Halfedge *)NULL) 
	PQmin += 1;
    answer.x = PQhash[PQmin].PQnext->vertex->coord.x;
    answer.y = PQhash[PQmin].PQnext->ystar;
    return (answer);
}

struct Halfedge *PQextractmin()
{
    struct Halfedge *curr;

    curr = PQhash[PQmin].PQnext;
    PQhash[PQmin].PQnext = curr->PQnext;
    PQcount -= 1;
    return(curr);
}

PQinitialize()
{
    int i; struct Point *s;

    PQcount = 0;
    PQmin = 0;
    PQhashsize = 4 * sqrt_nsites;
    PQhash = (struct Halfedge *) myalloc(PQhashsize * sizeof (struct Halfedge));
    for(i=0; i<PQhashsize; i+=1) 
	PQhash[i].PQnext = (struct Halfedge *)NULL;
}



int ntry, totalsearch;

ELinitialize()
{
    int i;

    freeinit(&hfl, sizeof **ELhash);
    ELhashsize = 2 * sqrt_nsites;
    ELhash = (struct Halfedge **) myalloc ( sizeof (*ELhash) * ELhashsize);
    for(i=0; i<ELhashsize; i +=1) 
	ELhash[i] = (struct Halfedge *)NULL;
    ELleftend = HEcreate( (struct Edge *)NULL, 0);
    ELrightend = HEcreate( (struct Edge *)NULL, 0);
    ELleftend->ELleft = (struct Halfedge *)NULL;
    ELleftend->ELright = ELrightend;
    ELrightend->ELleft = ELleftend;
    ELrightend->ELright = (struct Halfedge *)NULL;
    ELhash[0] = ELleftend;
    ELhash[ELhashsize-1] = ELrightend;
}

struct Halfedge *HEcreate(e, pm)
struct Edge *e;
int pm;
{
    struct Halfedge *answer;

    answer = (struct Halfedge *) getfree(&hfl);
    answer->ELedge = e;
    answer->ELpm = pm;
    answer->PQnext = (struct Halfedge *) NULL;
    answer->vertex = (struct Site *) NULL;
    answer->ELrefcnt = 0;
    return(answer);
}

ELinsert(lb, new)
struct	Halfedge *lb, *new;
{
    new->ELleft = lb;
    new->ELright = lb->ELright;
    (lb->ELright)->ELleft = new;
    lb->ELright = new;
}

/* Get entry from hash table, pruning any deleted nodes */
struct Halfedge *ELgethash(b)
int b;
{
    struct Halfedge *he;

    if(b<0 || b>=ELhashsize) 
	return((struct Halfedge *) NULL);
    he = ELhash[b]; 
    if (he == (struct Halfedge *)NULL || he->ELedge!=(struct Edge *)DELETED) 
	return (he);

/* Hash table points to deleted half edge.  Patch as necessary. */
    ELhash[b] = (struct Halfedge *) NULL;
    if ((he->ELrefcnt -= 1) == 0) makefree(he, &hfl);
	return ((struct Halfedge *) NULL);
}	

struct Halfedge *ELleftbnd(p)
struct Point *p;
{
    int i, bucket;
    struct Halfedge *he;

/* Use hash table to get close to desired halfedge */
    bucket = (p->x - xmin)/deltax * ELhashsize;
    if(bucket<0) 
	bucket =0;
    if(bucket>=ELhashsize) 
	bucket = ELhashsize - 1;
    he = ELgethash(bucket);
    if(he == (struct Halfedge *) NULL) {   
	for(i=1; 1 ; i += 1) {	
	    if ((he=ELgethash(bucket-i)) != (struct Halfedge *) NULL) 
		break;
	    if ((he=ELgethash(bucket+i)) != (struct Halfedge *) NULL) 
		break;
	}
	totalsearch += i;
    }
    ntry += 1;

/* Now search linear list of halfedges for the corect one */
    if (he==ELleftend  || (he != ELrightend && right_of(he,p))) {
	do {
	    he = he->ELright;
	} while (he!=ELrightend && right_of(he,p));
     	he = he->ELleft;
    } else 
	do {
	    he = he->ELleft;
	} while (he!=ELleftend && !right_of(he,p));

/* Update hash table and reference counts */
    if(bucket > 0 && bucket <ELhashsize-1) {	
	if(ELhash[bucket] != (struct Halfedge *) NULL) 
	    ELhash[bucket]->ELrefcnt -= 1;
	ELhash[bucket] = he;
	ELhash[bucket]->ELrefcnt += 1;
    }
    return (he);
}

/*
 *  This delete routine can't reclaim node, since pointers from hash
 *  table may be present.   
 */
ELdelete(he)
struct Halfedge *he;
{
    (he->ELleft)->ELright = he->ELright;
    (he->ELright)->ELleft = he->ELleft;
    he->ELedge = (struct Edge *)DELETED;
}

struct Halfedge *ELright(he)
struct Halfedge *he;
{
    return (he->ELright);
}

struct Halfedge *ELleft(he)
struct Halfedge *he;
{
    return (he->ELleft);
}

struct Site *leftreg(he)
struct Halfedge *he;
{
    if(he->ELedge == (struct Edge *)NULL) 
	return(bottomsite);
    return( he->ELpm == le ? he->ELedge->reg[le] : he->ELedge->reg[re]);
}

struct Site *rightreg(he)
struct Halfedge *he;
{
    if(he->ELedge == (struct Edge *)NULL) 
	return(bottomsite);
    return( he->ELpm == le ? he->ELedge->reg[re] : he->ELedge->reg[le]);
}


geominit()
{
    struct Edge e;
    float sn;

    freeinit (&efl, sizeof (struct Edge));
    nvertices = 0;
    nedges = 0;
    sn = nsites + 4;
    sqrt_nsites = sqrt(sn);
    deltay = ymax - ymin;
    deltax = xmax - xmin;
    siteidx = 0;
}

struct Edge *bisect(s1,s2)
struct	Site *s1,*s2;
{
    float dx,dy,adx,ady;
    struct Edge *newedge;

    newedge = (struct Edge *) getfree(&efl);

    newedge->reg[0] = s1;
    newedge->reg[1] = s2;
    ref(s1); 
    ref(s2);
    newedge->ep[0] = (struct Site *) NULL;
    newedge->ep[1] = (struct Site *) NULL;

    dx = s2->coord.x - s1->coord.x;
    dy = s2->coord.y - s1->coord.y;
    adx = dx>0 ? dx : -dx;
    ady = dy>0 ? dy : -dy;
    newedge->c = s1->coord.x * dx + s1->coord.y * dy + (dx*dx + dy*dy)*0.5;
    if (adx>ady) {	
	newedge->a = 1.0; 
	newedge->b = dy/dx; 
	newedge->c /= dx;
    } else {	
	newedge->b = 1.0; 
	newedge->a = dx/dy; 
	newedge->c /= dy;
    }
    newedge->edgenbr = nedges;
    out_bisector(newedge);
    nedges += 1;
    return(newedge);
}

struct Site *intersect(el1, el2, p)
struct Halfedge *el1, *el2;
struct Point *p;
{
    struct Edge *e1,*e2, *e;
    struct Halfedge *el;
    float d, xint, yint;
    int right_of_site;
    struct Site *v;

    e1 = el1->ELedge;
    e2 = el2->ELedge;
    if (e1 == (struct Edge*)NULL || e2 == (struct Edge*)NULL) 
	return ((struct Site *) NULL);
    if (e1->reg[1] == e2->reg[1]) 
	return ((struct Site *) NULL);
    d = e1->a * e2->b - e1->b * e2->a;
    if (-1.0e-10<d && d<1.0e-10) 
	return ((struct Site *) NULL);
    xint = (e1->c*e2->b - e2->c*e1->b)/d;
    yint = (e2->c*e1->a - e1->c*e2->a)/d;
    if ((e1->reg[1]->coord.y < e2->reg[1]->coord.y) ||
	    (e1->reg[1]->coord.y == e2->reg[1]->coord.y &&
		e1->reg[1]->coord.x < e2->reg[1]->coord.x) ) {	
	el = el1; 
	e = e1;
    } else {	
	el = el2; 
	e = e2;
    }
    right_of_site = xint >= e->reg[1]->coord.x;
    if ((right_of_site && el->ELpm == le) ||
	   (!right_of_site && el->ELpm == re)) 
	return ((struct Site *) NULL);
    v = (struct Site *) getfree(&sfl);
    v->refcnt = 0;
    v->coord.x = xint;
    v->coord.y = yint;
    return(v);
}

/* returns 1 if p is to right of halfedge e */
int right_of(el, p)
struct Halfedge *el;
struct Point *p;
{
    struct Edge *e;
    struct Site *topsite;
    int right_of_site, above, fast;
    float dxp, dyp, dxs, t1, t2, t3, yl;

    e = el->ELedge;
    topsite = e->reg[1];
    right_of_site = p->x > topsite->coord.x;
    if (right_of_site && el->ELpm == le) 
	return(1);
    if (!right_of_site && el->ELpm == re) 
	return (0);
    if (e->a == 1.0) {	
	dyp = p->y - topsite->coord.y;
	dxp = p->x - topsite->coord.x;
	fast = 0;
	if ((!right_of_site && e->b<0.0) || (right_of_site && e->b>=0.0) ) {	
	    above = dyp>= e->b*dxp;	
	    fast = above;
	} else {	
	    above = p->x + p->y*e->b > e->c;
	    if (e->b<0.0) 
		above = !above;
	    if (!above) 
		fast = 1;
	}
	if (!fast) {	
	    dxs = topsite->coord.x - (e->reg[0])->coord.x;
	    above = e->b * (dxp*dxp - dyp*dyp) <
				    dxs*dyp*(1.0+2.0*dxp/dxs + e->b*e->b);
	    if (e->b<0.0) 
		above = !above;
	}
    } else {	/*e->b==1.0 */
	yl = e->c - e->a*p->x;
	t1 = p->y - yl;
	t2 = p->x - topsite->coord.x;
	t3 = yl - topsite->coord.y;
	above = t1*t1 > t2*t2 + t3*t3;
    }
    return (el->ELpm==le ? above : !above);
}

void
v_endpoint(e, lr, s)
struct Edge *e;
int	lr;
struct Site *s;
{
    e->ep[lr] = s;
    ref(s);
    if (e->ep[re-lr]==(struct Site *)NULL) 
	return;
    out_ep(e);
    deref(e->reg[le]);
    deref(e->reg[re]);
    makefree(e, &efl);
}

float dist(s,t)
struct Site *s,*t;
{
    float dx,dy;
    dx = s->coord.x - t->coord.x;
    dy = s->coord.y - t->coord.y;
    return(sqrt(dx*dx + dy*dy));
}


int makevertex(v)
struct Site *v;
{
    v->sitenbr = nvertices;
    nvertices += 1;
    out_vertex(v);
}


deref(v)
struct	Site *v;
{
    v->refcnt -= 1;
    if (v->refcnt == 0) 
	makefree(v, &sfl);
}

ref(v)
struct Site *v;
{
    v->refcnt += 1;
}


freeinit(fl, size)
struct	Freelist *fl;
int	size;
{
    fl->head = (struct Freenode *) NULL;
    fl->nodesize = size;
}

char *getfree(fl)
struct	Freelist *fl;
{
    int i; struct Freenode *t;

    if(fl->head == (struct Freenode *) NULL) {	
	t = (struct Freenode *)myalloc(sqrt_nsites * fl->nodesize);
	for(i=0; i<sqrt_nsites; i+=1) 	
	    makefree((struct Freenode *)((char *)t+i*fl->nodesize), fl);
    }
    t = fl->head;
    fl->head = (fl->head)->nextfree;
    return((char *)t);
}

makefree(curr,fl)
struct Freenode *curr;
struct Freelist *fl;
{
    curr->nextfree = fl->head;
    fl->head = curr;
}

#define MAXALLOCS (8 * MAXSITES)
int total_alloc, num_allocs;
char *locations[MAXALLOCS];

char *myalloc(n)
unsigned n;
{
    char *t;
    if ((t=malloc(n)) == (char *) 0) {
	fprintf(stderr,"out of memory processing site %d (%d bytes in use)\n",
		siteidx, total_alloc);
    	exit(-1);
    }
    total_alloc += n;
    locations[num_allocs++] = t;
    return(t);
}

void 
myfreeall()
{
int k;

    for (k = 0; k < num_allocs; k++)
	free (locations [k]);
/*
    printf ("freed %d allocs (%d bytes)\n", num_allocs, total_alloc);
*/
    num_allocs = total_alloc = 0;
}


int triangulate = FALSE, sorted, plot, debug = FALSE, randmode, outputmode;
float xmin, xmax, ymin, ymax, deltax, deltay;

char execdir[80];

struct Site _sites[MAXSITES], *sites = _sites;

int nsites;
int siteidx;
int sqrt_nsites;
int nvertices;
struct Freelist sfl;
struct Site *bottomsite;

int nedges;
struct	Freelist efl;

struct Freelist	hfl;
struct Halfedge *ELleftend, *ELrightend;
int ELhashsize;
struct Halfedge **ELhash;

int PQhashsize;
struct Halfedge *PQhash;
int PQcount;
int PQmin;

/* sort sites on y, then x, coord */
int scomp(s1,s2)
struct Point *s1,*s2;
{
    if(s1->y < s2->y) return(-1);
    if(s1->y > s2->y) return(1);
    if(s1->x < s2->x) return(-1);
    if(s1->x > s2->x) return(1);
    return(0);
}

/* sort GLsites on y, then x, coord */
int GLscomp(s1,s2)
VERT *s1,*s2;
{
    if(s1->y < s2->y) return(-1);
    if(s1->y > s2->y) return(1);
    if(s1->x < s2->x) return(-1);
    if(s1->x > s2->x) return(1);
    return(0);
}

/* return a single in-storage site */
struct Site *
nextone()
{
    if(siteidx < nsites)
	return (&sites[siteidx++]);
    else
	return( (struct Site *)NULL);
}

sortGLsites()
{
float vx, vy;
int k;
char *a = NULL;
float randr(), dx, dy;

    qsort((char *)GLsites, nsites, sizeof (VERT), GLscomp);

    /* check: are there coincident points? */
    for (k = 1; k < nsites; k++)
	if ((GLsites[k-1].y == GLsites[k].y) && (GLsites[k-1].x == GLsites[k].x))   {
/*	    printf ("coincident sites at %d, %d!\n", k-1, k);*/
	    dx = GLsites[k].x * (1.0 / 4096.0);
	    dy = GLsites[k].y * (1.0 / 4096.0);
	    /* hack: jitter the last point randomly in its last significant digit */
	    GLsites[k-1].x += randr (-dx, dx);
	    GLsites[k-1].y += randr (-dy, dy);
	    }
}

/* implicit parameters: nsites, sqrt_nsites, xmin, xmax, ymin, ymax,
   deltax, deltay (can all be estimates).
   Performance suffers if they are wrong; better to make nsites,
   deltax, and deltay too big than too small.  (?) */
voronoi (nextsite)
struct Site *(*nextsite)();
{
struct Site *newsite, *bot, *top, *temp, *p;
struct Site *v;
struct Point newintstar;
int pm;
struct Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector;
struct Edge *e;

    myfreeall ();

    if (nsites <= 1)
	return;

    freeinit(&sfl, sizeof (struct Site));

    /* bboxinit();   copies static array into nsites */
    geominit();	    /* internal use of deltax, deltay  */

    PQinitialize();
    bottomsite = (*nextsite)();
    out_site(bottomsite);
    ELinitialize();

    newsite = (*nextsite)();
    while(1) {
	if(!PQempty()) newintstar = PQ_min();
/* new site is smallest */
	if (newsite != (struct Site *)NULL && (PQempty() 
		 || newsite->coord.y < newintstar.y
			 || (newsite->coord.y == newintstar.y 
			     && newsite->coord.x < newintstar.x))) {
	    out_site(newsite);
	    lbnd = ELleftbnd(&(newsite->coord));
	    rbnd = ELright(lbnd);
	    bot = rightreg(lbnd);
	    e = bisect(bot, newsite);
	    bisector = HEcreate(e, le);
	    ELinsert(lbnd, bisector);
	    if ((p = intersect(lbnd, bisector)) != (struct Site *) NULL) {
		PQdelete(lbnd);
		PQinsert(lbnd, p, dist(p,newsite));
	    }
	    lbnd = bisector;
	    bisector = HEcreate(e, re);
	    ELinsert(lbnd, bisector);
	    if ((p = intersect(bisector, rbnd)) != (struct Site *) NULL) {
		PQinsert(bisector, p, dist(p,newsite));	
	    }
	    newsite = (*nextsite)();	
	} else if (!PQempty()) {
	/* intersection is smallest */
	    lbnd = PQextractmin();
	    llbnd = ELleft(lbnd);
	    rbnd = ELright(lbnd);
	    rrbnd = ELright(rbnd);
	    bot = leftreg(lbnd);
	    top = rightreg(rbnd);
	    out_triple(bot, top, rightreg(lbnd));
	    v = lbnd->vertex;
	    makevertex(v);
	    v_endpoint(lbnd->ELedge,lbnd->ELpm,v);
	    v_endpoint(rbnd->ELedge,rbnd->ELpm,v);
	    ELdelete(lbnd); 
	    PQdelete(rbnd);
	    ELdelete(rbnd); 
	    pm = le;
	    if (bot->coord.y > top->coord.y) {	
		temp = bot; 
		bot = top; 
		top = temp; 
		pm = re;
	    }
	    e = bisect(bot, top);
	    bisector = HEcreate(e, pm);
	    ELinsert(llbnd, bisector);
	    v_endpoint(e, re-pm, v);
	    deref(v);
	    if((p = intersect(llbnd, bisector)) != (struct Site *) NULL) {
		PQdelete(llbnd);
		PQinsert(llbnd, p, dist(p,bot));
	    }
	    if((p = intersect(bisector, rrbnd)) != (struct Site *) NULL) 
		PQinsert(bisector, p, dist(p,bot));
	} else 
	    break;
    }

    for(lbnd=ELright(ELleftend); lbnd != ELrightend; lbnd=ELright(lbnd)) {
	e = lbnd->ELedge;
	out_ep(e);
        }
}


#define LBND 0
#define BBND 1
#define RBND 2
#define TBND 3

/* left, right, top bottom */
#define NBOUND(x,y) (x >= Pxmax ? RBND : (x <= Pxmin ? LBND : (y >= Pymax ? TBND : (y <= Pymin ? BBND : -1))))

#define BBOX_THERE() ((pxmin == txmin) && (pxmax == txmax) && (pymin == tymin) && (pymax && pymax))
#define BBOX_NOT_THERE() (!BBOX_THERE())

static float Pxmin, Pxmax, Pymin, Pymax;	/* floating point tolerance values */
static float pxmin = 1E10, pxmax = -1E10, pymin = 1E10, pymax = -1E10;
/* space for theta after x*y* */
static float xminymin[3], xminymax[3], xmaxymax[3], xmaxymin[3];
/* these point to the x*y* */
static VERTTHETA *corner[4];

static POINT v, v1, v2, v3, v4;

static int numverts = 0, numvedges = 0, numtris;

/* sites[i].coord.x,y defined and sorted in main() */

VERT   GLsites[MAXVERTS];
static VERT verts[MAXVERTS];
static EDGE vedges[MAXEDGES];
static TRI  tris[MAXTRIS];
static int  inverse[MAXVERTS];

float randr(low,high)
float low, high;
{
    float value;

    value = low + (random()/((1<<15)-1.0))*(high-low);
    return( value );
}

out_site(s)
struct Site *s;
{
}

out_bisector(e)
struct Edge *e;
{
}

out_ep(e)
struct Edge *e;
{
    if(!triangulate)
	clip_line(e);

    if(debug)	{
	printf("out_ep(): edge %d", e->edgenbr);
	printf(" ep %d", e->ep[le] != (struct Site *)NULL ? e->ep[le]->sitenbr : -1);
	printf(" ep %d", e->ep[re] != (struct Site *)NULL ? e->ep[re]->sitenbr : -1);
	printf(" reg %d", e->reg[le] != (struct Site *)NULL ? e->reg[le]->sitenbr : -1);
	printf(" reg %d\n", e->reg[re] != (struct Site *)NULL ? e->reg[re]->sitenbr : -1);
	}
}

out_vertex(v)
struct Site *v;
{
    if (!triangulate)	{
	verts[numverts].x = v->coord.x;
	verts[numverts].y = v->coord.y;
	if (numverts < MAXVERTS)
	    numverts++;
	else	{
	    fprintf (stderr, "\nvert list overflow!");
	    exit (-1);
	    }
	}

    if(debug)
	printf("vertex(%d) at %10.7f %10.7f\n", v->sitenbr, v->coord.x, v->coord.y);
}

out_triple(s1, s2, s3)
struct Site *s1, *s2, *s3;
{
    if(triangulate) {
	tris[numtris].v1 = (VERT *)&(sites[s1->sitenbr].coord);
	tris[numtris].v2 = (VERT *)&(sites[s2->sitenbr].coord);
	tris[numtris].v3 = (VERT *)&(sites[s3->sitenbr].coord);

	if (numtris < MAXTRIS)
	    numtris++;
	else	{
	    fprintf (stderr, "out_triple(): triangle list overflow!\n");
	    exit (-1);
	    }
	}
}

vdinit()
{
    numverts = numvedges = numtris = 0;
}

/* here, want to copy the gl sites INTO the voronoi array */
/* gl sites are written exactly once at beginning of time */
bboxinit()
{
int k;
float dx, dy, x, y;

    /* get tight bounding box */
    xmin =  1e10;
    xmax = -1e10;

    ymin =  1e10;
    ymax = -1e10;

    for(k = 0; k < nsites; k++) {	
	x = sites[k].coord.x = GLsites[k].x;
	y = sites[k].coord.y = GLsites[k].y;

        sites[k].refcnt = 0;
	sites[k].sitenbr = k;

	if (x < xmin) xmin = x;
	if (y < ymin) ymin = y;

	if (x > xmax) xmax = x;
	if (y > ymax) ymax = y;
	}

    /* NOW: xmin, ymin, xmax, ymax EXACT, as required by voronoi() */
    /* we shall fool with pxmin, pymin, pxmax, pymax, as used by clip() */

#define EPSILON 1.0

    /* compute 'loose' bounding box */
    dx = xmax - xmin;
    dx = MAX (dx, EPSILON);

    dy = ymax - ymin;
    dy = MAX (dy, EPSILON);

    pxmin = xmin - dx * 0.25;
    pymin = ymin - dy * 0.25;

    pxmax = xmax + dx * 0.25;
    pymax = ymax + dy * 0.25;

#ifdef STUPID_C_COMPILER
    printf ("/* xmin, ymin %10.7f %10.7f; xmax, ymax %10.7f %10.7f; */\n", xmin, ymin, xmax, ymax);
    printf ("/* pxmin, pymin %10.7f %10.7f; pxmax, pymax %10.7f %10.7f; crad %10.7f; */\n", pxmin, pymin, pxmax, pymax, cradius);
#endif
}

int clip_line(e)
struct Edge *e;
{
struct Site *s1, *s2;
struct Site *r1, *r2;
float x1,x2,y1,y2;

    if(e->a == 1.0 && e->b >= 0.0) {	
	s1 = e->ep[1];
	s2 = e->ep[0];

	r1 = e->reg[1];
	r2 = e->reg[0];
    } else {	
	s1 = e->ep[0];
	s2 = e->ep[1];

	r1 = e->reg[0];
	r2 = e->reg[1];
    }

    if(e->a == 1.0) {
	y1 = pymin;
	if (s1!=(struct Site *)NULL && s1->coord.y > pymin)
	    y1 = s1->coord.y;
	if(y1>pymax)
	    return;
	x1 = e->c - e->b * y1;
	y2 = pymax;
	if (s2!=(struct Site *)NULL && s2->coord.y < pymax)
	    y2 = s2->coord.y;
	if(y2<pymin)
	    return(0);
	x2 = e->c - e->b * y2;
	if ((x1> pxmax && x2>pxmax) || (x1<pxmin && x2<pxmin)) 
	    return;
	if(x1> pxmax) {
	    x1 = pxmax;
	    y1 = (e->c - x1)/e->b;
	}
	if(x1<pxmin) {
	    x1 = pxmin;
	    y1 = (e->c - x1)/e->b;
	}
	if(x2>pxmax) {
	    x2 = pxmax;
	    y2 = (e->c - x2)/e->b;
	}
	if(x2<pxmin) {
	    x2 = pxmin;
	    y2 = (e->c - x2)/e->b;
	}
    } else {
	x1 = pxmin;
	if (s1!=(struct Site *)NULL && s1->coord.x > pxmin) 
	    x1 = s1->coord.x;
	if(x1>pxmax)
	    return(0);
	y1 = e->c - e->a * x1;
	x2 = pxmax;
	if (s2!=(struct Site *)NULL && s2->coord.x < pxmax)
	    x2 = s2->coord.x;
	if(x2<pxmin)
	    return(0);
	y2 = e->c - e->a * x2;
	if ((y1> pymax && y2>pymax) || (y1<pymin && y2<pymin))
	    return(0);
	if(y1> pymax) {
	    y1 = pymax;
	    x1 = (e->c - y1)/e->a;
	}
	if(y1<pymin) {
	    y1 = pymin;
	    x1 = (e->c - y1)/e->a;
	}
	if(y2>pymax) {
	    y2 = pymax;
	    x2 = (e->c - y2)/e->a;
	}
	if(y2<pymin) {
	    y2 = pymin;
	    x2 = (e->c - y2)/e->a;
	}
    }

/*  fprintf (stderr, "clip_line(): edge %d is (%10.7f, %10.7f) to (%10.7f, %10.7f)\n", numvedges, x1, y1, x2, y2);*/

    if (!triangulate)	{
	vedges[numvedges].x1 = x1;
	vedges[numvedges].y1 = y1;

	vedges[numvedges].x2 = x2;
	vedges[numvedges].y2 = y2;

	vedges[numvedges].nbr1 = (r1 != NULL ? r1->sitenbr : -9998);
	vedges[numvedges].nbr2 = (r2 != NULL ? r2->sitenbr : -9999);

	if (r1 != NULL && r2 != NULL)	{
	    vedges[numvedges].xm = AVG (r1->coord.x, r2->coord.x);
	    vedges[numvedges].ym = AVG (r1->coord.y, r2->coord.y);
	    }

	if (debug)
	    printf ("clip_line puts edge induced by %d and %d\n", r1->sitenbr, r2->sitenbr);

	if (numvedges < MAXEDGES)
	    numvedges++;
	else	{
	    fprintf (stderr, "clip_line(): edge list overflow!\n");
	    exit (-1);
	    }
	}
}

/* load_vsites(): 
    accept the n voronoi sites (x_n, y_n)
    calculate and store the voronoi diagram over the n sites,
    clipping all infinite edges to bbox: [xmin, ymin, xmax, ymax].
    
    note: if (xmin,ymin,xmax,ymax are all == 0.0), OR
	if these do not enclose the data, a bounding box
	 will be computed over the input.

    returns:
	-1 if error
	 0 otherwise
*/
int
load_vsites (n, usites, uxmin, uymin, uxmax, uymax)
int n;
float usites[][2];	/* sites in x,y order */
float uxmin, uymin, uxmax, uymax;
{
int k, compute_bbox, sid, tid;
float dx, dy, x, y;

    if (n >= MAXSITES)	{
	fprintf (stderr, "load_vsites(): can't handle >= %d sites.\n", MAXSITES);
	return (-1);
	}

    compute_bbox = (uxmin == 0.0) && (uymin == 0.0) && (uxmax == 0.0) && (uymax == 0.0);

    /* copy the sites into GLsites and set global nsites */
    for (k = 0; k < n; k++) {
	GLsites[k].x = usites[k][0];
	GLsites[k].y = usites[k][1];
	GLsites[k].pid = k;
	}

    nsites = n;

    /* sort GL sites lexicographically by position */
    sortGLsites();

    /* copy sorted GLsites into voronoi alg sites */
    bboxinit();

    /* now, if user punted on bbox calculation, OR if user bbox does not truly enclose
       user data, we use our bbox (computed in initbbox).  otherwise we take the user's. */
    if (!(compute_bbox || uxmin > xmin || uymin > ymin || uxmax < xmax || uymax < ymax))	{
	pxmin = uxmin;
	pymin = uymin;

	pxmax = uxmax;
	pymax = uymax;
	}

    xminymax [0] = xminymin[0] = pxmin;
    xminymin [1] = xmaxymin[1] = pymin;

    xmaxymax [0] = xmaxymin[0] = pxmax;
    xmaxymax [1] = xminymax[1] = pymax;

    corner[0] = (VERTTHETA *)xminymin;
    corner[1] = (VERTTHETA *)xmaxymin;
    corner[2] = (VERTTHETA *)xmaxymax;
    corner[3] = (VERTTHETA *)xminymax;

    /* now: set the floating point tolerance values P*** to be 1 or 2 significant bits inside the p*** */
    /* be careful to use RELATIVE error; that is, to scale the tolerance to the bbox values themselves */
    /* now, if some user puts points way out in left field, our technique handles the ranges correctly */
    {
    float dx = (pxmax - pxmin) * (1.0 / (4096.0));      /* twelve binary digits out */
    float dy = (pymax - pymin) * (1.0 / (4096.0));      /* twelve binary digits out */

    Pxmin = pxmin + dx;
    Pxmax = pxmax - dx;

    Pymin = pymin + dy;
    Pymax = pymax - dy;
    }

    /* compute inverse of external->internal sid mapping */
    for (sid = 0; sid < nsites; sid++)
	inverse[sid] = GLsites[sid].pid;

    /* zero list lengths out */
    vdinit ();

    /* run the voronoi code, no triangulate */
    triangulate = FALSE;
    voronoi (nextone);

    /* RE-copy sorted GLsites into voronoi alg sites */
    bboxinit();

    /* run the voronoi code, do triangulate */
    triangulate = TRUE;
    voronoi (nextone);

    /* invert the integer id's in sites[], for find_dtriangles() */
    /* and restore the original vertex values (from GLsites)	 */
    for (sid = 0; sid < nsites; sid++)	{
        sites[sid].sitenbr = GLsites[sid].pid;
	sites[sid].coord.x = GLsites[sid].x;
	sites[sid].coord.y = GLsites[sid].y;
	}

    return (0);
}

static VERTTHETA vtlist[1024], slist[1024];
static int  vtnum;

int
vtcomp (vt1, vt2)
VERTTHETA *vt1, *vt2;
{
    if (vt1->theta < vt2->theta)
	return (-1);
    else if (vt1->theta > vt2->theta)
	return (1);
    else
	return (0);
}	

/*
    find_vregion(sid, plen, pverts)
	given a site id 'sid' from 0..nsites-1 inclusive, 
	returns the voronoi polygon associated with that site
	in the array 'pverts', and its length on call stack.

	the vertices are returned in counterclockwise order.

    returns:
	-1 if error condition
	plen > 2 [i.e., the # of verts in the voronoi polygon] otherwise
*/
int
find_vregion (vsid, pverts)
int vsid;
float pverts[][2];
{
int sid, b, k, vnum, bnd1, bnd2, bdiff, sleft, lag, lead, p1a, p1b, p2a, p2b;
float x1, y1, x2, y2, theta1, theta2, lasttheta, dtheta, h;

/* note that PUTV(v,sid,vid) has the side effect of incrementing vid */
#define PUTV(v,sid,vid)	{ pverts[vid][0] = (v)->x; pverts[vid][1] = (v)->y; vid++; }

    if (vsid < 0 || vsid >= nsites) {
	fprintf (stderr, "find_vregion(%d) called with illegal site id.\n", vsid);
	return (-1);
	}

    /* first thing is to convert user's 'virtual' site id to 
       a 'physical' site id, namely, an index into GLsites */
    for (sid = 0; sid < nsites; sid++)
	if (GLsites[sid].pid == vsid)
	    break;

    if (sid == nsites)	{
	fprintf (stderr, "find_vregion(%d) can't find requested site id.\n", vsid);
	return (-1);
	}

    for (k = 0; k < 4; k++)
	corner[k]->theta = atan2 (corner[k]->y - GLsites[sid].y, corner[k]->x - GLsites[sid].x);

    for (vtnum = 0, k = 0; k < numvedges; k++)
	if (vedges[k].nbr1 == sid || vedges[k].nbr2 == sid)   {
	    /* add both ends of the edge, and their thetas, to unsorted list (parent is edge k) */
	    slist[vtnum  ].e1 = slist[vtnum  ].e2 = k;
	    slist[vtnum  ].theta = atan2 (vedges[k].y1 - GLsites[sid].y, vedges[k].x1 - GLsites[sid].x);
	    slist[vtnum  ].x = vedges[k].x1;
	    slist[vtnum++].y = vedges[k].y1;

	    slist[vtnum  ].e1 = slist[vtnum  ].e2 = k;
            slist[vtnum  ].theta = atan2 (vedges[k].y2 - GLsites[sid].y, vedges[k].x2 - GLsites[sid].x);
	    slist[vtnum  ].x = vedges[k].x2;
	    slist[vtnum++].y = vedges[k].y2;
	    }

    /* now we have a list of vtnum thetas and vertices.  sort it on theta */
    qsort ((char *)slist, vtnum, sizeof (VERTTHETA), vtcomp);

    /* next, install the unique slist entries into vtlist */
    lag = lead = 0;
    vtlist[lag] = slist[lead];
    lasttheta = -10.0;

    while (lead < vtnum)	{
        if (fabs (slist[lead].theta - lasttheta) > 1E-4) {
	    lasttheta = slist[lead].theta;
	    vtlist[lag++] = slist[lead++];
	    }
	else    {
	    vtlist[lag-1].e2 = slist[lead++].e1;
	    }
	}

    vtnum = lag;
/*
    printf ("\n");
    for (k = 0; k < vtnum; k++)
        printf ("vtlist[%d]: x,y %f,%f; theta %f; edges %d %d\n",
		k, vtlist[k].x, vtlist[k].y, vtlist[k].theta, vtlist[k].e1, vtlist[k].e2);
*/

    for (vnum = 0, k = 0; k < vtnum; k++)	{
        if (fabs(vtlist[(k + vtnum - 1) % vtnum].theta - vtlist[k].theta) < 1E-4)   {
	    fprintf (stderr,"find_vregion(%d): vtlist %d, %d identical?\n", sid, (k + vtnum - 1) % vtnum, k);
	    return (-1);
	    }

        x1 = vtlist[(k + vtnum - 1) % vtnum].x;
        y1 = vtlist[(k + vtnum - 1) % vtnum].y;
	p1a = vtlist[(k + vtnum - 1) % vtnum].e1;
	p1b = vtlist[(k + vtnum - 1) % vtnum].e2;

	x2 = vtlist[k].x;
	y2 = vtlist[k].y;
	p2a = vtlist[k].e1;
	p2b = vtlist[k].e2;

	/* now: if the two vertices come from the same edge, output and continue */
	if (((p1a == p2a) || (p1a == p2b) || (p1b == p2a) || (p1b == p2b))  && (vtnum > 2))	{
	    PUTV (&vtlist[k], sid, vnum);
	    continue;
	    }

	/* otherwise must fill in missing corners between x1,y1 and x2,y2 */
	bnd1 = NBOUND(x1,y1);
	bnd2 = NBOUND(x2,y2);

	/* find number of CLOCKWISE steps around bbox */
	if (bnd1 >= 0 && bnd2 >= 0)
	    bdiff = ((bnd2 - bnd1) + 4) % 4;	/* from 0 to 3 */
	else
	    bdiff = 0;

	/* if they were on the same bounding box edge, output and continue */
	if (bdiff == 0) {		    
	    PUTV (&vtlist[k], sid, vnum);
	    continue;
	    }

	/* special case: exactly two vertices */
	if (vtnum == 2) {
	    theta1 = vtlist[0].theta;
	    theta2 = vtlist[1].theta;
	    dtheta = theta2 - theta1;    

	    /* add all corners OUTSIDE the angular range of the edge as seen from the site */
	    for (b = 0; b < 4; b++)	{
		if ((dtheta >= M_PI) && (corner[b]->theta > theta1) && (corner[b]->theta < theta2))
		    vtlist[vtnum++] = *corner[b];
		else if ((dtheta < M_PI) && ((corner[b]->theta < theta1) || (corner[b]->theta > theta2)))
		    vtlist[vtnum++] = *corner[b];
		}

	    /* resort the (small) vertex list by theta */
	    qsort ((char *)vtlist, vtnum, sizeof (VERTTHETA), vtcomp);

	    /* and output */
	    for (b = 0; b < vtnum; b++)
		PUTV (&vtlist[b], sid, vnum);

	    break;
	    }

	for (b = 0; b < bdiff; b++)
	    PUTV (corner[(bnd1 + b) % 4], sid, vnum);

	PUTV (&vtlist[k], sid, vnum);
	}	/* k 0 ..vtnum */

    return (vnum);
}

/*
    int
    find_dtriangles (**dtris)

    returns:
	-1 if error condition, *dtris == NULL
	o/wise, # of delaunay triangles, *dtris == array of TRIS (see voronoi.h)
*/
int
find_dtriangles (dtris)
TRI **dtris;
{
    if (numtris <= 0)	{
	*dtris = NULL;
	return (-1);
	}
    else    {
	*dtris = tris;
	return (numtris);
	}
}
@end
