/*
 * delaunayc.c
 *
 *  TRI = DELAUNAYC(X,Y) returns a set of triangles such that no data
 *	points are contained in any triangle's circumcircle. Each row of
 *	the M-by-3 matrix TRI defines one such triangle and contains
 *	indices into the vectors X and Y.
 *
 *  TRI = DELAUNAYC(X,Y,'sorted') assumes that the points X and Y are
 *  sorted first by Y and then by X, and that duplicate points have
 *  already been eliminated.
 *
 *  The Delaunay triangulation is used with: GRIDDATA (to
 *  interpolate scattered data), CONVHULL, VORONOI (to compute the
 *  VORONOI diagram), and is useful by itself to create a triangular
 *  grid for scattered data points.
 *
 *  The functions DSEARCH and TSEARCH search the triangulation to
 *  find nearest neighbor points or enclosing triangles, respectively.
 *	
 *  Based on sweepline Voronoi code by Steven Fortune (see copyright
 *  notice below).
 *	
 *   See also VORONOI, GRIDDATA, CONVHULL, DSEARCH, TSEARCH.
 *
 * This is a mex file for MATLAB.
 *
 * Clay M. Thompson 8-15-95
 * Copyright (c) 1984-98 by The MathWorks, Inc.
 * $Revision: 1.2 $ 
 */
 
/*
 * The author of this software is Steven Fortune.  Copyright (c) 1994 by AT&T
 * Bell Laboratories.
 * Permission to use, copy, modify, and distribute this software for any
 * purpose without fee is hereby granted, provided that this entire notice
 * is included in all copies of any software which is or includes a copy
 * or modification of this software and in all copies of the supporting
 * documentation for such software.
 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
 * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
 */

#include "mex.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>

#ifndef NULL
#define NULL 0
#endif
#define DELETED -2
#define le 0 /* Left edge */
#define re 1 /* Right edge */

/*
 * Structure definitions
 */

struct	Freenode	{
	struct	Freenode	*nextfree;
};

struct	Freelist	{
	struct	Freenode	*head;
	int		nodesize;
};

struct vPoint	{
	double x,y;
};

struct Site	{ /* structure used both for sites and for vertices */
	struct	vPoint	coord;
	int		sitenbr; /* Site number */
	int		refcnt;  /* Reference count */
};

struct Edge	{
	double	a,b,c;   			/* Equation for edge line a*x + b*y = c */
	struct	Site 	*ep[2];		/* End points of edge */
	struct	Site	*reg[2];	/* Original vertex points (left and right) */
	int		edgenbr;			/* Edge number */
};

struct Halfedge {
	struct 	Halfedge	*ELleft, *ELright;
	struct 	Edge	*ELedge;
	int		ELrefcnt;
	int		ELpm;
	struct	Site	*vertex;
	double	ystar;
	struct	Halfedge *PQnext;
};


/*
 * global variables
 */
 
double 	xmin, xmax, ymin, ymax, deltax, deltay;
int 	triangulate, sorted, debug, warned;
int 	ELhashsize;
int 	PQhashsize;
int 	PQcount;
int 	PQmin;
int 	nedges;
int		nsites;
int		siteidx;
int		sqrt_nsites;
int		nvertices;
struct	Freelist	hfl;
struct	Halfedge 	*ELleftend, *ELrightend;
struct	Halfedge 	*PQhash;
struct	Halfedge 	**ELhash;
struct	Freelist 	efl;
struct	Site		*sites;
struct 	Freelist 	sfl;
struct	Site		*bottomsite;

/* Buffer for triangle output */
double 	*ptriangles;  /* triangle data array (output array) */
int		maxtriangles; /* maximum number of triangles */
int		ntriangles;   /* number of triangles */

/*
 * Forward Declarations
 */

/* edgelist */
extern void ELdelete(struct Halfedge *he);
extern void ELinitialize(void);
extern void ELinsert(struct Halfedge *lb, struct Halfedge *new);
extern struct Halfedge *HEcreate(struct Edge *e, int pm);
extern struct Halfedge *ELleftbnd(struct vPoint *p);
extern struct Halfedge *ELright(struct Halfedge *he);
extern struct Halfedge *ELleft(struct Halfedge *he);
extern struct Site *rightreg(struct Halfedge *he);
extern struct Site *leftreg(struct Halfedge *he);

/* geometry */
extern void geominit(void);
extern void deref(struct	Site *v);
extern void ref(struct Site *v);
extern void endpoint(struct Edge *e, int lr, struct Site *s);
extern void makevertex(struct Site *v);
extern struct Site *intersect(struct Halfedge *el1, struct Halfedge *el2);
extern double dist(struct Site *s,struct Site *t);
extern struct Edge *bisect(struct Site *s1,struct Site *s2);
extern int right_of(struct Halfedge *el, struct vPoint *p);

