delaunay.c

this file is contain important facts aboute different programming langueges.

/*
 * 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.
 */

/*
 *  Routines modified to provide a procedure interface,
 *  callable from Fortran, protect namespace by hiding all externals
 *  except pm_delaunay(), free most dynamic memory after call,
 *  node numbering starting at 1 (rather than 0), some unused code deleted,
 *  and minor mods for compilation under an ANSI C compiler.
 *
 *  Dynamic memory allocated to free lists are kept around after each call,
 *  don't know if this will cause a memory leak or not.
 *
 *  John Faricelli, Compaq Computer Corporation
 */

/* Original file: main.c */

#
#include <stdio.h>
#include <stdlib.h> 
#include <math.h> 

#include "delaunay.h"

static struct Site *nextone();

/* Take care of Fortran namespace issues */
#if defined(vms) || defined(__vms)
/*  VMS does not need name change (all externals are upper case'd) */
#define FORTRAN_CALL           
#elif defined(_WIN32)
/*  Win32 name decoration */
#define FORTRAN_CALL   __stdcall
#define pm_delaunay    PM_DELAUNAY
#else
/*  Unix name decoration */
#define FORTRAN_CALL           
#define pm_delaunay pm_delaunay_
#endif

int FORTRAN_CALL pm_delaunay(struct Point *cord, struct Tri *nop,
                             int *npts, int *maxtri)
{
int c;
struct Site *(*next)();
struct Freenode *tmp;

sorted = 0; triangulate = 1; debug = 0;

mycord =  cord;
mynop  =  nop;
mynpts = *npts;
myntri = 0;
mymaxtri = *maxtri;

freeinit(&sfl, sizeof *sites);

loadsites();
next = nextone;

siteidx = 0;
geominit();

voronoi(triangulate, next);

free(sites);
/* We should also free the nodes on the other freelists, but this
   seems complicated to do correctly... */

/* return number of triangles */
return(myntri);

}

/* sort sites on y, then x, coord */
static 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);
}

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

/* read all sites, sort, and compute xmin, xmax, ymin, ymax */
static loadsites()
{
int i;
int scomp(const void *, const void *);

nsites=0;
sites = (struct Site *) myalloc(4000*sizeof *sites);

for(i=0; i<mynpts; i++)
{
        sites[nsites].coord.x = mycord[i].x;
        sites[nsites].coord.y = mycord[i].y;
	sites[nsites].sitenbr = nsites;
	sites[nsites].refcnt = 0;
	nsites += 1;
	if (nsites % 4000 == 0) 
		sites = (struct Site *) realloc(sites,(nsites+4000)*sizeof*sites);
};

qsort(sites, nsites, sizeof *sites, scomp);
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;
}

/* Original file: edgelist.c */

static int ntry, totalsearch;

static 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;
}


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


static 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 */
static 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);
}	

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


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

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


static 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]);
}

static 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]);
}

/* Original file: geometry.c */

static geominit()
{
struct Edge e;
float sn;

	freeinit(&efl, sizeof e);
	nvertices = 0;
	nedges = 0;
	sn = (float) nsites+4;
#if defined(vms) || defined(__vms__) || defined(ultrix) || defined(_WIN32)
	sqrt_nsites = (int) sqrt( (double) sn);
#else
	sqrt_nsites = (int) sqrtf(sn);
#endif
	deltay = ymax - ymin;
	deltax = xmax - xmin;
}


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


static 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 */
static 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;