PLplot  5.10.0
delaunay.c
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines