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