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