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