/* heap (Priority Queue) */
extern void PQinitialize(void);
extern void PQdelete(struct Halfedge *he);
extern void PQinsert(struct Halfedge *he, struct Site *v, double offset);
extern struct vPoint PQ_min(void);
extern struct Halfedge *PQextractmin(void);
extern int PQempty(void);

/* memory */
extern void freeinit(struct Freelist *fl,int size);
extern char *getfree(struct Freelist *fl);
extern void makefree(struct Freenode *curr,struct Freelist *fl);
extern char *myalloc(unsigned n);

/* output */
extern void out_bisector(struct Edge *e);
extern void out_ep(struct Edge *e);
extern void out_vertex(struct Site *v);
extern void out_triple(struct Site *s1,struct Site *s2,struct Site *s3);
extern void out_site(struct Site *s);

/* voronoi */
extern void voronoi(int triangulate, struct Site *(*nextsite)());

/*
 * Use mxMalloc to do garbage-collection of dynamically allocated memory
 * so that we don't leak memory 
 */
#define myalloc(n) mxCalloc(n,1)

/* main */
struct Site *nextone(void);

/* sort sites on y, then x, coord */
static int scomp(struct vPoint *s1,struct vPoint *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 
 *
 * Return the next input vertex point (as converted to Sites in mexFunction)
 */
struct Site *nextone(void)
{
	struct Site *s;
	if(siteidx < nsites)
	{	s = &sites[siteidx];
		siteidx += 1;
		return(s);
	}
	else	return( (struct Site *)NULL);
}


/* 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.  (?) */

void voronoi(int triangulate, struct Site *(*nextsite)())
{
#pragma unused(triangulate)
	struct Site *newsite, *bot, *top, *temp, *p;
	struct Site *v;
	struct vPoint newintstar;
	int pm;
	struct Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector;
	struct Edge *e;

	PQinitialize();
	bottomsite = (*nextsite)();
	out_site(bottomsite);
	ELinitialize();
	
	newsite = (*nextsite)();
	while(1)
	{

		if(!PQempty()) newintstar = PQ_min();
	
		if (newsite != (struct Site *)NULL 
		   && (PQempty() 
			 || newsite -> coord.y < newintstar.y
			 || (newsite->coord.y == newintstar.y 
				 && newsite->coord.x < newintstar.x)))
		{/* new site is smallest */
			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);
			endpoint(lbnd->ELedge,lbnd->ELpm,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);
			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);
		};
}

int ntry, totalsearch;

void ELinitialize(void)
{
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(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);
}


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

/* Get entry from hash table, pruning any deleted nodes */
static struct Halfedge *ELgethash(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((struct Freenode *) he, &hfl);
	return ((struct Halfedge *) NULL);
}	

