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