1 /****************************************************************************/
2 /* */
3 /* This file is part of CONCORDE */
4 /* */
5 /* (c) Copyright 1995--1999 by David Applegate, Robert Bixby, */
6 /* Vasek Chvatal, and William Cook */
7 /* */
8 /* Permission is granted for academic research use. For other uses, */
9 /* contact the authors for licensing options. */
10 /* */
11 /* Use at your own risk. We make no guarantees about the */
12 /* correctness or usefulness of this code. */
13 /* */
14 /****************************************************************************/
15
16 /****************************************************************************/
17 /* */
18 /* EDGESET OF DELAUNAY TRIANGULATION (L2-NORM) */
19 /* */
20 /* TSP CODE */
21 /* */
22 /* *** ONLY MINOR MODIFICATIONS TO STEVE FORTUNE'S SWEEP2 CODE *** */
23 /* *** SWEEP2 IS AVAILABLE FROM NETLIB *** */
24 /* */
25 /* *** SWEEP2 COMES WITH THE FOLLOWING NOTICE *** */
26 /* */
27 /* The author of this software is Steven Fortune. Copyright (c) 1994 by */
28 /* AT&T Bell Laboratories. Permission to use, copy, modify, and distribute */
29 /* this software for any purpose without fee is hereby granted, provided */
30 /* that this entire notice is included in all copies of any software which */
31 /* is or includes a copy or modification of this software and in all copies */
32 /* of the supporting documentation for such software. THIS SOFTWARE IS */
33 /* BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED WARRANTY. IN */
34 /* PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY REPRESENTATION OR */
35 /* WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR */
36 /* ITS FITNESS FOR ANY PARTICULAR PURPOSE. */
37 /* */
38 /* *** END OF SWEEP2 NOTICE *** */
39 /* */
40 /* */
41 /* Modifications by: Applegate, Bixby, Chvatal, and Cook */
42 /* Date: March 12, 1996 */
43 /* */
44 /* */
45 /* EXPORTED FUNCTIONS: */
46 /* */
47 /* int CCedgegen_delaunay (int ncount, CCdatagroup *dat, int wantlist, */
48 /* int *ecount, int **elist) */
49 /* RETURNS the edgeset of the (Euclidean-norm) Delaunay triangulation */
50 /* of the point set for norms of size CC_D2_NORM_SIZE (2-coordinates). */
51 /* -ncount is the number of nodes) */
52 /* -dat contains the info to generate edge lengths */
53 /* -wantlist is 1 if you want the function to return the edges. */
54 /* -ecount returns the number of edges if wantlist is 1 */
55 /* -elist returns the edges in end1 end2 format if wantlist is 1 */
56 /* NOTES: */
57 /* To get the triangles (rather than the edges), modify the */
58 /* function out_triple as in the comments. */
59 /* */
60 /* The code is NOT threadsafe and does not clean up its memory. */
61 /* */
62 /****************************************************************************/
63
64 #include "machdefs.h"
65 #include "util.h"
66 #include "delaunay.h"
67 #include "macrorus.h"
68
69 #define DELETED -2
70 #define le 0
71 #define re 1
72
73 struct Freenode {
74 struct Freenode *nextfree;
75 };
76
77 struct Freelist {
78 struct Freenode *head;
79 int nodesize;
80 };
81
82 struct Point {
83 double x, y;
84 };
85
86 /* structure used both for sites and for vertices */
87
88 struct Site {
89 struct Point coord;
90 int sitenbr;
91 int refcnt;
92 };
93
94
95 struct Edge {
96 double a, b, c;
97 struct Site *ep[2];
98 struct Site *reg[2];
99 int edgenbr;
100 };
101
102 struct Halfedge {
103 struct Halfedge *ELleft, *ELright;
104 struct Edge *ELedge;
105 int ELrefcnt;
106 char ELpm;
107 struct Site *vertex;
108 double ystar;
109 struct Halfedge *PQnext;
110 };
111
112 typedef struct intptr {
113 int this;
114 struct intptr *next;
115 } intptr;
116
117 typedef struct delaunaydat {
118 int nvertices;
119 int PQmin;
120 int PQcount;
121 int PQhashsize;
122 int ELhashsize;
123 int siteidx;
124 struct Halfedge *PQhash;
125 double xmin, xmax, ymin, ymax, deltax, deltay;
126 int nedges;
127 struct Freelist efl;
128 struct Freelist hfl;
129 struct Halfedge *ELleftend, *ELrightend;
130 struct Halfedge **ELhash;
131 struct Site *sites;
132 int nsites;
133 int sqrt_nsites;
134 struct Freelist sfl;
135 struct Site *bottomsite;
136 intptr **table;
137 CCptrworld intptr_world;
138 int total_alloc;
139 } delaunaydat;
140
141
142 CC_PTRWORLD_LIST_ROUTINES (intptr, int, intptralloc, intptr_bulk_alloc,
143 intptrfree, intptr_listadd, intptr_listfree)
144 CC_PTRWORLD_LEAKS_ROUTINE (intptr, intptr_check_leaks, this, int)
145
146
147 static int
148 put_in_table (delaunaydat *dd, int i, int j, int *added);
149
150
151
152 static void
153 PQinsert (delaunaydat *dd, struct Halfedge *he, struct Site *v,
154 double offset),
155 PQdelete (delaunaydat *dd, struct Halfedge *he),
156 PQinitialize (delaunaydat *dd),
157 ELinitialize (delaunaydat *dd),
158 ELinsert (struct Halfedge *lb, struct Halfedge *new),
159 ELdelete (struct Halfedge *he),
160 freeinit (struct Freelist *fl, int size),
161 makefree (struct Freenode *curr, struct Freelist *fl),
162 geominit (delaunaydat *dd),
163 endpoint (delaunaydat *dd, struct Edge *e, int lr, struct Site *s),
164 deref (delaunaydat *dd, struct Site *v),
165 ref (struct Site *v),
166 makevertex (delaunaydat *dd, struct Site *v);
167
168 static char
169 *getfree (delaunaydat *dd, struct Freelist *fl),
170 *vor_myalloc (delaunaydat *dd, unsigned n);
171
172 static int
173 out_triple (delaunaydat *dd, struct Site *s1, struct Site *s2,
174 struct Site *s3, int *ntotal),
175 scomp (const void *s1_, const void *s2_),
176 set_up_sites (delaunaydat *dd, int ncount, double *x, double *y),
177 PQbucket (delaunaydat *dd, struct Halfedge *he),
178 PQempty (delaunaydat *dd),
179 right_of (struct Halfedge *el, struct Point *p);
180
181 static double
182 dist (struct Site *s, struct Site *t);
183
184 static struct Site
185 *nextsite (delaunaydat *dd),
186 *leftreg (delaunaydat *dd, struct Halfedge *he),
187 *rightreg (delaunaydat *dd, struct Halfedge *he),
188 *intersect (delaunaydat *dd, struct Halfedge *el1, struct Halfedge *el2);
189
190 static struct Edge
191 *bisect (delaunaydat *dd, struct Site *s1, struct Site *s2);
192
193 static struct Point
194 PQ_min (delaunaydat *dd);
195
196 static struct Halfedge
197 *PQextractmin (delaunaydat *dd),
198 *HEcreate (delaunaydat *dd, struct Edge *e, int pm),
199 *ELgethash (delaunaydat *dd, int b),
200 *ELleftbnd (delaunaydat *dd, struct Point *p),
201 *ELright (struct Halfedge *he),
202 *ELleft (struct Halfedge *he);
203
204
205
CCedgegen_delaunay(int ncount,CCdatagroup * dat,int wantlist,int * ecount,int ** elist)206 int CCedgegen_delaunay (int ncount, CCdatagroup *dat, int wantlist,
207 int *ecount, int **elist)
208 {
209 struct Site *newsite, *bot, *top, *temp, *p;
210 struct Site *v;
211 struct Point newintstar;
212 int i, pm;
213 struct Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector;
214 struct Edge *e;
215 int ntotal = 0;
216 int rval = 0;
217 intptr *ip, *ipnext;
218 int onlist, total;
219 int norm;
220 delaunaydat dd;
221
222 CCptrworld_init (&dd.intptr_world);
223 dd.total_alloc = 0;
224
225 dd.siteidx = 0;
226
227 if (wantlist) {
228 *ecount = 0;
229 *elist = (int *) NULL;
230 }
231
232 CCutil_dat_getnorm (dat, &norm);
233 if ((norm & CC_NORM_SIZE_BITS) != CC_D2_NORM_SIZE) {
234 printf ("Cannot compute Delaunay triangulation with norm %d\n",
235 norm);
236 fflush (stdout);
237 return 0;
238 }
239
240 dd.table = CC_SAFE_MALLOC (ncount, intptr *);
241 if (!dd.table) {
242 rval = 1;
243 goto CLEANUP;
244 }
245 for (i = 0; i < ncount; i++)
246 dd.table[i] = (intptr *) NULL;
247
248 freeinit (&dd.sfl, sizeof *dd.sites);
249 if (set_up_sites (&dd, ncount, dat->x, dat->y)) {
250 fprintf (stderr, "set_up_sites failed\n");
251 rval = 1;
252 goto CLEANUP;
253 }
254 geominit (&dd);
255
256 PQinitialize (&dd);
257 dd.bottomsite = nextsite (&dd);
258 ELinitialize (&dd);
259
260 newsite = nextsite (&dd);
261 while (1) {
262 if (!PQempty (&dd))
263 newintstar = PQ_min (&dd);
264
265 if (newsite != (struct Site *) NULL
266 && (PQempty (&dd)
267 || newsite->coord.y < newintstar.y
268 || (newsite->coord.y == newintstar.y
269 && newsite->coord.x < newintstar.x))) {
270 /* new site is smallest */
271 lbnd = ELleftbnd (&dd, &(newsite->coord));
272 rbnd = ELright (lbnd);
273 bot = rightreg (&dd, lbnd);
274 e = bisect (&dd, bot, newsite);
275 bisector = HEcreate (&dd, e, le);
276 ELinsert (lbnd, bisector);
277 if ((p = intersect (&dd, lbnd, bisector)) != (struct Site *) NULL) {
278 PQdelete (&dd, lbnd);
279 PQinsert (&dd, lbnd, p, dist (p, newsite));
280 }
281 lbnd = bisector;
282 bisector = HEcreate (&dd, e, re);
283 ELinsert (lbnd, bisector);
284 if ((p = intersect (&dd, bisector, rbnd)) != (struct Site *) NULL) {
285 PQinsert (&dd, bisector, p, dist (p, newsite));
286 }
287 newsite = nextsite (&dd);
288 } else if (!PQempty (&dd)) {
289 /* intersection is smallest */
290 lbnd = PQextractmin (&dd);
291 llbnd = ELleft (lbnd);
292 rbnd = ELright (lbnd);
293 rrbnd = ELright (rbnd);
294 bot = leftreg (&dd, lbnd);
295 top = rightreg (&dd, rbnd);
296 if (wantlist)
297 out_triple (&dd, bot, top, rightreg (&dd, lbnd), &ntotal);
298 v = lbnd->vertex;
299 makevertex (&dd, v);
300 endpoint (&dd, lbnd->ELedge, lbnd->ELpm, v);
301 endpoint (&dd, rbnd->ELedge, rbnd->ELpm, v);
302 ELdelete (lbnd);
303 PQdelete (&dd, rbnd);
304 ELdelete (rbnd);
305 pm = le;
306 if (bot->coord.y > top->coord.y) {
307 temp = bot;
308 bot = top;
309 top = temp;
310 pm = re;
311 }
312 e = bisect (&dd, bot, top);
313 bisector = HEcreate (&dd, e, pm);
314 ELinsert (llbnd, bisector);
315 endpoint (&dd, e, re - pm, v);
316 deref (&dd, v);
317 if ((p = intersect (&dd, llbnd, bisector)) != (struct Site *) NULL) {
318 PQdelete (&dd, llbnd);
319 PQinsert (&dd, llbnd, p, dist (p, bot));
320 }
321 if ((p = intersect (&dd, bisector, rrbnd)) != (struct Site *) NULL) {
322 PQinsert (&dd, bisector, p, dist (p, bot));
323 }
324 } else
325 break;
326 }
327
328 if (wantlist) {
329 int j = 0;
330 *elist = CC_SAFE_MALLOC (2 * ntotal, int);
331 if (!(*elist)) {
332 rval = 1;
333 goto CLEANUP;
334 }
335 *ecount = ntotal;
336 for (i = 0; i < ncount; i++) {
337 for (ip = dd.table[i]; ip; ip = ipnext) {
338 ipnext = ip->next;
339 (*elist)[j++] = i;
340 (*elist)[j++] = ip->this;
341 intptrfree (&dd.intptr_world, ip);
342 }
343 dd.table[i] = (intptr *) NULL;
344 }
345 } else {
346 for (i = 0; i < ncount; i++) {
347 intptr_listfree (&dd.intptr_world, dd.table[i]);
348 dd.table[i] = (intptr *) NULL;
349 }
350 }
351
352 if (intptr_check_leaks (&dd.intptr_world, &total, &onlist)) {
353 fprintf (stderr, "WARNING: %d outstanding intptrs in delaunay\n",
354 total - onlist);
355 }
356
357 CLEANUP:
358
359 CCptrworld_delete (&dd.intptr_world);
360 if (dd.table)
361 CC_FREE (dd.table, intptr *);
362 if (dd.sites)
363 CC_FREE (dd.sites, struct Site);
364
365 return rval;
366 }
367
put_in_table(delaunaydat * dd,int i,int j,int * added)368 static int put_in_table (delaunaydat *dd, int i, int j, int *added)
369 {
370 intptr *ip;
371
372 if (j < i) {
373 int temp;
374 CC_SWAP(i, j, temp);
375 }
376
377 for (ip = dd->table[i]; ip; ip = ip->next) {
378 if (ip->this == j) {
379 *added = 0;
380 return 0;
381 }
382 }
383 if (intptr_listadd (&dd->table[i], j, &dd->intptr_world)) {
384 *added = 0;
385 return 1;
386 }
387 *added = 1;
388 return 0;
389 }
390
391 /* sort sites on y, then x, coord */
392
scomp(const void * s1_,const void * s2_)393 static int scomp (const void *s1_, const void *s2_)
394 {
395 const struct Point *s1,*s2;
396 s1=(const struct Point*) s1_;
397 s2=(const struct Point*) s2_;
398
399 if (s1->y < s2->y)
400 return (-1);
401 if (s1->y > s2->y)
402 return (1);
403 if (s1->x < s2->x)
404 return (-1);
405 if (s1->x > s2->x)
406 return (1);
407 return (0);
408 }
409
set_up_sites(delaunaydat * dd,int ncount,double * x,double * y)410 static int set_up_sites (delaunaydat *dd, int ncount, double *x, double *y)
411 {
412 int i;
413
414 dd->nsites = ncount;
415 dd->sites = CC_SAFE_MALLOC (ncount, struct Site);
416 if (!dd->sites) {
417 fprintf (stderr, "out of memory in set_up_sites\n");
418 fflush (stdout);
419 return 1;
420 }
421
422 for (i = 0; i < ncount; i++) {
423 dd->sites[i].coord.x = x[i];
424 dd->sites[i].coord.y = y[i];
425 dd->sites[i].sitenbr = i;
426 dd->sites[i].refcnt = 0;
427 }
428 qsort ((void *)dd->sites, (size_t)dd->nsites, sizeof *dd->sites, scomp);
429 dd->xmin = dd->sites[0].coord.x;
430 dd->xmax = dd->sites[0].coord.x;
431 for (i = 1; i < dd->nsites; i += 1) {
432 if (dd->sites[i].coord.x < dd->xmin)
433 dd->xmin = dd->sites[i].coord.x;
434 if (dd->sites[i].coord.x > dd->xmax)
435 dd->xmax = dd->sites[i].coord.x;
436 }
437 dd->ymin = dd->sites[0].coord.y;
438 dd->ymax = dd->sites[dd->nsites - 1].coord.y;
439
440 return 0;
441 }
442
out_triple(delaunaydat * dd,struct Site * s1,struct Site * s2,struct Site * s3,int * ntotal)443 static int out_triple (delaunaydat *dd, struct Site *s1, struct Site *s2,
444 struct Site *s3, int *ntotal)
445 {
446 int added;
447 /*
448 to just print the triangles, use:
449
450 printf ("%d %d %d\n", s1->sitenbr, s2->sitenbr, s3->sitenbr);
451 */
452
453 if (put_in_table (dd, s1->sitenbr, s2->sitenbr, &added)) {
454 fprintf (stderr, "put_in_table failed\n");
455 return 1;
456 }
457 *ntotal += added;
458 if (put_in_table (dd, s2->sitenbr, s3->sitenbr, &added)) {
459 fprintf (stderr, "put_in_table failed\n");
460 return 1;
461 }
462 *ntotal += added;
463 if (put_in_table (dd, s3->sitenbr, s1->sitenbr, &added)) {
464 fprintf (stderr, "put_in_table failed\n");
465 return 1;
466 }
467 *ntotal += added;
468 return 0;
469 }
470
nextsite(delaunaydat * dd)471 static struct Site *nextsite (delaunaydat *dd)
472 {
473 struct Site *s;
474
475 if (dd->siteidx < dd->nsites) {
476 s = &dd->sites[dd->siteidx];
477 dd->siteidx += 1;
478 return (s);
479 } else
480 return ((struct Site *) NULL);
481 }
482
483
484 /******************* edgelist.c ***********************/
485
ELinitialize(delaunaydat * dd)486 static void ELinitialize (delaunaydat *dd)
487 {
488 int i;
489
490 freeinit (&dd->hfl, sizeof **dd->ELhash);
491 dd->ELhashsize = 2 * dd->sqrt_nsites;
492 dd->ELhash = (struct Halfedge **) vor_myalloc (dd, (unsigned) (sizeof (*dd->ELhash) * dd->ELhashsize));
493 for (i = 0; i < dd->ELhashsize; i += 1)
494 dd->ELhash[i] = (struct Halfedge *) NULL;
495 dd->ELleftend = HEcreate (dd, (struct Edge *) NULL, 0);
496 dd->ELrightend = HEcreate (dd, (struct Edge *) NULL, 0);
497 dd->ELleftend->ELleft = (struct Halfedge *) NULL;
498 dd->ELleftend->ELright = dd->ELrightend;
499 dd->ELrightend->ELleft = dd->ELleftend;
500 dd->ELrightend->ELright = (struct Halfedge *) NULL;
501 dd->ELhash[0] = dd->ELleftend;
502 dd->ELhash[dd->ELhashsize - 1] = dd->ELrightend;
503 }
504
HEcreate(delaunaydat * dd,struct Edge * e,int pm)505 static struct Halfedge *HEcreate (delaunaydat *dd, struct Edge *e, int pm)
506 {
507 struct Halfedge *answer;
508
509 answer = (struct Halfedge *) getfree (dd, &dd->hfl);
510 answer->ELedge = e;
511 answer->ELpm = pm;
512 answer->PQnext = (struct Halfedge *) NULL;
513 answer->vertex = (struct Site *) NULL;
514 answer->ELrefcnt = 0;
515 return (answer);
516 }
517
ELinsert(struct Halfedge * lb,struct Halfedge * new)518 static void ELinsert (struct Halfedge *lb, struct Halfedge *new)
519 {
520 new->ELleft = lb;
521 new->ELright = lb->ELright;
522 (lb->ELright)->ELleft = new;
523 lb->ELright = new;
524 }
525
526 /* Get entry from hash table, pruning any deleted nodes */
527
ELgethash(delaunaydat * dd,int b)528 static struct Halfedge *ELgethash (delaunaydat *dd, int b)
529 {
530 struct Halfedge *he;
531
532 if (b < 0 || b >= dd->ELhashsize)
533 return ((struct Halfedge *) NULL);
534 he = dd->ELhash[b];
535 if (he == (struct Halfedge *) NULL ||
536 he->ELedge != (struct Edge *) (size_t) DELETED)
537 return (he);
538
539 /* Hash table points to deleted half edge. Patch as necessary. */
540 dd->ELhash[b] = (struct Halfedge *) NULL;
541 if ((he->ELrefcnt -= 1) == 0)
542 makefree ((struct Freenode *) he, &dd->hfl);
543 return ((struct Halfedge *) NULL);
544 }
545
ELleftbnd(delaunaydat * dd,struct Point * p)546 static struct Halfedge *ELleftbnd (delaunaydat *dd, struct Point *p)
547 {
548 int i, bucket;
549 struct Halfedge *he;
550
551 /* Use hash table to get close to desired halfedge */
552 bucket = (p->x - dd->xmin) / dd->deltax * dd->ELhashsize;
553 if (bucket < 0)
554 bucket = 0;
555 if (bucket >= dd->ELhashsize)
556 bucket = dd->ELhashsize - 1;
557 he = ELgethash (dd, bucket);
558 if (he == (struct Halfedge *) NULL) {
559 for (i = 1; 1; i += 1) {
560 if ((he = ELgethash (dd, bucket - i)) != (struct Halfedge *) NULL)
561 break;
562 if ((he = ELgethash (dd, bucket + i)) != (struct Halfedge *) NULL)
563 break;
564 }
565 }
566 /* Now search linear list of halfedges for the corect one */
567 if (he == dd->ELleftend || (he != dd->ELrightend && right_of (he, p))) {
568 do {
569 he = he->ELright;
570 } while (he != dd->ELrightend && right_of (he, p));
571 he = he->ELleft;
572 } else
573 do {
574 he = he->ELleft;
575 } while (he != dd->ELleftend && !right_of (he, p));
576
577 /* Update hash table and reference counts */
578 if (bucket > 0 && bucket < dd->ELhashsize - 1) {
579 if (dd->ELhash[bucket] != (struct Halfedge *) NULL)
580 dd->ELhash[bucket]->ELrefcnt -= 1;
581 dd->ELhash[bucket] = he;
582 dd->ELhash[bucket]->ELrefcnt += 1;
583 };
584 return (he);
585 }
586
587
588 /* This delete routine can't reclaim node, since pointers from hash table may
589 * be present. */
590
ELdelete(struct Halfedge * he)591 static void ELdelete (struct Halfedge *he)
592 {
593 (he->ELleft)->ELright = he->ELright;
594 (he->ELright)->ELleft = he->ELleft;
595 he->ELedge = (struct Edge *) (size_t) DELETED;
596 }
597
598
ELright(struct Halfedge * he)599 static struct Halfedge *ELright (struct Halfedge *he)
600 {
601 return (he->ELright);
602 }
603
ELleft(struct Halfedge * he)604 static struct Halfedge *ELleft (struct Halfedge *he)
605 {
606 return (he->ELleft);
607 }
608
609
leftreg(delaunaydat * dd,struct Halfedge * he)610 static struct Site *leftreg (delaunaydat *dd, struct Halfedge *he)
611 {
612 if (he->ELedge == (struct Edge *) NULL)
613 return (dd->bottomsite);
614 return (he->ELpm == le ?
615 he->ELedge->reg[le] : he->ELedge->reg[re]);
616 }
617
rightreg(delaunaydat * dd,struct Halfedge * he)618 static struct Site *rightreg (delaunaydat *dd, struct Halfedge *he)
619 {
620 if (he->ELedge == (struct Edge *) NULL)
621 return (dd->bottomsite);
622 return (he->ELpm == le ?
623 he->ELedge->reg[re] : he->ELedge->reg[le]);
624 }
625
626
627 /***************** geometry.c ********************/
628
geominit(delaunaydat * dd)629 static void geominit (delaunaydat *dd)
630 {
631 struct Edge e;
632 double sn;
633
634 freeinit (&dd->efl, sizeof e);
635 dd->nvertices = 0;
636 dd->nedges = 0;
637 sn = dd->nsites + 4;
638 dd->sqrt_nsites = sqrt (sn);
639 dd->deltay = dd->ymax - dd->ymin;
640 dd->deltax = dd->xmax - dd->xmin;
641 }
642
bisect(delaunaydat * dd,struct Site * s1,struct Site * s2)643 static struct Edge *bisect (delaunaydat *dd, struct Site *s1, struct Site *s2)
644 {
645 double dx, dy, adx, ady;
646 struct Edge *newedge;
647
648 newedge = (struct Edge *) getfree (dd, &dd->efl);
649
650 newedge->reg[0] = s1;
651 newedge->reg[1] = s2;
652 ref (s1);
653 ref (s2);
654 newedge->ep[0] = (struct Site *) NULL;
655 newedge->ep[1] = (struct Site *) NULL;
656
657 dx = s2->coord.x - s1->coord.x;
658 dy = s2->coord.y - s1->coord.y;
659 adx = dx > 0 ? dx : -dx;
660 ady = dy > 0 ? dy : -dy;
661 newedge->c = s1->coord.x * dx + s1->coord.y * dy +
662 (dx * dx + dy * dy) * 0.5;
663 if (adx > ady) {
664 newedge->a = 1.0;
665 newedge->b = dy / dx;
666 newedge->c /= dx;
667 } else {
668 newedge->b = 1.0;
669 newedge->a = dx / dy;
670 newedge->c /= dy;
671 };
672
673 newedge->edgenbr = dd->nedges;
674 dd->nedges += 1;
675 return (newedge);
676 }
677
intersect(delaunaydat * dd,struct Halfedge * el1,struct Halfedge * el2)678 static struct Site *intersect (delaunaydat *dd, struct Halfedge *el1,
679 struct Halfedge *el2)
680 {
681 struct Edge *e1, *e2, *e;
682 struct Halfedge *el;
683 double d, xint, yint;
684 int right_of_site;
685 struct Site *v;
686
687 e1 = el1->ELedge;
688 e2 = el2->ELedge;
689 if (e1 == (struct Edge *) NULL || e2 == (struct Edge *) NULL)
690 return ((struct Site *) NULL);
691 if (e1->reg[1] == e2->reg[1])
692 return ((struct Site *) NULL);
693
694 d = e1->a * e2->b - e1->b * e2->a;
695 if (-1.0e-10 < d && d < 1.0e-10)
696 return ((struct Site *) NULL);
697
698 xint = (e1->c * e2->b - e2->c * e1->b) / d;
699 yint = (e2->c * e1->a - e1->c * e2->a) / d;
700
701 if ((e1->reg[1]->coord.y < e2->reg[1]->coord.y) ||
702 (e1->reg[1]->coord.y == e2->reg[1]->coord.y &&
703 e1->reg[1]->coord.x < e2->reg[1]->coord.x)) {
704 el = el1;
705 e = e1;
706 } else {
707 el = el2;
708 e = e2;
709 };
710 right_of_site = xint >= e->reg[1]->coord.x;
711 if ((right_of_site && el->ELpm == le) ||
712 (!right_of_site && el->ELpm == re))
713 return ((struct Site *) NULL);
714
715 v = (struct Site *) getfree (dd, &dd->sfl);
716 v->refcnt = 0;
717 v->coord.x = xint;
718 v->coord.y = yint;
719 return (v);
720 }
721
722 /* returns 1 if p is to right of halfedge e */
723
right_of(struct Halfedge * el,struct Point * p)724 static int right_of (struct Halfedge *el, struct Point *p)
725 {
726 struct Edge *e;
727 struct Site *topsite;
728 int right_of_site, above, fast;
729 double dxp, dyp, dxs, t1, t2, t3, yl;
730
731 e = el->ELedge;
732 topsite = e->reg[1];
733 right_of_site = p->x > topsite->coord.x;
734 if (right_of_site && el->ELpm == le)
735 return (1);
736 if (!right_of_site && el->ELpm == re)
737 return (0);
738
739 if (e->a == 1.0) {
740 dyp = p->y - topsite->coord.y;
741 dxp = p->x - topsite->coord.x;
742 fast = 0;
743 if ((!right_of_site & (e->b < 0.0)) |
744 (right_of_site & (e->b >= 0.0))) {
745 above = dyp >= e->b * dxp;
746 fast = above;
747 } else {
748 above = p->x + p->y * e->b > e->c;
749 if (e->b < 0.0)
750 above = !above;
751 if (!above)
752 fast = 1;
753 };
754 if (!fast) {
755 dxs = topsite->coord.x - (e->reg[0])->coord.x;
756 above = e->b * (dxp * dxp - dyp * dyp) <
757 dxs * dyp * (1.0 + 2.0 * dxp / dxs + e->b * e->b);
758 if (e->b < 0.0)
759 above = !above;
760 };
761 } else { /* e->b==1.0 */
762 yl = e->c - e->a * p->x;
763 t1 = p->y - yl;
764 t2 = p->x - topsite->coord.x;
765 t3 = yl - topsite->coord.y;
766 above = t1 * t1 > t2 * t2 + t3 * t3;
767 };
768 return (el->ELpm == le ? above : !above);
769 }
770
771
endpoint(delaunaydat * dd,struct Edge * e,int lr,struct Site * s)772 static void endpoint (delaunaydat *dd, struct Edge *e, int lr, struct Site *s)
773 {
774 e->ep[lr] = s;
775 ref (s);
776 if (e->ep[re - lr] == (struct Site *) NULL)
777 return;
778 deref (dd, e->reg[le]);
779 deref (dd, e->reg[re]);
780 makefree ((struct Freenode *) e, &dd->efl);
781 }
782
783
dist(struct Site * s,struct Site * t)784 static double dist (struct Site *s, struct Site *t)
785 {
786 double dx, dy;
787
788 dx = s->coord.x - t->coord.x;
789 dy = s->coord.y - t->coord.y;
790 return (sqrt (dx * dx + dy * dy));
791 }
792
793
makevertex(delaunaydat * dd,struct Site * v)794 static void makevertex (delaunaydat *dd, struct Site *v)
795 {
796 v->sitenbr = dd->nvertices;
797 dd->nvertices += 1;
798 }
799
800
deref(delaunaydat * dd,struct Site * v)801 static void deref (delaunaydat *dd, struct Site *v)
802 {
803 v->refcnt -= 1;
804 if (v->refcnt == 0)
805 makefree ((struct Freenode *) v, &dd->sfl);
806 }
807
ref(struct Site * v)808 static void ref (struct Site *v)
809 {
810 v->refcnt += 1;
811 }
812
813
814 /*********************** heap.c *******************/
815
PQinsert(delaunaydat * dd,struct Halfedge * he,struct Site * v,double offset)816 static void PQinsert (delaunaydat *dd, struct Halfedge *he, struct Site *v,
817 double offset)
818 {
819 struct Halfedge *last, *next;
820
821 he->vertex = v;
822 ref (v);
823 he->ystar = v->coord.y + offset;
824 last = &dd->PQhash[PQbucket (dd, he)];
825 while ((next = last->PQnext) != (struct Halfedge *) NULL &&
826 (he->ystar > next->ystar ||
827 (he->ystar == next->ystar && v->coord.x > next->vertex->coord.x))) {
828 last = next;
829 }
830 he->PQnext = last->PQnext;
831 last->PQnext = he;
832 dd->PQcount += 1;
833 }
834
PQdelete(delaunaydat * dd,struct Halfedge * he)835 static void PQdelete (delaunaydat *dd, struct Halfedge *he)
836 {
837 struct Halfedge *last;
838
839 if (he->vertex != (struct Site *) NULL) {
840 last = &dd->PQhash[PQbucket (dd, he)];
841 while (last->PQnext != he)
842 last = last->PQnext;
843 last->PQnext = he->PQnext;
844 dd->PQcount -= 1;
845 deref (dd, he->vertex);
846 he->vertex = (struct Site *) NULL;
847 }
848 }
849
PQbucket(delaunaydat * dd,struct Halfedge * he)850 static int PQbucket (delaunaydat *dd, struct Halfedge *he)
851 {
852 int bucket;
853
854 bucket = (he->ystar - dd->ymin) / dd->deltay * dd->PQhashsize;
855 if (bucket < 0)
856 bucket = 0;
857 if (bucket >= dd->PQhashsize)
858 bucket = dd->PQhashsize - 1;
859 if (bucket < dd->PQmin)
860 dd->PQmin = bucket;
861 return (bucket);
862 }
863
864
865
PQempty(delaunaydat * dd)866 static int PQempty (delaunaydat *dd)
867 {
868 return (dd->PQcount == 0);
869 }
870
871
PQ_min(delaunaydat * dd)872 static struct Point PQ_min (delaunaydat *dd)
873 {
874 struct Point answer;
875
876 while (dd->PQhash[dd->PQmin].PQnext == (struct Halfedge *) NULL) {
877 dd->PQmin += 1;
878 }
879 answer.x = dd->PQhash[dd->PQmin].PQnext->vertex->coord.x;
880 answer.y = dd->PQhash[dd->PQmin].PQnext->ystar;
881 return (answer);
882 }
883
PQextractmin(delaunaydat * dd)884 static struct Halfedge *PQextractmin (delaunaydat *dd)
885 {
886 struct Halfedge *curr;
887
888 curr = dd->PQhash[dd->PQmin].PQnext;
889 dd->PQhash[dd->PQmin].PQnext = curr->PQnext;
890 dd->PQcount -= 1;
891 return (curr);
892 }
893
894
PQinitialize(delaunaydat * dd)895 static void PQinitialize (delaunaydat *dd)
896 {
897 int i;
898
899 dd->PQcount = 0;
900 dd->PQmin = 0;
901 dd->PQhashsize = 4 * dd->sqrt_nsites;
902 dd->PQhash = (struct Halfedge *) vor_myalloc (dd, (unsigned) (dd->PQhashsize * sizeof (*dd->PQhash)));
903 for (i = 0; i < dd->PQhashsize; i += 1)
904 dd->PQhash[i].PQnext = (struct Halfedge *) NULL;
905 }
906
907
908 /***************** memory.c *****************/
909
freeinit(struct Freelist * fl,int size)910 static void freeinit (struct Freelist *fl, int size)
911 {
912 fl->head = (struct Freenode *) NULL;
913 fl->nodesize = size;
914 }
915
getfree(delaunaydat * dd,struct Freelist * fl)916 static char *getfree (delaunaydat *dd, struct Freelist *fl)
917 {
918 int i;
919 struct Freenode *t;
920
921 if (fl->head == (struct Freenode *) NULL) {
922 t = (struct Freenode *) vor_myalloc (dd, (unsigned) (dd->sqrt_nsites * fl->nodesize));
923 for (i = 0; i < dd->sqrt_nsites; i += 1)
924 makefree ((struct Freenode *) ((char *) t + i * fl->nodesize), fl);
925 };
926 t = fl->head;
927 fl->head = (fl->head)->nextfree;
928 return ((char *) t);
929 }
930
makefree(struct Freenode * curr,struct Freelist * fl)931 static void makefree (struct Freenode *curr, struct Freelist *fl)
932 {
933 curr->nextfree = fl->head;
934 fl->head = curr;
935 }
936
vor_myalloc(delaunaydat * dd,unsigned n)937 static char *vor_myalloc (delaunaydat *dd, unsigned n)
938 {
939 char *t;
940
941 if ((t = malloc (n)) == (char *) 0) {
942 fprintf (stderr, "Out of memory processing %d (%d bytes in use)\n",
943 dd->siteidx, dd->total_alloc);
944 exit (1);
945 };
946 dd->total_alloc += n;
947 return (t);
948 }
949