PLplot
5.10.0
|
00001 //-------------------------------------------------------------------------- 00002 // 00003 // File: delaunay.c 00004 // 00005 // Created: 04/08/2000 00006 // 00007 // Author: Pavel Sakov 00008 // CSIRO Marine Research 00009 // 00010 // Purpose: Delaunay triangulation - a wrapper to triangulate() 00011 // 00012 // Description: None 00013 // 00014 // Revisions: 10/06/2003 PS: delaunay_build(); delaunay_destroy(); 00015 // struct delaunay: from now on, only shallow copy of the 00016 // input data is contained in struct delaunay. This saves 00017 // memory and is consistent with libcsa. 00018 // 00019 // Modified: Joao Cardoso, 4/2/2003 00020 // Adapted for use with Qhull instead of "triangle". 00021 // Andrew Ross 20/10/2008 00022 // Fix bug in delaunay_circles_find() when checking 00023 // whether a circle has been found. 00024 // 00025 //-------------------------------------------------------------------------- 00026 00027 #define USE_QHULL 00028 00029 #include <stdlib.h> 00030 #include <stdio.h> 00031 #include <assert.h> 00032 #include <math.h> 00033 #include <string.h> 00034 #include <limits.h> 00035 #include <float.h> 00036 #ifdef USE_QHULL 00037 #include <qhull/qhull_a.h> 00038 #else 00039 #include "triangle.h" 00040 #endif 00041 #include "istack.h" 00042 #include "nan.h" 00043 #include "delaunay.h" 00044 00045 int circle_build( circle* c, point* p0, point* p1, point* p2 ); 00046 int circle_contains( circle* c, point* p ); 00047 int delaunay_xytoi( delaunay* d, point* p, int id ); 00048 void delaunay_circles_find( delaunay* d, point* p, int* n, int** out ); 00049 00050 #ifdef USE_QHULL 00051 static int cw( delaunay *d, triangle *t ); 00052 #endif 00053 00054 #ifndef USE_QHULL 00055 static void tio_init( struct triangulateio* tio ) 00056 { 00057 tio->pointlist = NULL; 00058 tio->pointattributelist = NULL; 00059 tio->pointmarkerlist = NULL; 00060 tio->numberofpoints = 0; 00061 tio->numberofpointattributes = 0; 00062 tio->trianglelist = NULL; 00063 tio->triangleattributelist = NULL; 00064 tio->trianglearealist = NULL; 00065 tio->neighborlist = NULL; 00066 tio->numberoftriangles = 0; 00067 tio->numberofcorners = 0; 00068 tio->numberoftriangleattributes = 0; 00069 tio->segmentlist = 0; 00070 tio->segmentmarkerlist = NULL; 00071 tio->numberofsegments = 0; 00072 tio->holelist = NULL; 00073 tio->numberofholes = 0; 00074 tio->regionlist = NULL; 00075 tio->numberofregions = 0; 00076 tio->edgelist = NULL; 00077 tio->edgemarkerlist = NULL; 00078 tio->normlist = NULL; 00079 tio->numberofedges = 0; 00080 } 00081 00082 static void tio_destroy( struct triangulateio* tio ) 00083 { 00084 if ( tio->pointlist != NULL ) 00085 free( tio->pointlist ); 00086 if ( tio->pointattributelist != NULL ) 00087 free( tio->pointattributelist ); 00088 if ( tio->pointmarkerlist != NULL ) 00089 free( tio->pointmarkerlist ); 00090 if ( tio->trianglelist != NULL ) 00091 free( tio->trianglelist ); 00092 if ( tio->triangleattributelist != NULL ) 00093 free( tio->triangleattributelist ); 00094 if ( tio->trianglearealist != NULL ) 00095 free( tio->trianglearealist ); 00096 if ( tio->neighborlist != NULL ) 00097 free( tio->neighborlist ); 00098 if ( tio->segmentlist != NULL ) 00099 free( tio->segmentlist ); 00100 if ( tio->segmentmarkerlist != NULL ) 00101 free( tio->segmentmarkerlist ); 00102 if ( tio->holelist != NULL ) 00103 free( tio->holelist ); 00104 if ( tio->regionlist != NULL ) 00105 free( tio->regionlist ); 00106 if ( tio->edgelist != NULL ) 00107 free( tio->edgelist ); 00108 if ( tio->edgemarkerlist != NULL ) 00109 free( tio->edgemarkerlist ); 00110 if ( tio->normlist != NULL ) 00111 free( tio->normlist ); 00112 } 00113 00114 static delaunay* delaunay_create() 00115 { 00116 delaunay* d = malloc( sizeof ( delaunay ) ); 00117 00118 d->npoints = 0; 00119 d->points = NULL; 00120 d->xmin = DBL_MAX; 00121 d->xmax = -DBL_MAX; 00122 d->ymin = DBL_MAX; 00123 d->ymax = -DBL_MAX; 00124 d->ntriangles = 0; 00125 d->triangles = NULL; 00126 d->circles = NULL; 00127 d->neighbours = NULL; 00128 d->n_point_triangles = NULL; 00129 d->point_triangles = NULL; 00130 d->nedges = 0; 00131 d->edges = NULL; 00132 d->flags = NULL; 00133 d->first_id = -1; 00134 d->t_in = NULL; 00135 d->t_out = NULL; 00136 00137 return d; 00138 } 00139 00140 static void tio2delaunay( struct triangulateio* tio_out, delaunay* d ) 00141 { 00142 int i, j; 00143 00144 // 00145 // I assume that all input points appear in tio_out in the same order as 00146 // they were written to tio_in. I have seen no exceptions so far, even 00147 // if duplicate points were presented. Just in case, let us make a couple 00148 // of checks. 00149 // 00150 assert( tio_out->numberofpoints == d->npoints ); 00151 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 ); 00152 00153 for ( i = 0, j = 0; i < d->npoints; ++i ) 00154 { 00155 point* p = &d->points[i]; 00156 00157 if ( p->x < d->xmin ) 00158 d->xmin = p->x; 00159 if ( p->x > d->xmax ) 00160 d->xmax = p->x; 00161 if ( p->y < d->ymin ) 00162 d->ymin = p->y; 00163 if ( p->y > d->ymax ) 00164 d->ymax = p->y; 00165 } 00166 if ( nn_verbose ) 00167 { 00168 fprintf( stderr, "input:\n" ); 00169 for ( i = 0, j = 0; i < d->npoints; ++i ) 00170 { 00171 point* p = &d->points[i]; 00172 00173 fprintf( stderr, " %d: %15.7g %15.7g %15.7g\n", i, p->x, p->y, p->z ); 00174 } 00175 } 00176 00177 d->ntriangles = tio_out->numberoftriangles; 00178 if ( d->ntriangles > 0 ) 00179 { 00180 d->triangles = malloc( d->ntriangles * sizeof ( triangle ) ); 00181 d->neighbours = malloc( d->ntriangles * sizeof ( triangle_neighbours ) ); 00182 d->circles = malloc( d->ntriangles * sizeof ( circle ) ); 00183 d->n_point_triangles = calloc( d->npoints, sizeof ( int ) ); 00184 d->point_triangles = malloc( d->npoints * sizeof ( int* ) ); 00185 d->flags = calloc( d->ntriangles, sizeof ( int ) ); 00186 } 00187 00188 if ( nn_verbose ) 00189 fprintf( stderr, "triangles:\n" ); 00190 for ( i = 0; i < d->ntriangles; ++i ) 00191 { 00192 int offset = i * 3; 00193 triangle * t = &d->triangles[i]; 00194 triangle_neighbours* n = &d->neighbours[i]; 00195 circle * c = &d->circles[i]; 00196 00197 t->vids[0] = tio_out->trianglelist[offset]; 00198 t->vids[1] = tio_out->trianglelist[offset + 1]; 00199 t->vids[2] = tio_out->trianglelist[offset + 2]; 00200 00201 n->tids[0] = tio_out->neighborlist[offset]; 00202 n->tids[1] = tio_out->neighborlist[offset + 1]; 00203 n->tids[2] = tio_out->neighborlist[offset + 2]; 00204 00205 circle_build( c, &d->points[t->vids[0]], &d->points[t->vids[1]], &d->points[t->vids[2]] ); 00206 00207 if ( nn_verbose ) 00208 fprintf( stderr, " %d: (%d,%d,%d)\n", i, t->vids[0], t->vids[1], t->vids[2] ); 00209 } 00210 00211 for ( i = 0; i < d->ntriangles; ++i ) 00212 { 00213 triangle* t = &d->triangles[i]; 00214 00215 for ( j = 0; j < 3; ++j ) 00216 d->n_point_triangles[t->vids[j]]++; 00217 } 00218 if ( d->ntriangles > 0 ) 00219 { 00220 for ( i = 0; i < d->npoints; ++i ) 00221 { 00222 if ( d->n_point_triangles[i] > 0 ) 00223 d->point_triangles[i] = malloc( d->n_point_triangles[i] * sizeof ( int ) ); 00224 else 00225 d->point_triangles[i] = NULL; 00226 d->n_point_triangles[i] = 0; 00227 } 00228 } 00229 for ( i = 0; i < d->ntriangles; ++i ) 00230 { 00231 triangle* t = &d->triangles[i]; 00232 00233 for ( j = 0; j < 3; ++j ) 00234 { 00235 int vid = t->vids[j]; 00236 00237 d->point_triangles[vid][d->n_point_triangles[vid]] = i; 00238 d->n_point_triangles[vid]++; 00239 } 00240 } 00241 00242 if ( tio_out->edgelist != NULL ) 00243 { 00244 d->nedges = tio_out->numberofedges; 00245 d->edges = malloc( d->nedges * 2 * sizeof ( int ) ); 00246 memcpy( d->edges, tio_out->edgelist, d->nedges * 2 * sizeof ( int ) ); 00247 } 00248 } 00249 #endif 00250 00251 // Builds Delaunay triangulation of the given array of points. 00252 // 00253 // @param np Number of points 00254 // @param points Array of points [np] (input) 00255 // @param ns Number of forced segments 00256 // @param segments Array of (forced) segment endpoint indices [2*ns] 00257 // @param nh Number of holes 00258 // @param holes Array of hole (x,y) coordinates [2*nh] 00259 // @return Delaunay triangulation structure with triangulation results 00260 // 00261 delaunay* delaunay_build( int np, point points[], int ns, int segments[], int nh, double holes[] ) 00262 #ifndef USE_QHULL 00263 { 00264 delaunay * d = delaunay_create(); 00265 struct triangulateio tio_in; 00266 struct triangulateio tio_out; 00267 char cmd[64] = "eznC"; 00268 int i, j; 00269 00270 assert( sizeof ( REAL ) == sizeof ( double ) ); 00271 00272 tio_init( &tio_in ); 00273 00274 if ( np == 0 ) 00275 { 00276 free( d ); 00277 return NULL; 00278 } 00279 00280 tio_in.pointlist = malloc( np * 2 * sizeof ( double ) ); 00281 tio_in.numberofpoints = np; 00282 for ( i = 0, j = 0; i < np; ++i ) 00283 { 00284 tio_in.pointlist[j++] = points[i].x; 00285 tio_in.pointlist[j++] = points[i].y; 00286 } 00287 00288 if ( ns > 0 ) 00289 { 00290 tio_in.segmentlist = malloc( ns * 2 * sizeof ( int ) ); 00291 tio_in.numberofsegments = ns; 00292 memcpy( tio_in.segmentlist, segments, ns * 2 * sizeof ( int ) ); 00293 } 00294 00295 if ( nh > 0 ) 00296 { 00297 tio_in.holelist = malloc( nh * 2 * sizeof ( double ) ); 00298 tio_in.numberofholes = nh; 00299 memcpy( tio_in.holelist, holes, nh * 2 * sizeof ( double ) ); 00300 } 00301 00302 tio_init( &tio_out ); 00303 00304 if ( !nn_verbose ) 00305 strcat( cmd, "Q" ); 00306 else if ( nn_verbose > 1 ) 00307 strcat( cmd, "VV" ); 00308 if ( ns != 0 ) 00309 strcat( cmd, "p" ); 00310 00311 if ( nn_verbose ) 00312 fflush( stderr ); 00313 00314 // 00315 // climax 00316 // 00317 triangulate( cmd, &tio_in, &tio_out, NULL ); 00318 00319 if ( nn_verbose ) 00320 fflush( stderr ); 00321 00322 d->npoints = np; 00323 d->points = points; 00324 00325 tio2delaunay( &tio_out, d ); 00326 00327 tio_destroy( &tio_in ); 00328 tio_destroy( &tio_out ); 00329 00330 return d; 00331 } 00332 #else // USE_QHULL 00333 { 00334 delaunay* d = malloc( sizeof ( delaunay ) ); 00335 00336 coordT *qpoints; // array of coordinates for each point 00337 boolT ismalloc = False; // True if qhull should free points 00338 char flags[64] = "qhull d Qbb Qt"; // option flags for qhull 00339 facetT *facet, *neighbor, **neighborp; // variables to walk through facets 00340 vertexT *vertex, **vertexp; // variables to walk through vertex 00341 00342 int curlong, totlong; // memory remaining after qh_memfreeshort 00343 FILE *outfile = stdout; 00344 FILE *errfile = stderr; // error messages from qhull code 00345 00346 int i, j; 00347 int exitcode; 00348 int dim, ntriangles; 00349 int numfacets, numsimplicial, numridges, totneighbors, numcoplanars, numtricoplanars; 00350 00351 (void) segments; // Cast to void to suppress compiler warnings about unused parameters 00352 (void) holes; 00353 00354 dim = 2; 00355 00356 assert( sizeof ( realT ) == sizeof ( double ) ); // Qhull was compiled with doubles? 00357 00358 if ( np == 0 || ns > 0 || nh > 0 ) 00359 { 00360 fprintf( stderr, "segments=%d holes=%d\n, aborting Qhull implementation, use 'triangle' instead.\n", ns, nh ); 00361 free( d ); 00362 return NULL; 00363 } 00364 00365 qpoints = (coordT *) malloc( (size_t) ( np * ( dim + 1 ) ) * sizeof ( coordT ) ); 00366 00367 for ( i = 0; i < np; i++ ) 00368 { 00369 qpoints[i * dim] = points[i].x; 00370 qpoints[i * dim + 1] = points[i].y; 00371 } 00372 00373 if ( !nn_verbose ) 00374 outfile = NULL; 00375 if ( nn_verbose ) 00376 strcat( flags, " s" ); 00377 if ( nn_verbose > 1 ) 00378 strcat( flags, " Ts" ); 00379 00380 if ( nn_verbose ) 00381 fflush( stderr ); 00382 00383 // 00384 // climax 00385 // 00386 00387 exitcode = qh_new_qhull( dim, np, qpoints, ismalloc, 00388 flags, outfile, errfile ); 00389 00390 if ( !exitcode ) 00391 { 00392 if ( nn_verbose ) 00393 fflush( stderr ); 00394 00395 d->xmin = DBL_MAX; 00396 d->xmax = -DBL_MAX; 00397 d->ymin = DBL_MAX; 00398 d->ymax = -DBL_MAX; 00399 00400 d->npoints = np; 00401 d->points = malloc( (size_t) np * sizeof ( point ) ); 00402 for ( i = 0; i < np; ++i ) 00403 { 00404 point* p = &d->points[i]; 00405 00406 p->x = points[i].x; 00407 p->y = points[i].y; 00408 p->z = points[i].z; 00409 00410 if ( p->x < d->xmin ) 00411 d->xmin = p->x; 00412 if ( p->x > d->xmax ) 00413 d->xmax = p->x; 00414 if ( p->y < d->ymin ) 00415 d->ymin = p->y; 00416 if ( p->y > d->ymax ) 00417 d->ymax = p->y; 00418 } 00419 00420 if ( nn_verbose ) 00421 { 00422 fprintf( stderr, "input:\n" ); 00423 for ( i = 0; i < np; ++i ) 00424 { 00425 point* p = &d->points[i]; 00426 00427 fprintf( stderr, " %d: %15.7g %15.7g %15.7g\n", 00428 i, p->x, p->y, p->z ); 00429 } 00430 } 00431 00432 qh_findgood_all( qh facet_list ); 00433 qh_countfacets( qh facet_list, NULL, !qh_ALL, &numfacets, 00434 &numsimplicial, &totneighbors, &numridges, 00435 &numcoplanars, &numtricoplanars ); 00436 00437 ntriangles = 0; 00438 FORALLfacets { 00439 if ( !facet->upperdelaunay && facet->simplicial ) 00440 ntriangles++; 00441 } 00442 00443 d->ntriangles = ntriangles; 00444 d->triangles = malloc( (size_t) d->ntriangles * sizeof ( triangle ) ); 00445 d->neighbours = malloc( (size_t) d->ntriangles * sizeof ( triangle_neighbours ) ); 00446 d->circles = malloc( (size_t) d->ntriangles * sizeof ( circle ) ); 00447 00448 if ( nn_verbose ) 00449 fprintf( stderr, "triangles:\tneighbors:\n" ); 00450 00451 i = 0; 00452 FORALLfacets { 00453 if ( !facet->upperdelaunay && facet->simplicial ) 00454 { 00455 triangle * t = &d->triangles[i]; 00456 triangle_neighbours* n = &d->neighbours[i]; 00457 circle * c = &d->circles[i]; 00458 00459 j = 0; 00460 FOREACHvertex_( facet->vertices ) 00461 t->vids[j++] = qh_pointid( vertex->point ); 00462 00463 j = 0; 00464 FOREACHneighbor_( facet ) 00465 n->tids[j++] = ( neighbor->visitid > 0 ) ? (int) neighbor->visitid - 1 : -1; 00466 00467 // Put triangle vertices in counterclockwise order, as 00468 // 'triangle' do. 00469 // The same needs to be done with the neighbors. 00470 // 00471 // The following works, i.e., it seems that Qhull maintains a 00472 // relationship between the vertices and the neighbors 00473 // triangles, but that is not said anywhere, so if this stop 00474 // working in a future Qhull release, you know what you have 00475 // to do, reorder the neighbors. 00476 // 00477 00478 if ( cw( d, t ) ) 00479 { 00480 int tmp = t->vids[1]; 00481 t->vids[1] = t->vids[2]; 00482 t->vids[2] = tmp; 00483 00484 tmp = n->tids[1]; 00485 n->tids[1] = n->tids[2]; 00486 n->tids[2] = tmp; 00487 } 00488 00489 circle_build( c, &d->points[t->vids[0]], &d->points[t->vids[1]], 00490 &d->points[t->vids[2]] ); 00491 00492 if ( nn_verbose ) 00493 fprintf( stderr, " %d: (%d,%d,%d)\t(%d,%d,%d)\n", 00494 i, t->vids[0], t->vids[1], t->vids[2], n->tids[0], 00495 n->tids[1], n->tids[2] ); 00496 00497 i++; 00498 } 00499 } 00500 00501 d->flags = calloc( (size_t) ( d->ntriangles ), sizeof ( int ) ); 00502 00503 d->n_point_triangles = calloc( (size_t) ( d->npoints ), sizeof ( int ) ); 00504 for ( i = 0; i < d->ntriangles; ++i ) 00505 { 00506 triangle* t = &d->triangles[i]; 00507 00508 for ( j = 0; j < 3; ++j ) 00509 d->n_point_triangles[t->vids[j]]++; 00510 } 00511 d->point_triangles = malloc( (size_t) ( d->npoints ) * sizeof ( int* ) ); 00512 for ( i = 0; i < d->npoints; ++i ) 00513 { 00514 if ( d->n_point_triangles[i] > 0 ) 00515 d->point_triangles[i] = malloc( (size_t) ( d->n_point_triangles[i] ) * sizeof ( int ) ); 00516 else 00517 d->point_triangles[i] = NULL; 00518 d->n_point_triangles[i] = 0; 00519 } 00520 for ( i = 0; i < d->ntriangles; ++i ) 00521 { 00522 triangle* t = &d->triangles[i]; 00523 00524 for ( j = 0; j < 3; ++j ) 00525 { 00526 int vid = t->vids[j]; 00527 00528 d->point_triangles[vid][d->n_point_triangles[vid]] = i; 00529 d->n_point_triangles[vid]++; 00530 } 00531 } 00532 00533 d->nedges = 0; 00534 d->edges = NULL; 00535 00536 d->t_in = NULL; 00537 d->t_out = NULL; 00538 d->first_id = -1; 00539 } 00540 else 00541 { 00542 free( d ); 00543 d = NULL; 00544 } 00545 00546 free( qpoints ); 00547 qh_freeqhull( !qh_ALL ); // free long memory 00548 qh_memfreeshort( &curlong, &totlong ); // free short memory and memory allocator 00549 if ( curlong || totlong ) 00550 fprintf( errfile, 00551 "qhull: did not free %d bytes of long memory (%d pieces)\n", 00552 totlong, curlong ); 00553 00554 return d; 00555 } 00556 00557 // returns 1 if a,b,c are clockwise ordered 00558 static int cw( delaunay *d, triangle *t ) 00559 { 00560 point* pa = &d->points[t->vids[0]]; 00561 point* pb = &d->points[t->vids[1]]; 00562 point* pc = &d->points[t->vids[2]]; 00563 00564 return ( ( pb->x - pa->x ) * ( pc->y - pa->y ) < 00565 ( pc->x - pa->x ) * ( pb->y - pa->y ) ); 00566 } 00567 00568 #endif 00569 00570 // Releases memory engaged in the Delaunay triangulation structure. 00571 // 00572 // @param d Structure to be destroyed 00573 // 00574 void delaunay_destroy( delaunay* d ) 00575 { 00576 if ( d == NULL ) 00577 return; 00578 00579 if ( d->point_triangles != NULL ) 00580 { 00581 int i; 00582 00583 for ( i = 0; i < d->npoints; ++i ) 00584 if ( d->point_triangles[i] != NULL ) 00585 free( d->point_triangles[i] ); 00586 free( d->point_triangles ); 00587 } 00588 if ( d->nedges > 0 ) 00589 free( d->edges ); 00590 #ifdef USE_QHULL 00591 // This is a shallow copy if we're not using qhull so we don't 00592 // need to free it 00593 if ( d->points != NULL ) 00594 free( d->points ); 00595 #endif 00596 if ( d->n_point_triangles != NULL ) 00597 free( d->n_point_triangles ); 00598 if ( d->flags != NULL ) 00599 free( d->flags ); 00600 if ( d->circles != NULL ) 00601 free( d->circles ); 00602 if ( d->neighbours != NULL ) 00603 free( d->neighbours ); 00604 if ( d->triangles != NULL ) 00605 free( d->triangles ); 00606 if ( d->t_in != NULL ) 00607 istack_destroy( d->t_in ); 00608 if ( d->t_out != NULL ) 00609 istack_destroy( d->t_out ); 00610 free( d ); 00611 } 00612 00613 // Returns whether the point p is on the right side of the vector (p0, p1). 00614 // 00615 static int on_right_side( point* p, point* p0, point* p1 ) 00616 { 00617 return ( p1->x - p->x ) * ( p0->y - p->y ) > ( p0->x - p->x ) * ( p1->y - p->y ); 00618 } 00619 00620 // Finds triangle specified point belongs to (if any). 00621 // 00622 // @param d Delaunay triangulation 00623 // @param p Point to be mapped 00624 // @param seed Triangle index to start with 00625 // @return Triangle id if successful, -1 otherwhile 00626 // 00627 int delaunay_xytoi( delaunay* d, point* p, int id ) 00628 { 00629 triangle* t; 00630 int i; 00631 00632 if ( p->x < d->xmin || p->x > d->xmax || p->y < d->ymin || p->y > d->ymax ) 00633 return -1; 00634 00635 if ( id < 0 || id > d->ntriangles ) 00636 id = 0; 00637 t = &d->triangles[id]; 00638 do 00639 { 00640 for ( i = 0; i < 3; ++i ) 00641 { 00642 int i1 = ( i + 1 ) % 3; 00643 00644 if ( on_right_side( p, &d->points[t->vids[i]], &d->points[t->vids[i1]] ) ) 00645 { 00646 id = d->neighbours[id].tids[( i + 2 ) % 3]; 00647 if ( id < 0 ) 00648 return id; 00649 t = &d->triangles[id]; 00650 break; 00651 } 00652 } 00653 } while ( i < 3 ); 00654 00655 return id; 00656 } 00657 00658 // Finds all tricircles specified point belongs to. 00659 // 00660 // @param d Delaunay triangulation 00661 // @param p Point to be mapped 00662 // @param n Pointer to the number of tricircles within `d' containing `p' 00663 // (output) 00664 // @param out Pointer to an array of indices of the corresponding triangles 00665 // [n] (output) 00666 // 00667 // There is a standard search procedure involving search through triangle 00668 // neighbours (not through vertex neighbours). It must be a bit faster due to 00669 // the smaller number of triangle neighbours (3 per triangle) but can fail 00670 // for a point outside convex hall. 00671 // 00672 // We may wish to modify this procedure in future: first check if the point 00673 // is inside the convex hall, and depending on that use one of the two 00674 // search algorithms. It not 100% clear though whether this will lead to a 00675 // substantial speed gains because of the check on convex hall involved. 00676 // 00677 void delaunay_circles_find( delaunay* d, point* p, int* n, int** out ) 00678 { 00679 int i; 00680 00681 if ( d->t_in == NULL ) 00682 { 00683 d->t_in = istack_create(); 00684 d->t_out = istack_create(); 00685 } 00686 00687 // 00688 // It is important to have a reasonable seed here. If the last search 00689 // was successful -- start with the last found tricircle, otherwhile (i) 00690 // try to find a triangle containing (x,y); if fails then (ii) check 00691 // tricircles from the last search; if fails then (iii) make linear 00692 // search through all tricircles 00693 // 00694 if ( d->first_id < 0 || !circle_contains( &d->circles[d->first_id], p ) ) 00695 { 00696 // 00697 // if any triangle contains (x,y) -- start with this triangle 00698 // 00699 d->first_id = delaunay_xytoi( d, p, d->first_id ); 00700 00701 // 00702 // if no triangle contains (x,y), there still is a chance that it is 00703 // inside some of circumcircles 00704 // 00705 if ( d->first_id < 0 ) 00706 { 00707 int nn = d->t_out->n; 00708 int tid = -1; 00709 00710 // 00711 // first check results of the last search 00712 // 00713 for ( i = 0; i < nn; ++i ) 00714 { 00715 tid = d->t_out->v[i]; 00716 if ( circle_contains( &d->circles[tid], p ) ) 00717 break; 00718 } 00719 // 00720 // if unsuccessful, search through all circles 00721 // 00722 if ( tid < 0 || i == nn ) 00723 { 00724 double nt = d->ntriangles; 00725 00726 for ( tid = 0; tid < nt; ++tid ) 00727 { 00728 if ( circle_contains( &d->circles[tid], p ) ) 00729 break; 00730 } 00731 if ( tid == nt ) 00732 { 00733 istack_reset( d->t_out ); 00734 *n = 0; 00735 *out = NULL; 00736 return; // failed 00737 } 00738 } 00739 d->first_id = tid; 00740 } 00741 } 00742 00743 istack_reset( d->t_in ); 00744 istack_reset( d->t_out ); 00745 00746 istack_push( d->t_in, d->first_id ); 00747 d->flags[d->first_id] = 1; 00748 00749 // 00750 // main cycle 00751 // 00752 while ( d->t_in->n > 0 ) 00753 { 00754 int tid = istack_pop( d->t_in ); 00755 triangle* t = &d->triangles[tid]; 00756 00757 if ( circle_contains( &d->circles[tid], p ) ) 00758 { 00759 istack_push( d->t_out, tid ); 00760 for ( i = 0; i < 3; ++i ) 00761 { 00762 int vid = t->vids[i]; 00763 int nt = d->n_point_triangles[vid]; 00764 int j; 00765 00766 for ( j = 0; j < nt; ++j ) 00767 { 00768 int ntid = d->point_triangles[vid][j]; 00769 00770 if ( d->flags[ntid] == 0 ) 00771 { 00772 istack_push( d->t_in, ntid ); 00773 d->flags[ntid] = 1; 00774 } 00775 } 00776 } 00777 } 00778 } 00779 00780 *n = d->t_out->n; 00781 *out = d->t_out->v; 00782 }