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