37 #include <qhull/qhull_a.h>
55 static void tio_init(
struct triangulateio* tio )
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;
70 tio->segmentmarkerlist = NULL;
71 tio->numberofsegments = 0;
73 tio->numberofholes = 0;
74 tio->regionlist = NULL;
75 tio->numberofregions = 0;
77 tio->edgemarkerlist = NULL;
79 tio->numberofedges = 0;
82 static void tio_destroy(
struct triangulateio* tio )
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 );
140 static void tio2delaunay(
struct triangulateio* tio_out,
delaunay* d )
150 assert( tio_out->numberofpoints == d->
npoints );
153 for ( i = 0, j = 0; i < d->
npoints; ++i )
157 if ( p->
x < d->
xmin )
159 if ( p->
x > d->
xmax )
161 if ( p->
y < d->
ymin )
163 if ( p->
y > d->
ymax )
168 fprintf( stderr,
"input:\n" );
169 for ( i = 0, j = 0; i < d->
npoints; ++i )
173 fprintf( stderr,
" %d: %15.7g %15.7g %15.7g\n", i, p->
x, p->
y, p->
z );
189 fprintf( stderr,
"triangles:\n" );
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];
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];
208 fprintf( stderr,
" %d: (%d,%d,%d)\n", i, t->
vids[0], t->
vids[1], t->
vids[2] );
215 for ( j = 0; j < 3; ++j )
220 for ( i = 0; i < d->
npoints; ++i )
233 for ( j = 0; j < 3; ++j )
235 int vid = t->
vids[j];
242 if ( tio_out->edgelist != NULL )
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 ) );
265 struct triangulateio tio_in;
266 struct triangulateio tio_out;
267 char cmd[64] =
"eznC";
270 assert(
sizeof ( REAL ) ==
sizeof (
double ) );
280 tio_in.pointlist = malloc( np * 2 *
sizeof (
double ) );
281 tio_in.numberofpoints = np;
282 for ( i = 0, j = 0; i < np; ++i )
284 tio_in.pointlist[j++] = points[i].
x;
285 tio_in.pointlist[j++] = points[i].
y;
290 tio_in.segmentlist = malloc( ns * 2 *
sizeof (
int ) );
291 tio_in.numberofsegments = ns;
292 memcpy( tio_in.segmentlist, segments, ns * 2 * sizeof (
int ) );
297 tio_in.holelist = malloc( nh * 2 *
sizeof (
double ) );
298 tio_in.numberofholes = nh;
299 memcpy( tio_in.holelist, holes, nh * 2 * sizeof (
double ) );
302 tio_init( &tio_out );
317 triangulate( cmd, &tio_in, &tio_out, NULL );
325 tio2delaunay( &tio_out, d );
327 tio_destroy( &tio_in );
328 tio_destroy( &tio_out );
337 boolT ismalloc = False;
338 char flags[64] =
"qhull d Qbb Qt";
339 facetT *facet, *neighbor, **neighborp;
340 vertexT *vertex, **vertexp;
342 int curlong, totlong;
343 FILE *outfile = stdout;
344 FILE *errfile = stderr;
349 int numfacets, numsimplicial, numridges, totneighbors, numcoplanars, numtricoplanars;
356 assert(
sizeof ( realT ) ==
sizeof (
double ) );
358 if ( np == 0 || ns > 0 || nh > 0 )
360 fprintf( stderr,
"segments=%d holes=%d\n, aborting Qhull implementation, use 'triangle' instead.\n", ns, nh );
365 qpoints = (coordT *) malloc( (
size_t) ( np * ( dim + 1 ) ) *
sizeof ( coordT ) );
367 for ( i = 0; i < np; i++ )
369 qpoints[i * dim] = points[i].
x;
370 qpoints[i * dim + 1] = points[i].
y;
376 strcat( flags,
" s" );
378 strcat( flags,
" Ts" );
387 exitcode = qh_new_qhull( dim, np, qpoints, ismalloc,
388 flags, outfile, errfile );
401 d->
points = malloc( (
size_t) np *
sizeof (
point ) );
402 for ( i = 0; i < np; ++i )
410 if ( p->
x < d->
xmin )
412 if ( p->
x > d->
xmax )
414 if ( p->
y < d->
ymin )
416 if ( p->
y > d->
ymax )
422 fprintf( stderr,
"input:\n" );
423 for ( i = 0; i < np; ++i )
427 fprintf( stderr,
" %d: %15.7g %15.7g %15.7g\n",
428 i, p->
x, p->
y, p->
z );
432 qh_findgood_all( qh facet_list );
433 qh_countfacets( qh facet_list, NULL, !qh_ALL, &numfacets,
434 &numsimplicial, &totneighbors, &numridges,
435 &numcoplanars, &numtricoplanars );
439 if ( !facet->upperdelaunay && facet->simplicial )
449 fprintf( stderr,
"triangles:\tneighbors:\n" );
453 if ( !facet->upperdelaunay && facet->simplicial )
460 FOREACHvertex_( facet->vertices )
461 t->vids[j++] = qh_pointid( vertex->
point );
464 FOREACHneighbor_( facet )
465 n->tids[j++] = ( neighbor->visitid > 0 ) ? (
int) neighbor->visitid - 1 : -1;
480 int tmp = t->vids[1];
481 t->vids[1] = t->vids[2];
493 fprintf( stderr,
" %d: (%d,%d,%d)\t(%d,%d,%d)\n",
508 for ( j = 0; j < 3; ++j )
512 for ( i = 0; i < d->
npoints; ++i )
524 for ( j = 0; j < 3; ++j )
526 int vid = t->
vids[j];
547 qh_freeqhull( !qh_ALL );
548 qh_memfreeshort( &curlong, &totlong );
549 if ( curlong || totlong )
551 "qhull: did not free %d bytes of long memory (%d pieces)\n",
564 return ( ( pb->
x - pa->
x ) * ( pc->
y - pa->
y ) <
565 ( pc->
x - pa->
x ) * ( pb->
y - pa->
y ) );
583 for ( i = 0; i < d->
npoints; ++i )
598 if ( d->
flags != NULL )
606 if ( d->
t_in != NULL )
608 if ( d->
t_out != NULL )
617 return ( p1->
x - p->
x ) * ( p0->
y - p->
y ) > ( p0->
x - p->
x ) * ( p1->
y - p->
y );
640 for ( i = 0; i < 3; ++i )
642 int i1 = ( i + 1 ) % 3;
681 if ( d->
t_in == NULL )
713 for ( i = 0; i < nn; ++i )
722 if ( tid < 0 || i == nn )
726 for ( tid = 0; tid < nt; ++tid )
752 while ( d->
t_in->
n > 0 )
760 for ( i = 0; i < 3; ++i )
762 int vid = t->
vids[i];
766 for ( j = 0; j < nt; ++j )
770 if ( d->
flags[ntid] == 0 )
int circle_build(circle *c, point *p0, point *p1, point *p2)
triangle_neighbours * neighbours
def cmd
Now do the PLplot API.
void istack_push(istack *s, int v)
static int on_right_side(point *p, point *p0, point *p1)
void delaunay_destroy(delaunay *d)
static int cw(delaunay *d, triangle *t)
delaunay * delaunay_build(int np, point points[], int ns, int segments[], int nh, double holes[])
void istack_destroy(istack *s)
int delaunay_xytoi(delaunay *d, point *p, int id)
void istack_reset(istack *s)
int circle_contains(circle *c, point *p)
int istack_pop(istack *s)
void delaunay_circles_find(delaunay *d, point *p, int *n, int **out)