Graphviz  2.41.20170921.2350
delaunay.c
Go to the documentation of this file.
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 
28 static 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 
44 static 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  */
55 typedef struct {
56  GtsVertex v;
57  int idx;
58 } GVertex;
59 
60 typedef struct {
61  GtsVertexClass parent_class;
62 } GVertexClass;
63 
64 static 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 
85 typedef struct {
86  GtsFace v;
87  int idx;
88 } GFace;
89 
90 typedef struct {
91  GtsFaceClass parent_class;
92 } GFaceClass;
93 
94 static 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  */
118 static void
119 destroy (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  */
144 static GtsSurface*
145 tri(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 
235 typedef struct {
236  int n;
237  v_data *delaunay;
238 } estats;
239 
240 static 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 
249 static void
250 edgeStats (GtsSurface* s, estats* sp)
251 {
252  gts_surface_foreach_edge (s, (GtsFunc) cnt_edge, sp);
253 }
254 
255 static 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 
264 v_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 
300 typedef struct {
301  int n;
302  int* edges;
303 } estate;
304 
305 static 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  */
318 static double* _vals;
319 typedef int (*qsort_cmpf) (const void *, const void *);
320 
321 static int
322 vcmp (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  */
344 int *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 
398 static void cntFace (GFace* fp, int* ip)
399 {
400  fp->idx = *ip;
401  *ip += 1;
402 }
403 
404 typedef struct {
405  GtsSurface* s;
406  int* faces;
407  int* neigh;
408 } fstate;
409 
410 typedef struct {
411  int nneigh;
412  int* neigh;
413 } ninfo;
414 
415 static void addNeighbor (GFace* f, ninfo* es)
416 {
417  es->neigh[es->nneigh] = f->idx;
418  es->nneigh++;
419 }
420 
421 static 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 
441 static 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  */
460 surface_t*
461 mkSurface (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  */
512 int*
513 get_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 
534 void
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 
547 int*
548 get_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
617 int*
618 delaunay_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 
692 v_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 
729 surface_t*
730 mkSurface (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 }
736 void
738 {
739  agerr (AGERR, "freeSurface not yet implemented using Triangle library\n");
740  assert (0);
741 }
742 #else
743 static char* err = "Graphviz built without any triangulation library\n";
744 int* get_triangles (double *x, int n, int* tris)
745 {
746  agerr(AGERR, "get_triangles: %s\n", err);
747  return 0;
748 }
749 v_data *delaunay_triangulation(double *x, double *y, int n)
750 {
751  agerr(AGERR, "delaunay_triangulation: %s\n", err);
752  return 0;
753 }
754 int *delaunay_tri(double *x, double *y, int n, int* nedges)
755 {
756  agerr(AGERR, "delaunay_tri: %s\n", err);
757  return 0;
758 }
759 surface_t*
760 mkSurface (double *x, double *y, int n, int* segs, int nsegs)
761 {
762  agerr(AGERR, "mkSurface: %s\n", err);
763  return 0;
764 }
765 void
767 {
768  agerr (AGERR, "freeSurface: %s\n", err);
769 }
770 #endif
771 
772 static 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 
784 v_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 
897 void 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 
908 void 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 
Definition: cgraph.h:388
v_data * delaunay_triangulation(double *x, double *y, int n)
Definition: delaunay.c:749
#define FREE
Definition: PriorityQueue.c:23
#define assert(x)
Definition: cghdr.h:47
int * faces
Definition: delaunay.h:23
int * neigh
Definition: delaunay.h:24
int agerr(agerrlevel_t level, const char *fmt,...)
Definition: agerror.c:141
void add_edge(edgelist *list, Agedge_t *e)
Definition: edgelist.c:65
struct _tri tri
void remove_edge(edgelist *list, Agedge_t *e)
Definition: edgelist.c:73
int
Definition: grammar.c:1264
int nfaces
Definition: delaunay.h:22
int(* qsort_cmpf)(const void *, const void *)
Definition: types.h:45
#define static
Definition: dtlist.c:171
Definition: grammar.c:79
int * delaunay_tri(double *x, double *y, int n, int *nedges)
Definition: delaunay.c:754
Agraph_t * graph(char *name)
Definition: gv.cpp:38
int * get_triangles(double *x, int n, int *tris)
Definition: delaunay.c:744
void freeSurface(surface_t *s)
Definition: delaunay.c:766
#define NULL
Definition: logic.h:39
#define GNEW(t)
Definition: memory.h:37
void freeGraphData(vtx_data *graph)
Definition: delaunay.c:908
surface_t * mkSurface(double *x, double *y, int n, int *segs, int nsegs)
Definition: delaunay.c:760
int nedges
Definition: delaunay.h:20
int * edges
Definition: delaunay.h:21
float * ewgts
Definition: sparsegraph.h:76
Definition: cdt.h:99
int * edges
Definition: sparsegraph.h:75
int nedges
Definition: sparsegraph.h:74
#define N_GNEW(n, t)
Definition: agxbuf.c:20
#define FALSE
Definition: cgraph.h:35
void freeGraph(v_data *graph)
Definition: delaunay.c:897
v_data * UG_graph(double *x, double *y, int n, int accurate_computation)
Definition: delaunay.c:784
#define TRUE
Definition: cgraph.h:38