This is VoronoiViewPart.m in view mode; [Download] [Up]
#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
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.