Graphviz 13.1.3~dev.20250829.0113
Loading...
Searching...
No Matches
delaunay.c
Go to the documentation of this file.
1/*************************************************************************
2 * Copyright (c) 2011 AT&T Intellectual Property
3 * All rights reserved. This program and the accompanying materials
4 * are made available under the terms of the Eclipse Public License v1.0
5 * which accompanies this distribution, and is available at
6 * https://www.eclipse.org/legal/epl-v10.html
7 *
8 * Contributors: Details at https://graphviz.org
9 *************************************************************************/
10
11#include "config.h"
12
13#include <stdbool.h>
14#include <stddef.h>
15#include <stdio.h>
16#include <stdlib.h>
17#include <string.h>
18#include <math.h>
19#include <cgraph/cgraph.h> /* for agerr() and friends */
20#include <neatogen/delaunay.h>
21#include <util/alloc.h>
22#include <util/sort.h>
23
24#ifdef HAVE_GTS
25#include <gts.h>
26
27static int triangle_is_hole(void *triangle, void *ignored) {
28 GtsTriangle *t = triangle;
29 (void)ignored;
30
31 GtsEdge *e1, *e2, *e3;
32 GtsVertex *v1, *v2, *v3;
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 return TRUE;
40
41 return FALSE;
42}
43
44static unsigned delaunay_remove_holes(GtsSurface *surface) {
45 return gts_surface_foreach_face_remove(surface, triangle_is_hole, NULL);
46}
47
48/* Derived classes for vertices and faces so we can assign integer ids
49 * to make it easier to identify them. In particular, this allows the
50 * segments and faces to refer to vertices using the order in which
51 * they were passed in.
52 */
53typedef struct {
54 GtsVertex v;
55 int idx;
56} GVertex;
57
59static GVertex *downcast(GtsVertex *base) {
60 return (GVertex *)((char *)base - offsetof(GVertex, v));
61}
62
64static GtsVertex *upcast(GVertex *derived) {
65 return &derived->v;
66}
67
68typedef struct {
69 GtsVertexClass parent_class;
70} GVertexClass;
71
72static GtsVertexClass *g_vertex_class(void) {
73 static GVertexClass *klass = NULL;
74
75 if (klass == NULL) {
76 GtsObjectClassInfo vertex_info = {
77 .name = "GVertex",
78 .object_size = sizeof(GVertex),
79 .class_size = sizeof(GVertexClass),
80 };
81 klass = gts_object_class_new(GTS_OBJECT_CLASS(gts_vertex_class()),
82 &vertex_info);
83 }
84
85 return &klass->parent_class;
86}
87
88typedef struct {
89 GtsFace v;
90 int idx;
91} GFace;
92
93typedef struct {
94 GtsFaceClass parent_class;
95} GFaceClass;
96
97static GtsFaceClass *g_face_class(void) {
98 static GFaceClass *klass = NULL;
99
100 if (klass == NULL) {
101 GtsObjectClassInfo face_info = {
102 .name = "GFace",
103 .object_size = sizeof(GFace),
104 .class_size = sizeof(GFaceClass),
105 };
106 klass = gts_object_class_new(GTS_OBJECT_CLASS(gts_face_class()),
107 &face_info);
108 }
109
110 return &klass->parent_class;
111}
112
113/* destroy:
114 * Destroy each edge using v, then destroy v
115 */
116static void
117destroy (GtsVertex* v)
118{
119 GSList * i;
120
121 i = v->segments;
122 while (i) {
123 GSList * next = i->next;
124 gts_object_destroy (i->data);
125 i = next;
126 }
127 g_assert (v->segments == NULL);
128 gts_object_destroy(GTS_OBJECT(v));
129}
130
131/* tri:
132 * Main entry point to using GTS for triangulation.
133 * Input is npt points with x and y coordinates stored either separately
134 * in x[] and y[] (sepArr != 0) or consecutively in x[] (sepArr == 0).
135 * Optionally, the input can include nsegs line segments, whose endpoint
136 * indices are supplied in segs[2*i] and segs[2*i+1] yielding a constrained
137 * triangulation.
138 *
139 * The return value is the corresponding gts surface, which can be queries for
140 * the triangles and line segments composing the triangulation.
141 */
142static GtsSurface*
143tri(double *x, double *y, int npt, int *segs, int nsegs, int sepArr)
144{
145 int i;
146 GtsSurface *surface;
147 GVertex **vertices = gv_calloc(npt, sizeof(GVertex *));
148 GtsEdge **edges = gv_calloc(nsegs, sizeof(GtsEdge *));
149 GSList *list = NULL;
150 GtsVertex *v1, *v2, *v3;
151 GtsTriangle *t;
152 GtsVertexClass *vcl = g_vertex_class();
153 GtsEdgeClass *ecl = GTS_EDGE_CLASS (gts_constraint_class ());
154
155 if (sepArr) {
156 for (i = 0; i < npt; i++) {
157 GVertex *p = downcast(gts_vertex_new(vcl, x[i], y[i], 0));
158 p->idx = i;
159 vertices[i] = p;
160 }
161 }
162 else {
163 for (i = 0; i < npt; i++) {
164 GVertex *p = downcast(gts_vertex_new(vcl, x[2*i], x[2*i+1], 0));
165 p->idx = i;
166 vertices[i] = p;
167 }
168 }
169
170 /* N.B. Edges need to be created here, presumably before the
171 * the vertices are added to the face. In particular, they cannot
172 * be added created and added vi gts_delaunay_add_constraint() below.
173 */
174 for (i = 0; i < nsegs; i++) {
175 edges[i] = gts_edge_new(ecl,
176 upcast(vertices[segs[2 * i]]),
177 upcast(vertices[segs[2 * i + 1]]));
178 }
179
180 for (i = 0; i < npt; i++)
181 list = g_slist_prepend(list, vertices[i]);
182 t = gts_triangle_enclosing(gts_triangle_class(), list, 100.);
183 g_slist_free(list);
184
185 gts_triangle_vertices(t, &v1, &v2, &v3);
186
187 surface = gts_surface_new(gts_surface_class(), g_face_class(),
188 gts_edge_class(), gts_vertex_class());
189 gts_surface_add_face(surface, gts_face_new(gts_face_class(),
190 t->e1, t->e2, t->e3));
191
192 for (i = 0; i < npt; i++) {
193 GtsVertex *v4 = upcast(vertices[i]);
194 GtsVertex *v = gts_delaunay_add_vertex(surface, v4, NULL);
195
196 /* if v != NULL, it is a previously added pt with the same
197 * coordinates as v4, in which case we replace v4 with v
198 */
199 if (v && v4 != v) {
200 gts_vertex_replace(v4, v);
201 }
202 }
203
204 for (i = 0; i < nsegs; i++) {
205 gts_delaunay_add_constraint(surface,GTS_CONSTRAINT(edges[i]));
206 }
207
208 /* destroy enclosing triangle */
209 gts_allow_floating_vertices = TRUE;
210 gts_allow_floating_edges = TRUE;
211 destroy(v1);
212 destroy(v2);
213 destroy(v3);
214 gts_allow_floating_edges = FALSE;
215 gts_allow_floating_vertices = FALSE;
216
217 if (nsegs)
218 delaunay_remove_holes(surface);
219
220 free (edges);
221 free(vertices);
222 return surface;
223}
224
225typedef struct {
226 int n;
227 v_data *delaunay;
228} estats;
229
230static int cnt_edge(void *edge, void *stats) {
231 GtsSegment *e = edge;
232 estats *sp = stats;
233
234 sp->n++;
235 if (sp->delaunay) {
236 sp->delaunay[downcast(e->v1)->idx].nedges++;
237 sp->delaunay[downcast(e->v2)->idx].nedges++;
238 }
239
240 return 0;
241}
242
243static void
244edgeStats (GtsSurface* s, estats* sp)
245{
246 gts_surface_foreach_edge(s, cnt_edge, sp);
247}
248
249static int add_edge(void *edge, void *data) {
250 GtsSegment *e = edge;
251 v_data *delaunay = data;
252
253 int source = downcast(e->v1)->idx;
254 int dest = downcast(e->v2)->idx;
255
256 delaunay[source].edges[delaunay[source].nedges++] = dest;
257 delaunay[dest].edges[delaunay[dest].nedges++] = source;
258
259 return 0;
260}
261
262static v_data *delaunay_triangulation(double *x, double *y, int n) {
263 GtsSurface* s = tri(x, y, n, NULL, 0, 1);
264 int i, nedges;
265 estats stats;
266
267 if (!s) return NULL;
268
269 v_data *delaunay = gv_calloc(n, sizeof(v_data));
270
271 for (i = 0; i < n; i++) {
272 delaunay[i].ewgts = NULL;
273 delaunay[i].nedges = 1;
274 }
275
276 stats.n = 0;
277 stats.delaunay = delaunay;
278 edgeStats (s, &stats);
279 nedges = stats.n;
280 int *edges = gv_calloc(2 * nedges + n, sizeof(int));
281
282 for (i = 0; i < n; i++) {
283 delaunay[i].edges = edges;
284 edges += delaunay[i].nedges;
285 delaunay[i].edges[0] = i;
286 delaunay[i].nedges = 1;
287 }
288 gts_surface_foreach_edge(s, add_edge, delaunay);
289
290 gts_object_destroy (GTS_OBJECT (s));
291
292 return delaunay;
293}
294
295typedef struct {
296 int n;
297 int* edges;
298} estate;
299
300static int addEdge(void *edge, void *state) {
301 GtsSegment *e = edge;
302 estate *es = state;
303
304 int source = downcast(e->v1)->idx;
305 int dest = downcast(e->v2)->idx;
306
307 es->edges[2 * es->n] = source;
308 es->edges[2 * es->n + 1] = dest;
309 es->n += 1;
310
311 return 0;
312}
313
314static int vcmp(const void *x, const void *y, void *values) {
315 const int *a = x;
316 const int *b = y;
317 const double *_vals = values;
318 double va = _vals[*a];
319 double vb = _vals[*b];
320
321 if (va < vb) return -1;
322 if (va > vb) return 1;
323 return 0;
324}
325
326/* delaunay_tri:
327 * Given n points whose coordinates are in the x[] and y[]
328 * arrays, compute a Delaunay triangulation of the points.
329 * The number of edges in the triangulation is returned in pnedges.
330 * The return value itself is an array e of 2*(*pnedges) integers,
331 * with edge i having points whose indices are e[2*i] and e[2*i+1].
332 *
333 * If the points are collinear, GTS fails with 0 edges.
334 * In this case, we sort the points by x coordinates (or y coordinates
335 * if the points form a vertical line). We then return a "triangulation"
336 * consisting of the n-1 pairs of adjacent points.
337 */
338int *delaunay_tri(double *x, double *y, int n, int* pnedges)
339{
340 GtsSurface* s = tri(x, y, n, NULL, 0, 1);
341 int nedges;
342 int* edges;
343 estats stats;
344 estate state;
345
346 if (!s) return NULL;
347
348 stats.n = 0;
349 stats.delaunay = NULL;
350 edgeStats (s, &stats);
351 *pnedges = nedges = stats.n;
352
353 if (nedges) {
354 edges = gv_calloc(2 * nedges, sizeof(int));
355 state.n = 0;
356 state.edges = edges;
357 gts_surface_foreach_edge(s, addEdge, &state);
358 }
359 else {
360 int* vs = gv_calloc(n, sizeof(int));
361 int* ip;
362 int i, hd, tl;
363
364 *pnedges = nedges = n-1;
365 ip = edges = gv_calloc(2 * nedges, sizeof(int));
366
367 for (i = 0; i < n; i++)
368 vs[i] = i;
369
370 gv_sort(vs, n, sizeof(int), vcmp, x[0] == x[1] /* vertical line? */ ? y : x);
371
372 tl = vs[0];
373 for (i = 1; i < n; i++) {
374 hd = vs[i];
375 *ip++ = tl;
376 *ip++ = hd;
377 tl = hd;
378 }
379
380 free (vs);
381 }
382
383 gts_object_destroy (GTS_OBJECT (s));
384
385 return edges;
386}
387
388static int cntFace(void *face, void *data) {
389 GFace *fp = face;
390 int *ip = data;
391
392 fp->idx = *ip;
393 *ip += 1;
394
395 return 0;
396}
397
398typedef struct {
399 GtsSurface* s;
400 int* faces;
401 int* neigh;
402} fstate;
403
404typedef struct {
405 int nneigh;
406 int* neigh;
407} ninfo;
408
409static int addNeighbor(void *face, void *ni) {
410 GFace *f = face;
411 ninfo *es = ni;
412
413 es->neigh[es->nneigh] = f->idx;
414 es->nneigh++;
415
416 return 0;
417}
418
419static int addFace(void *face, void *state) {
420 GFace *f = face;
421 fstate *es = state;
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++ = downcast(v1)->idx;
431 *ip++ = downcast(v2)->idx;
432 *ip++ = downcast(v3)->idx;
433
434 ni.nneigh = 0;
435 ni.neigh = neigh;
436 gts_face_foreach_neighbor((GtsFace*)f, 0, addNeighbor, &ni);
437 for (i = ni.nneigh; i < 3; i++)
438 neigh[i] = -1;
439
440 return 0;
441}
442
443static int addTri(void *face, void *state) {
444 GFace *f = face;
445 fstate *es = state;
446
447 int myid = f->idx;
448 int* ip = es->faces + 3*myid;
449 GtsVertex *v1, *v2, *v3;
450
451 gts_triangle_vertices (&f->v.triangle, &v1, &v2, &v3);
452 *ip++ = downcast(v1)->idx;
453 *ip++ = downcast(v2)->idx;
454 *ip++ = downcast(v3)->idx;
455
456 return 0;
457}
458
459/* mkSurface:
460 * Given n points whose coordinates are in x[] and y[], and nsegs line
461 * segments whose end point indices are given in segs, return a surface
462 * corresponding the constrained Delaunay triangulation.
463 * The surface records the line segments, the triangles, and the neighboring
464 * triangles.
465 */
466surface_t*
467mkSurface (double *x, double *y, int n, int* segs, int nsegs)
468{
469 GtsSurface* s = tri(x, y, n, segs, nsegs, 1);
470 estats stats;
471 estate state;
472 fstate statf;
473 int nfaces = 0;
474
475 if (!s) return NULL;
476
477 surface_t *sf = gv_alloc(sizeof(surface_t));
478 stats.n = 0;
479 stats.delaunay = NULL;
480 edgeStats (s, &stats);
481 nsegs = stats.n;
482 segs = gv_calloc(2 * nsegs, sizeof(int));
483
484 state.n = 0;
485 state.edges = segs;
486 gts_surface_foreach_edge(s, addEdge, &state);
487
488 gts_surface_foreach_face(s, cntFace, &nfaces);
489
490 int *faces = gv_calloc(3 * nfaces, sizeof(int));
491 int *neigh = gv_calloc(3 * nfaces, sizeof(int));
492
493 statf.faces = faces;
494 statf.neigh = neigh;
495 gts_surface_foreach_face(s, addFace, &statf);
496
497 sf->nedges = nsegs;
498 sf->edges = segs;
499 sf->nfaces = nfaces;
500 sf->faces = faces;
501 sf->neigh = neigh;
502
503 gts_object_destroy (GTS_OBJECT (s));
504
505 return sf;
506}
507
508/* get_triangles:
509 * Given n points whose coordinates are stored as (x[2*i],x[2*i+1]),
510 * compute a Delaunay triangulation of the points.
511 * The number of triangles in the triangulation is returned in tris.
512 * The return value t is an array of 3*(*tris) integers,
513 * with triangle i having points whose indices are t[3*i], t[3*i+1] and t[3*i+2].
514 */
515int*
516get_triangles (double *x, int n, int* tris)
517{
518 GtsSurface* s;
519 int nfaces = 0;
520 fstate statf;
521
522 if (n <= 2) return NULL;
523
524 s = tri(x, NULL, n, NULL, 0, 0);
525 if (!s) return NULL;
526
527 gts_surface_foreach_face(s, cntFace, &nfaces);
528 statf.faces = gv_calloc(3 * nfaces, sizeof(int));
529 gts_surface_foreach_face(s, addTri, &statf);
530
531 gts_object_destroy (GTS_OBJECT (s));
532
533 *tris = nfaces;
534 return statf.faces;
535}
536
537void
539{
540 free (s->edges);
541 free (s->faces);
542 free (s->neigh);
543 free(s);
544}
545#else
546static char* err = "Graphviz built without any triangulation library\n";
547int* get_triangles (double *x, int n, int* tris)
548{
549 (void)x;
550 (void)n;
551 (void)tris;
552
553 agerrorf("get_triangles: %s\n", err);
554 return 0;
555}
556static v_data *delaunay_triangulation(double *x, double *y, int n) {
557 (void)x;
558 (void)y;
559 (void)n;
560
561 agerrorf("delaunay_triangulation: %s\n", err);
562 return 0;
563}
564int *delaunay_tri(double *x, double *y, int n, int* nedges)
565{
566 (void)x;
567 (void)y;
568 (void)n;
569 (void)nedges;
570
571 agerrorf("delaunay_tri: %s\n", err);
572 return 0;
573}
574surface_t*
575mkSurface (double *x, double *y, int n, int* segs, int nsegs)
576{
577 (void)x;
578 (void)y;
579 (void)n;
580 (void)segs;
581 (void)nsegs;
582
583 agerrorf("mkSurface: %s\n", err);
584 return 0;
585}
586void
588{
589 (void)s;
590
591 agerrorf("freeSurface: %s\n", err);
592}
593#endif
594
595static void remove_edge(v_data * graph, int source, int dest)
596{
597 int i;
598 for (i = 1; i < graph[source].nedges; i++) {
599 if (graph[source].edges[i] == dest) {
600 graph[source].edges[i] = graph[source].edges[--graph[source].nedges];
601 break;
602 }
603 }
604}
605
606v_data *UG_graph(double *x, double *y, int n) {
607 v_data *delaunay;
608 int i;
609 double dist_ij, dist_ik, dist_jk, x_i, y_i, x_j, y_j;
610 int j, k, neighbor_j, neighbor_k;
611
612 if (n == 2) {
613 int *edges = gv_calloc(4, sizeof(int));
614 delaunay = gv_calloc(n, sizeof(v_data));
615 delaunay[0].ewgts = NULL;
616 delaunay[0].edges = edges;
617 delaunay[0].nedges = 2;
618 delaunay[0].edges[0] = 0;
619 delaunay[0].edges[1] = 1;
620 delaunay[1].edges = edges + 2;
621 delaunay[1].ewgts = NULL;
622 delaunay[1].nedges = 2;
623 delaunay[1].edges[0] = 1;
624 delaunay[1].edges[1] = 0;
625 return delaunay;
626 } else if (n == 1) {
627 int *edges = gv_calloc(1, sizeof(int));
628 delaunay = gv_calloc(n, sizeof(v_data));
629 delaunay[0].ewgts = NULL;
630 delaunay[0].edges = edges;
631 delaunay[0].nedges = 1;
632 delaunay[0].edges[0] = 0;
633 return delaunay;
634 }
635
636 delaunay = delaunay_triangulation(x, y, n);
637
638 // 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)
639 for (i = 0; i < n; i++) {
640 x_i = x[i];
641 y_i = y[i];
642 for (j = 1; j < delaunay[i].nedges;) {
643 neighbor_j = delaunay[i].edges[j];
644 x_j = x[neighbor_j];
645 y_j = y[neighbor_j];
646 dist_ij = (x_j - x_i) * (x_j - x_i) + (y_j - y_i) * (y_j - y_i);
647 // now look at i'th neighbors to see whether there is a node in the "forbidden region"
648 // we will also go through neighbor_j's neighbors when we traverse the edge from its other side
649 bool removed = false;
650 for (k = 1; k < delaunay[i].nedges && !removed; k++) {
651 neighbor_k = delaunay[i].edges[k];
652 dist_ik = (x[neighbor_k] - x_i) * (x[neighbor_k] - x_i) +
653 (y[neighbor_k] - y_i) * (y[neighbor_k] - y_i);
654 if (dist_ik < dist_ij) {
655 dist_jk = (x[neighbor_k] - x_j) * (x[neighbor_k] - x_j) +
656 (y[neighbor_k] - y_j) * (y[neighbor_k] - y_j);
657 if (dist_jk < dist_ij) {
658 // remove the edge beteween i and neighbor j
659 delaunay[i].edges[j] = delaunay[i].edges[--delaunay[i].nedges];
660 remove_edge(delaunay, neighbor_j, i);
661 removed = true;
662 }
663 }
664 }
665 if (!removed) {
666 j++;
667 }
668 }
669 }
670 return delaunay;
671}
672
674{
675 if (graph != NULL) {
676 free(graph[0].edges);
677 free(graph[0].ewgts);
678 free(graph);
679 }
680}
681
683{
684 if (graph != NULL) {
685 free(graph[0].edges);
686 free(graph[0].ewgts);
687#ifdef DIGCOLA
688 free(graph[0].edists);
689#endif
690 free(graph);
691 }
692}
693
Memory allocation wrappers that exit on failure.
static void * gv_calloc(size_t nmemb, size_t size)
Definition alloc.h:26
static void * gv_alloc(size_t size)
Definition alloc.h:47
abstract graph C library, Cgraph API
int * get_triangles(double *x, int n, int *tris)
Definition delaunay.c:547
void freeGraph(v_data *graph)
Definition delaunay.c:673
v_data * UG_graph(double *x, double *y, int n)
Definition delaunay.c:606
void freeGraphData(vtx_data *graph)
Definition delaunay.c:682
static char * err
Definition delaunay.c:546
static v_data * delaunay_triangulation(double *x, double *y, int n)
Definition delaunay.c:556
surface_t * mkSurface(double *x, double *y, int n, int *segs, int nsegs)
Definition delaunay.c:575
static void remove_edge(v_data *graph, int source, int dest)
Definition delaunay.c:595
int * delaunay_tri(double *x, double *y, int n, int *nedges)
Definition delaunay.c:564
void freeSurface(surface_t *s)
Definition delaunay.c:587
void add_edge(edgelist *list, Agedge_t *e)
Definition edgelist.c:57
void free(void *)
edge
Definition gmlparse.y:244
node NULL
Definition grammar.y:181
void agerrorf(const char *fmt,...)
Definition agerror.c:165
Agraph_t * graph(char *name)
Definition gv.cpp:30
static void addEdge(edge_t *de, edge_t *e)
Definition layout.c:349
struct _tri tri
static tri * addTri(int i, int j, tri *oldp)
static int nedges
total no. of edges used in routing
Definition routespl.c:32
qsort with carried along context
static void gv_sort(void *base, size_t nmemb, size_t size, int(*compar)(const void *, const void *, void *), void *arg)
qsort with an extra state parameter, ala qsort_r
Definition sort.h:24
Definition legal.c:50
int nedges
Definition delaunay.h:22
int * edges
Definition delaunay.h:23
int * faces
Definition delaunay.h:25
int nfaces
Definition delaunay.h:24
int * neigh
Definition delaunay.h:26
int nedges
Definition sparsegraph.h:28
int * edges
Definition sparsegraph.h:29
float * ewgts
Definition sparsegraph.h:30
Definition grammar.c:90