1 /******************************************************************************
2 *
3 * File: delaunay.c
4 *
5 * Created: 04/08/2000
6 *
7 * Author: Pavel Sakov
8 * CSIRO Marine Research
9 *
10 * Purpose: Delaunay triangulation - a wrapper to triangulate()
11 *
12 * Description: None
13 *
14 * Revisions: 10/06/2003 PS: delaunay_build(); delaunay_destroy();
15 * struct delaunay: from now on, only shallow copy of the
16 * input data is contained in struct delaunay. This saves
17 * memory and is consistent with libcsa.
18 * 30/10/2007 PS: added delaunay_addflag() and
19 * delaunay_resetflags(); modified delaunay_circles_find()
20 * to reset the flags to 0 on return. This is very important
21 * for large datasets, many thanks to John Gerschwitz,
22 * Petroleum Geo-Services, for identifying the problem.
23 *
24 *****************************************************************************/
25
26 #define ANSI_DECLARATORS /* for triangle.h */
27
28 #include <stdlib.h>
29 #include <stdio.h>
30 #include <assert.h>
31 #include <math.h>
32 #include <string.h>
33 #include <limits.h>
34 #include <float.h>
35 #include "istack.h"
36 #include "nan.h"
37 #include "delaunay.h"
38 #include "nn.h"
39 #include "nn_internal.h"
40
41 /*
42 * This parameter is used in search of tricircles containing a given point:
43 * if there are no more triangles than N_SEARCH_TURNON
44 * do linear search
45 * else
46 * do more complicated stuff
47 */
48 #define N_SEARCH_TURNON 20
49 #define N_FLAGS_TURNON 1000
50 #define N_FLAGS_INC 100
51
delaunay_create()52 static delaunay* delaunay_create()
53 {
54 delaunay* d = malloc(sizeof(delaunay));
55
56 d->npoints = 0;
57 d->points = NULL;
58 d->xmin = DBL_MAX;
59 d->xmax = -DBL_MAX;
60 d->ymin = DBL_MAX;
61 d->ymax = -DBL_MAX;
62 d->ntriangles = 0;
63 d->triangles = NULL;
64 d->circles = NULL;
65 d->neighbours = NULL;
66 d->n_point_triangles = NULL;
67 d->point_triangles = NULL;
68 d->nedges = 0;
69 d->edges = NULL;
70 d->flags = NULL;
71 d->first_id = -1;
72 d->t_in = NULL;
73 d->t_out = NULL;
74 d->nflags = 0;
75 d->nflagsallocated = 0;
76 d->flagids = NULL;
77
78 return d;
79 }
80
81 //---------------------------------------------------------
82 #ifndef USE_QHULL
83
84 #include "triangle.h"
85
tio_init(struct triangulateio * tio)86 static void tio_init(struct triangulateio* tio)
87 {
88 tio->pointlist = NULL;
89 tio->pointattributelist = NULL;
90 tio->pointmarkerlist = NULL;
91 tio->numberofpoints = 0;
92 tio->numberofpointattributes = 0;
93 tio->trianglelist = NULL;
94 tio->triangleattributelist = NULL;
95 tio->trianglearealist = NULL;
96 tio->neighborlist = NULL;
97 tio->numberoftriangles = 0;
98 tio->numberofcorners = 0;
99 tio->numberoftriangleattributes = 0;
100 tio->segmentlist = 0;
101 tio->segmentmarkerlist = NULL;
102 tio->numberofsegments = 0;
103 tio->holelist = NULL;
104 tio->numberofholes = 0;
105 tio->regionlist = NULL;
106 tio->numberofregions = 0;
107 tio->edgelist = NULL;
108 tio->edgemarkerlist = NULL;
109 tio->normlist = NULL;
110 tio->numberofedges = 0;
111 }
112
tio_destroy(struct triangulateio * tio)113 static void tio_destroy(struct triangulateio* tio)
114 {
115 if (tio->pointlist != NULL)
116 free(tio->pointlist);
117 if (tio->pointattributelist != NULL)
118 free(tio->pointattributelist);
119 if (tio->pointmarkerlist != NULL)
120 free(tio->pointmarkerlist);
121 if (tio->trianglelist != NULL)
122 free(tio->trianglelist);
123 if (tio->triangleattributelist != NULL)
124 free(tio->triangleattributelist);
125 if (tio->trianglearealist != NULL)
126 free(tio->trianglearealist);
127 if (tio->neighborlist != NULL)
128 free(tio->neighborlist);
129 if (tio->segmentlist != NULL)
130 free(tio->segmentlist);
131 if (tio->segmentmarkerlist != NULL)
132 free(tio->segmentmarkerlist);
133 if (tio->holelist != NULL)
134 free(tio->holelist);
135 if (tio->regionlist != NULL)
136 free(tio->regionlist);
137 if (tio->edgelist != NULL)
138 free(tio->edgelist);
139 if (tio->edgemarkerlist != NULL)
140 free(tio->edgemarkerlist);
141 if (tio->normlist != NULL)
142 free(tio->normlist);
143 }
144
tio2delaunay(struct triangulateio * tio_out,delaunay * d)145 static void tio2delaunay(struct triangulateio* tio_out, delaunay* d)
146 {
147 int i, j;
148
149 /*
150 * I assume that all input points appear in tio_out in the same order as
151 * they were written to tio_in. I have seen no exceptions so far, even
152 * if duplicate points were presented. Just in case, let us make a couple
153 * of checks.
154 */
155 assert(tio_out->numberofpoints == d->npoints);
156 assert(tio_out->pointlist[2 * d->npoints - 2] == d->points[d->npoints - 1].x && tio_out->pointlist[2 * d->npoints - 1] == d->points[d->npoints - 1].y);
157
158 for (i = 0, j = 0; i < d->npoints; ++i) {
159 point* p = &d->points[i];
160
161 if (p->x < d->xmin)
162 d->xmin = p->x;
163 if (p->x > d->xmax)
164 d->xmax = p->x;
165 if (p->y < d->ymin)
166 d->ymin = p->y;
167 if (p->y > d->ymax)
168 d->ymax = p->y;
169 }
170 if (nn_verbose) {
171 fprintf(stderr, "input:\n");
172 for (i = 0, j = 0; i < d->npoints; ++i) {
173 point* p = &d->points[i];
174
175 fprintf(stderr, " %d: %15.7g %15.7g %15.7g\n", i, p->x, p->y, p->z);
176 }
177 }
178
179 d->ntriangles = tio_out->numberoftriangles;
180 if (d->ntriangles > 0) {
181 d->triangles = malloc(d->ntriangles * sizeof(triangle));
182 d->neighbours = malloc(d->ntriangles * sizeof(triangle_neighbours));
183 d->circles = malloc(d->ntriangles * sizeof(circle));
184 d->n_point_triangles = calloc(d->npoints, sizeof(int));
185 d->point_triangles = malloc(d->npoints * sizeof(int*));
186 d->flags = calloc(d->ntriangles, sizeof(int));
187 }
188
189 if (nn_verbose)
190 fprintf(stderr, "triangles:\n");
191 for (i = 0; i < d->ntriangles; ++i) {
192 int offset = i * 3;
193 triangle* t = &d->triangles[i];
194 triangle_neighbours* n = &d->neighbours[i];
195 circle* c = &d->circles[i];
196 int status;
197
198 t->vids[0] = tio_out->trianglelist[offset];
199 t->vids[1] = tio_out->trianglelist[offset + 1];
200 t->vids[2] = tio_out->trianglelist[offset + 2];
201
202 n->tids[0] = tio_out->neighborlist[offset];
203 n->tids[1] = tio_out->neighborlist[offset + 1];
204 n->tids[2] = tio_out->neighborlist[offset + 2];
205
206 status = circle_build1(c, &d->points[t->vids[0]], &d->points[t->vids[1]], &d->points[t->vids[2]]);
207 assert(status);
208
209 if (nn_verbose)
210 fprintf(stderr, " %d: (%d,%d,%d)\n", i, t->vids[0], t->vids[1], t->vids[2]);
211 }
212
213 for (i = 0; i < d->ntriangles; ++i) {
214 triangle* t = &d->triangles[i];
215
216 for (j = 0; j < 3; ++j)
217 d->n_point_triangles[t->vids[j]]++;
218 }
219 if (d->ntriangles > 0) {
220 for (i = 0; i < d->npoints; ++i) {
221 if (d->n_point_triangles[i] > 0)
222 d->point_triangles[i] = malloc(d->n_point_triangles[i] * sizeof(int));
223 else
224 d->point_triangles[i] = NULL;
225 d->n_point_triangles[i] = 0;
226 }
227 }
228 for (i = 0; i < d->ntriangles; ++i) {
229 triangle* t = &d->triangles[i];
230
231 for (j = 0; j < 3; ++j) {
232 int vid = t->vids[j];
233
234 d->point_triangles[vid][d->n_point_triangles[vid]] = i;
235 d->n_point_triangles[vid]++;
236 }
237 }
238
239 if (tio_out->edgelist != NULL) {
240 d->nedges = tio_out->numberofedges;
241 d->edges = malloc(d->nedges * 2 * sizeof(int));
242 memcpy(d->edges, tio_out->edgelist, d->nedges * 2 * sizeof(int));
243 }
244 }
245
246 /* Builds Delaunay triangulation of the given array of points.
247 *
248 * @param np Number of points
249 * @param points Array of points [np] (input)
250 * @param ns Number of forced segments
251 * @param segments Array of (forced) segment endpoint indices [2*ns]
252 * @param nh Number of holes
253 * @param holes Array of hole (x,y) coordinates [2*nh]
254 * @return Delaunay triangulation structure with triangulation results
255 */
delaunay_build(int np,point points[],int ns,int segments[],int nh,double holes[])256 delaunay* delaunay_build(int np, point points[], int ns, int segments[], int nh, double holes[])
257 {
258 delaunay* d = delaunay_create();
259 struct triangulateio tio_in;
260 struct triangulateio tio_out;
261 char cmd[64] = "eznC";
262 int i, j;
263
264 assert(sizeof(REAL) == sizeof(double));
265
266 tio_init(&tio_in);
267
268 if (np == 0) {
269 free(d);
270 return NULL;
271 }
272
273 tio_in.pointlist = malloc(np * 2 * sizeof(double));
274 tio_in.numberofpoints = np;
275 for (i = 0, j = 0; i < np; ++i) {
276 tio_in.pointlist[j++] = points[i].x;
277 tio_in.pointlist[j++] = points[i].y;
278 }
279
280 if (ns > 0) {
281 tio_in.segmentlist = malloc(ns * 2 * sizeof(int));
282 tio_in.numberofsegments = ns;
283 memcpy(tio_in.segmentlist, segments, ns * 2 * sizeof(int));
284 }
285
286 if (nh > 0) {
287 tio_in.holelist = malloc(nh * 2 * sizeof(double));
288 tio_in.numberofholes = nh;
289 memcpy(tio_in.holelist, holes, nh * 2 * sizeof(double));
290 }
291
292 tio_init(&tio_out);
293
294 if (!nn_verbose)
295 strcat(cmd, "Q");
296 else if (nn_verbose > 1)
297 strcat(cmd, "VV");
298 if (ns != 0)
299 strcat(cmd, "p");
300
301 if (nn_verbose)
302 fflush(stderr);
303
304 /*
305 * climax
306 */
307 triangulate(cmd, &tio_in, &tio_out, NULL);
308
309 if (nn_verbose)
310 fflush(stderr);
311
312 d->npoints = np;
313 d->points = points;
314
315 tio2delaunay(&tio_out, d);
316
317 tio_destroy(&tio_in);
318 tio_destroy(&tio_out);
319
320 return d;
321 }
322
323 //---------------------------------------------------------
324 #else /* USE_QHULL */
325
326 #ifdef HAVE_LIBQHULL_R_QHULL_RA_H
327 #include <libqhull_r/qhull_ra.h>
328 #elif HAVE_LIBQHULL_QHULL_A_H
329 #include <libqhull/qhull_a.h>
330 #else
331 #ifdef HAVE_QHULL_QHULL_A_H
332 #include <qhull/qhull_a.h>
333 #else
334 #error Failed to include qhull_a.h. Is Qhull installed?
335 #endif
336 #endif
337
338 /* returns 1 if a,b,c are clockwise ordered */
cw(delaunay * d,triangle * t)339 static int cw(delaunay *d, triangle *t)
340 {
341 point* pa = &d->points[t->vids[0]];
342 point* pb = &d->points[t->vids[1]];
343 point* pc = &d->points[t->vids[2]];
344
345 return ((pb->x - pa->x)*(pc->y - pa->y) < (pc->x - pa->x)*(pb->y - pa->y));
346 }
347
delaunay_build(int np,point points[],int ns,int segments[],int nh,double holes[])348 delaunay* delaunay_build(int np, point points[], int ns, int segments[], int nh, double holes[])
349 {
350 delaunay* d = delaunay_create();
351
352 coordT *qpoints; /* array of coordinates for each point */
353 boolT ismalloc = False; /* True if qhull should free points */
354 char flags[64] = "qhull d Qbb Qt"; /* option flags for qhull */
355 facetT *facet,*neighbor,**neighborp; /* variables to walk through facets */
356 vertexT *vertex, **vertexp; /* variables to walk through vertex */
357
358 int curlong, totlong; /* memory remaining after qh_memfreeshort */
359 FILE *outfile = stdout;
360 FILE *errfile = stderr; /* error messages from qhull code */
361
362 int i, j;
363 int exitcode;
364 int dim, ntriangles;
365 int numfacets, numsimplicial, numridges, totneighbors, numcoplanars, numtricoplanars;
366
367 dim = 2;
368
369 assert(sizeof(realT) == sizeof(double)); /* Qhull was compiled with doubles? */
370
371 if (np == 0 || ns > 0 || nh > 0) {
372 fprintf(stderr, "segments=%d holes=%d\n, aborting Qhull implementation, use 'triangle' instead.\n", ns, nh);
373 free(d);
374 return NULL;
375 }
376
377 qpoints = (coordT *) malloc(np * (dim+1) * sizeof(coordT));
378
379 for (i=0; i<np; i++) {
380 qpoints[i*dim] = points[i].x;
381 qpoints[i*dim+1] = points[i].y;
382 }
383
384 if (!nn_verbose)
385 outfile = NULL;
386 if (nn_verbose)
387 strcat(flags, " s");
388 if (nn_verbose > 1)
389 strcat(flags, " Ts");
390
391 if (nn_verbose)
392 fflush(stderr);
393
394 /*
395 * climax
396 */
397
398 #ifdef HAVE_LIBQHULL_R_QHULL_RA_H
399 qhT qh_qh;
400 qhT *qh= &qh_qh;
401
402 QHULL_LIB_CHECK
403
404 qh_zero(qh, errfile);
405
406 exitcode = qh_new_qhull (qh, dim, np, qpoints, ismalloc,
407 flags, outfile, errfile);
408 #else
409 exitcode = qh_new_qhull (dim, np, qpoints, ismalloc,
410 flags, outfile, errfile);
411 #endif
412
413 if(!exitcode) {
414
415 if (nn_verbose)
416 fflush(stderr);
417
418 d->xmin = DBL_MAX;
419 d->xmax = -DBL_MAX;
420 d->ymin = DBL_MAX;
421 d->ymax = -DBL_MAX;
422
423 d->npoints = np;
424 d->points = malloc(np * sizeof(point));
425 for (i = 0; i < np; ++i) {
426 point* p = &d->points[i];
427
428 p->x = points[i].x;
429 p->y = points[i].y;
430 p->z = points[i].z;
431
432 if (p->x < d->xmin)
433 d->xmin = p->x;
434 if (p->x > d->xmax)
435 d->xmax = p->x;
436 if (p->y < d->ymin)
437 d->ymin = p->y;
438 if (p->y > d->ymax)
439 d->ymax = p->y;
440 }
441
442 if (nn_verbose) {
443 fprintf(stderr, "input:\n");
444 for (i = 0; i < np; ++i) {
445 point* p = &d->points[i];
446
447 fprintf(stderr, " %d: %15.7g %15.7g %15.7g\n",
448 i, p->x, p->y, p->z);
449 }
450 }
451
452 #ifdef HAVE_LIBQHULL_R_QHULL_RA_H
453 qh_findgood_all (qh, qh->facet_list);
454 qh_countfacets (qh, qh->facet_list, NULL, !qh_ALL, &numfacets,
455 &numsimplicial, &totneighbors, &numridges,
456 &numcoplanars, &numtricoplanars);
457 #else
458 qh_findgood_all (qh facet_list);
459 qh_countfacets (qh facet_list, NULL, !qh_ALL, &numfacets,
460 &numsimplicial, &totneighbors, &numridges,
461 &numcoplanars, &numtricoplanars);
462 #endif
463
464 ntriangles = 0;
465 FORALLfacets {
466 if (!facet->upperdelaunay && facet->simplicial)
467 ntriangles++;
468 }
469
470 d->ntriangles = ntriangles;
471 d->triangles = malloc(d->ntriangles * sizeof(triangle));
472 d->neighbours = malloc(d->ntriangles * sizeof(triangle_neighbours));
473 d->circles = malloc(d->ntriangles * sizeof(circle));
474
475 if (nn_verbose)
476 fprintf(stderr, "triangles:\tneighbors:\n");
477
478 i = 0;
479 FORALLfacets {
480 if (!facet->upperdelaunay && facet->simplicial) {
481 triangle* t = &d->triangles[i];
482 triangle_neighbours* n = &d->neighbours[i];
483 circle* c = &d->circles[i];
484
485 j = 0;
486 FOREACHvertex_(facet->vertices)
487
488 #ifdef HAVE_LIBQHULL_R_QHULL_RA_H
489 t->vids[j++] = qh_pointid(qh, vertex->point);
490 #else
491 t->vids[j++] = qh_pointid(vertex->point);
492 #endif
493
494 j = 0;
495 FOREACHneighbor_(facet)
496 n->tids[j++] = neighbor->visitid ? neighbor->visitid - 1 : - 1;
497
498 /* Put triangle vertices in counterclockwise order, as
499 * 'triangle' do.
500 * The same needs to be done with the neighbors.
501 *
502 * The following works, i.e., it seems that Qhull maintains a
503 * relationship between the vertices and the neighbors
504 * triangles, but that is not said anywhere, so if this stop
505 * working in a future Qhull release, you know what you have
506 * to do, reorder the neighbors.
507 */
508
509 if(cw(d, t)) {
510 int tmp = t->vids[1];
511 t->vids[1] = t->vids[2];
512 t->vids[2] = tmp;
513
514 tmp = n->tids[1];
515 n->tids[1] = n->tids[2];
516 n->tids[2] = tmp;
517 }
518
519 circle_build1(c, &d->points[t->vids[0]], &d->points[t->vids[1]],
520 &d->points[t->vids[2]]);
521
522 if (nn_verbose)
523 fprintf(stderr, " %d: (%d,%d,%d)\t(%d,%d,%d)\n",
524 i, t->vids[0], t->vids[1], t->vids[2], n->tids[0],
525 n->tids[1], n->tids[2]);
526
527 i++;
528 }
529 }
530
531 d->flags = calloc(d->ntriangles, sizeof(int));
532
533 d->n_point_triangles = calloc(d->npoints, sizeof(int));
534 for (i = 0; i < d->ntriangles; ++i) {
535 triangle* t = &d->triangles[i];
536
537 for (j = 0; j < 3; ++j)
538 d->n_point_triangles[t->vids[j]]++;
539 }
540 d->point_triangles = malloc(d->npoints * sizeof(int*));
541 for (i = 0; i < d->npoints; ++i) {
542 if (d->n_point_triangles[i] > 0)
543 d->point_triangles[i] = malloc(d->n_point_triangles[i] * sizeof(int));
544 else
545 d->point_triangles[i] = NULL;
546 d->n_point_triangles[i] = 0;
547 }
548 for (i = 0; i < d->ntriangles; ++i) {
549 triangle* t = &d->triangles[i];
550
551 for (j = 0; j < 3; ++j) {
552 int vid = t->vids[j];
553
554 d->point_triangles[vid][d->n_point_triangles[vid]] = i;
555 d->n_point_triangles[vid]++;
556 }
557 }
558
559 d->nedges = 0;
560 d->edges = NULL;
561
562 d->t_in = NULL;
563 d->t_out = NULL;
564 d->first_id = -1;
565
566 } else {
567 free(d);
568 d = NULL;
569 }
570
571 free(qpoints);
572
573 #ifdef HAVE_LIBQHULL_R_QHULL_RA_H
574 qh_freeqhull(qh, !qh_ALL); /* free long memory */
575 qh_memfreeshort (qh, &curlong, &totlong); /* free short memory and memory allocator */
576 #else
577 qh_freeqhull(!qh_ALL); /* free long memory */
578 qh_memfreeshort (&curlong, &totlong); /* free short memory and memory allocator */
579 #endif
580
581 if (curlong || totlong)
582 fprintf (errfile,
583 "qhull: did not free %d bytes of long memory (%d pieces)\n",
584 totlong, curlong);
585
586 return d;
587 }
588 #endif /* USE_QHULL */
589
590
591 /* Destroys Delaunay triangulation.
592 *
593 * @param d Structure to be destroyed
594 */
delaunay_destroy(delaunay * d)595 void delaunay_destroy(delaunay* d)
596 {
597 if (d == NULL)
598 return;
599
600 if (d->point_triangles != NULL) {
601 int i;
602
603 for (i = 0; i < d->npoints; ++i)
604 if (d->point_triangles[i] != NULL)
605 free(d->point_triangles[i]);
606 free(d->point_triangles);
607 }
608 if (d->nedges > 0)
609 free(d->edges);
610 #ifdef USE_QHULL
611 /* This is a shallow copy if we're not using qhull so we don't
612 * need to free it */
613 if (d->points != NULL)
614 free(d->points);
615 #endif
616 if (d->n_point_triangles != NULL)
617 free(d->n_point_triangles);
618 if (d->flags != NULL)
619 free(d->flags);
620 if (d->circles != NULL)
621 free(d->circles);
622 if (d->neighbours != NULL)
623 free(d->neighbours);
624 if (d->triangles != NULL)
625 free(d->triangles);
626 if (d->t_in != NULL)
627 istack_destroy(d->t_in);
628 if (d->t_out != NULL)
629 istack_destroy(d->t_out);
630 if (d->flagids != NULL)
631 free(d->flagids);
632 free(d);
633 }
634
635 /* Returns whether the point p is on the right side of the vector (p0, p1).
636 */
onrightside(point * p,point * p0,point * p1)637 static int onrightside(point* p, point* p0, point* p1)
638 {
639 return (p1->x - p->x) * (p0->y - p->y) > (p0->x - p->x) * (p1->y - p->y);
640 }
641
642 /* Finds triangle specified point belongs to (if any).
643 *
644 * @param d Delaunay triangulation
645 * @param p Point to be mapped
646 * @param seed Triangle index to start with
647 * @return Triangle id if successful, -1 otherwhile
648 */
delaunay_xytoi(delaunay * d,point * p,int id)649 int delaunay_xytoi(delaunay* d, point* p, int id)
650 {
651 triangle* t;
652 int i;
653
654 if (p->x < d->xmin || p->x > d->xmax || p->y < d->ymin || p->y > d->ymax)
655 return -1;
656
657 if (id < 0 || id > d->ntriangles)
658 id = 0;
659 t = &d->triangles[id];
660 do {
661 for (i = 0; i < 3; ++i) {
662 int i1 = (i + 1) % 3;
663
664 if (onrightside(p, &d->points[t->vids[i]], &d->points[t->vids[i1]])) {
665 id = d->neighbours[id].tids[(i + 2) % 3];
666 if (id < 0)
667 return id;
668 t = &d->triangles[id];
669 break;
670 }
671 }
672 } while (i < 3);
673
674 return id;
675 }
676
delaunay_addflag(delaunay * d,int i)677 static void delaunay_addflag(delaunay* d, int i)
678 {
679 if (d->nflags == d->nflagsallocated) {
680 d->nflagsallocated += N_FLAGS_INC;
681 d->flagids = realloc(d->flagids, d->nflagsallocated * sizeof(int));
682 }
683 d->flagids[d->nflags] = i;
684 d->nflags++;
685 }
686
delaunay_resetflags(delaunay * d)687 static void delaunay_resetflags(delaunay* d)
688 {
689 int i;
690
691 for (i = 0; i < d->nflags; ++i)
692 d->flags[d->flagids[i]] = 0;
693 d->nflags = 0;
694 }
695
696 /* Finds all tricircles specified point belongs to.
697 *
698 * @param d Delaunay triangulation
699 * @param p Point to be mapped
700 * @param n Pointer to the number of tricircles within `d' containing `p'
701 * (output)
702 * @param out Pointer to an array of indices of the corresponding triangles
703 * [n] (output)
704 *
705 * There is a standard search procedure involving search through triangle
706 * neighbours (not through vertex neighbours). It must be a bit faster due to
707 * the smaller number of triangle neighbours (3 per triangle) but may fail
708 * for a point outside convex hall.
709 *
710 * We may wish to modify this procedure in future: first check if the point
711 * is inside the convex hall, and depending on that use one of the two
712 * search algorithms. It not 100% clear though whether this will lead to a
713 * substantial speed gains because of the check on convex hall involved.
714 */
delaunay_circles_find(delaunay * d,point * p,int * n,int ** out)715 void delaunay_circles_find(delaunay* d, point* p, int* n, int** out)
716 {
717 /*
718 * This flag was introduced as a hack to handle some degenerate cases. It
719 * is set to 1 only if the triangle associated with the first circle is
720 * already known to contain the point. In this case the circle is assumed
721 * to contain the point without a check. In my practice this turned
722 * useful in some cases when point p coincided with one of the vertices
723 * of a thin triangle.
724 */
725 int contains = 0;
726 int i;
727
728 if (d->t_in == NULL) {
729 d->t_in = istack_create();
730 d->t_out = istack_create();
731 }
732
733 /*
734 * if there are only a few data points, do linear search
735 */
736 if (d->ntriangles <= N_SEARCH_TURNON) {
737 istack_reset(d->t_out);
738
739 for (i = 0; i < d->ntriangles; ++i) {
740 if (circle_contains(&d->circles[i], p)) {
741 istack_push(d->t_out, i);
742 }
743 }
744
745 *n = d->t_out->n;
746 *out = d->t_out->v;
747
748 return;
749 }
750 /*
751 * otherwise, do a more complicated stuff
752 */
753
754 /*
755 * It is important to have a reasonable seed here. If the last search
756 * was successful -- start with the last found tricircle, otherwhile (i)
757 * try to find a triangle containing p; if fails then (ii) check
758 * tricircles from the last search; if fails then (iii) make linear
759 * search through all tricircles
760 */
761 if (d->first_id < 0 || !circle_contains(&d->circles[d->first_id], p)) {
762 /*
763 * if any triangle contains p -- start with this triangle
764 */
765 d->first_id = delaunay_xytoi(d, p, d->first_id);
766 contains = (d->first_id >= 0);
767
768 /*
769 * if no triangle contains p, there still is a chance that it is
770 * inside some of circumcircles
771 */
772 if (d->first_id < 0) {
773 int nn = d->t_out->n;
774 int tid = -1;
775
776 /*
777 * first check results of the last search
778 */
779 for (i = 0; i < nn; ++i) {
780 tid = d->t_out->v[i];
781 if (circle_contains(&d->circles[tid], p))
782 break;
783 }
784 /*
785 * if unsuccessful, search through all circles
786 */
787 if (tid < 0 || i == nn) {
788 double nt = d->ntriangles;
789
790 for (tid = 0; tid < nt; ++tid) {
791 if (circle_contains(&d->circles[tid], p))
792 break;
793 }
794 if (tid == nt) {
795 istack_reset(d->t_out);
796 *n = 0;
797 *out = NULL;
798 return; /* failed */
799 }
800 }
801 d->first_id = tid;
802 }
803 }
804
805 istack_reset(d->t_in);
806 istack_reset(d->t_out);
807
808 istack_push(d->t_in, d->first_id);
809 d->flags[d->first_id] = 1;
810 delaunay_addflag(d, d->first_id);
811
812 /*
813 * main cycle
814 */
815 while (d->t_in->n > 0) {
816 int tid = istack_pop(d->t_in);
817 triangle* t = &d->triangles[tid];
818
819 if (contains || circle_contains(&d->circles[tid], p)) {
820 istack_push(d->t_out, tid);
821 for (i = 0; i < 3; ++i) {
822 int vid = t->vids[i];
823 int nt = d->n_point_triangles[vid];
824 int j;
825
826 for (j = 0; j < nt; ++j) {
827 int ntid = d->point_triangles[vid][j];
828
829 if (d->flags[ntid] == 0) {
830 istack_push(d->t_in, ntid);
831 d->flags[ntid] = 1;
832 delaunay_addflag(d, ntid);
833 }
834 }
835 }
836 }
837 contains = 0;
838 }
839
840 *n = d->t_out->n;
841 *out = d->t_out->v;
842 delaunay_resetflags(d);
843 }
844