1/* $Id$ $Revision$ */
2/* vim:set shiftwidth=4 ts=8: */
3
4/*************************************************************************
5 * Copyright (c) 2011 AT&T Intellectual Property
6 * All rights reserved. This program and the accompanying materials
7 * are made available under the terms of the Eclipse Public License v1.0
8 * which accompanies this distribution, and is available at
9 * http://www.eclipse.org/legal/epl-v10.html
10 *
11 * Contributors: See CVS logs. Details at http://www.graphviz.org/
12 *************************************************************************/
13
14#include "config.h"
15
16#include <stdio.h>
17#include <stdlib.h>
18#include <string.h>
19#include <math.h>
20#include "cgraph.h" /* for agerr() and friends */
21#include "delaunay.h"
22#include "memory.h"
23#include "logic.h"
24
25#if HAVE_GTS
26#include <gts.h>
27
28static gboolean triangle_is_hole(GtsTriangle * t)
29{
30 GtsEdge *e1, *e2, *e3;
31 GtsVertex *v1, *v2, *v3;
32 gboolean ret;
33
34 gts_triangle_vertices_edges(t, NULL, &v1, &v2, &v3, &e1, &e2, &e3);
35
36 if ((GTS_IS_CONSTRAINT(e1) && GTS_SEGMENT(e1)->v1 != v1) ||
37 (GTS_IS_CONSTRAINT(e2) && GTS_SEGMENT(e2)->v1 != v2) ||
38 (GTS_IS_CONSTRAINT(e3) && GTS_SEGMENT(e3)->v1 != v3))
39 ret = TRUE;
40 else ret = FALSE;
41 return ret;
42}
43
44static guint delaunay_remove_holes(GtsSurface * surface)
45{
46 return gts_surface_foreach_face_remove(surface,
47 (GtsFunc) triangle_is_hole, NULL);
48}
49
50/* Derived classes for vertices and faces so we can assign integer ids
51 * to make it easier to identify them. In particular, this allows the
52 * segments and faces to refer to vertices using the order in which
53 * they were passed in.
54 */
55typedef struct {
56 GtsVertex v;
57 int idx;
58} GVertex;
59
60typedef struct {
61 GtsVertexClass parent_class;
62} GVertexClass;
63
64static GVertexClass *g_vertex_class(void)
65{
66 static GVertexClass *klass = NULL;
67
68 if (klass == NULL) {
69 GtsObjectClassInfo vertex_info = {
70 "GVertex",
71 sizeof(GVertex),
72 sizeof(GVertexClass),
73 (GtsObjectClassInitFunc) NULL,
74 (GtsObjectInitFunc) NULL,
75 (GtsArgSetFunc) NULL,
76 (GtsArgGetFunc) NULL
77 };
78 klass = gts_object_class_new(GTS_OBJECT_CLASS(gts_vertex_class()),
79 &vertex_info);
80 }
81
82 return klass;
83}
84
85typedef struct {
86 GtsFace v;
87 int idx;
88} GFace;
89
90typedef struct {
91 GtsFaceClass parent_class;
92} GFaceClass;
93
94static GFaceClass *g_face_class(void)
95{
96 static GFaceClass *klass = NULL;
97
98 if (klass == NULL) {
99 GtsObjectClassInfo face_info = {
100 "GFace",
101 sizeof(GFace),
102 sizeof(GFaceClass),
103 (GtsObjectClassInitFunc) NULL,
104 (GtsObjectInitFunc) NULL,
105 (GtsArgSetFunc) NULL,
106 (GtsArgGetFunc) NULL
107 };
108 klass = gts_object_class_new(GTS_OBJECT_CLASS(gts_face_class()),
109 &face_info);
110 }
111
112 return klass;
113}
114
115/* destroy:
116 * Destroy each edge using v, then destroy v
117 */
118static void
119destroy (GtsVertex* v)
120{
121 GSList * i;
122
123 i = v->segments;
124 while (i) {
125 GSList * next = i->next;
126 gts_object_destroy (i->data);
127 i = next;
128 }
129 g_assert (v->segments == NULL);
130 gts_object_destroy(GTS_OBJECT(v));
131}
132
133/* tri:
134 * Main entry point to using GTS for triangulation.
135 * Input is npt points with x and y coordinates stored either separately
136 * in x[] and y[] (sepArr != 0) or consecutively in x[] (sepArr == 0).
137 * Optionally, the input can include nsegs line segments, whose endpoint
138 * indices are supplied in segs[2*i] and segs[2*i+1] yielding a constrained
139 * triangulation.
140 *
141 * The return value is the corresponding gts surface, which can be queries for
142 * the triangles and line segments composing the triangulation.
143 */
144static GtsSurface*
145tri(double *x, double *y, int npt, int *segs, int nsegs, int sepArr)
146{
147 int i;
148 GtsSurface *surface;
149 GVertex **vertices = N_GNEW(npt, GVertex *);
150 GtsEdge **edges = N_GNEW(nsegs, GtsEdge*);
151 GSList *list = NULL;
152 GtsVertex *v1, *v2, *v3;
153 GtsTriangle *t;
154 GtsVertexClass *vcl = (GtsVertexClass *) g_vertex_class();
155 GtsEdgeClass *ecl = GTS_EDGE_CLASS (gts_constraint_class ());
156
157 if (sepArr) {
158 for (i = 0; i < npt; i++) {
159 GVertex *p = (GVertex *) gts_vertex_new(vcl, x[i], y[i], 0);
160 p->idx = i;
161 vertices[i] = p;
162 }
163 }
164 else {
165 for (i = 0; i < npt; i++) {
166 GVertex *p = (GVertex *) gts_vertex_new(vcl, x[2*i], x[2*i+1], 0);
167 p->idx = i;
168 vertices[i] = p;
169 }
170 }
171
172 /* N.B. Edges need to be created here, presumably before the
173 * the vertices are added to the face. In particular, they cannot
174 * be added created and added vi gts_delaunay_add_constraint() below.
175 */
176 for (i = 0; i < nsegs; i++) {
177 edges[i] = gts_edge_new(ecl,
178 (GtsVertex *) (vertices[ segs[ 2 * i]]),
179 (GtsVertex *) (vertices[ segs[ 2 * i + 1]]));
180 }
181
182 for (i = 0; i < npt; i++)
183 list = g_slist_prepend(list, vertices[i]);
184 t = gts_triangle_enclosing(gts_triangle_class(), list, 100.);
185 g_slist_free(list);
186
187 gts_triangle_vertices(t, &v1, &v2, &v3);
188
189 surface = gts_surface_new(gts_surface_class(),
190 (GtsFaceClass *) g_face_class(),
191 gts_edge_class(),
192 gts_vertex_class());
193 gts_surface_add_face(surface, gts_face_new(gts_face_class(),
194 t->e1, t->e2, t->e3));
195
196 for (i = 0; i < npt; i++) {
197 GtsVertex *v1 = (GtsVertex *) vertices[i];
198 GtsVertex *v = gts_delaunay_add_vertex(surface, v1, NULL);
199
200 /* if v != NULL, it is a previously added pt with the same
201 * coordinates as v1, in which case we replace v1 with v
202 */
203 if (v) {
204 /* agerr (AGWARN, "Duplicate point %d %d\n", i, ((GVertex*)v)->idx); */
205 gts_vertex_replace (v1, v);
206 }
207 }
208
209 for (i = 0; i < nsegs; i++) {
210 gts_delaunay_add_constraint(surface,GTS_CONSTRAINT(edges[i]));
211 }
212
213 /* destroy enclosing triangle */
214 gts_allow_floating_vertices = TRUE;
215 gts_allow_floating_edges = TRUE;
216/*
217 gts_object_destroy(GTS_OBJECT(v1));
218 gts_object_destroy(GTS_OBJECT(v2));
219 gts_object_destroy(GTS_OBJECT(v3));
220*/
221 destroy(v1);
222 destroy(v2);
223 destroy(v3);
224 gts_allow_floating_edges = FALSE;
225 gts_allow_floating_vertices = FALSE;
226
227 if (nsegs)
228 delaunay_remove_holes(surface);
229
230 free (edges);
231 free(vertices);
232 return surface;
233}
234
235typedef struct {
236 int n;
237 v_data *delaunay;
238} estats;
239
240static void cnt_edge (GtsSegment * e, estats* sp)
241{
242 sp->n++;
243 if (sp->delaunay) {
244 sp->delaunay[((GVertex*)(e->v1))->idx].nedges++;
245 sp->delaunay[((GVertex*)(e->v2))->idx].nedges++;
246 }
247}
248
249static void
250edgeStats (GtsSurface* s, estats* sp)
251{
252 gts_surface_foreach_edge (s, (GtsFunc) cnt_edge, sp);
253}
254
255static void add_edge (GtsSegment * e, v_data* delaunay)
256{
257 int source = ((GVertex*)(e->v1))->idx;
258 int dest = ((GVertex*)(e->v2))->idx;
259
260 delaunay[source].edges[delaunay[source].nedges++] = dest;
261 delaunay[dest].edges[delaunay[dest].nedges++] = source;
262}
263
264v_data *delaunay_triangulation(double *x, double *y, int n)
265{
266 v_data *delaunay;
267 GtsSurface* s = tri(x, y, n, NULL, 0, 1);
268 int i, nedges;
269 int* edges;
270 estats stats;
271
272 if (!s) return NULL;
273
274 delaunay = N_GNEW(n, v_data);
275
276 for (i = 0; i < n; i++) {
277 delaunay[i].ewgts = NULL;
278 delaunay[i].nedges = 1;
279 }
280
281 stats.n = 0;
282 stats.delaunay = delaunay;
283 edgeStats (s, &stats);
284 nedges = stats.n;
285 edges = N_GNEW(2 * nedges + n, int);
286
287 for (i = 0; i < n; i++) {
288 delaunay[i].edges = edges;
289 edges += delaunay[i].nedges;
290 delaunay[i].edges[0] = i;
291 delaunay[i].nedges = 1;
292 }
293 gts_surface_foreach_edge (s, (GtsFunc) add_edge, delaunay);
294
295 gts_object_destroy (GTS_OBJECT (s));
296
297 return delaunay;
298}
299
300typedef struct {
301 int n;
302 int* edges;
303} estate;
304
305static void addEdge (GtsSegment * e, estate* es)
306{
307 int source = ((GVertex*)(e->v1))->idx;
308 int dest = ((GVertex*)(e->v2))->idx;
309
310 es->edges[2*(es->n)] = source;
311 es->edges[2*(es->n)+1] = dest;
312 es->n += 1;
313}
314
315/* If qsort_r ever becomes standardized, this should be used
316 * instead of having a global variable.
317 */
318static double* _vals;
319typedef int (*qsort_cmpf) (const void *, const void *);
320
321static int
322vcmp (int* a, int* b)
323{
324 double va = _vals[*a];
325 double vb = _vals[*b];
326
327 if (va < vb) return -1;
328 else if (va > vb) return 1;
329 else return 0;
330}
331
332/* delaunay_tri:
333 * Given n points whose coordinates are in the x[] and y[]
334 * arrays, compute a Delaunay triangulation of the points.
335 * The number of edges in the triangulation is returned in pnedges.
336 * The return value itself is an array e of 2*(*pnedges) integers,
337 * with edge i having points whose indices are e[2*i] and e[2*i+1].
338 *
339 * If the points are collinear, GTS fails with 0 edges.
340 * In this case, we sort the points by x coordinates (or y coordinates
341 * if the points form a vertical line). We then return a "triangulation"
342 * consisting of the n-1 pairs of adjacent points.
343 */
344int *delaunay_tri(double *x, double *y, int n, int* pnedges)
345{
346 GtsSurface* s = tri(x, y, n, NULL, 0, 1);
347 int nedges;
348 int* edges;
349 estats stats;
350 estate state;
351
352 if (!s) return NULL;
353
354 stats.n = 0;
355 stats.delaunay = NULL;
356 edgeStats (s, &stats);
357 *pnedges = nedges = stats.n;
358
359 if (nedges) {
360 edges = N_GNEW(2 * nedges, int);
361 state.n = 0;
362 state.edges = edges;
363 gts_surface_foreach_edge (s, (GtsFunc) addEdge, &state);
364 }
365 else {
366 int* vs = N_GNEW(n, int);
367 int* ip;
368 int i, hd, tl;
369
370 *pnedges = nedges = n-1;
371 ip = edges = N_GNEW(2 * nedges, int);
372
373 for (i = 0; i < n; i++)
374 vs[i] = i;
375
376 if (x[0] == x[1]) /* vertical line */
377 _vals = y;
378 else
379 _vals = x;
380 qsort (vs, n, sizeof(int), (qsort_cmpf)vcmp);
381
382 tl = vs[0];
383 for (i = 1; i < n; i++) {
384 hd = vs[i];
385 *ip++ = tl;
386 *ip++ = hd;
387 tl = hd;
388 }
389
390 free (vs);
391 }
392
393 gts_object_destroy (GTS_OBJECT (s));
394
395 return edges;
396}
397
398static void cntFace (GFace* fp, int* ip)
399{
400 fp->idx = *ip;
401 *ip += 1;
402}
403
404typedef struct {
405 GtsSurface* s;
406 int* faces;
407 int* neigh;
408} fstate;
409
410typedef struct {
411 int nneigh;
412 int* neigh;
413} ninfo;
414
415static void addNeighbor (GFace* f, ninfo* es)
416{
417 es->neigh[es->nneigh] = f->idx;
418 es->nneigh++;
419}
420
421static void addFace (GFace* f, fstate* es)
422{
423 int i, myid = f->idx;
424 int* ip = es->faces + 3*myid;
425 int* neigh = es->neigh + 3*myid;
426 ninfo ni;
427 GtsVertex *v1, *v2, *v3;
428
429 gts_triangle_vertices (&f->v.triangle, &v1, &v2, &v3);
430 *ip++ = ((GVertex*)(v1))->idx;
431 *ip++ = ((GVertex*)(v2))->idx;
432 *ip++ = ((GVertex*)(v3))->idx;
433
434 ni.nneigh = 0;
435 ni.neigh = neigh;
436 gts_face_foreach_neighbor ((GtsFace*)f, 0, (GtsFunc) addNeighbor, &ni);
437 for (i = ni.nneigh; i < 3; i++)
438 neigh[i] = -1;
439}
440
441static void addTri (GFace* f, fstate* es)
442{
443 int myid = f->idx;
444 int* ip = es->faces + 3*myid;
445 GtsVertex *v1, *v2, *v3;
446
447 gts_triangle_vertices (&f->v.triangle, &v1, &v2, &v3);
448 *ip++ = ((GVertex*)(v1))->idx;
449 *ip++ = ((GVertex*)(v2))->idx;
450 *ip++ = ((GVertex*)(v3))->idx;
451}
452
453/* mkSurface:
454 * Given n points whose coordinates are in x[] and y[], and nsegs line
455 * segments whose end point indices are given in segs, return a surface
456 * corresponding the constrained Delaunay triangulation.
457 * The surface records the line segments, the triangles, and the neighboring
458 * triangles.
459 */
460surface_t*
461mkSurface (double *x, double *y, int n, int* segs, int nsegs)
462{
463 GtsSurface* s = tri(x, y, n, segs, nsegs, 1);
464 estats stats;
465 estate state;
466 fstate statf;
467 surface_t* sf;
468 int nfaces = 0;
469 int* faces;
470 int* neigh;
471
472 if (!s) return NULL;
473
474 sf = GNEW(surface_t);
475 stats.n = 0;
476 stats.delaunay = NULL;
477 edgeStats (s, &stats);
478 nsegs = stats.n;
479 segs = N_GNEW(2 * nsegs, int);
480
481 state.n = 0;
482 state.edges = segs;
483 gts_surface_foreach_edge (s, (GtsFunc) addEdge, &state);
484
485 gts_surface_foreach_face (s, (GtsFunc) cntFace, &nfaces);
486
487 faces = N_GNEW(3 * nfaces, int);
488 neigh = N_GNEW(3 * nfaces, int);
489
490 statf.faces = faces;
491 statf.neigh = neigh;
492 gts_surface_foreach_face (s, (GtsFunc) addFace, &statf);
493
494 sf->nedges = nsegs;
495 sf->edges = segs;
496 sf->nfaces = nfaces;
497 sf->faces = faces;
498 sf->neigh = neigh;
499
500 gts_object_destroy (GTS_OBJECT (s));
501
502 return sf;
503}
504
505/* get_triangles:
506 * Given n points whose coordinates are stored as (x[2*i],x[2*i+1]),
507 * compute a Delaunay triangulation of the points.
508 * The number of triangles in the triangulation is returned in tris.
509 * The return value t is an array of 3*(*tris) integers,
510 * with triangle i having points whose indices are t[3*i], t[3*i+1] and t[3*i+2].
511 */
512int*
513get_triangles (double *x, int n, int* tris)
514{
515 GtsSurface* s;
516 int nfaces = 0;
517 fstate statf;
518
519 if (n <= 2) return NULL;
520
521 s = tri(x, NULL, n, NULL, 0, 0);
522 if (!s) return NULL;
523
524 gts_surface_foreach_face (s, (GtsFunc) cntFace, &nfaces);
525 statf.faces = N_GNEW(3 * nfaces, int);
526 gts_surface_foreach_face (s, (GtsFunc) addTri, &statf);
527
528 gts_object_destroy (GTS_OBJECT (s));
529
530 *tris = nfaces;
531 return statf.faces;
532}
533
534void
535freeSurface (surface_t* s)
536{
537 free (s->edges);
538 free (s->faces);
539 free (s->neigh);
540}
541#elif HAVE_TRIANGLE
542#define TRILIBRARY
543#include "triangle.c"
544#include "assert.h"
545#include "general.h"
546
547int*
548get_triangles (double *x, int n, int* tris)
549{
550 struct triangulateio in, mid, vorout;
551 int i;
552
553 if (n <= 2) return NULL;
554
555 in.numberofpoints = n;
556 in.numberofpointattributes = 0;
557 in.pointlist = (REAL *) N_GNEW(in.numberofpoints * 2, REAL);
558
559 for (i = 0; i < n; i++){
560 in.pointlist[i*2] = x[i*2];
561 in.pointlist[i*2 + 1] = x[i*2 + 1];
562 }
563 in.pointattributelist = NULL;
564 in.pointmarkerlist = NULL;
565 in.numberofsegments = 0;
566 in.numberofholes = 0;
567 in.numberofregions = 0;
568 in.regionlist = NULL;
569 mid.pointlist = (REAL *) NULL; /* Not needed if -N switch used. */
570 mid.pointattributelist = (REAL *) NULL;
571 mid.pointmarkerlist = (int *) NULL; /* Not needed if -N or -B switch used. */
572 mid.trianglelist = (int *) NULL; /* Not needed if -E switch used. */
573 mid.triangleattributelist = (REAL *) NULL;
574 mid.neighborlist = (int *) NULL; /* Needed only if -n switch used. */
575 mid.segmentlist = (int *) NULL;
576 mid.segmentmarkerlist = (int *) NULL;
577 mid.edgelist = (int *) NULL; /* Needed only if -e switch used. */
578 mid.edgemarkerlist = (int *) NULL; /* Needed if -e used and -B not used. */
579 vorout.pointlist = (REAL *) NULL; /* Needed only if -v switch used. */
580 vorout.pointattributelist = (REAL *) NULL;
581 vorout.edgelist = (int *) NULL; /* Needed only if -v switch used. */
582 vorout.normlist = (REAL *) NULL; /* Needed only if -v switch used. */
583
584 /* Triangulate the points. Switches are chosen to read and write a */
585 /* PSLG (p), preserve the convex hull (c), number everything from */
586 /* zero (z), assign a regional attribute to each element (A), and */
587 /* produce an edge list (e), a Voronoi diagram (v), and a triangle */
588 /* neighbor list (n). */
589
590 triangulate("Qenv", &in, &mid, &vorout);
591 assert (mid.numberofcorners == 3);
592
593 *tris = mid.numberoftriangles;
594
595 FREE(in.pointlist);
596 FREE(in.pointattributelist);
597 FREE(in.pointmarkerlist);
598 FREE(in.regionlist);
599 FREE(mid.pointlist);
600 FREE(mid.pointattributelist);
601 FREE(mid.pointmarkerlist);
602 FREE(mid.triangleattributelist);
603 FREE(mid.neighborlist);
604 FREE(mid.segmentlist);
605 FREE(mid.segmentmarkerlist);
606 FREE(mid.edgelist);
607 FREE(mid.edgemarkerlist);
608 FREE(vorout.pointlist);
609 FREE(vorout.pointattributelist);
610 FREE(vorout.edgelist);
611 FREE(vorout.normlist);
612
613 return mid.trianglelist;
614}
615
616// maybe it should be replaced by RNG - relative neighborhood graph, or by GG - gabriel graph
617int*
618delaunay_tri (double *x, double *y, int n, int* nedges)
619{
620 struct triangulateio in, out;
621 int i;
622
623 in.pointlist = N_GNEW(2 * n, REAL);
624 for (i = 0; i < n; i++) {
625 in.pointlist[2 * i] = x[i];
626 in.pointlist[2 * i + 1] = y[i];
627 }
628
629 in.pointattributelist = NULL;
630 in.pointmarkerlist = NULL;
631 in.numberofpoints = n;
632 in.numberofpointattributes = 0;
633 in.trianglearealist = NULL;
634 in.triangleattributelist = NULL;
635 in.numberoftriangleattributes = 0;
636 in.neighborlist = NULL;
637 in.segmentlist = NULL;
638 in.segmentmarkerlist = NULL;
639 in.holelist = NULL;
640 in.numberofholes = 0;
641 in.regionlist = NULL;
642 in.edgelist = NULL;
643 in.edgemarkerlist = NULL;
644 in.normlist = NULL;
645
646 out.pointattributelist = NULL;
647 out.pointmarkerlist = NULL;
648 out.numberofpoints = n;
649 out.numberofpointattributes = 0;
650 out.trianglearealist = NULL;
651 out.triangleattributelist = NULL;
652 out.numberoftriangleattributes = 0;
653 out.neighborlist = NULL;
654 out.segmentlist = NULL;
655 out.segmentmarkerlist = NULL;
656 out.holelist = NULL;
657 out.numberofholes = 0;
658 out.regionlist = NULL;
659 out.edgelist = NULL;
660 out.edgemarkerlist = NULL;
661 out.normlist = NULL;
662
663 triangulate("zQNEeB", &in, &out, NULL);
664
665 *nedges = out.numberofedges;
666 free (in.pointlist);
667 free (in.pointattributelist);
668 free (in.pointmarkerlist);
669 free (in.trianglearealist);
670 free (in.triangleattributelist);
671 free (in.neighborlist);
672 free (in.segmentlist);
673 free (in.segmentmarkerlist);
674 free (in.holelist);
675 free (in.regionlist);
676 free (in.edgemarkerlist);
677 free (in.normlist);
678 free (out.pointattributelist);
679 free (out.pointmarkerlist);
680 free (out.trianglearealist);
681 free (out.triangleattributelist);
682 free (out.neighborlist);
683 free (out.segmentlist);
684 free (out.segmentmarkerlist);
685 free (out.holelist);
686 free (out.regionlist);
687 free (out.edgemarkerlist);
688 free (out.normlist);
689 return out.edgelist;
690}
691
692v_data *delaunay_triangulation(double *x, double *y, int n)
693{
694 v_data *delaunay;
695 int nedges;
696 int *edges;
697 int source, dest;
698 int* edgelist = delaunay_tri (x, y, n, &nedges);
699 int i;
700
701 delaunay = N_GNEW(n, v_data);
702 edges = N_GNEW(2 * nedges + n, int);
703
704 for (i = 0; i < n; i++) {
705 delaunay[i].ewgts = NULL;
706 delaunay[i].nedges = 1;
707 }
708
709 for (i = 0; i < 2 * nedges; i++)
710 delaunay[edgelist[i]].nedges++;
711
712 for (i = 0; i < n; i++) {
713 delaunay[i].edges = edges;
714 edges += delaunay[i].nedges;
715 delaunay[i].edges[0] = i;
716 delaunay[i].nedges = 1;
717 }
718 for (i = 0; i < nedges; i++) {
719 source = edgelist[2 * i];
720 dest = edgelist[2 * i + 1];
721 delaunay[source].edges[delaunay[source].nedges++] = dest;
722 delaunay[dest].edges[delaunay[dest].nedges++] = source;
723 }
724
725 free(edgelist);
726 return delaunay;
727}
728
729surface_t*
730mkSurface (double *x, double *y, int n, int* segs, int nsegs)
731{
732 agerr (AGERR, "mkSurface not yet implemented using Triangle library\n");
733 assert (0);
734 return 0;
735}
736void
737freeSurface (surface_t* s)
738{
739 agerr (AGERR, "freeSurface not yet implemented using Triangle library\n");
740 assert (0);
741}
742#else
743static char* err = "Graphviz built without any triangulation library\n";
744int* get_triangles (double *x, int n, int* tris)
745{
746 agerr(AGERR, "get_triangles: %s\n", err);
747 return 0;
748}
749v_data *delaunay_triangulation(double *x, double *y, int n)
750{
751 agerr(AGERR, "delaunay_triangulation: %s\n", err);
752 return 0;
753}
754int *delaunay_tri(double *x, double *y, int n, int* nedges)
755{
756 agerr(AGERR, "delaunay_tri: %s\n", err);
757 return 0;
758}
759surface_t*
760mkSurface (double *x, double *y, int n, int* segs, int nsegs)
761{
762 agerr(AGERR, "mkSurface: %s\n", err);
763 return 0;
764}
765void
766freeSurface (surface_t* s)
767{
768 agerr (AGERR, "freeSurface: %s\n", err);
769}
770#endif
771
772static void remove_edge(v_data * graph, int source, int dest)
773{
774 int i;
775 for (i = 1; i < graph[source].nedges; i++) {
776 if (graph[source].edges[i] == dest) {
777 graph[source].edges[i] =
778 graph[source].edges[--graph[source].nedges];
779 break;
780 }
781 }
782}
783
784v_data *UG_graph(double *x, double *y, int n, int accurate_computation)
785{
786 v_data *delaunay;
787 int i;
788 double dist_ij, dist_ik, dist_jk, x_i, y_i, x_j, y_j;
789 int j, k, neighbor_j, neighbor_k;
790 int removed;
791
792 if (n == 2) {
793 int *edges = N_GNEW(4, int);
794 delaunay = N_GNEW(n, v_data);
795 delaunay[0].ewgts = NULL;
796 delaunay[0].edges = edges;
797 delaunay[0].nedges = 2;
798 delaunay[0].edges[0] = 0;
799 delaunay[0].edges[1] = 1;
800 delaunay[1].edges = edges + 2;
801 delaunay[1].ewgts = NULL;
802 delaunay[1].nedges = 2;
803 delaunay[1].edges[0] = 1;
804 delaunay[1].edges[1] = 0;
805 return delaunay;
806 } else if (n == 1) {
807 int *edges = N_GNEW(1, int);
808 delaunay = N_GNEW(n, v_data);
809 delaunay[0].ewgts = NULL;
810 delaunay[0].edges = edges;
811 delaunay[0].nedges = 1;
812 delaunay[0].edges[0] = 0;
813 return delaunay;
814 }
815
816 delaunay = delaunay_triangulation(x, y, n);
817
818 if (accurate_computation) {
819 for (i = 0; i < n; i++) {
820 x_i = x[i];
821 y_i = y[i];
822 for (j = 1; j < delaunay[i].nedges;) {
823 neighbor_j = delaunay[i].edges[j];
824 if (neighbor_j < i) {
825 j++;
826 continue;
827 }
828 x_j = x[neighbor_j];
829 y_j = y[neighbor_j];
830 dist_ij =
831 (x_j - x_i) * (x_j - x_i) + (y_j - y_i) * (y_j - y_i);
832 removed = FALSE;
833 for (k = 0; k < n && !removed; k++) {
834 dist_ik =
835 (x[k] - x_i) * (x[k] - x_i) + (y[k] -
836 y_i) * (y[k] - y_i);
837 if (dist_ik < dist_ij) {
838 dist_jk =
839 (x[k] - x_j) * (x[k] - x_j) + (y[k] -
840 y_j) * (y[k] -
841 y_j);
842 if (dist_jk < dist_ij) {
843 // remove the edge beteween i and neighbor j
844 delaunay[i].edges[j] =
845 delaunay[i].edges[--delaunay[i].nedges];
846 remove_edge(delaunay, neighbor_j, i);
847 removed = TRUE;
848 }
849 }
850 }
851 if (!removed) {
852 j++;
853 }
854 }
855 }
856 } else {
857 // remove all edges v-u if there is w, neighbor of u or v, that is closer to both u and v than dist(u,v)
858 for (i = 0; i < n; i++) {
859 x_i = x[i];
860 y_i = y[i];
861 for (j = 1; j < delaunay[i].nedges;) {
862 neighbor_j = delaunay[i].edges[j];
863 x_j = x[neighbor_j];
864 y_j = y[neighbor_j];
865 dist_ij =
866 (x_j - x_i) * (x_j - x_i) + (y_j - y_i) * (y_j - y_i);
867 // now look at i'th neighbors to see whether there is a node in the "forbidden region"
868 // we will also go through neighbor_j's neighbors when we traverse the edge from its other side
869 removed = FALSE;
870 for (k = 1; k < delaunay[i].nedges && !removed; k++) {
871 neighbor_k = delaunay[i].edges[k];
872 dist_ik =
873 (x[neighbor_k] - x_i) * (x[neighbor_k] - x_i) +
874 (y[neighbor_k] - y_i) * (y[neighbor_k] - y_i);
875 if (dist_ik < dist_ij) {
876 dist_jk =
877 (x[neighbor_k] - x_j) * (x[neighbor_k] - x_j) +
878 (y[neighbor_k] - y_j) * (y[neighbor_k] - y_j);
879 if (dist_jk < dist_ij) {
880 // remove the edge beteween i and neighbor j
881 delaunay[i].edges[j] =
882 delaunay[i].edges[--delaunay[i].nedges];
883 remove_edge(delaunay, neighbor_j, i);
884 removed = TRUE;
885 }
886 }
887 }
888 if (!removed) {
889 j++;
890 }
891 }
892 }
893 }
894 return delaunay;
895}
896
897void freeGraph (v_data * graph)
898{
899 if (graph != NULL) {
900 if (graph[0].edges != NULL)
901 free(graph[0].edges);
902 if (graph[0].ewgts != NULL)
903 free(graph[0].ewgts);
904 free(graph);
905 }
906}
907
908void freeGraphData(vtx_data * graph)
909{
910 if (graph != NULL) {
911 if (graph[0].edges != NULL)
912 free(graph[0].edges);
913 if (graph[0].ewgts != NULL)
914 free(graph[0].ewgts);
915#ifdef USE_STYLES
916 if (graph[0].styles != NULL)
917 free(graph[0].styles);
918#endif
919#ifdef DIGCOLA
920 if (graph[0].edists != NULL)
921 free(graph[0].edists);
922#endif
923 free(graph);
924 }
925}
926
927