GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-8828094037
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
46  *s4
47  ** s5 is set to first segment and s6 to second
48  ** neighbours are taken as crossing each other only if overlap
49  ** returns: 1 found
50  ** -1 found overlap
51  ** 0 not found
52  */
53 static int find_cross(struct line_pnts *Points, int s1, int s2, int s3, int s4,
54  int *s5, int *s6)
55 {
56  int i, j, ret;
57  double *x, *y;
58 
59  G_debug(5, "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 = dig_test_for_intersection(x[i], y[i], x[i + 1], y[i + 1],
71  x[j], y[j], x[j + 1], y[j + 1]);
72  if (ret == 1 && ((i - j) > 1 || (i - j) < -1)) {
73  *s5 = i;
74  *s6 = j;
75  G_debug(5, " intersection: s5 = %d, s6 = %d", *s5, *s6);
76  return 1;
77  }
78  if (ret == -1) {
79  *s5 = i;
80  *s6 = j;
81  G_debug(5, " overlap: s5 = %d, s6 = %d", *s5, *s6);
82  return -1;
83  }
84  }
85  }
86  G_debug(5, " no intersection");
87  return 0;
88 }
89 
90 /* point_in_buf - test if point px,py is in d buffer of Points
91  ** returns: 1 in buffer
92  ** 0 not in buffer
93  */
94 static int point_in_buf(struct line_pnts *Points, double px, double py,
95  double d)
96 {
97  int i, np;
98  double sd;
99 
100  np = Points->n_points;
101  d *= d;
102  for (i = 0; i < np - 1; i++) {
103  sd = dig_distance2_point_to_line(px, py, 0, Points->x[i], Points->y[i],
104  0, Points->x[i + 1], Points->y[i + 1],
105  0, 0, NULL, NULL, NULL, NULL, NULL);
106  if (sd <= d) {
107  return 1;
108  }
109  }
110  return 0;
111 }
112 
113 /* clean_parallel - clean parallel line created by parallel_line:
114  ** - looking for loops and if loop doesn't contain any other loop
115  ** and centroid of loop is in buffer removes this loop (repeated)
116  ** - optionally removes all end points in buffer
117  * parameters:
118  * Points - parallel line
119  * origPoints - original line
120  * d - offset
121  * rm_end - remove end points in buffer
122  ** note1: on some lines (multiply selfcrossing; lines with end points
123  ** in buffer of line other; some shapes of ends ) may create nosense
124  ** note2: this function is stupid and slow, somebody more clever
125  ** than I am should write paralle_line + clean_parallel
126  ** better; RB March 2000
127  */
128 static void clean_parallel(struct line_pnts *Points,
129  struct line_pnts *origPoints, double d, int rm_end)
130 {
131  int i, j, np, npn, sa, sb;
132  int sa_max = 0;
133  int first = 0, current, last, lcount;
134  double *x, *y, px, py, ix, iy;
135  static struct line_pnts *sPoints = NULL;
136 
137  G_debug(4, "clean_parallel(): npoints = %d, d = %f, rm_end = %d",
138  Points->n_points, d, rm_end);
139 
140  x = Points->x;
141  y = Points->y;
142  np = Points->n_points;
143 
144  if (sPoints == NULL)
145  sPoints = Vect_new_line_struct();
146 
147  Vect_reset_line(sPoints);
148 
149  npn = 1;
150 
151  /* remove loops */
152  while (first < np - 2) {
153  /* find first loop which doesn't contain any other loop */
154  current = first;
155  last = Points->n_points - 2;
156  lcount = 0;
157  while (find_cross(Points, current, last - 1, current + 1, last, &sa,
158  &sb) != 0) {
159  if (lcount == 0) {
160  first = sa;
161  } /* move first forward */
162 
163  current = sa + 1;
164  last = sb;
165  lcount++;
166  G_debug(5, " current = %d, last = %d, lcount = %d", current, last,
167  lcount);
168  }
169  if (lcount == 0) {
170  break;
171  } /* loop not found */
172 
173  /* ensure sa is monotonically increasing, so npn doesn't reset low */
174  if (sa > sa_max)
175  sa_max = sa;
176  if (sa < sa_max)
177  break;
178 
179  /* remove loop if in buffer */
180  if ((sb - sa) == 1) { /* neighbouring lines overlap */
181  j = sb + 1;
182  npn = sa + 1;
183  }
184  else {
185  Vect_reset_line(sPoints);
186  dig_find_intersection(x[sa], y[sa], x[sa + 1], y[sa + 1], x[sb],
187  y[sb], x[sb + 1], y[sb + 1], &ix, &iy);
188  Vect_append_point(sPoints, ix, iy, 0);
189  for (i = sa + 1; i < sb + 1; i++) { /* create loop polygon */
190  Vect_append_point(sPoints, x[i], y[i], 0);
191  }
192  Vect_find_poly_centroid(sPoints, &px, &py);
193  if (point_in_buf(origPoints, px, py, d)) { /* is loop in buffer ? */
194  npn = sa + 1;
195  x[npn] = ix;
196  y[npn] = iy;
197  j = sb + 1;
198  npn++;
199  if (lcount == 0) {
200  first = sb;
201  }
202  }
203  else { /* loop is not in buffer */
204  first = sb;
205  continue;
206  }
207  }
208 
209  for (i = j; i < Points->n_points; i++) { /* move points down */
210  x[npn] = x[i];
211  y[npn] = y[i];
212  npn++;
213  }
214  Points->n_points = npn;
215  }
216 
217  if (rm_end) {
218  /* remove points from start in buffer */
219  j = 0;
220  for (i = 0; i < Points->n_points - 1; i++) {
221  px = (x[i] + x[i + 1]) / 2;
222  py = (y[i] + y[i + 1]) / 2;
223  if (point_in_buf(origPoints, x[i], y[i], d * 0.9999) &&
224  point_in_buf(origPoints, px, py, d * 0.9999)) {
225  j++;
226  }
227  else {
228  break;
229  }
230  }
231  if (j > 0) {
232  npn = 0;
233  for (i = j; i < Points->n_points; i++) {
234  x[npn] = x[i];
235  y[npn] = y[i];
236  npn++;
237  }
238  Points->n_points = npn;
239  }
240  /* remove points from end in buffer */
241  j = 0;
242  for (i = Points->n_points - 1; i >= 1; i--) {
243  px = (x[i] + x[i - 1]) / 2;
244  py = (y[i] + y[i - 1]) / 2;
245  if (point_in_buf(origPoints, x[i], y[i], d * 0.9999) &&
246  point_in_buf(origPoints, px, py, d * 0.9999)) {
247  j++;
248  }
249  else {
250  break;
251  }
252  }
253  if (j > 0) {
254  Points->n_points -= j;
255  }
256  }
257 }
258 
259 /* parallel_line - remove duplicate points from input line and
260  * creates new parallel line in 'd' offset distance;
261  * 'tol' is tolerance between arc and polyline;
262  * this function doesn't care about created loops;
263  *
264  * New line is written to existing nPoints structure.
265  */
266 static void parallel_line(struct line_pnts *Points, double d, double tol,
267  struct line_pnts *nPoints)
268 {
269  int i, j, np, na, side;
270  double *x, *y, nx, ny, tx, ty, vx, vy, ux, uy, wx, wy;
271  double atol, atol2, a, av, aw;
272 
273  G_debug(4, "parallel_line()");
274 
275  Vect_reset_line(nPoints);
276 
277  Vect_line_prune(Points);
278  np = Points->n_points;
279  x = Points->x;
280  y = Points->y;
281 
282  if (np == 0)
283  return;
284 
285  if (np == 1) {
286  Vect_append_point(nPoints, x[0], y[0],
287  0); /* ? OK, should make circle for points ? */
288  return;
289  }
290 
291  if (d == 0) {
292  Vect_copy_xyz_to_pnts(nPoints, x, y, NULL, np);
293  return;
294  }
295 
296  side = (int)(d / fabs(d));
297  atol = 2 * acos(1 - tol / fabs(d));
298 
299  for (i = 0; i < np - 1; i++) {
300  vect(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);
301  vx = ty * d;
302  vy = -tx * d;
303 
304  nx = x[i] + vx;
305  ny = y[i] + vy;
306  Vect_append_point(nPoints, nx, ny, 0);
307 
308  nx = x[i + 1] + vx;
309  ny = y[i + 1] + vy;
310  Vect_append_point(nPoints, nx, ny, 0);
311 
312  if (i <
313  np - 2) { /* use polyline instead of arc between line segments */
314  vect(x[i + 1], y[i + 1], x[i + 2], y[i + 2], &ux, &uy);
315  wx = uy * d;
316  wy = -ux * d;
317  av = atan2(vy, vx);
318  aw = atan2(wy, wx);
319  a = (aw - av) * side;
320  if (a < 0)
321  a += 2 * PI;
322 
323  /* TODO: a <= PI can probably fail because of representation error
324  */
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
348  segments \param rm_end remove end points falling into distance \param[out]
349  OutPoints output line
350 
351  \return
352  */
353 void Vect_line_parallel(struct line_pnts *InPoints, double distance,
354  double tolerance, int rm_end,
355  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
379  segments \param[out] OutPoints output line
380  */
381 void Vect_line_buffer(const struct line_pnts *InPoints, double distance,
382  double tolerance, struct line_pnts *OutPoints)
383 {
384  double dangle;
385  int side, npoints;
386  static struct line_pnts *Points = NULL;
387  static struct line_pnts *PPoints = NULL;
388 
389  distance = fabs(distance);
390 
391  dangle = 2 * acos(1 - tolerance / fabs(distance)); /* angle step */
392 
393  if (Points == NULL)
394  Points = Vect_new_line_struct();
395 
396  if (PPoints == NULL)
397  PPoints = Vect_new_line_struct();
398 
399  /* Copy and prune input */
400  Vect_reset_line(Points);
401  Vect_append_points(Points, InPoints, GV_FORWARD);
402  Vect_line_prune(Points);
403 
404  Vect_reset_line(OutPoints);
405 
406  npoints = Points->n_points;
407  if (npoints <= 0) {
408  return;
409  }
410  else if (npoints == 1) { /* make a circle */
411  double angle, x, y;
412 
413  for (angle = 0; angle < 2 * PI; angle += dangle) {
414  x = Points->x[0] + distance * cos(angle);
415  y = Points->y[0] + distance * sin(angle);
416  Vect_append_point(OutPoints, x, y, 0);
417  }
418  /* Close polygon */
419  Vect_append_point(OutPoints, OutPoints->x[0], OutPoints->y[0], 0);
420  }
421  else { /* 2 and more points */
422  for (side = 0; side < 2; side++) {
423  double angle, sangle;
424  double lx1, ly1, lx2, ly2;
425  double x, y, nx, ny, sx, sy, ex, ey;
426 
427  /* Parallel on one side */
428  if (side == 0) {
429  Vect_line_parallel(Points, distance, tolerance, 0, PPoints);
430  Vect_append_points(OutPoints, PPoints, GV_FORWARD);
431  }
432  else {
433  Vect_line_parallel(Points, -distance, tolerance, 0, PPoints);
434  Vect_append_points(OutPoints, PPoints, GV_BACKWARD);
435  }
436 
437  /* Arc at the end */
438  /* 2 points at theend of original line */
439  if (side == 0) {
440  lx1 = Points->x[npoints - 2];
441  ly1 = Points->y[npoints - 2];
442  lx2 = Points->x[npoints - 1];
443  ly2 = Points->y[npoints - 1];
444  }
445  else {
446  lx1 = Points->x[1];
447  ly1 = Points->y[1];
448  lx2 = Points->x[0];
449  ly2 = Points->y[0];
450  }
451 
452  /* normalized vector */
453  vect(lx1, ly1, lx2, ly2, &nx, &ny);
454 
455  /* starting point */
456  sangle = atan2(-nx, ny); /* starting angle */
457  sx = lx2 + ny * distance;
458  sy = ly2 - nx * distance;
459 
460  /* end point */
461  ex = lx2 - ny * distance;
462  ey = ly2 + nx * distance;
463 
464  Vect_append_point(OutPoints, sx, sy, 0);
465 
466  /* arc */
467  for (angle = dangle; angle < PI; angle += dangle) {
468  x = lx2 + distance * cos(sangle + angle);
469  y = ly2 + distance * sin(sangle + angle);
470  Vect_append_point(OutPoints, x, y, 0);
471  }
472 
473  Vect_append_point(OutPoints, ex, ey, 0);
474  }
475 
476  /* Close polygon */
477  Vect_append_point(OutPoints, OutPoints->x[0], OutPoints->y[0], 0);
478  }
479  Vect_line_prune(OutPoints);
480 
481  return;
482 }
#define LENGTH(DX, DY)
Definition: buffer.c:24
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:353
#define PI
Definition: buffer.c:25
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:381
#define NULL
Definition: ccmath.h:32
int G_debug(int, const char *,...) __attribute__((format(printf
int Vect_find_poly_centroid(const struct line_pnts *, double *, double *)
Get centroid of polygon.
Definition: Vlib/poly.c:355
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
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
void Vect_reset_line(struct line_pnts *)
Reset line.
Definition: line.c:129
int Vect_line_prune(struct line_pnts *)
Remove duplicate points, i.e. zero length segments.
Definition: line.c:279
int Vect_append_point(struct line_pnts *, double, double, double)
Appends one point to the end of a line.
Definition: line.c:148
int Vect_append_points(struct line_pnts *, const struct line_pnts *, int)
Appends points to the end of a line.
Definition: line.c:335
#define GV_FORWARD
Line direction indicator forward/backward.
Definition: dig_defines.h:179
#define GV_BACKWARD
Definition: dig_defines.h:180
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:181
int dig_test_for_intersection(double, double, double, double, double, double, double, double)
Definition: linecros.c:57
double l
Definition: r_raster.c:39
Feature geometry info - coordinates.
Definition: dig_structs.h:1651
double * y
Array of Y coordinates.
Definition: dig_structs.h:1659
double * x
Array of X coordinates.
Definition: dig_structs.h:1655
int n_points
Number of points.
Definition: dig_structs.h:1667
#define x