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