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