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