struct Halfedge *ELleftbnd(struct vPoint *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 correct 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.   */
void ELdelete(struct Halfedge *he)
{
	(he -> ELleft) -> ELright = he -> ELright;
	(he -> ELright) -> ELleft = he -> ELleft;
	he -> ELedge = (struct Edge *)DELETED;
}


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

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


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


void out_bisector(struct Edge *e)
{
if(!triangulate & !debug)
	printf("l %f %f %f", e->a, e->b, e->c);
if(debug)
	printf("line(%d) %g x + %g y = %g, bisecting %d %d\n", e->edgenbr,
	    e->a, e->b, e->c, e->reg[le]->sitenbr, e->reg[re]->sitenbr);
}


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

void out_vertex(struct Site *v)
{
if(!triangulate  & !debug)
	printf ("v %f %f\n", v->coord.x, v->coord.y);
if(debug)
	printf("vertex(%d) at %f %f\n", v->sitenbr, v->coord.x, v->coord.y);
}


void out_site(struct Site *s)
{
if(!triangulate & !debug)
	printf("s %f %f\n", s->coord.x, s->coord.y);
if(debug)
	printf("site (%d) at %f %f\n", s->sitenbr, s->coord.x, s->coord.y);
}


void out_triple(struct Site *s1,struct Site *s2,struct Site *s3)
/* Fill in triangle buffer */
{
if(triangulate)
{	/* Fill triangles buffer (as rows) */
	if (ntriangles < maxtriangles)
	{	*ptriangles = s1->sitenbr+1;
		*(ptriangles+maxtriangles) = s2->sitenbr+1; 
		*(ptriangles+2*maxtriangles) = s3->sitenbr+1;
		++ptriangles;
		++ntriangles;
		/* printf("%d %d %d\n", s1->sitenbr, s2->sitenbr, s3->sitenbr); */
	}
	else mexErrMsgTxt("Not enough triangles allocated.");
}
if(debug)
	printf("circle through left = %d right = %d bottom = %d\n", 
		s1->sitenbr, s2->sitenbr, s3->sitenbr);
}


void geominit(void)
{
	int sn;

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


struct Edge *bisect(struct Site *s1,struct Site *s2)
{
	double 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(struct Halfedge *el1, struct Halfedge *el2)
{
	struct	Edge *e1,*e2, *e;
	struct  Halfedge *el;
	double 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) {
		if (!warned) {
			mexWarnMsgTxt("Degenerate case: collinear points. Possible incorrect triangulation.");
			warned = 1;
		}
		if (d==0) d = mxGetEps();
		if (debug) {
			printf("Degenerate case (d = %g) between Site[%d] Site[%d] Site[%d] Site[%d]\n",d,
					e1->reg[0]->sitenbr, e1->reg[1]->sitenbr, 
					e2->reg[0]->sitenbr, e2->reg[1]->sitenbr);
		}
		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(struct Halfedge *el, struct vPoint *p)
{
	struct Edge *e;
	struct Site *topsite;
	int right_of_site, above, fast;
	double 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 + 2*dxp*dyp + dxs*dyp*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 endpoint(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((struct Freenode *) e, &efl);
}


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


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


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

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


static int PQbucket(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);
}

void PQinsert(struct Halfedge *he, struct Site *v, double 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;
}

void PQdelete(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 PQempty(void)
{
	return(PQcount==0);
}


struct vPoint PQ_min(void)
{
struct vPoint 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(void)
{
struct Halfedge *curr;

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


void PQinitialize(void)
{
int i;

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


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

char *getfree(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);
}



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

#if 0
int total_alloc;
char *myalloc(unsigned n)
{
	char *t;
	if ((t=malloc(n)) == (char *) 0) mexErrMsgTxt("Out of memory.");
	total_alloc += n;
	return(t);
}
#endif

void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
{
    int 	i,n;
    double 	*px, *py;
    int     	buflen = 7; /* long enough to fit 'sorted\0' */
    char    	buf[7];
	
    triangulate = 1; /* Do triangulation for now */
    debug = 0;
    warned = 0;
	
    if ((nrhs != 2)&&(nrhs != 3)) {
	mexErrMsgTxt("Wrong number of input parameters.");
    }

    if (nlhs > 1) {
	mexErrMsgTxt("Too many output arguments.");
    }
		
    n = mxGetM(prhs[0])*mxGetN(prhs[0]);
    if (n != mxGetM(prhs[1])*mxGetN(prhs[1])) {
	mexErrMsgTxt("X and Y must be the same size.");
    }

    if (nrhs == 2)
    {
        sorted = 0;
    }
    else /* nrhs == 3 */
    {
        if (mxIsChar(prhs[2]) && !mxGetString(prhs[2],buf,buflen) &&
            !strcmp(buf,"sorted")) {
	    sorted = 1;
	} else {
            mexErrMsgTxt("Third input parameter must be the string \'sorted\'");
        }
    }

    /* Initialize output triangle buffer */
    ntriangles = 0;

    /* If the inputs to the function are both empty */
    if (n == 0) {
	plhs[0] = mxCreateDoubleMatrix(0,3,mxREAL);
	return;
    }

    maxtriangles = 2*n+100; /* Rule of thumb: there are 2 triangles per 
			     *  vertex point (plus 100 just to be safe) */
    plhs[0] = mxCreateDoubleMatrix(maxtriangles,3,mxREAL);
    ptriangles = mxGetPr(plhs[0]);
	
    freeinit(&sfl, sizeof *sites);
    sites = (struct Site *) myalloc(n*sizeof *sites);
    for(nsites=0, px = mxGetPr(prhs[0]), py = mxGetPr(prhs[1]);
	nsites<n;
	++nsites,++px,++py)
      {
	  sites[nsites].coord.x = *px;
	  sites[nsites].coord.y = *py;
	  sites[nsites].sitenbr = nsites;
	  sites[nsites].refcnt = 0;
      };	

    if (!sorted)
      {
	  qsort(sites, nsites, sizeof *sites, (int (*)(const void*,const void *))scomp);
      }

    for(i=1; i<nsites; i++)
      {
	  if ((sites[i].coord.y == sites[i-1].coord.y) &&
	      (sites[i].coord.x == sites[i-1].coord.x))
            mexErrMsgTxt("x-y data must be unique");
      }
    
    xmin=sites[0].coord.x; 
    xmax=sites[0].coord.x;
    for(i=1; i<nsites; i+=1)
      {
	  if(sites[i].coord.x < xmin) xmin = sites[i].coord.x;
	  if(sites[i].coord.x > xmax) xmax = sites[i].coord.x;
      };	
    ymin = sites[0].coord.y;
    ymax = sites[nsites-1].coord.y;
    
    siteidx = 0;
    geominit();
	
    voronoi(triangulate, nextone);

    /* Shift triangle data in memory */
    if (ntriangles != maxtriangles)
      {	ptriangles = mxGetPr(plhs[0]);
      memmove(ptriangles+  ntriangles, ptriangles+  maxtriangles, ntriangles*sizeof(double));
      memmove(ptriangles+2*ntriangles, ptriangles+2*maxtriangles, ntriangles*sizeof(double));
      };
	
    mxSetM(plhs[0],ntriangles);
	
} /* mexFunction() */