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