GRASS GIS 7 Programmer's Manual  7.5.svn(2018)-r72255
buffer.c
Go to the documentation of this file.
1 /*!
2  \file lib/vector/Vlib/buffer.c
3
4  \brief Vector library - nearest, adjust, parallel lines
5
6  Higher level functions for reading/writing/manipulating vectors.
7
8  See buffer2.c for replacement.
9
10  (C) 2001-2009 by the GRASS Development Team
11
12  This program is free software under the
13  GNU General Public License (>=v2).
14  Read the file COPYING that comes with GRASS
15  for details.
16
18  */
19
20 #include <stdlib.h>
21 #include <math.h>
22 #include <grass/vector.h>
23
24 #define LENGTH(DX, DY) ( sqrt( (DX*DX)+(DY*DY) ) )
25 #define PI M_PI
26
27 /* vector() calculates normalized vector form two points */
28 static void vect(double x1, double y1, double x2, double y2, double *x,
29  double *y)
30 {
31  double dx, dy, l;
32
33  dx = x2 - x1;
34  dy = y2 - y1;
35  l = LENGTH(dx, dy);
36  if (l == 0) {
37  /* assume that dx == dy == 0, which should give (NaN,NaN) */
38  /* without this, very small dx or dy could result in Infinity */
39  dx = dy = 0;
40  }
41  *x = dx / l;
42  *y = dy / l;
43 }
44
45 /* find_cross find first crossing between segments from s1 to s2 and from s3 to s4
46  ** s5 is set to first segment and s6 to second
47  ** neighbours are taken as crossing each other only if overlap
48  ** returns: 1 found
49  ** -1 found overlap
51  */
52 static int find_cross(struct line_pnts *Points, int s1, int s2, int s3,
53  int s4, int *s5, int *s6)
54 {
55  int i, j, ret;
56  double *x, *y;
57
58  G_debug(5,
59  "find_cross(): npoints = %d, s1 = %d, s2 = %d, s3 = %d, s4 = %d",
60  Points->n_points, s1, s2, s3, s4);
61
62  x = Points->x;
63  y = Points->y;
64
65  for (i = s1; i <= s2; i++) {
66  for (j = s3; j <= s4; j++) {
67  if (j == i) {
68  continue;
69  }
70  ret =
71  dig_test_for_intersection(x[i], y[i], x[i + 1], y[i + 1],
72  x[j], y[j], x[j + 1], y[j + 1]);
73  if (ret == 1 && ((i - j) > 1 || (i - j) < -1)) {
74  *s5 = i;
75  *s6 = j;
76  G_debug(5, " intersection: s5 = %d, s6 = %d", *s5, *s6);
77  return 1;
78  }
79  if (ret == -1) {
80  *s5 = i;
81  *s6 = j;
82  G_debug(5, " overlap: s5 = %d, s6 = %d", *s5, *s6);
83  return -1;
84  }
85  }
86  }
87  G_debug(5, " no intersection");
88  return 0;
89 }
90
91 /* point_in_buf - test if point px,py is in d buffer of Points
92  ** returns: 1 in buffer
93  ** 0 not in buffer
94  */
95 static int point_in_buf(struct line_pnts *Points, double px, double py,
96  double d)
97 {
98  int i, np;
99  double sd;
100
101  np = Points->n_points;
102  d *= d;
103  for (i = 0; i < np - 1; i++) {
104  sd = dig_distance2_point_to_line(px, py, 0,
105  Points->x[i], Points->y[i], 0,
106  Points->x[i + 1], Points->y[i + 1],
107  0, 0, NULL, NULL, NULL, NULL, NULL);
108  if (sd <= d) {
109  return 1;
110  }
111  }
112  return 0;
113 }
114
115 /* clean_parallel - clean parallel line created by parallel_line:
116  ** - looking for loops and if loop doesn't contain any other loop
117  ** and centroid of loop is in buffer removes this loop (repeated)
118  ** - optionally removes all end points in buffer
119  * parameters:
120  * Points - parallel line
121  * origPoints - original line
122  * d - offset
123  * rm_end - remove end points in buffer
124  ** note1: on some lines (multiply selfcrossing; lines with end points
125  ** in buffer of line other; some shapes of ends ) may create nosense
126  ** note2: this function is stupid and slow, somebody more clever
127  ** than I am should write paralle_line + clean_parallel
128  ** better; RB March 2000
129  */
130 static void clean_parallel(struct line_pnts *Points,
131  struct line_pnts *origPoints, double d, int rm_end)
132 {
133  int i, j, np, npn, sa, sb;
134  int sa_max = 0;
135  int first = 0, current, last, lcount;
136  double *x, *y, px, py, ix, iy;
137  static struct line_pnts *sPoints = NULL;
138
139  G_debug(4, "clean_parallel(): npoints = %d, d = %f, rm_end = %d",
140  Points->n_points, d, rm_end);
141
142  x = Points->x;
143  y = Points->y;
144  np = Points->n_points;
145
146  if (sPoints == NULL)
147  sPoints = Vect_new_line_struct();
148
149  Vect_reset_line(sPoints);
150
151  npn = 1;
152
153  /* remove loops */
154  while (first < np - 2) {
155  /* find first loop which doesn't contain any other loop */
156  current = first;
157  last = Points->n_points - 2;
158  lcount = 0;
159  while (find_cross
160  (Points, current, last - 1, current + 1, last, &sa,
161  &sb) != 0) {
162  if (lcount == 0) {
163  first = sa;
164  } /* move first forward */
165
166  current = sa + 1;
167  last = sb;
168  lcount++;
169  G_debug(5, " current = %d, last = %d, lcount = %d", current,
170  last, lcount);
171  }
172  if (lcount == 0) {
173  break;
175
176  /* ensure sa is monotonically increasing, so npn doesn't reset low */
177  if (sa > sa_max)
178  sa_max = sa;
179  if (sa < sa_max)
180  break;
181
182  /* remove loop if in buffer */
183  if ((sb - sa) == 1) { /* neighbouring lines overlap */
184  j = sb + 1;
185  npn = sa + 1;
186  }
187  else {
188  Vect_reset_line(sPoints);
189  dig_find_intersection(x[sa], y[sa], x[sa + 1], y[sa + 1], x[sb],
190  y[sb], x[sb + 1], y[sb + 1], &ix, &iy);
191  Vect_append_point(sPoints, ix, iy, 0);
192  for (i = sa + 1; i < sb + 1; i++) { /* create loop polygon */
193  Vect_append_point(sPoints, x[i], y[i], 0);
194  }
195  Vect_find_poly_centroid(sPoints, &px, &py);
196  if (point_in_buf(origPoints, px, py, d)) { /* is loop in buffer ? */
197  npn = sa + 1;
198  x[npn] = ix;
199  y[npn] = iy;
200  j = sb + 1;
201  npn++;
202  if (lcount == 0) {
203  first = sb;
204  }
205  }
206  else { /* loop is not in buffer */
207  first = sb;
208  continue;
209  }
210  }
211
212  for (i = j; i < Points->n_points; i++) { /* move points down */
213  x[npn] = x[i];
214  y[npn] = y[i];
215  npn++;
216  }
217  Points->n_points = npn;
218  }
219
220  if (rm_end) {
221  /* remove points from start in buffer */
222  j = 0;
223  for (i = 0; i < Points->n_points - 1; i++) {
224  px = (x[i] + x[i + 1]) / 2;
225  py = (y[i] + y[i + 1]) / 2;
226  if (point_in_buf(origPoints, x[i], y[i], d * 0.9999)
227  && point_in_buf(origPoints, px, py, d * 0.9999)) {
228  j++;
229  }
230  else {
231  break;
232  }
233  }
234  if (j > 0) {
235  npn = 0;
236  for (i = j; i < Points->n_points; i++) {
237  x[npn] = x[i];
238  y[npn] = y[i];
239  npn++;
240  }
241  Points->n_points = npn;
242  }
243  /* remove points from end in buffer */
244  j = 0;
245  for (i = Points->n_points - 1; i >= 1; i--) {
246  px = (x[i] + x[i - 1]) / 2;
247  py = (y[i] + y[i - 1]) / 2;
248  if (point_in_buf(origPoints, x[i], y[i], d * 0.9999)
249  && point_in_buf(origPoints, px, py, d * 0.9999)) {
250  j++;
251  }
252  else {
253  break;
254  }
255  }
256  if (j > 0) {
257  Points->n_points -= j;
258  }
259  }
260 }
261
262 /* parallel_line - remove duplicate points from input line and
263  * creates new parallel line in 'd' offset distance;
264  * 'tol' is tolerance between arc and polyline;
265  * this function doesn't care about created loops;
266  *
267  * New line is written to existing nPoints structure.
268  */
269 static void parallel_line(struct line_pnts *Points, double d, double tol,
270  struct line_pnts *nPoints)
271 {
272  int i, j, np, na, side;
273  double *x, *y, nx, ny, tx, ty, vx, vy, ux, uy, wx, wy;
274  double atol, atol2, a, av, aw;
275
276  G_debug(4, "parallel_line()");
277
278  Vect_reset_line(nPoints);
279
280  Vect_line_prune(Points);
281  np = Points->n_points;
282  x = Points->x;
283  y = Points->y;
284
285  if (np == 0)
286  return;
287
288  if (np == 1) {
289  Vect_append_point(nPoints, x[0], y[0], 0); /* ? OK, should make circle for points ? */
290  return;
291  }
292
293  if (d == 0) {
294  Vect_copy_xyz_to_pnts(nPoints, x, y, NULL, np);
295  return;
296  }
297
298  side = (int)(d / fabs(d));
299  atol = 2 * acos(1 - tol / fabs(d));
300
301  for (i = 0; i < np - 1; i++) {
302  vect(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);
303  vx = ty * d;
304  vy = -tx * d;
305
306  nx = x[i] + vx;
307  ny = y[i] + vy;
308  Vect_append_point(nPoints, nx, ny, 0);
309
310  nx = x[i + 1] + vx;
311  ny = y[i + 1] + vy;
312  Vect_append_point(nPoints, nx, ny, 0);
313
314  if (i < np - 2) { /* use polyline instead of arc between line segments */
315  vect(x[i + 1], y[i + 1], x[i + 2], y[i + 2], &ux, &uy);
316  wx = uy * d;
317  wy = -ux * d;
318  av = atan2(vy, vx);
319  aw = atan2(wy, wx);
320  a = (aw - av) * side;
321  if (a < 0)
322  a += 2 * PI;
323
324  /* TODO: a <= PI can probably fail because of representation error */
325  if (a <= PI && a > atol) {
326  na = (int)(a / atol);
327  atol2 = a / (na + 1) * side;
328  for (j = 0; j < na; j++) {
329  av += atol2;
330  nx = x[i + 1] + fabs(d) * cos(av);
331  ny = y[i + 1] + fabs(d) * sin(av);
332  Vect_append_point(nPoints, nx, ny, 0);
333  }
334  }
335  }
336  }
337  Vect_line_prune(nPoints);
338 }
339
340 /*!
341  \brief Create parallel line
342
343  This function is replaced by Vect_line_parallel2().
344
345  \param InPoints input line
346  \param distance create parallel line in distance
347  \param tolerance maximum distance between theoretical arc and polygon segments
348  \param rm_end remove end points falling into distance
349  \param[out] OutPoints output line
350
351  \return
352  */
353 void
354 Vect_line_parallel(struct line_pnts *InPoints, double distance,
355  double tolerance, int rm_end, struct line_pnts *OutPoints)
356 {
357  G_debug(4,
358  "Vect_line_parallel(): npoints = %d, distance = %f, tolerance = %f",
359  InPoints->n_points, distance, tolerance);
360
361  parallel_line(InPoints, distance, tolerance, OutPoints);
362
363  clean_parallel(OutPoints, InPoints, distance, rm_end);
364
365  return;
366 }
367
368 /*!
369  \brief Create buffer around the line line.
370
371  This function is replaced by Vect_line_buffer().
372
373  Buffer is closed counter clockwise polygon. Warning: output line
374  may contain loops!
375
376  \param InPoints input line
377  \param distance create buffer in distance
378  \param tolerance maximum distance between theoretical arc and polygon segments
379  \param[out] OutPoints output line
380  */
381 void
382 Vect_line_buffer(const struct line_pnts *InPoints, double distance,
383  double tolerance, struct line_pnts *OutPoints)
384 {
385  double dangle;
386  int side, npoints;
387  static struct line_pnts *Points = NULL;
388  static struct line_pnts *PPoints = NULL;
389
390  distance = fabs(distance);
391
392  dangle = 2 * acos(1 - tolerance / fabs(distance)); /* angle step */
393
394  if (Points == NULL)
395  Points = Vect_new_line_struct();
396
397  if (PPoints == NULL)
398  PPoints = Vect_new_line_struct();
399
400  /* Copy and prune input */
401  Vect_reset_line(Points);
402  Vect_append_points(Points, InPoints, GV_FORWARD);
403  Vect_line_prune(Points);
404
405  Vect_reset_line(OutPoints);
406
407  npoints = Points->n_points;
408  if (npoints <= 0) {
409  return;
410  }
411  else if (npoints == 1) { /* make a circle */
412  double angle, x, y;
413
414  for (angle = 0; angle < 2 * PI; angle += dangle) {
415  x = Points->x[0] + distance * cos(angle);
416  y = Points->y[0] + distance * sin(angle);
417  Vect_append_point(OutPoints, x, y, 0);
418  }
419  /* Close polygon */
420  Vect_append_point(OutPoints, OutPoints->x[0], OutPoints->y[0], 0);
421  }
422  else { /* 2 and more points */
423  for (side = 0; side < 2; side++) {
424  double angle, sangle;
425  double lx1, ly1, lx2, ly2;
426  double x, y, nx, ny, sx, sy, ex, ey;
427
428  /* Parallel on one side */
429  if (side == 0) {
430  Vect_line_parallel(Points, distance, tolerance, 0, PPoints);
431  Vect_append_points(OutPoints, PPoints, GV_FORWARD);
432  }
433  else {
434  Vect_line_parallel(Points, -distance, tolerance, 0, PPoints);
435  Vect_append_points(OutPoints, PPoints, GV_BACKWARD);
436  }
437
438  /* Arc at the end */
439  /* 2 points at theend of original line */
440  if (side == 0) {
441  lx1 = Points->x[npoints - 2];
442  ly1 = Points->y[npoints - 2];
443  lx2 = Points->x[npoints - 1];
444  ly2 = Points->y[npoints - 1];
445  }
446  else {
447  lx1 = Points->x[1];
448  ly1 = Points->y[1];
449  lx2 = Points->x[0];
450  ly2 = Points->y[0];
451  }
452
453  /* normalized vector */
454  vect(lx1, ly1, lx2, ly2, &nx, &ny);
455
456  /* starting point */
457  sangle = atan2(-nx, ny); /* starting angle */
458  sx = lx2 + ny * distance;
459  sy = ly2 - nx * distance;
460
461  /* end point */
462  ex = lx2 - ny * distance;
463  ey = ly2 + nx * distance;
464
465  Vect_append_point(OutPoints, sx, sy, 0);
466
467  /* arc */
468  for (angle = dangle; angle < PI; angle += dangle) {
469  x = lx2 + distance * cos(sangle + angle);
470  y = ly2 + distance * sin(sangle + angle);
471  Vect_append_point(OutPoints, x, y, 0);
472  }
473
474  Vect_append_point(OutPoints, ex, ey, 0);
475  }
476
477  /* Close polygon */
478  Vect_append_point(OutPoints, OutPoints->x[0], OutPoints->y[0], 0);
479  }
480  Vect_line_prune(OutPoints);
481
482  return;
483 }
#define GV_FORWARD
Line direction indicator forward/backward.
Definition: dig_defines.h:178
void Vect_line_buffer(const struct line_pnts *InPoints, double distance, double tolerance, struct line_pnts *OutPoints)
Create buffer around the line line.
Definition: buffer.c:382
#define LENGTH(DX, DY)
Definition: buffer.c:24
int first
Definition: rbtree.h:99
#define GV_BACKWARD
Definition: dig_defines.h:179
struct line_pnts * Vect_new_line_struct()
Creates and initializes a line_pnts structure.
Definition: line.c:45
int n_points
Number of points.
Definition: dig_structs.h:1692
int Vect_copy_xyz_to_pnts(struct line_pnts *Points, const double *x, const double *y, const double *z, int n)
Copy points from array to line_pnts structure.
Definition: line.c:99
#define NULL
Definition: ccmath.h:32
#define PI
Definition: buffer.c:25
#define x
int Vect_append_point(struct line_pnts *Points, double x, double y, double z)
Appends one point to the end of a line.
Definition: line.c:149
int Vect_append_points(struct line_pnts *Points, const struct line_pnts *APoints, int direction)
Appends points to the end of a line.
Definition: line.c:337
double * x
Array of X coordinates.
Definition: dig_structs.h:1680
double l
Definition: r_raster.c:39
Feature geometry info - coordinates.
Definition: dig_structs.h:1675
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
double dig_distance2_point_to_line(double, double, double, double, double, double, double, double, double, int, double *, double *, double *, double *, int *)
int dig_find_intersection(double, double, double, double, double, double, double, double, double *, double *)
Definition: linecros.c:185
int dig_test_for_intersection(double, double, double, double, double, double, double, double)
Definition: linecros.c:58
int Vect_find_poly_centroid(const struct line_pnts *points, double *cent_x, double *cent_y)
Get centroid of polygon.
Definition: Vlib/poly.c:363
double * y
Array of Y coordinates.
Definition: dig_structs.h:1684
int Vect_line_prune(struct line_pnts *Points)
Remove duplicate points, i.e. zero length segments.
Definition: line.c:281
void Vect_line_parallel(struct line_pnts *InPoints, double distance, double tolerance, int rm_end, struct line_pnts *OutPoints)
Create parallel line.
Definition: buffer.c:354
void Vect_reset_line(struct line_pnts *Points)
Reset line.
Definition: line.c:130
int
Reads the categories file for map name in mapset and stores the categories in the pcats structure...