GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
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 
17  \author Radim Blazek
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
50  ** 0 not found
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;
174  } /* loop not found */
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 }
int Vect_append_points(struct line_pnts *, const struct line_pnts *, int)
Appends points to the end of a line.
Definition: line.c:337
#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
int n_points
Number of points.
Definition: dig_structs.h:1692
int Vect_copy_xyz_to_pnts(struct line_pnts *, const double *, const double *, const double *, int)
Copy points from array to line_pnts structure.
Definition: line.c:99
int Vect_line_prune(struct line_pnts *)
Remove duplicate points, i.e. zero length segments.
Definition: line.c:281
#define NULL
Definition: ccmath.h:32
#define PI
Definition: buffer.c:25
#define x
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
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
double dig_distance2_point_to_line(double, double, double, double, double, double, double, double, double, int, double *, double *, double *, double *, int *)
int Vect_append_point(struct line_pnts *, double, double, double)
Appends one point to the end of a line.
Definition: line.c:149
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
double * y
Array of Y coordinates.
Definition: dig_structs.h:1684
int Vect_find_poly_centroid(const struct line_pnts *, double *, double *)
Get centroid of polygon.
Definition: Vlib/poly.c:363
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
int G_debug(int, const char *,...) __attribute__((format(printf
void Vect_reset_line(struct line_pnts *)
Reset line.
Definition: line.c:130