GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-112dd97adf
dgraph.c
Go to the documentation of this file.
1 /*!
2  \file lib/vector/Vlib/dgraph.c
3 
4  \brief Vector library - intersection (lower level functions)
5 
6  Higher level functions for reading/writing/manipulating vectors.
7 
8  (C) 2008-2009 by the GRASS Development Team
9 
10  This program is free software under the GNU General Public License
11  (>=v2). Read the file COPYING that comes with GRASS for details.
12 
13  \author Rewritten by Rosen Matev (Google Summer of Code 2008)
14  */
15 
16 #include <stdlib.h>
17 #include <string.h>
18 #include <math.h>
19 #include <grass/gis.h>
20 #include <grass/vector.h>
21 #include <grass/glocale.h>
22 #include "dgraph.h"
23 #include "e_intersect.h"
24 
25 #define LENGTH(DX, DY) (sqrt((DX * DX) + (DY * DY)))
26 #define PI M_PI
27 
28 struct intersection_point {
29  double x;
30  double y;
31  int group; /* IPs with very similar dist will be in the same group */
32 };
33 
34 struct seg_intersection {
35  int with; /* second segment */
36  int ip; /* index of the IP */
37  double dist; /* distance from first point of first segment to intersection
38  point (IP) */
39 };
40 
41 struct seg_intersection_list {
42  int count;
43  int allocated;
44  struct seg_intersection *a;
45 };
46 
47 struct seg_intersections {
48  int ipcount;
49  int ipallocated;
50  struct intersection_point *ip;
51  int ilcount;
52  struct seg_intersection_list *il;
53  int group_count;
54 };
55 
56 struct seg_intersections *create_si_struct(int segments_count)
57 {
58  struct seg_intersections *si;
59 
60  int i;
61 
62  si = G_malloc(sizeof(struct seg_intersections));
63  si->ipcount = 0;
64  si->ipallocated = segments_count + 16;
65  si->ip = G_malloc((si->ipallocated) * sizeof(struct intersection_point));
66  si->ilcount = segments_count;
67  si->il = G_malloc(segments_count * sizeof(struct seg_intersection_list));
68  for (i = 0; i < segments_count; i++) {
69  si->il[i].count = 0;
70  si->il[i].allocated = 0;
71  si->il[i].a = NULL;
72  }
73 
74  return si;
75 }
76 
77 void destroy_si_struct(struct seg_intersections *si)
78 {
79  int i;
80 
81  for (i = 0; i < si->ilcount; i++)
82  G_free(si->il[i].a);
83  G_free(si->il);
84  G_free(si->ip);
85  G_free(si);
86 
87  return;
88 }
89 
90 /* internal use */
91 void add_ipoint1(struct seg_intersection_list *il, int with, double dist,
92  int ip)
93 {
94  struct seg_intersection *s;
95 
96  if (il->count == il->allocated) {
97  il->allocated += 4;
98  il->a =
99  G_realloc(il->a, (il->allocated) * sizeof(struct seg_intersection));
100  }
101  s = &(il->a[il->count]);
102  s->with = with;
103  s->ip = ip;
104  s->dist = dist;
105  il->count++;
106 
107  return;
108 }
109 
110 /* adds intersection point to the structure */
111 void add_ipoint(const struct line_pnts *Points, int first_seg, int second_seg,
112  double x, double y, struct seg_intersections *si)
113 {
114  struct intersection_point *t;
115 
116  int ip;
117 
118  G_debug(4, "add_ipoint()");
119 
120  if (si->ipcount == si->ipallocated) {
121  si->ipallocated += 16;
122  si->ip = G_realloc(si->ip, (si->ipallocated) *
123  sizeof(struct intersection_point));
124  }
125  ip = si->ipcount;
126  t = &(si->ip[ip]);
127  t->x = x;
128  t->y = y;
129  t->group = -1;
130  si->ipcount++;
131 
132  add_ipoint1(&(si->il[first_seg]), second_seg,
133  LENGTH((Points->x[first_seg] - x), (Points->y[first_seg] - y)),
134  ip);
135  if (second_seg >= 0)
136  add_ipoint1(
137  &(si->il[second_seg]), first_seg,
138  LENGTH((Points->x[second_seg] - x), (Points->y[second_seg] - y)),
139  ip);
140 }
141 
142 void sort_intersection_list(struct seg_intersection_list *il)
143 {
144  int n, i, j, min;
145  struct seg_intersection t;
146 
147  G_debug(4, "sort_intersection_list()");
148 
149  n = il->count;
150  G_debug(4, " n=%d", n);
151  for (i = 0; i < n - 1; i++) {
152  min = i;
153  for (j = i + 1; j < n; j++) {
154  if (il->a[j].dist < il->a[min].dist) {
155  min = j;
156  }
157  }
158  if (min != i) {
159  t = il->a[i];
160  il->a[i] = il->a[min];
161  il->a[min] = t;
162  }
163  }
164 
165  return;
166 }
167 
168 int compare(const void *a, const void *b)
169 {
170  struct intersection_point *aa, *bb;
171 
172  aa = *((struct intersection_point **)a);
173  bb = *((struct intersection_point **)b);
174 
175  if (aa->x < bb->x)
176  return -1;
177  else if (aa->x > bb->x)
178  return 1;
179  else
180  return (aa->y < bb->y) ? -1 : ((aa->y > bb->y) ? 1 : 0);
181 }
182 
183 /* O(Points->n_points) time */
184 double get_epsilon(struct line_pnts *Points)
185 {
186  int i, np;
187 
188  double min, t;
189 
190  double *x, *y;
191 
192  np = Points->n_points;
193  x = Points->x;
194  y = Points->y;
195 
196  min = MAX(fabs(x[1] - x[0]), fabs(y[1] - y[0]));
197  for (i = 1; i <= np - 2; i++) {
198  t = MAX(fabs(x[i + 1] - x[i]), fabs(y[i + 1] - y[i]));
199  if ((t > 0) && (t < min)) {
200  min = t;
201  }
202  }
203 
204  /* ??? is 0.001 ok ??? */
205  return min * 0.000001;
206 }
207 
208 /* currently O(n*n); future implementation O(nlogn) */
209 struct seg_intersections *find_all_intersections(const struct line_pnts *Points)
210 {
211  int i, j, np;
212  int group, t;
213  int looped;
214 
215  /* double EPSILON = 0.00000001; */
216  double EPSILON = GRASS_EPSILON;
217  double *x, *y;
218  double x1, y1, x2, y2;
219  int res;
220 
221  /*int res2
222  double x1_, y1_, x2_, y2_, z1_, z2_; */
223  struct seg_intersections *si;
224  struct seg_intersection_list *il;
225  struct intersection_point **sorted;
226 
227  G_debug(3, "find_all_intersections()");
228 
229  np = Points->n_points;
230  x = Points->x;
231  y = Points->y;
232 
233  si = create_si_struct(np - 1);
234 
235  looped = ((x[0] == x[np - 1]) && (y[0] == y[np - 1]));
236  G_debug(3, " looped=%d", looped);
237 
238  G_debug(3, " finding intersections...");
239  for (i = 0; i < np - 1; i++) {
240  for (j = i + 1; j < np - 1; j++) {
241  G_debug(4, " checking %d-%d %d-%d", i, i + 1, j, j + 1);
242  /*res = segment_intersection_2d_e(x[i], y[i], x[i+1], y[i+1], x[j],
243  * y[j], x[j+1], y[j+1], &x1, &y1, &x2, &y2); */
244  res = segment_intersection_2d(x[i], y[i], x[i + 1], y[i + 1], x[j],
245  y[j], x[j + 1], y[j + 1], &x1, &y1,
246  &x2, &y2);
247 
248  /* res2 = segment_intersection_2d_e(x[i], y[i], x[i+1],
249  y[i+1], x[j], y[j], x[j+1], y[j+1], &x1_, &y1_, &x2_, &y2_); if
250  ((res != res2) || ((res != 0) && (x1!=x1_ || y1!=y1_)) ) {
251  G_debug(1, "exact=%d orig=%d", res, res2);
252  segment_intersection_2d_test(x[i], y[i], x[i+1], y[i+1], x[j],
253  y[j], x[j+1], y[j+1], &x1, &y1, &x2, &y2);
254  }
255  */
256  G_debug(4, " intersection type = %d", res);
257  if (res == 1) {
258  add_ipoint(Points, i, j, x1, y1, si);
259  }
260  else if ((res >= 2) && (res <= 5)) {
261  add_ipoint(Points, i, j, x1, y1, si);
262  add_ipoint(Points, i, j, x2, y2, si);
263  }
264  }
265  }
266  if (!looped) {
267  /* these are not really intersection points */
268  add_ipoint(Points, 0, -1, Points->x[0], Points->y[0], si);
269  add_ipoint(Points, np - 2, -1, Points->x[np - 1], Points->y[np - 1],
270  si);
271  }
272  G_debug(3, " finding intersections...done");
273 
274  G_debug(3, " postprocessing...");
275  if (si->ipallocated > si->ipcount) {
276  si->ipallocated = si->ipcount;
277  si->ip = G_realloc(si->ip,
278  (si->ipcount) * sizeof(struct intersection_point));
279  }
280  for (i = 0; i < si->ilcount; i++) {
281  il = &(si->il[i]);
282  if (il->allocated > il->count) {
283  il->allocated = il->count;
284  il->a =
285  G_realloc(il->a, (il->count) * sizeof(struct seg_intersection));
286  }
287 
288  if (il->count > 0) {
290  /* is it ok to use qsort here ? */
291  }
292  }
293 
294  /* si->ip will not be reallocated any more so we can use pointers */
295  sorted = G_malloc((si->ipcount) * sizeof(struct intersection_point *));
296  for (i = 0; i < si->ipcount; i++)
297  sorted[i] = &(si->ip[i]);
298 
299  qsort(sorted, si->ipcount, sizeof(struct intersection_point *), compare);
300 
301  /* assign groups */
302  group = 0; /* next available group number */
303  for (i = 0; i < si->ipcount; i++) {
304 
305  t = group;
306  for (j = i - 1; j >= 0; j--) {
307  if (!FEQUAL(sorted[j]->x, sorted[i]->x, EPSILON))
308  /* if (!almost_equal(sorted[j]->x, sorted[i]->x, 16))
309  */
310  break;
311  if (FEQUAL(sorted[j]->y, sorted[i]->y, EPSILON)) {
312  /* if (almost_equal(sorted[j]->y, sorted[i]->y, 16))
313  * { */
314  t = sorted[j]->group;
315  break;
316  }
317  }
318  G_debug(4, " group=%d, ip=%d", t,
319  (int)(sorted[i] - &(si->ip[0])));
320  sorted[i]->group = t;
321  if (t == group)
322  group++;
323  }
324  si->group_count = group;
325 
326  G_debug(3, " postprocessing...done");
327 
328  /* output contents of si */
329  for (i = 0; i < si->ilcount; i++) {
330  G_debug(4, "%d-%d :", i, i + 1);
331  for (j = 0; j < si->il[i].count; j++) {
332  G_debug(4, " %d-%d, group=%d", si->il[i].a[j].with,
333  si->il[i].a[j].with + 1, si->ip[si->il[i].a[j].ip].group);
334  G_debug(4, " dist=%.18f", si->il[i].a[j].dist);
335  G_debug(4, " x=%.18f, y=%.18f",
336  si->ip[si->il[i].a[j].ip].x, si->ip[si->il[i].a[j].ip].y);
337  }
338  }
339  G_free(sorted);
340 
341  return si;
342 }
343 
344 /* create's graph with n vertices and allocates memory for e edges */
345 /* trying to add more than e edges, produces fatal error */
346 struct planar_graph *pg_create_struct(int n, int e)
347 {
348  struct planar_graph *pg;
349 
350  pg = G_malloc(sizeof(struct planar_graph));
351  pg->vcount = n;
352  pg->v = G_malloc(n * sizeof(struct pg_vertex));
353  memset(pg->v, 0, n * sizeof(struct pg_vertex));
354  pg->ecount = 0;
355  pg->eallocated = MAX(e, 0);
356  pg->e = NULL;
357  pg->e = G_malloc(e * sizeof(struct pg_edge));
358 
359  return pg;
360 }
361 
363 {
364  int i;
365 
366  for (i = 0; i < pg->vcount; i++) {
367  G_free(pg->v[i].edges);
368  G_free(pg->v[i].angles);
369  }
370  G_free(pg->v);
371  G_free(pg->e);
372  G_free(pg);
373 }
374 
375 /* v1 and v2 must be valid */
376 int pg_existsedge(struct planar_graph *pg, int v1, int v2)
377 {
378  struct pg_vertex *v;
379  struct pg_edge *e;
380  int i, ecount;
381 
382  if (pg->v[v1].ecount <= pg->v[v2].ecount)
383  v = &(pg->v[v1]);
384  else
385  v = &(pg->v[v2]);
386 
387  ecount = v->ecount;
388  for (i = 0; i < ecount; i++) {
389  e = v->edges[i];
390  if (((e->v1 == v1) && (e->v2 == v2)) ||
391  ((e->v1 == v2) && (e->v2 == v1)))
392  return 1;
393  }
394 
395  return 0;
396 }
397 
398 /* for internal use */
399 void pg_addedge1(struct pg_vertex *v, struct pg_edge *e)
400 {
401  if (v->ecount == v->eallocated) {
402  v->eallocated += 4;
403  v->edges =
404  G_realloc(v->edges, (v->eallocated) * sizeof(struct pg_edge *));
405  }
406  v->edges[v->ecount] = e;
407  v->ecount++;
408 }
409 
410 void pg_addedge(struct planar_graph *pg, int v1, int v2)
411 {
412  struct pg_edge *e;
413 
414  G_debug(4, "pg_addedge(), v1=%d, v2=%d", v1, v2);
415 
416  if ((v1 == v2) || (v1 < 0) || (v1 >= pg->vcount) || (v2 < 0) ||
417  (v2 >= pg->vcount)) {
418  G_fatal_error(" pg_addedge(), v1 and/or v2 is invalid");
419  return;
420  }
421 
422  if (pg_existsedge(pg, v1, v2))
423  return;
424 
425  if (pg->ecount == pg->eallocated) {
426  G_fatal_error(_("Trying to add more edges to the planar_graph "
427  "than the initial allocation size allows"));
428  }
429  e = &(pg->e[pg->ecount]);
430  e->v1 = v1;
431  e->v2 = v2;
432  e->winding_left =
433  0; /* winding is undefined if the corresponding side is not visited */
434  e->winding_right = 0;
435  e->visited_left = 0;
436  e->visited_right = 0;
437  pg->ecount++;
438  pg_addedge1(&(pg->v[v1]), e);
439  pg_addedge1(&(pg->v[v2]), e);
440 
441  return;
442 }
443 
444 struct planar_graph *pg_create(const struct line_pnts *Points)
445 {
446  struct seg_intersections *si;
447  struct planar_graph *pg;
448  struct intersection_point *ip;
449  struct pg_vertex *vert;
450  struct pg_edge *edge;
451 
452  int i, j, t, v;
453 
454  G_debug(3, "pg_create()");
455 
456  si = find_all_intersections(Points);
457  pg = pg_create_struct(si->group_count, 2 * (si->ipcount));
458 
459  /* set vertices info */
460  for (i = 0; i < si->ipcount; i++) {
461  ip = &(si->ip[i]);
462  t = ip->group;
463  pg->v[t].x = ip->x;
464  pg->v[t].y = ip->y;
465  }
466 
467  /* add all edges */
468  for (i = 0; i < si->ilcount; i++) {
469  v = si->ip[si->il[i].a[0].ip].group;
470  for (j = 1; j < si->il[i].count; j++) {
471  t = si->ip[si->il[i].a[j].ip].group;
472  if (t != v) {
473  pg_addedge(pg, v, t); /* edge direction is v ---> t */
474  v = t;
475  }
476  }
477  }
478 
479  /* precalculate angles with 0x */
480  for (i = 0; i < pg->vcount; i++) {
481  vert = &(pg->v[i]);
482  vert->angles = G_malloc((vert->ecount) * sizeof(double));
483  for (j = 0; j < vert->ecount; j++) {
484  edge = vert->edges[j];
485  t = (edge->v1 != i) ? (edge->v1) : (edge->v2);
486  vert->angles[j] = atan2(pg->v[t].y - vert->y, pg->v[t].x - vert->x);
487  }
488  }
489 
490  destroy_si_struct(si);
491  /*
492  I'm not sure if shrinking of the allocated memory always preserves it's
493  physical place. That's why I don't want to do this: if (pg->ecount <
494  pg->eallocated) { pg->eallocated = pg->ecount; pg->e = G_realloc(pg->e,
495  (pg->ecount)*sizeof(struct pg_edge));
496  }
497  */
498 
499  /* very time consuming */
500  /*
501  for (i = 0; i < pg->vcount; i++) {
502  if (pg->v[i].ecount < pg->v[i].eallocated) {
503  pg->v[i].eallocated = pg->v[i].ecount;
504  pg->v[i].edges = G_realloc(pg->v[i].edges,
505  (pg->v[i].ecount)*sizeof(struct pg_edges));
506  }
507  }
508  */
509 
510  /* output pg */
511  for (i = 0; i < pg->vcount; i++) {
512  G_debug(4, " vertex %d (%g, %g)", i, pg->v[i].x, pg->v[i].y);
513  for (j = 0; j < pg->v[i].ecount; j++) {
514  G_debug(4, " edge %d-%d", pg->v[i].edges[j]->v1,
515  pg->v[i].edges[j]->v2);
516  }
517  }
518 
519  return pg;
520 }
#define EPSILON
Definition: blas_level_2.c:26
#define NULL
Definition: ccmath.h:32
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:150
#define G_realloc(p, n)
Definition: defs/gis.h:96
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
#define G_malloc(n)
Definition: defs/gis.h:94
int G_debug(int, const char *,...) __attribute__((format(printf
void add_ipoint1(struct seg_intersection_list *il, int with, double dist, int ip)
Definition: dgraph.c:91
struct planar_graph * pg_create_struct(int n, int e)
Definition: dgraph.c:346
#define LENGTH(DX, DY)
Definition: dgraph.c:25
struct seg_intersections * create_si_struct(int segments_count)
Definition: dgraph.c:56
void add_ipoint(const struct line_pnts *Points, int first_seg, int second_seg, double x, double y, struct seg_intersections *si)
Definition: dgraph.c:111
int pg_existsedge(struct planar_graph *pg, int v1, int v2)
Definition: dgraph.c:376
void sort_intersection_list(struct seg_intersection_list *il)
Definition: dgraph.c:142
void pg_addedge(struct planar_graph *pg, int v1, int v2)
Definition: dgraph.c:410
double get_epsilon(struct line_pnts *Points)
Definition: dgraph.c:184
struct planar_graph * pg_create(const struct line_pnts *Points)
Definition: dgraph.c:444
void destroy_si_struct(struct seg_intersections *si)
Definition: dgraph.c:77
int compare(const void *a, const void *b)
Definition: dgraph.c:168
void pg_destroy_struct(struct planar_graph *pg)
Definition: dgraph.c:362
struct seg_intersections * find_all_intersections(const struct line_pnts *Points)
Definition: dgraph.c:209
void pg_addedge1(struct pg_vertex *v, struct pg_edge *e)
Definition: dgraph.c:399
#define min(x, y)
Definition: draw2.c:29
int segment_intersection_2d(double ax1, double ay1, double ax2, double ay2, double bx1, double by1, double bx2, double by2, double *x1, double *y1, double *x2, double *y2)
Definition: e_intersect.c:692
#define FEQUAL(X, Y, TOL)
Definition: e_intersect.h:5
#define GRASS_EPSILON
Definition: gis.h:173
#define MAX(a, b)
Definition: gis.h:149
#define _(str)
Definition: glocale.h:10
int count
double b
Definition: r_raster.c:39
double t
Definition: r_raster.c:39
Feature geometry info - coordinates.
Definition: dig_structs.h:1651
double * y
Array of Y coordinates.
Definition: dig_structs.h:1659
double * x
Array of X coordinates.
Definition: dig_structs.h:1655
int n_points
Number of points.
Definition: dig_structs.h:1667
Definition: dgraph.h:6
char winding_left
Definition: dgraph.h:11
int v1
Definition: dgraph.h:7
int v2
Definition: dgraph.h:8
char winding_right
Definition: dgraph.h:12
char visited_right
Definition: dgraph.h:10
char visited_left
Definition: dgraph.h:9
int ecount
Definition: dgraph.h:18
struct pg_edge ** edges
Definition: dgraph.h:20
double * angles
Definition: dgraph.h:21
double x
Definition: dgraph.h:16
int eallocated
Definition: dgraph.h:19
double y
Definition: dgraph.h:17
int ecount
Definition: dgraph.h:27
struct pg_edge * e
Definition: dgraph.h:29
int eallocated
Definition: dgraph.h:28
struct pg_vertex * v
Definition: dgraph.h:26
int vcount
Definition: dgraph.h:25
#define x