GRASS Programmer's Manual  6.5.svn(2014)-r66266
buffer.c
Go to the documentation of this file.
1
21 /* BUGS
22  *
23  * buffering these shapes is unresolved for buffer distances which
24  * should create a closed loop in the inside
25  *
26  * ----------- ------------
27  * | \ /
28  * |
29  * |
30  * | / \
31  * ----------- ------------
32  *
33  *
34  * -----------
35  * | |
36  * |
37  * | |
38  * -----------
39  *
40  *
41  * for certain buffer distances, undefined behaviour for
42  *
43  * ---------
44  * | |
45  * |
46  * --------------
47  * |
48  * --------------
49  * |
50  * | |
51  * ---------
52  *
53  * MM April 2011
54  * */
55
56 #include <stdlib.h>
57 #include <math.h>
58 #include <grass/Vect.h>
59 #include <grass/gis.h>
60
61 #define LENGTH(DX, DY) ( sqrt( (DX*DX)+(DY*DY) ) )
62 #define PI M_PI
63 #define D_MULT 0.99999999 /* distance multiplier for point_in_buf() */
64
65 /* vector() calculates normalized vector from two points */
66 static void vect(double x1, double y1, double x2, double y2, double *x,
67  double *y)
68 {
69  double dx, dy, l;
70
71  dx = x2 - x1;
72  dy = y2 - y1;
73  l = LENGTH(dx, dy);
74  if (l == 0) {
75  /* assume that dx == dy == 0, which should give (NaN,NaN) */
76  /* without this, very small dx or dy could result in Infinity */
77  dx = dy = 0;
78  }
79  *x = dx / l;
80  *y = dy / l;
81 }
82
83 /* find_cross find first crossing between segments from s1 to s2 and from s3 to s4
84  ** s5 is set to first segment and s6 to second
85  ** neighbours are taken as crossing each other only if overlap
86  ** returns: 1 found
87  ** -1 found overlap
89  */
90 static int find_cross(struct line_pnts *Points, int s1, int s2, int s3,
91  int s4, int *s5, int *s6)
92 {
93  int i, j, jstart, np, ret;
94  double *x, *y;
95
96  G_debug(5,
97  "find_cross(): npoints = %d, s1 = %d, s2 = %d, s3 = %d, s4 = %d",
98  Points->n_points, s1, s2, s3, s4);
99
100  x = Points->x;
101  y = Points->y;
102  np = Points->n_points;
103
104  for (i = s1; i <= s2; i++) {
105  jstart = i + 1 < s3 ? s3 : i + 1;
106  for (j = jstart; j <= s4; j++) {
107  if (j == i) {
108  continue;
109  }
110  ret =
111  dig_test_for_intersection(x[i], y[i], x[i + 1], y[i + 1],
112  x[j], y[j], x[j + 1], y[j + 1]);
113  /* j should always be > i */
114  if (ret == 1 && ((i - j) > 1 || (i - j) < -1)) {
115  *s5 = i;
116  *s6 = j;
117  G_debug(5, " intersection: s5 = %d, s6 = %d", *s5, *s6);
118  return 1;
119  }
120  if (ret == -1) {
121  *s5 = i;
122  *s6 = j;
123  G_debug(5, " overlap: s5 = %d, s6 = %d", *s5, *s6);
124  return -1;
125  }
126  }
127  }
128  G_debug(5, "find_cross() -> no intersection");
129  return 0;
130 }
131
132 /* find_cross2 find first crossing between segments from s1 to s2 and from s3 to s4
133  ** proceed from s1 to s2 and from s4 to s3, i.e. forward/backward
134  ** s5 is set to first segment and s6 to second
135  ** neighbours are taken as crossing each other only if overlap
136  ** returns: 1 found
137  ** -1 found overlap
139  */
140 static int find_cross_reverse(struct line_pnts *Points, int s1, int s2, int s3,
141  int s4, int *s5, int *s6)
142 {
143  int i, j, jend, np, ret;
144  double *x, *y;
145
146  G_debug(5,
147  "find_cross_reverse(): npoints = %d, s1 = %d, s2 = %d, s3 = %d, s4 = %d",
148  Points->n_points, s1, s2, s3, s4);
149
150  x = Points->x;
151  y = Points->y;
152  np = Points->n_points;
153
154  for (i = s1; i <= s2; i++) {
155  jend = i + 1 < s3 ? s3 : i + 1;
156  for (j = s4 - 1; j >= jend; j--) {
157  if (j == i) {
158  continue;
159  }
160  ret =
161  dig_test_for_intersection(x[i], y[i], x[i + 1], y[i + 1],
162  x[j], y[j], x[j + 1], y[j + 1]);
163  /* j should always be >= i */
164  if (ret == 1 && ((i - j) > 1 || (i - j) < -1)) {
165  *s5 = i;
166  *s6 = j;
167  G_debug(5, " intersection: s5 = %d, s6 = %d", *s5, *s6);
168  return 1;
169  }
170  if (ret == -1) {
171  *s5 = i;
172  *s6 = j;
173  G_debug(5, " overlap: s5 = %d, s6 = %d", *s5, *s6);
174  return -1;
175  }
176  }
177  }
178  G_debug(5, "find_cross_reverse() -> no intersection");
179  return 0;
180 }
181
182 /* find_cross_from_start find first crossing of segment s1 to s1 + 1 with all
183  ** segments from s1 + 1 to Points->n_points - 2
184  ** ix is set to East crossing and iy to North crossing
185  ** s2 is set to the intersecting segment
186  ** neighbours are taken as crossing each other only if overlap
187  ** returns: 1 found
189  */
190 static int find_cross_from_start(struct line_pnts *Points, int s1, int *s2,
191  double *ix, double *iy)
192 {
193  int i, np, ret, found = 0;
194  double *x, *y, min_dist, new_dist, new_x, new_y, dx, dy;
195
196  G_debug(5, "find_cross_from_start(): npoints = %d, s1 = %d",
197  Points->n_points, s1);
198
199  x = Points->x;
200  y = Points->y;
201  np = Points->n_points;
202
203  min_dist = PORT_DOUBLE_MAX;
204
205  for (i = np - 2; i > s1; i--) {
206  ret =
207  dig_test_for_intersection(x[s1], y[s1], x[s1 + 1], y[s1 + 1],
208  x[i], y[i], x[i + 1], y[i + 1]);
209
210  if (ret) {
211  dig_find_intersection(x[s1], y[s1], x[s1 + 1], y[s1 + 1], x[i],
212  y[i], x[i + 1], y[i + 1], &new_x, &new_y);
213  if ((new_x != x[s1] || new_y != y[s1]) &&
214  (new_x != x[s1 + 1] || new_y != y[s1 + 1])) {
215  found = 1;
216  dx = x[s1] - new_x;
217  dy = y[s1] - new_y;
218  new_dist = LENGTH(dx, dy);
219  if (min_dist > new_dist) {
220  min_dist = new_dist;
221  *ix = new_x;
222  *iy = new_y;
223  *s2 = i;
224  }
225  }
226  }
227  }
228
229  G_debug(5, "find_cross_from_start(): intersection %s",
231  return found;
232 }
233
234 /* find_cross_from_end find first crossing of segment s1 to s1 + 1 with all
235  ** segments from 0 to s1 - 1
236  ** ix is set to East crossing and iy to North crossing
237  ** s2 is set to the intersecting segment
238  ** neighbours are taken as crossing each other only if overlap
239  ** returns: 1 found
241  */
242 static int find_cross_from_end(struct line_pnts *Points, int s1, int *s2,
243  double *ix, double *iy)
244 {
245  int i, np, ret, found = 0;
246  double *x, *y, min_dist, new_dist, new_x, new_y, dx, dy;
247
248  G_debug(5, "find_cross_from_end(): npoints = %d, s1 = %d",
249  Points->n_points, s1);
250
251  x = Points->x;
252  y = Points->y;
253  np = Points->n_points;
254
255  min_dist = PORT_DOUBLE_MAX;
256
257  for (i = 0; i < s1; i++) {
258  ret =
259  dig_test_for_intersection(x[s1], y[s1], x[s1 + 1], y[s1 + 1],
260  x[i], y[i], x[i + 1], y[i + 1]);
261
262  if (ret) {
263  dig_find_intersection(x[s1], y[s1], x[s1 + 1], y[s1 + 1], x[i],
264  y[i], x[i + 1], y[i + 1], &new_x, &new_y);
265  if ((new_x != x[s1] || new_y != y[s1]) &&
266  (new_x != x[s1 + 1] || new_y != y[s1 + 1])) {
267  found = 1;
268  dx = x[s1 + 1] - new_x;
269  dy = y[s1 + 1] - new_y;
270  new_dist = LENGTH(dx, dy);
271  if (min_dist > new_dist) {
272  min_dist = new_dist;
273  *ix = new_x;
274  *iy = new_y;
275  *s2 = i;
276  }
277  }
278  }
279  }
280
281  G_debug(5, "find_cross_from_end(): intersection %s",
283  return found;
284 }
285
286 /* point_in_buf - test if point px,py is in d buffer of Points
287  ** returns: 1 in buffer
288  ** 0 not in buffer
289  */
290 static int point_in_buf(struct line_pnts *Points, double px, double py,
291  double d)
292 {
293  int i, np;
294  double sd;
295
296  np = Points->n_points;
297  /* d must be the squared distance */
298  for (i = 0; i < np - 1; i++) {
299  sd = dig_distance2_point_to_line(px, py, 0,
300  Points->x[i], Points->y[i], 0,
301  Points->x[i + 1], Points->y[i + 1],
302  0, 0, NULL, NULL, NULL, NULL, NULL);
303  if (sd <= d) {
304  return 1;
305  }
306  }
307  return 0;
308 }
309
310 /* clean_parallel - clean parallel line created by parallel_line:
311  ** - looking for loops and if loop doesn't contain any other loop
312  ** and centroid of loop is in buffer removes this loop (repeated)
313  ** - optionally removes all end points in buffer
314  ** - optionally closes output line if input line is looped
315  * parameters:
316  * Points - parallel line
317  * origPoints - original line
318  * d - offset
319  * rm_end - remove end points in buffer
320  ** note1: on some lines (multiple selfcrossing; lines with end points
321  ** in buffer of other line; some shapes of ends ) may create nosense
322  ** note2: this function is stupid and slow, somebody more clever
323  ** than I am should write paralle_line + clean_parallel
324  ** better; RB March 2000
325  ** speed a bit improved, more thorough cleaning MM Feb 2011
326  */
327 static void clean_parallel(struct line_pnts *Points,
328  struct line_pnts *origPoints,
329  double d, double tol, int rm_end, int close_loop)
330 {
331  int i, j, np, npn, sa, sb, ret, side;
332  int sa_max = 0;
333  int first = 0, current, last, lcount;
334  double *x, *y, px, py, ix, iy;
335  static struct line_pnts *sPoints = NULL;
336  double d2 = d * d * D_MULT;
337
338  G_debug(4, "clean_parallel(): npoints = %d, d = %f, rm_end = %d",
339  Points->n_points, d, rm_end);
340
341  Vect_line_prune(Points);
342  x = Points->x;
343  y = Points->y;
344  np = Points->n_points;
345
346  if (np < 3)
347  return;
348
349  /* duplicate consecutive points cause problems for
350  * dig_test_for_intersection() and dig_find_intersection()
351  * -> output line must always be pruned */
352
353  if (rm_end) {
354  int inserted;
355
356  /* remove points from end in buffer */
357  j = inserted = 0;
358  for (i = Points->n_points - 1; i >= 1; i--) {
359  x = Points->x;
360  y = Points->y;
361  if (rm_end < 2) {
362  px = (x[i] + x[i - 1]) / 2;
363  py = (y[i] + y[i - 1]) / 2;
364  }
365  else {
366  /* this is for the last arc created by Vect_line_buffer() */
367  px = x[i - 1];
368  py = y[i - 1];
369  }
370  if (point_in_buf(origPoints, x[i], y[i], d2)) {
371
372  /* check for intersection of this segment */
373  /* must not be the first or last point of this segment */
374  current = i - 1;
375  if (find_cross_from_end(Points, current, &sa, &ix, &iy) != 0) {
376  /* insert intersection point into this segment */
377  /* and adjust i and Points->n_points */
378  int k = 0;
379
380  if ((ix == x[sa] && iy == y[sa]) ||
381  (ix == x[sa + 1] && iy == y[sa + 1])) {
382
383  /* do not insert ix,iy into sa segment */
384  Vect_append_point(Points, Points->x[Points->n_points - 1],
385  Points->y[Points->n_points - 1], 0);
386
387  for (k = Points->n_points - 2; k > current + 1; k--) {
388  Points->x[k] = Points->x[k - 1];
389  Points->y[k] = Points->y[k - 1];
390  }
391  Points->x[current + 1] = ix;
392  Points->y[current + 1] = iy;
393  i += 2;
394  }
395  else {
396  Vect_append_point(Points, Points->x[Points->n_points - 1],
397  Points->y[Points->n_points - 1], 0);
398  Vect_append_point(Points, Points->x[Points->n_points - 1],
399  Points->y[Points->n_points - 1], 0);
400
401  for (k = Points->n_points - 2; k > current + 2; k--) {
402  Points->x[k] = Points->x[k - 2];
403  Points->y[k] = Points->y[k - 2];
404  }
405  Points->x[current + 2] = ix;
406  Points->y[current + 2] = iy;
407  for (k = current + 1; k > sa + 1; k--) {
408  Points->x[k] = Points->x[k - 1];
409  Points->y[k] = Points->y[k - 1];
410  }
411  Points->x[sa + 1] = ix;
412  Points->y[sa + 1] = iy;
413  i += 3;
414  }
415
416  inserted++;
417  }
418  else if (point_in_buf(origPoints, px, py, d2)) {
419  j++;
420  if (inserted)
421  close_loop = 0;
422  }
423  else {
424  break;
425  }
426  }
427  else {
428  break;
429  }
430  }
431  if (j > 0) {
432  Points->n_points -= j;
433  }
434
435  /* remove points from start in buffer */
436  j = inserted = 0;
437  for (i = 0; i < Points->n_points - 1; i++) {
438  x = Points->x;
439  y = Points->y;
440  px = (x[i] + x[i + 1]) / 2;
441  py = (y[i] + y[i + 1]) / 2;
442  if (point_in_buf(origPoints, x[i], y[i], d2)) {
443
444  /* check for intersection of this segment */
445  /* must not be the first or last point of this segment */
446  current = i;
447  if (find_cross_from_start(Points, current, &sa, &ix, &iy) != 0) {
448  /* insert intersection point into this segment */
449  /* and adjust i and Points->n_points */
450  int k = 0;
451
452  if ((ix == x[sa] && iy == y[sa]) ||
453  (ix == x[sa + 1] && iy == y[sa + 1])) {
454
455  /* do not insert ix,iy into sa segment */
456  Vect_append_point(Points, Points->x[Points->n_points - 1],
457  Points->y[Points->n_points - 1], 0);
458
459  for (k = Points->n_points - 2; k > current + 1; k--) {
460  Points->x[k] = Points->x[k - 1];
461  Points->y[k] = Points->y[k - 1];
462  }
463  Points->x[current + 1] = ix;
464  Points->y[current + 1] = iy;
465  i--;
466  }
467  else {
468  Vect_append_point(Points, Points->x[Points->n_points - 1],
469  Points->y[Points->n_points - 1], 0);
470  Vect_append_point(Points, Points->x[Points->n_points - 1],
471  Points->y[Points->n_points - 1], 0);
472
473  for (k = Points->n_points - 2; k > sa + 2; k--) {
474  Points->x[k] = Points->x[k - 2];
475  Points->y[k] = Points->y[k - 2];
476  }
477  Points->x[sa + 2] = ix;
478  Points->y[sa + 2] = iy;
479  for (k = sa + 1; k > current + 1; k--) {
480  Points->x[k] = Points->x[k - 1];
481  Points->y[k] = Points->y[k - 1];
482  }
483  Points->x[current + 1] = ix;
484  Points->y[current + 1] = iy;
485  i--;
486  }
487  inserted++;
488  }
489  else if (point_in_buf(origPoints, px, py, d2)) {
490  j++;
491  if (inserted)
492  close_loop = 0;
493  }
494  else {
495  break;
496  }
497  }
498  else {
499  break;
500  }
501  }
502  x = Points->x;
503  y = Points->y;
504  if (j > 0) {
505  npn = 0;
506  for (i = j; i < Points->n_points; i++) {
507  x[npn] = x[i];
508  y[npn] = y[i];
509  npn++;
510  }
511  Points->n_points = npn;
512  }
513
514  /* test for intersection of end segments */
515  /* need to do this here, otherwise remove loops below will
516  * remove everything but the end dangles */
517  np = Points->n_points;
518  if (np > 3 && (x[0] != x[np - 1] || y[0] != y[np - 1])) {
519  i = 0;
520  j = np - 2;
521  ret =
522  dig_test_for_intersection(x[i], y[i], x[i + 1], y[i + 1],
523  x[j], y[j], x[j + 1], y[j + 1]);
524
525  if (ret == -1) {
526  /* overlap */
527  x[np - 1] = x[0];
528  y[np - 1] = y[0];
529  }
530  else if (ret == 1) {
531  dig_find_intersection(x[i], y[i], x[i + 1], y[i + 1], x[j],
532  y[j], x[j + 1], y[j + 1], &ix, &iy);
533  x[0] = ix;
534  y[0] = iy;
535  x[np - 1] = ix;
536  y[np - 1] = iy;
537  }
538  }
539  Vect_line_prune(Points);
540  }
541
542  if (sPoints == NULL)
543  sPoints = Vect_new_line_struct();
544
545  Vect_reset_line(sPoints);
546
547  np = Points->n_points;
548  npn = 1;
549
550  side = (int)(d / fabs(d));
551
552  /* remove loops */
553  while (np > 3 && first < np - 2) {
554  /* find first loop which doesn't contain any other loop */
555  current = first;
556  G_debug(5, "current: %d", current);
557  last = Points->n_points - 2;
558  lcount = 0;
559  sa = current;
560  while (find_cross
561  (Points, current, last - 1, current + 1, last, &sa, &sb) != 0) {
562
563  if (lcount == 0)
564  first = sa; /* move first forward */
565
566  current = sa + 1;
567  last = sb;
568  lcount++;
569  G_debug(5, " current = %d, last = %d, lcount = %d", current,
570  last, lcount);
571  }
572  if (lcount == 0) {
573  break;
575
576  /* ensure sa is monotonically increasing, so npn doesn't reset low */
577  /* disabled, otherwise nested loops are missed
578  if (sa > sa_max)
579  sa_max = sa;
580  if (sa < sa_max)
581  break;
582  */
583  G_debug(4, "sa: %d", sa);
584
585  /* remove loop if in buffer */
586  if ((sb - sa) == 1) { /* neighbouring lines overlap */
587  j = sb + 1;
588  npn = sa + 1;
589  /* first = sa; */
590  }
591  else {
592  double area_size;
593  int in_buffer = 0;
594
595  Vect_reset_line(sPoints);
596  dig_find_intersection(x[sa], y[sa], x[sa + 1], y[sa + 1], x[sb],
597  y[sb], x[sb + 1], y[sb + 1], &ix, &iy);
598
599  if (ix == x[0] && iy == y[0] && ix == x[np - 1] && iy == y[np - 1])
600  break;
601
602  Vect_append_point(sPoints, ix, iy, 0);
603  for (i = sa + 1; i < sb + 1; i++) { /* create loop polygon */
604  Vect_append_point(sPoints, x[i], y[i], 0);
605  }
606  /* close polygon */
607  Vect_append_point(sPoints, ix, iy, 0);
608  Vect_line_prune(sPoints);
609
610  /* is loop in buffer ? */
611  dig_find_area_poly(sPoints, &area_size);
612  in_buffer = area_size * side < 0;
613  if (!in_buffer && area_size) {
614  /* Vect_find_poly_centroid() may produce coords outside polygon */
615  Vect_find_poly_centroid(sPoints, &px, &py);
616  in_buffer = point_in_buf(origPoints, px, py, d2) != 0;
617  }
618
619  if (in_buffer) {
620  npn = sa + 1;
621  /* do not create zero-length segments */
622  if (ix != x[sa] && iy != y[sa]) {
623  x[npn] = ix;
624  y[npn] = iy;
625  npn++;
626  }
627  if (ix != x[sb + 1] && iy != y[sb + 1]) {
628  j = sb + 1;
629  }
630  else {
631  j = sb + 2;
632  }
633
634  /* lcount == 0 can not happen
635  if (lcount == 0) {
636  first = sa;
637  }
638  */
639  }
640  else { /* loop is not in buffer */
641  first = sb;
642  continue;
643  }
644  }
645
646  for (i = j; i < Points->n_points; i++) { /* move points down */
647  x[npn] = x[i];
648  y[npn] = y[i];
649  npn++;
650  }
651  Points->n_points = np = npn;
652  }
653
654  Vect_line_prune(Points);
655
656  /* looped input line ? */
657  if (origPoints->x[0] == origPoints->x[origPoints->n_points - 1] &&
658  origPoints->y[0] == origPoints->y[origPoints->n_points - 1] &&
659  close_loop) {
660  double tx, ty, vx, vy, ux, uy, wx, wy, a, aw, av;
661  int sa, sb, side;
662  double atol, atol2;
663  double *xorg, *yorg;
664
665  atol = 2 * acos(1 - tol / fabs(d));
666
667  side = (int)(d / fabs(d));
668
669  xorg = origPoints->x;
670  yorg = origPoints->y;
671
672  np = origPoints->n_points;
673
674  /* last segment */
675  vect(xorg[np - 2], yorg[np - 2], xorg[np - 1], yorg[np - 1], &tx, &ty);
676  vx = ty * d;
677  vy = -tx * d;
678  /* first segment */
679  vect(xorg[0], yorg[0], xorg[1], yorg[1], &ux, &uy);
680  wx = uy * d;
681  wy = -ux * d;
682  av = atan2(vy, vx);
683  aw = atan2(wy, wx);
684  a = (aw - av) * side;
685  if (a < 0)
686  a += 2 * PI;
687
688  if (a > PI) {
689  G_debug(5, "search intersection");
690  /* there can be two intersections !!! */
691  /* a > PI is probably handled by remove ends above */
692  if (find_cross_reverse(Points, 0, Points->n_points - 2, 1, Points->n_points - 1, &sa,
693  &sb) != 0) {
694  G_debug(5, "found intersection");
695  dig_find_intersection(x[sa], y[sa], x[sa + 1], y[sa + 1], x[sb],
696  y[sb], x[sb + 1], y[sb + 1], &ix, &iy);
697  x[0] = ix;
698  y[0] = iy;
699  npn = 1;
700  /* move points down */
701  for (i = sa + 1; i <= sb; i++) {
702  x[npn] = x[i];
703  y[npn] = y[i];
704  npn++;
705  }
706  Points->n_points = npn;
707  }
708  }
709  /* TODO: a <= PI can probably fail because of representation error */
710  else {
711  if (a <= PI && a > atol) {
712  int na;
713  double nx, ny;
714
715  /* OK to close parallel line ? */
716  na = (int)(a / atol);
717  atol2 = a / (na + 1) * side;
718  for (j = 0; j < na; j++) {
719  av += atol2;
720  nx = xorg[0] + fabs(d) * cos(av);
721  ny = yorg[0] + fabs(d) * sin(av);
722  Vect_append_point(Points, nx, ny, 0);
723  }
724  }
725  }
726  /* always close for looped input line */
727  Vect_append_point(Points, Points->x[0], Points->y[0], 0);
728
729  Vect_line_prune(Points);
730  }
731
732 }
733
734 /* parallel_line - remove duplicate points from input line and
735  * creates new parallel line in 'd' offset distance;
736  * 'tol' is tolerance between arc and polyline;
737  * this function doesn't care about created loops;
738  *
739  * New line is written to existing nPoints structure.
740  */
741 static void parallel_line(struct line_pnts *Points, double d, double tol,
742  struct line_pnts *nPoints)
743 {
744  int i, j, np, na, side;
745  double *x, *y, nx, ny, tx, ty, vx, vy, ux, uy, wx, wy;
746  double atol, atol2, a, av, aw;
747
748  G_debug(4, "parallel_line()");
749
750  Vect_reset_line(nPoints);
751
752  np = Points->n_points;
753  x = Points->x;
754  y = Points->y;
755
756  if (np == 0)
757  return;
758
759  if (np == 1) {
760  Vect_append_point(nPoints, x[0], y[0], 0); /* ? OK, should make circle for points ? */
761  return;
762  }
763
764  if (d == 0) {
765  Vect_copy_xyz_to_pnts(nPoints, x, y, NULL, np);
766  return;
767  }
768
769  side = (int)(d / fabs(d));
770  atol = 2 * acos(1 - tol / fabs(d));
771
772  for (i = 0; i < np - 1; i++) {
773  vect(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);
774  vx = ty * d;
775  vy = -tx * d;
776
777  nx = x[i] + vx;
778  ny = y[i] + vy;
779  Vect_append_point(nPoints, nx, ny, 0);
780
781  nx = x[i + 1] + vx;
782  ny = y[i + 1] + vy;
783  Vect_append_point(nPoints, nx, ny, 0);
784
785  if (i < np - 2) { /* use polyline instead of arc between line segments */
786  vect(x[i + 1], y[i + 1], x[i + 2], y[i + 2], &ux, &uy);
787  wx = uy * d;
788  wy = -ux * d;
789  av = atan2(vy, vx);
790  aw = atan2(wy, wx);
791  a = (aw - av) * side;
792  if (a < 0)
793  a += 2 * PI;
794
795  /* TODO: a <= PI can probably fail because of representation error */
796  if (a <= PI && a > atol) {
797  na = (int)(a / atol);
798  atol2 = a / (na + 1) * side;
799  for (j = 0; j < na; j++) {
800  av += atol2;
801  nx = x[i + 1] + fabs(d) * cos(av);
802  ny = y[i + 1] + fabs(d) * sin(av);
803  Vect_append_point(nPoints, nx, ny, 0);
804  }
805  }
806  }
807  }
808  Vect_line_prune(nPoints);
809 }
810
822 void
823 Vect_line_parallel(struct line_pnts *InPoints, double distance,
824  double tolerance, int rm_end, struct line_pnts *OutPoints)
825 {
826  G_debug(4,
827  "Vect_line_parallel(): npoints = %d, distance = %f, tolerance = %f",
828  InPoints->n_points, distance, tolerance);
829
830  parallel_line(InPoints, distance, tolerance, OutPoints);
831  G_debug(4, "%d outpoints", OutPoints->n_points);
832
833  clean_parallel(OutPoints, InPoints, distance, tolerance, rm_end, 1);
834  G_debug(4, "%d outpoints after cleaning", OutPoints->n_points);
835
836  return;
837 }
838
850 void
851 Vect_line_buffer(struct line_pnts *InPoints, double distance,
852  double tolerance, struct line_pnts *OutPoints)
853 {
854  double dangle;
855  int side, npoints;
856  static struct line_pnts *Points = NULL;
857  static struct line_pnts *PPoints = NULL;
858  double d2 = distance * distance * D_MULT;
859
860  distance = fabs(distance);
861
862  dangle = 2 * acos(1 - tolerance / fabs(distance)); /* angle step */
863
864  if (Points == NULL)
865  Points = Vect_new_line_struct();
866
867  if (PPoints == NULL)
868  PPoints = Vect_new_line_struct();
869
870  /* Copy input */
871  Vect_reset_line(Points);
872  Vect_append_points(Points, InPoints, GV_FORWARD);
873
874  Vect_reset_line(OutPoints);
875
876  npoints = Points->n_points;
877  if (npoints <= 0) {
878  return;
879  }
880  else if (npoints == 1) { /* make a circle */
881  double angle, x, y;
882
883  for (angle = 0; angle < 2 * PI; angle += dangle) {
884  x = Points->x[0] + distance * cos(angle);
885  y = Points->y[0] + distance * sin(angle);
886  Vect_append_point(OutPoints, x, y, 0);
887  }
888  /* Close polygon */
889  Vect_append_point(OutPoints, OutPoints->x[0], OutPoints->y[0], 0);
890  }
891  else { /* 2 and more points */
892  for (side = 0; side < 2; side++) {
893  double angle, sangle;
894  double lx1, ly1, lx2, ly2;
895  double x, y, nx, ny, sx, sy, ex, ey;
896
897  /* Parallel on one side */
898  if (side == 0) {
899  /* Vect_line_parallel(Points, distance, tolerance, 0, PPoints); */
900
901  /* just create parallel line, clean later */
902  parallel_line(InPoints, distance, tolerance, PPoints);
903  G_debug(4, "%d outpoints", PPoints->n_points);
904
905  Vect_append_points(OutPoints, PPoints, GV_FORWARD);
906  }
907  else {
908  /* Vect_line_parallel(Points, -distance, tolerance, 0, PPoints); */
909
910  parallel_line(InPoints, -distance, tolerance, PPoints);
911  G_debug(4, "%d outpoints", PPoints->n_points);
912
913  Vect_append_points(OutPoints, PPoints, GV_BACKWARD);
914  }
915
916  /* Arc at the end */
917  /* 2 points at the end of original line */
918  if (side == 0) {
919  lx1 = Points->x[npoints - 2];
920  ly1 = Points->y[npoints - 2];
921  lx2 = Points->x[npoints - 1];
922  ly2 = Points->y[npoints - 1];
923  }
924  else {
925  lx1 = Points->x[1];
926  ly1 = Points->y[1];
927  lx2 = Points->x[0];
928  ly2 = Points->y[0];
929  }
930
931  /* normalized vector */
932  vect(lx1, ly1, lx2, ly2, &nx, &ny);
933
934  /* starting point */
935  sangle = atan2(-nx, ny); /* starting angle */
936  sx = lx2 + ny * distance;
937  sy = ly2 - nx * distance;
938
939  /* end point */
940  ex = lx2 - ny * distance;
941  ey = ly2 + nx * distance;
942
943  Vect_append_point(OutPoints, sx, sy, 0);
944
945  /* arc */
946  for (angle = dangle; angle < PI; angle += dangle) {
947  x = lx2 + distance * cos(sangle + angle);
948  y = ly2 + distance * sin(sangle + angle);
949
950  Vect_append_point(OutPoints, x, y, 0);
951  }
952
953  Vect_append_point(OutPoints, ex, ey, 0);
954  }
955
956  /* clean up arcs at the end of input lines */
957  Vect_line_prune(OutPoints);
958  if (OutPoints->x[0] == OutPoints->x[OutPoints->n_points - 1] &&
959  OutPoints->y[0] == OutPoints->y[OutPoints->n_points - 1])
960  OutPoints->n_points--;
961
962  clean_parallel(OutPoints, InPoints, distance, tolerance, 2, 0);
963
964  /* Close polygon */
965  if (OutPoints->x[0] != OutPoints->x[OutPoints->n_points - 1] ||
966  OutPoints->y[0] != OutPoints->y[OutPoints->n_points - 1]) {
967
968  if (point_in_buf(InPoints, OutPoints->x[0], OutPoints->y[0], d2)) {
969  OutPoints->x[0] = OutPoints->x[OutPoints->n_points - 1];
970  OutPoints->y[0] = OutPoints->y[OutPoints->n_points - 1];
971  }
972  else if (point_in_buf(InPoints, OutPoints->x[OutPoints->n_points - 1],
973  OutPoints->y[OutPoints->n_points - 1], d2)) {
974  OutPoints->x[OutPoints->n_points - 1] = OutPoints->x[0];
975  OutPoints->y[OutPoints->n_points - 1] = OutPoints->y[0];
976  }
977  else
978  Vect_append_point(OutPoints, OutPoints->x[0], OutPoints->y[0], 0);
979  }
980  }
981
982  return;
983 }
int l
#define LENGTH(DX, DY)
Definition: buffer.c:61
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:62
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:851
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
Definition: g3dcolor.c:48
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
#define D_MULT
Definition: buffer.c:63
return NULL
Definition: dbfopen.c:1394
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
int dig_find_area_poly(struct line_pnts *Points, double *totalarea)
Definition: diglib/poly.c:93
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:823