GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
buffer2.c
Go to the documentation of this file.
1 /*!
2  \file lib/vector/Vlib/buffer2.c
3
4  \brief Vector library - nearest, adjust, parallel lines
5
6  Higher level functions for reading/writing/manipulating vectors.
7
8  (C) 2001-2009 by the GRASS Development Team
9
10  This program is free software under the
11  GNU General Public License (>=v2).
12  Read the file COPYING that comes with GRASS
13  for details.
14
15  \author Original author Radim Blazek (see buffer.c)
16  \author Rewritten by Rosen Matev (Google Summer of Code 2008)
17  */
18
19 #include <stdlib.h>
20 #include <math.h>
21 #include <grass/vector.h>
22 #include <grass/glocale.h>
23
24 #include "dgraph.h"
25
26 #define LENGTH(DX, DY) (sqrt((DX*DX)+(DY*DY)))
27 #ifndef MIN
28 #define MIN(X,Y) ((X<Y)?X:Y)
29 #endif
30 #ifndef MAX
31 #define MAX(X,Y) ((X>Y)?X:Y)
32 #endif
33 #define PI M_PI
34 #define RIGHT_SIDE 1
35 #define LEFT_SIDE -1
36 #define LOOPED_LINE 1
37 #define NON_LOOPED_LINE 0
38
39 /* norm_vector() calculates normalized vector form two points */
40 static void norm_vector(double x1, double y1, double x2, double y2, double *x,
41  double *y)
42 {
43  double dx, dy, l;
44
45  dx = x2 - x1;
46  dy = y2 - y1;
47  if ((dx == 0) && (dy == 0)) {
48  /* assume that dx == dy == 0, which should give (NaN,NaN) */
49  /* without this, very small dx or dy could result in Infinity */
50  *x = 0;
51  *y = 0;
52  return;
53  }
54  l = LENGTH(dx, dy);
55  *x = dx / l;
56  *y = dy / l;
57
58  return;
59 }
60
61 static void rotate_vector(double x, double y, double cosa, double sina,
62  double *nx, double *ny)
63 {
64  *nx = x * cosa - y * sina;
65  *ny = x * sina + y * cosa;
66
67  return;
68 }
69
70 /*
71  * (x,y) should be normalized vector for common transforms; This func transforms (x,y) to a vector corresponding to da, db, dalpha params
72  * dalpha is in radians
73  */
74 static void elliptic_transform(double x, double y, double da, double db,
75  double dalpha, double *nx, double *ny)
76 {
77  double cosa = cos(dalpha);
78  double sina = sin(dalpha);
79
80  /* double cc = cosa*cosa;
81  double ss = sina*sina;
82  double t = (da-db)*sina*cosa;
83
84  *nx = (da*cc + db*ss)*x + t*y;
85  *ny = (da*ss + db*cc)*y + t*x;
86  return; */
87
88  double va, vb;
89
90  va = (x * cosa + y * sina) * da;
91  vb = (x * (-sina) + y * cosa) * db;
92  *nx = va * cosa + vb * (-sina);
93  *ny = va * sina + vb * cosa;
94
95  return;
96 }
97
98 /*
99  * vect(x,y) must be normalized
100  * gives the tangent point of the tangent to ellpise(da,db,dalpha) parallel to vect(x,y)
101  * dalpha is in radians
102  * ellipse center is in (0,0)
103  */
104 static void elliptic_tangent(double x, double y, double da, double db,
105  double dalpha, double *px, double *py)
106 {
107  double cosa = cos(dalpha);
108  double sina = sin(dalpha);
109  double u, v, len;
110
111  /* rotate (x,y) -dalpha radians */
112  rotate_vector(x, y, cosa, -sina, &x, &y);
113  /*u = (x + da*y/db)/2;
114  v = (y - db*x/da)/2; */
115  u = da * da * y;
116  v = -db * db * x;
117  len = da * db / sqrt(da * da * v * v + db * db * u * u);
118  u *= len;
119  v *= len;
120  rotate_vector(u, v, cosa, sina, px, py);
121
122  return;
123 }
124
125
126 /*
127  * !!! This is not line in GRASS' sense. See http://en.wikipedia.org/wiki/Line_%28mathematics%29
128  */
129 static void line_coefficients(double x1, double y1, double x2, double y2,
130  double *a, double *b, double *c)
131 {
132  *a = y2 - y1;
133  *b = x1 - x2;
134  *c = x2 * y1 - x1 * y2;
135
136  return;
137 }
138
139 /*
140  * Finds intersection of two straight lines. Returns 0 if the lines are parallel, 1 if they cross,
141  * 2 if they are the same line.
142  * !!!!!!!!!!!!!!!! FIX THIS TOLLERANCE CONSTANTS BAD (and UGLY) CODE !!!!!!!!!
143  */
144 static int line_intersection(double a1, double b1, double c1, double a2,
145  double b2, double c2, double *x, double *y)
146 {
147  double d;
148
149  if (fabs(a2 * b1 - a1 * b2) == 0) {
150  if (fabs(a2 * c1 - a1 * c2) == 0)
151  return 2;
152  else
153  return 0;
154  }
155  else {
156  d = a1 * b2 - a2 * b1;
157  *x = (b1 * c2 - b2 * c1) / d;
158  *y = (c1 * a2 - c2 * a1) / d;
159  return 1;
160  }
161 }
162
163 static double angular_tolerance(double tol, double da, double db)
164 {
165  double a = MAX(da, db);
166
167  if (tol > a)
168  tol = a;
169
170  return 2 * acos(1 - tol / a);
171 }
172
173 /*
174  * This function generates parallel line (with loops, but not like the old ones).
175  * It is not to be used directly for creating buffers.
176  * + added elliptical buffers/par.lines support
177  *
178  * dalpha - direction of elliptical buffer major axis in degrees
179  * da - distance along major axis
180  * db: distance along minor (perp.) axis
181  * side: side >= 0 - right side, side < 0 - left side
182  * when (da == db) we have plain distances (old case)
183  * round - 1 for round corners, 0 for sharp corners. (tol is used only if round == 1)
184  */
185 static void parallel_line(struct line_pnts *Points, double da, double db,
186  double dalpha, int side, int round, int caps, int looped,
187  double tol, struct line_pnts *nPoints)
188 {
189  int i, j, res, np;
190  double *x, *y;
191  double tx, ty, vx, vy, wx, wy, nx, ny, mx, my, rx, ry;
192  double vx1, vy1, wx1, wy1;
193  double a0, b0, c0, a1, b1, c1;
194  double phi1, phi2, delta_phi;
195  double nsegments, angular_tol, angular_step;
196  int inner_corner, turns360;
197
198  G_debug(3, "parallel_line()");
199
200  if (looped && 0) {
201  /* start point != end point */
202  return;
203  }
204
205  Vect_reset_line(nPoints);
206
207  if (looped) {
208  Vect_append_point(Points, Points->x[1], Points->y[1], Points->z[1]);
209  }
210  np = Points->n_points;
211  x = Points->x;
212  y = Points->y;
213
214  if ((np == 0) || (np == 1))
215  return;
216
217  if ((da == 0) || (db == 0)) {
218  Vect_copy_xyz_to_pnts(nPoints, x, y, NULL, np);
219  return;
220  }
221
222  side = (side >= 0) ? (1) : (-1); /* normalize variable */
223  dalpha *= PI / 180; /* convert dalpha from degrees to radians */
224  angular_tol = angular_tolerance(tol, da, db);
225
226  for (i = 0; i < np - 1; i++) {
227  /* save the old values */
228  a0 = a1;
229  b0 = b1;
230  c0 = c1;
231  wx = vx;
232  wy = vy;
233
234
235  norm_vector(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);
236  if ((tx == 0) && (ty == 0))
237  continue;
238
239  elliptic_tangent(side * tx, side * ty, da, db, dalpha, &vx, &vy);
240
241  nx = x[i] + vx;
242  ny = y[i] + vy;
243
244  mx = x[i + 1] + vx;
245  my = y[i + 1] + vy;
246
247  line_coefficients(nx, ny, mx, my, &a1, &b1, &c1);
248
249  if (i == 0) {
250  if (!looped)
251  Vect_append_point(nPoints, nx, ny, 0);
252  continue;
253  }
254
255  delta_phi = atan2(ty, tx) - atan2(y[i] - y[i - 1], x[i] - x[i - 1]);
256  if (delta_phi > PI)
257  delta_phi -= 2 * PI;
258  else if (delta_phi <= -PI)
259  delta_phi += 2 * PI;
260  /* now delta_phi is in [-pi;pi] */
261  turns360 = (fabs(fabs(delta_phi) - PI) < 1e-15);
262  inner_corner = (side * delta_phi <= 0) && (!turns360);
263
264  if ((turns360) && (!(caps && round))) {
265  if (caps) {
266  norm_vector(0, 0, vx, vy, &tx, &ty);
267  elliptic_tangent(side * tx, side * ty, da, db, dalpha, &tx,
268  &ty);
269  }
270  else {
271  tx = 0;
272  ty = 0;
273  }
274  Vect_append_point(nPoints, x[i] + wx + tx, y[i] + wy + ty, 0);
275  Vect_append_point(nPoints, nx + tx, ny + ty, 0); /* nx == x[i] + vx, ny == y[i] + vy */
276  }
277  else if ((!round) || inner_corner) {
278  res = line_intersection(a0, b0, c0, a1, b1, c1, &rx, &ry);
279  /* if (res == 0) {
280  G_debug(4, "a0=%.18f, b0=%.18f, c0=%.18f, a1=%.18f, b1=%.18f, c1=%.18f", a0, b0, c0, a1, b1, c1);
281  G_fatal_error("Two consecutive line segments are parallel, but not on one straight line! This should never happen.");
282  return;
283  } */
284  if (res == 1) {
285  if (!round)
286  Vect_append_point(nPoints, rx, ry, 0);
287  else {
288  /* d = dig_distance2_point_to_line(rx, ry, 0, x[i-1], y[i-1], 0, x[i], y[i], 0,
289  0, NULL, NULL, NULL, NULL, NULL);
290  if ( */
291  Vect_append_point(nPoints, rx, ry, 0);
292  }
293  }
294  }
295  else {
296  /* we should draw elliptical arc for outside corner */
297
298  /* inverse transforms */
299  elliptic_transform(wx, wy, 1 / da, 1 / db, dalpha, &wx1, &wy1);
300  elliptic_transform(vx, vy, 1 / da, 1 / db, dalpha, &vx1, &vy1);
301
302  phi1 = atan2(wy1, wx1);
303  phi2 = atan2(vy1, vx1);
304  delta_phi = side * (phi2 - phi1);
305
306  /* make delta_phi in [0, 2pi] */
307  if (delta_phi < 0)
308  delta_phi += 2 * PI;
309
310  nsegments = (int)(delta_phi / angular_tol) + 1;
311  angular_step = side * (delta_phi / nsegments);
312
313  for (j = 0; j <= nsegments; j++) {
314  elliptic_transform(cos(phi1), sin(phi1), da, db, dalpha, &tx,
315  &ty);
316  Vect_append_point(nPoints, x[i] + tx, y[i] + ty, 0);
317  phi1 += angular_step;
318  }
319  }
320
321  if ((!looped) && (i == np - 2)) {
322  Vect_append_point(nPoints, mx, my, 0);
323  }
324  }
325
326  if (looped) {
327  Vect_append_point(nPoints, nPoints->x[0], nPoints->y[0],
328  nPoints->z[0]);
329  }
330
331  Vect_line_prune(nPoints);
332
333  if (looped) {
334  Vect_line_delete_point(Points, Points->n_points - 1);
335  }
336 }
337
338 /* input line must be looped */
339 static void convolution_line(struct line_pnts *Points, double da, double db,
340  double dalpha, int side, int round, int caps,
341  double tol, struct line_pnts *nPoints)
342 {
343  int i, j, res, np;
344  double *x, *y;
345  double tx, ty, vx, vy, wx, wy, nx, ny, mx, my, rx, ry;
346  double vx1, vy1, wx1, wy1;
347  double a0, b0, c0, a1, b1, c1;
348  double phi1, phi2, delta_phi;
349  double nsegments, angular_tol, angular_step;
350  double angle0, angle1;
351  int inner_corner, turns360;
352
353  G_debug(3, "convolution_line() side = %d", side);
354
355  np = Points->n_points;
356  x = Points->x;
357  y = Points->y;
358  if ((np == 0) || (np == 1))
359  return;
360  if ((x[0] != x[np - 1]) || (y[0] != y[np - 1])) {
361  G_fatal_error(_("Line is not looped"));
362  return;
363  }
364
365  Vect_reset_line(nPoints);
366
367  if ((da == 0) || (db == 0)) {
368  Vect_copy_xyz_to_pnts(nPoints, x, y, NULL, np);
369  return;
370  }
371
372  side = (side >= 0) ? (1) : (-1); /* normalize variable */
373  dalpha *= PI / 180; /* convert dalpha from degrees to radians */
374  angular_tol = angular_tolerance(tol, da, db);
375
376  i = np - 2;
377  norm_vector(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);
378  elliptic_tangent(side * tx, side * ty, da, db, dalpha, &vx, &vy);
379  angle1 = atan2(ty, tx);
380  nx = x[i] + vx;
381  ny = y[i] + vy;
382  mx = x[i + 1] + vx;
383  my = y[i + 1] + vy;
384  if (!round)
385  line_coefficients(nx, ny, mx, my, &a1, &b1, &c1);
386
387  for (i = 0; i <= np - 2; i++) {
388  G_debug(4, "point %d, segment %d-%d", i, i, i + 1);
389  /* save the old values */
390  if (!round) {
391  a0 = a1;
392  b0 = b1;
393  c0 = c1;
394  }
395  wx = vx;
396  wy = vy;
397  angle0 = angle1;
398
399  norm_vector(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);
400  if ((tx == 0) && (ty == 0))
401  continue;
402  elliptic_tangent(side * tx, side * ty, da, db, dalpha, &vx, &vy);
403  angle1 = atan2(ty, tx);
404  nx = x[i] + vx;
405  ny = y[i] + vy;
406  mx = x[i + 1] + vx;
407  my = y[i + 1] + vy;
408  if (!round)
409  line_coefficients(nx, ny, mx, my, &a1, &b1, &c1);
410
411
412  delta_phi = angle1 - angle0;
413  if (delta_phi > PI)
414  delta_phi -= 2 * PI;
415  else if (delta_phi <= -PI)
416  delta_phi += 2 * PI;
417  /* now delta_phi is in [-pi;pi] */
418  turns360 = (fabs(fabs(delta_phi) - PI) < 1e-15);
419  inner_corner = (side * delta_phi <= 0) && (!turns360);
420
421
422  /* if <line turns 360> and (<caps> and <not round>) */
423  if (turns360 && caps && (!round)) {
424  norm_vector(0, 0, vx, vy, &tx, &ty);
425  elliptic_tangent(side * tx, side * ty, da, db, dalpha, &tx, &ty);
426  Vect_append_point(nPoints, x[i] + wx + tx, y[i] + wy + ty, 0);
427  G_debug(4, " append point (c) x=%.16f y=%.16f", x[i] + wx + tx,
428  y[i] + wy + ty);
429  Vect_append_point(nPoints, nx + tx, ny + ty, 0); /* nx == x[i] + vx, ny == y[i] + vy */
430  G_debug(4, " append point (c) x=%.16f y=%.16f", nx + tx, ny + ty);
431  }
432
433  if ((!turns360) && (!round) && (!inner_corner)) {
434  res = line_intersection(a0, b0, c0, a1, b1, c1, &rx, &ry);
435  if (res == 1) {
436  Vect_append_point(nPoints, rx, ry, 0);
437  G_debug(4, " append point (o) x=%.16f y=%.16f", rx, ry);
438  }
439  else if (res == 2) {
440  /* no need to append point in this case */
441  }
442  else
443  G_fatal_error(_("Unexpected result of line_intersection() res = %d"),
444  res);
445  }
446
447  if (round && (!inner_corner) && (!turns360 || caps)) {
448  /* we should draw elliptical arc for outside corner */
449
450  /* inverse transforms */
451  elliptic_transform(wx, wy, 1 / da, 1 / db, dalpha, &wx1, &wy1);
452  elliptic_transform(vx, vy, 1 / da, 1 / db, dalpha, &vx1, &vy1);
453
454  phi1 = atan2(wy1, wx1);
455  phi2 = atan2(vy1, vx1);
456  delta_phi = side * (phi2 - phi1);
457
458  /* make delta_phi in [0, 2pi] */
459  if (delta_phi < 0)
460  delta_phi += 2 * PI;
461
462  nsegments = (int)(delta_phi / angular_tol) + 1;
463  angular_step = side * (delta_phi / nsegments);
464
465  phi1 += angular_step;
466  for (j = 1; j <= nsegments - 1; j++) {
467  elliptic_transform(cos(phi1), sin(phi1), da, db, dalpha, &tx,
468  &ty);
469  Vect_append_point(nPoints, x[i] + tx, y[i] + ty, 0);
470  G_debug(4, " append point (r) x=%.16f y=%.16f", x[i] + tx,
471  y[i] + ty);
472  phi1 += angular_step;
473  }
474  }
475
476  Vect_append_point(nPoints, nx, ny, 0);
477  G_debug(4, " append point (s) x=%.16f y=%.16f", nx, ny);
478  Vect_append_point(nPoints, mx, my, 0);
479  G_debug(4, " append point (s) x=%.16f y=%.16f", mx, my);
480  }
481
482  /* close the output line */
483  Vect_append_point(nPoints, nPoints->x[0], nPoints->y[0], nPoints->z[0]);
484  Vect_line_prune ( nPoints );
485 }
486
487 /*
488  * side: side >= 0 - extracts contour on right side of edge, side < 0 - extracts contour on left side of edge
489  * if the extracted contour is the outer contour, it is returned in ccw order
490  * else if it is inner contour, it is returned in cw order
491  */
492 static void extract_contour(struct planar_graph *pg, struct pg_edge *first,
493  int side, int winding, int stop_at_line_end,
494  struct line_pnts *nPoints)
495 {
496  int j;
497  int v; /* current vertex number */
498  int v0;
499  int eside; /* side of the current edge */
500  double eangle; /* current edge angle with Ox (according to the current direction) */
501  struct pg_vertex *vert; /* current vertex */
502  struct pg_vertex *vert0; /* last vertex */
503  struct pg_edge *edge; /* current edge; must be edge of vert */
504
505  /* int cs; *//* on which side are we turning along the contour */
506  /* we will always turn right and don't need that one */
507  double opt_angle, tangle;
508  int opt_j, opt_side, opt_flag;
509
510  G_debug(3, "extract_contour(): v1=%d, v2=%d, side=%d, stop_at_line_end=%d",
511  first->v1, first->v2, side, stop_at_line_end);
512
513  Vect_reset_line(nPoints);
514
515  edge = first;
516  if (side >= 0) {
517  eside = 1;
518  v0 = edge->v1;
519  v = edge->v2;
520  }
521  else {
522  eside = -1;
523  v0 = edge->v2;
524  v = edge->v1;
525  }
526  vert0 = &(pg->v[v0]);
527  vert = &(pg->v[v]);
528  eangle = atan2(vert->y - vert0->y, vert->x - vert0->x);
529
530  while (1) {
531  Vect_append_point(nPoints, vert0->x, vert0->y, 0);
532  G_debug(4, "ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d", v0,
533  v, eside, edge->v1, edge->v2);
534  G_debug(4, "ec: append point x=%.18f y=%.18f", vert0->x, vert0->y);
535
536  /* mark current edge as visited on the appropriate side */
537  if (eside == 1) {
538  edge->visited_right = 1;
539  edge->winding_right = winding;
540  }
541  else {
542  edge->visited_left = 1;
543  edge->winding_left = winding;
544  }
545
546  opt_flag = 1;
547  for (j = 0; j < vert->ecount; j++) {
548  /* exclude current edge */
549  if (vert->edges[j] != edge) {
550  tangle = vert->angles[j] - eangle;
551  if (tangle < -PI)
552  tangle += 2 * PI;
553  else if (tangle > PI)
554  tangle -= 2 * PI;
555  /* now tangle is in (-PI, PI) */
556
557  if (opt_flag || (tangle < opt_angle)) {
558  opt_j = j;
559  opt_side = (vert->edges[j]->v1 == v) ? (1) : (-1);
560  opt_angle = tangle;
561  opt_flag = 0;
562  }
563  }
564  }
565
566  /*
567  G_debug(4, "ec: opt: side=%d opt_flag=%d opt_angle=%.18f opt_j=%d opt_step=%d",
568  side, opt_flag, opt_angle, opt_j, opt_step);
569  */
570
571  /* if line end is reached (no other edges at curr vertex) */
572  if (opt_flag) {
573  if (stop_at_line_end) {
574  G_debug(3, " end has been reached, will stop here");
575  break;
576  }
577  else {
578  opt_j = 0; /* the only edge of vert is vert->edges[0] */
579  opt_side = -eside; /* go to the other side of the current edge */
580  G_debug(3, " end has been reached, turning around");
581  }
582  }
583
584  /* break condition */
585  if ((vert->edges[opt_j] == first) && (opt_side == side))
586  break;
587  if (opt_side == 1) {
588  if (vert->edges[opt_j]->visited_right) {
589  G_warning(_("Next edge was visited (right) but it is not the first one !!! breaking loop"));
590  G_debug(4,
591  "ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d",
592  v, (edge->v1 == v) ? (edge->v2) : (edge->v1),
593  opt_side, vert->edges[opt_j]->v1,
594  vert->edges[opt_j]->v2);
595  break;
596  }
597  }
598  else {
599  if (vert->edges[opt_j]->visited_left) {
600  G_warning(_("Next edge was visited (left) but it is not the first one !!! breaking loop"));
601  G_debug(4,
602  "ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d",
603  v, (edge->v1 == v) ? (edge->v2) : (edge->v1),
604  opt_side, vert->edges[opt_j]->v1,
605  vert->edges[opt_j]->v2);
606  break;
607  }
608  }
609
610  edge = vert->edges[opt_j];
611  eside = opt_side;
612  v0 = v;
613  v = (edge->v1 == v) ? (edge->v2) : (edge->v1);
614  vert0 = vert;
615  vert = &(pg->v[v]);
616  eangle = vert0->angles[opt_j];
617  }
618  Vect_append_point(nPoints, vert->x, vert->y, 0);
619  Vect_line_prune(nPoints);
620  G_debug(4, "ec: append point x=%.18f y=%.18f", vert->x, vert->y);
621
622  return;
623 }
624
625 /*
626  * This function extracts the outer contour of a (self crossing) line.
627  * It can generate left/right contour if none of the line ends are in a loop.
628  * If one or both of them is in a loop, then there's only one contour
629  *
630  * side: side > 0 - right contour, side < 0 - left contour, side = 0 - outer contour
631  * if side != 0 and there's only one contour, the function returns it
632  *
633  * TODO: Implement side != 0 feature;
634  */
635 static void extract_outer_contour(struct planar_graph *pg, int side,
636  struct line_pnts *nPoints)
637 {
638  int i;
639  int flag;
640  int v;
641  struct pg_vertex *vert;
642  struct pg_edge *edge;
643  double min_x, min_angle;
644
645  G_debug(3, "extract_outer_contour()");
646
647  if (side != 0) {
648  G_fatal_error(_("side != 0 feature not implemented"));
649  return;
650  }
651
652  /* find a line segment which is on the outer contour */
653  flag = 1;
654  for (i = 0; i < pg->vcount; i++) {
655  if (flag || (pg->v[i].x < min_x)) {
656  v = i;
657  min_x = pg->v[i].x;
658  flag = 0;
659  }
660  }
661  vert = &(pg->v[v]);
662
663  flag = 1;
664  for (i = 0; i < vert->ecount; i++) {
665  if (flag || (vert->angles[i] < min_angle)) {
666  edge = vert->edges[i];
667  min_angle = vert->angles[i];
668  flag = 0;
669  }
670  }
671
672  /* the winding on the outer contour is 0 */
673  extract_contour(pg, edge, (edge->v1 == v) ? RIGHT_SIDE : LEFT_SIDE, 0, 0,
674  nPoints);
675
676  return;
677 }
678
679 /*
680  * Extracts contours which are not visited.
681  * IMPORTANT: the outer contour must be visited (you should call extract_outer_contour() to do that),
682  * so that extract_inner_contour() doesn't return it
683  *
684  * returns: 0 when there are no more inner contours; otherwise, 1
685  */
686 static int extract_inner_contour(struct planar_graph *pg, int *winding,
687  struct line_pnts *nPoints)
688 {
689  int i, w;
690  struct pg_edge *edge;
691
692  G_debug(3, "extract_inner_contour()");
693
694  for (i = 0; i < pg->ecount; i++) {
695  edge = &(pg->e[i]);
696  if (edge->visited_left) {
697  if (!(pg->e[i].visited_right)) {
698  w = edge->winding_left - 1;
699  extract_contour(pg, &(pg->e[i]), RIGHT_SIDE, w, 0, nPoints);
700  *winding = w;
701  return 1;
702  }
703  }
704  else {
705  if (pg->e[i].visited_right) {
706  w = edge->winding_right + 1;
707  extract_contour(pg, &(pg->e[i]), LEFT_SIDE, w, 0, nPoints);
708  *winding = w;
709  return 1;
710  }
711  }
712  }
713
714  return 0;
715 }
716
717 /* point_in_buf - test if point px,py is in d buffer of Points
718  ** dalpha is in degrees
719  ** returns: 1 in buffer
720  ** 0 not in buffer
721  */
722 static int point_in_buf(struct line_pnts *Points, double px, double py, double da,
723  double db, double dalpha)
724 {
725  int i, np;
726  double cx, cy;
727  double delta, delta_k, k;
728  double vx, vy, wx, wy, mx, my, nx, ny;
729  double len, tx, ty, d, da2;
730
731  G_debug(3, "point_in_buf()");
732
733  dalpha *= PI / 180; /* convert dalpha from degrees to radians */
734
735  np = Points->n_points;
736  da2 = da * da;
737  for (i = 0; i < np - 1; i++) {
738  vx = Points->x[i];
739  vy = Points->y[i];
740  wx = Points->x[i + 1];
741  wy = Points->y[i + 1];
742
743  if (da != db) {
744  mx = wx - vx;
745  my = wy - vy;
746  len = LENGTH(mx, my);
747  elliptic_tangent(mx / len, my / len, da, db, dalpha, &cx, &cy);
748
749  delta = mx * cy - my * cx;
750  delta_k = (px - vx) * cy - (py - vy) * cx;
751  k = delta_k / delta;
752  /* G_debug(4, "k = %g, k1 = %g", k, (mx * (px - vx) + my * (py - vy)) / (mx * mx + my * my)); */
753  if (k <= 0) {
754  nx = vx;
755  ny = vy;
756  }
757  else if (k >= 1) {
758  nx = wx;
759  ny = wy;
760  }
761  else {
762  nx = vx + k * mx;
763  ny = vy + k * my;
764  }
765
766  /* inverse transform */
767  elliptic_transform(px - nx, py - ny, 1 / da, 1 / db, dalpha, &tx,
768  &ty);
769
770  d = dig_distance2_point_to_line(nx + tx, ny + ty, 0, vx, vy, 0,
771  wx, wy, 0, 0, NULL, NULL, NULL,
772  NULL, NULL);
773
774  /* G_debug(4, "sqrt(d)*da = %g, len' = %g, olen = %g", sqrt(d)*da, da*LENGTH(tx,ty), LENGTH((px-nx),(py-ny))); */
775  if (d <= 1) {
776  /* G_debug(1, "d=%g", d); */
777  return 1;
778  }
779  }
780  else {
781  d = dig_distance2_point_to_line(px, py, 0, vx, vy, 0, wx, wy, 0,
782  0, NULL, NULL, NULL, NULL, NULL);
783  /* G_debug(4, "sqrt(d) = %g", sqrt(d)); */
784  if (d <= da2) {
785  return 1;
786  }
787  }
788  }
789
790  return 0;
791 }
792
793 /* returns 0 for ccw, non-zero for cw
794  */
795 static int get_polygon_orientation(const double *x, const double *y, int n)
796 {
797  double x1, y1, x2, y2;
798  double area;
799
800  x2 = x[n - 1];
801  y2 = y[n - 1];
802
803  area = 0;
804  while (--n >= 0) {
805  x1 = x2;
806  y1 = y2;
807
808  x2 = *x++;
809  y2 = *y++;
810
811  area += (y2 + y1) * (x2 - x1);
812  }
813
814  return (area > 0);
815 }
816
817 /* internal */
818 static void add_line_to_array(struct line_pnts *Points,
819  struct line_pnts ***arrPoints, int *count,
820  int *allocated, int more)
821 {
822  if (*allocated == *count) {
823  *allocated += more;
824  *arrPoints =
825  G_realloc(*arrPoints, (*allocated) * sizeof(struct line_pnts *));
826  }
827  (*arrPoints)[*count] = Points;
828  (*count)++;
829
830  return;
831 }
832
833 static void destroy_lines_array(struct line_pnts **arr, int count)
834 {
835  int i;
836
837  for (i = 0; i < count; i++)
838  Vect_destroy_line_struct(arr[i]);
839  G_free(arr);
840 }
841
842 /* area_outer and area_isles[i] must be closed non self-intersecting lines
843  side: 0 - auto, 1 - right, -1 left
844  */
845 static void buffer_lines(struct line_pnts *area_outer, struct line_pnts **area_isles,
846  int isles_count, int side, double da, double db,
847  double dalpha, int round, int caps, double tol,
848  struct line_pnts **oPoints, struct line_pnts ***iPoints,
849  int *inner_count)
850 {
851  struct planar_graph *pg2;
852  struct line_pnts *sPoints, *cPoints;
853  struct line_pnts **arrPoints;
854  int i, count = 0;
855  int res, winding;
856  int auto_side;
857  int more = 8;
858  int allocated = 0;
859  double px, py;
860
861  G_debug(3, "buffer_lines()");
862
863  auto_side = (side == 0);
864
865  /* initializations */
866  sPoints = Vect_new_line_struct();
867  cPoints = Vect_new_line_struct();
868  arrPoints = NULL;
869
870  /* outer contour */
871  G_debug(3, " processing outer contour");
872  *oPoints = Vect_new_line_struct();
873  if (auto_side)
874  side =
875  get_polygon_orientation(area_outer->x, area_outer->y,
876  area_outer->n_points -
877  1) ? LEFT_SIDE : RIGHT_SIDE;
878  convolution_line(area_outer, da, db, dalpha, side, round, caps, tol,
879  sPoints);
880  pg2 = pg_create(sPoints);
881  extract_outer_contour(pg2, 0, *oPoints);
882  res = extract_inner_contour(pg2, &winding, cPoints);
883  while (res != 0) {
884  if (winding == 0) {
885  int check_poly = 1;
886  double area_size;
887
888  dig_find_area_poly(cPoints, &area_size);
889  if (area_size == 0) {
890  G_warning(_("zero area size"));
891  check_poly = 0;
892  }
893  if (cPoints->x[0] != cPoints->x[cPoints->n_points - 1] ||
894  cPoints->y[0] != cPoints->y[cPoints->n_points - 1]) {
895
896  G_warning(_("Line was not closed"));
897  check_poly = 0;
898  }
899
900  if (check_poly && !Vect_point_in_poly(cPoints->x[0], cPoints->y[0], area_outer)) {
901  if (Vect_get_point_in_poly(cPoints, &px, &py) == 0) {
902  if (!point_in_buf(area_outer, px, py, da, db, dalpha)) {
904  more);
905  cPoints = Vect_new_line_struct();
906  }
907  }
908  else {
909  G_warning(_("Vect_get_point_in_poly() failed"));
910  }
911  }
912  }
913  res = extract_inner_contour(pg2, &winding, cPoints);
914  }
915  pg_destroy_struct(pg2);
916
917  /* inner contours */
918  G_debug(3, " processing inner contours");
919  for (i = 0; i < isles_count; i++) {
920  if (auto_side)
921  side =
922  get_polygon_orientation(area_isles[i]->x, area_isles[i]->y,
923  area_isles[i]->n_points -
924  1) ? RIGHT_SIDE : LEFT_SIDE;
925  convolution_line(area_isles[i], da, db, dalpha, side, round, caps,
926  tol, sPoints);
927  pg2 = pg_create(sPoints);
928  extract_outer_contour(pg2, 0, cPoints);
929  res = extract_inner_contour(pg2, &winding, cPoints);
930  while (res != 0) {
931  if (winding == -1) {
932  int check_poly = 1;
933  double area_size;
934
935  dig_find_area_poly(cPoints, &area_size);
936  if (area_size == 0) {
937  G_warning(_("zero area size"));
938  check_poly = 0;
939  }
940  if (cPoints->x[0] != cPoints->x[cPoints->n_points - 1] ||
941  cPoints->y[0] != cPoints->y[cPoints->n_points - 1]) {
942
943  G_warning(_("Line was not closed"));
944  check_poly = 0;
945  }
946
947  /* we need to check if the area is in the buffer.
948  I've simplfied convolution_line(), so that it runs faster,
949  however that leads to occasional problems */
950  if (check_poly && Vect_point_in_poly
951  (cPoints->x[0], cPoints->y[0], area_isles[i])) {
952  if (Vect_get_point_in_poly(cPoints, &px, &py) == 0) {
953  if (!point_in_buf(area_isles[i], px, py, da, db, dalpha)) {
955  &allocated, more);
956  cPoints = Vect_new_line_struct();
957  }
958  }
959  else {
960  G_warning(_("Vect_get_point_in_poly() failed"));
961  }
962  }
963  }
964  res = extract_inner_contour(pg2, &winding, cPoints);
965  }
966  pg_destroy_struct(pg2);
967  }
968
969  arrPoints = G_realloc(arrPoints, count * sizeof(struct line_pnts *));
970  *inner_count = count;
971  *iPoints = arrPoints;
972
973  Vect_destroy_line_struct(sPoints);
974  Vect_destroy_line_struct(cPoints);
975
976  G_debug(3, "buffer_lines() ... done");
977
978  return;
979 }
980
981
982 /*!
983  \brief Creates buffer around line.
984
986
987  Shape of buffer endings is managed by two parameters - round and cap.
988  Setting round=1, cap=1 gives "classical" buffer, while
989  round=0, cap=1 gives square end, but cap=0 – butt.
990  See v.buffer manual or SVG stroke-linecap for examples.
991
992  To get "classical" buffer, set db equal to da, and dalpha to 0.
993
994  \param Points input line geometry
995  \param da distance along major axis
996  \param db distance along minor axis
997  \param dalpha angle between 0x and major axis
998  \param round make corners round (0 - square, not 0 - round)
999  \param caps add caps at line ends (0 - butt, not 0 - caps)
1000  \param tol maximum distance between theoretical arc and output segments
1001  \param[out] oPoints output polygon outer border (ccw order)
1002  \param[out] iPoints array of output polygon's holes (cw order)
1003  \param[out] inner_count number of holes
1004  */
1005 void Vect_line_buffer2(const struct line_pnts *Points, double da, double db,
1006  double dalpha, int round, int caps, double tol,
1007  struct line_pnts **oPoints,
1008  struct line_pnts ***iPoints, int *inner_count)
1009 {
1010  struct planar_graph *pg;
1011  struct line_pnts *tPoints, *outer;
1012  struct line_pnts **isles;
1013  int isles_count = 0;
1014  int res, winding;
1015  int more = 8;
1016  int isles_allocated = 0;
1017
1018  G_debug(2, "Vect_line_buffer()");
1019
1020  Vect_line_prune((struct line_pnts *)Points);
1021
1022  if (Points->n_points == 1)
1023  return Vect_point_buffer2(Points->x[0], Points->y[0], da, db,
1024  dalpha, round, tol, oPoints);
1025
1026  /* initializations */
1027  tPoints = Vect_new_line_struct();
1028  isles = NULL;
1029  pg = pg_create(Points);
1030
1031  /* outer contour */
1032  outer = Vect_new_line_struct();
1033  extract_outer_contour(pg, 0, outer);
1034
1035  /* inner contours */
1036  res = extract_inner_contour(pg, &winding, tPoints);
1037  while (res != 0) {
1039  more);
1040  tPoints = Vect_new_line_struct();
1041  res = extract_inner_contour(pg, &winding, tPoints);
1042  }
1043
1044  buffer_lines(outer, isles, isles_count, RIGHT_SIDE, da, db, dalpha, round,
1045  caps, tol, oPoints, iPoints, inner_count);
1046
1047  Vect_destroy_line_struct(tPoints);
1048  Vect_destroy_line_struct(outer);
1049  destroy_lines_array(isles, isles_count);
1050  pg_destroy_struct(pg);
1051
1052  return;
1053 }
1054
1055 /*!
1056  \brief Creates buffer around area.
1057
1058  \param Map vector map
1059  \param area area id
1060  \param da distance along major axis
1061  \param db distance along minor axis
1062  \param dalpha angle between 0x and major axis
1063  \param round make corners round
1064  \param caps add caps at line ends
1065  \param tol maximum distance between theoretical arc and output segments
1066  \param[out] oPoints output polygon outer border (ccw order)
1067  \param[out] inner_count number of holes
1068  \param[out] iPoints array of output polygon's holes (cw order)
1069  */
1070 void Vect_area_buffer2(const struct Map_info *Map, int area, double da, double db,
1071  double dalpha, int round, int caps, double tol,
1072  struct line_pnts **oPoints,
1073  struct line_pnts ***iPoints, int *inner_count)
1074 {
1075  struct line_pnts *tPoints, *outer;
1076  struct line_pnts **isles;
1077  int isles_count = 0, n_isles;
1078  int i, isle;
1079  int more = 8;
1080  int isles_allocated = 0;
1081
1082  G_debug(2, "Vect_area_buffer()");
1083
1084  /* initializations */
1085  tPoints = Vect_new_line_struct();
1086  n_isles = Vect_get_area_num_isles(Map, area);
1087  isles_allocated = n_isles;
1088  isles = G_malloc(isles_allocated * sizeof(struct line_pnts *));
1089
1090  /* outer contour */
1091  outer = Vect_new_line_struct();
1092  Vect_get_area_points(Map, area, outer);
1093  /* does not work with zero length line segments */
1094  Vect_line_prune(outer);
1095
1096  /* inner contours */
1097  for (i = 0; i < n_isles; i++) {
1098  isle = Vect_get_area_isle(Map, area, i);
1099  Vect_get_isle_points(Map, isle, tPoints);
1100
1101  /* Check if the isle is big enough */
1102  /*
1103  if (Vect_line_length(tPoints) < 2*PI*max)
1104  continue;
1105  */
1106  /* does not work with zero length line segments */
1107  Vect_line_prune(tPoints);
1109  more);
1110  tPoints = Vect_new_line_struct();
1111  }
1112
1113  buffer_lines(outer, isles, isles_count, 0, da, db, dalpha, round, caps,
1114  tol, oPoints, iPoints, inner_count);
1115
1116  Vect_destroy_line_struct(tPoints);
1117  Vect_destroy_line_struct(outer);
1118  destroy_lines_array(isles, isles_count);
1119
1120  return;
1121 }
1122
1123 /*!
1124  \brief Creates buffer around the point (px, py).
1125
1126  \param px input point x-coordinate
1127  \param py input point y-coordinate
1128  \param da distance along major axis
1129  \param db distance along minor axis
1130  \param dalpha angle between 0x and major axis
1131  \param round make corners round
1132  \param tol maximum distance between theoretical arc and output segments
1133  \param[out] oPoints output polygon outer border (ccw order)
1134  */
1135 void Vect_point_buffer2(double px, double py, double da, double db,
1136  double dalpha, int round, double tol,
1137  struct line_pnts **oPoints)
1138 {
1139  double tx, ty;
1140  double angular_tol, angular_step, phi1;
1141  int j, nsegments;
1142
1143  G_debug(2, "Vect_point_buffer()");
1144
1145  *oPoints = Vect_new_line_struct();
1146
1147  dalpha *= PI / 180; /* convert dalpha from degrees to radians */
1148
1149  if (round || (!round)) {
1150  angular_tol = angular_tolerance(tol, da, db);
1151
1152  nsegments = (int)(2 * PI / angular_tol) + 1;
1153  angular_step = 2 * PI / nsegments;
1154
1155  phi1 = 0;
1156  for (j = 0; j < nsegments; j++) {
1157  elliptic_transform(cos(phi1), sin(phi1), da, db, dalpha, &tx,
1158  &ty);
1159  Vect_append_point(*oPoints, px + tx, py + ty, 0);
1160  phi1 += angular_step;
1161  }
1162  }
1163  else {
1164
1165  }
1166
1167  /* close the output line */
1168  Vect_append_point(*oPoints, (*oPoints)->x[0], (*oPoints)->y[0],
1169  (*oPoints)->z[0]);
1170
1171  return;
1172 }
1173
1174
1175 /*
1176  \brief Create parallel line
1177
1179
1180  \param InPoints input line geometry
1181  \param da distance along major axis
1182  \param da distance along minor axis
1183  \param dalpha angle between 0x and major axis
1184  \param round make corners round
1185  \param tol maximum distance between theoretical arc and output segments
1186  \param[out] OutPoints output line
1187  */
1188 void Vect_line_parallel2(struct line_pnts *InPoints, double da, double db,
1189  double dalpha, int side, int round, double tol,
1190  struct line_pnts *OutPoints)
1191 {
1192  G_debug(2, "Vect_line_parallel(): npoints = %d, da = %f, "
1193  "db = %f, dalpha = %f, side = %d, round_corners = %d, tol = %f",
1194  InPoints->n_points, da, db, dalpha, side, round, tol);
1195
1196  parallel_line(InPoints, da, db, dalpha, side, round, 1, NON_LOOPED_LINE,
1197  tol, OutPoints);
1198
1199  /* if (!loops)
1200  clean_parallel(OutPoints, InPoints, distance, rm_end);
1201  */
1202
1203  return;
1204 }
#define G_malloc(n)
Definition: defs/gis.h:112
#define PI
Definition: buffer2.c:33
int Vect_point_in_poly(double, double, const struct line_pnts *)
Determines if a point (X,Y) is inside a polygon.
Definition: Vlib/poly.c:826
double x
Definition: dgraph.h:16
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
char winding_right
Definition: dgraph.h:12
void pg_destroy_struct(struct planar_graph *pg)
Definition: dgraph.c:368
#define MAX(X, Y)
Definition: buffer2.c:31
int Vect_get_area_isle(const struct Map_info *, int, int)
Returns isle id for area.
int n_points
Number of points.
Definition: dig_structs.h:1692
char visited_left
Definition: dgraph.h:9
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
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:149
int count
int Vect_get_area_num_isles(const struct Map_info *, int)
Returns number of isles for given area.
int ecount
Definition: dgraph.h:27
int Vect_line_prune(struct line_pnts *)
Remove duplicate points, i.e. zero length segments.
Definition: line.c:281
struct pg_edge * e
Definition: dgraph.h:29
#define NULL
Definition: ccmath.h:32
#define x
double * x
Array of X coordinates.
Definition: dig_structs.h:1680
int vcount
Definition: dgraph.h:25
double l
Definition: r_raster.c:39
Feature geometry info - coordinates.
Definition: dig_structs.h:1675
struct planar_graph * pg_create(const struct line_pnts *Points)
Definition: dgraph.c:449
double * angles
Definition: dgraph.h:21
double b
Definition: r_raster.c:39
int Vect_get_point_in_poly(const struct line_pnts *, double *, double *)
Get point inside polygon.
Definition: Vlib/poly.c:211
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
struct pg_vertex * v
Definition: dgraph.h:26
struct pg_edge ** edges
Definition: dgraph.h:20
char visited_right
Definition: dgraph.h:10
double dig_distance2_point_to_line(double, double, double, double, double, double, double, double, double, int, double *, double *, double *, double *, int *)
int Vect_get_isle_points(const struct Map_info *, int, struct line_pnts *)
Returns polygon array of points for given isle.
int Vect_append_point(struct line_pnts *, double, double, double)
Appends one point to the end of a line.
Definition: line.c:149
char winding_left
Definition: dgraph.h:11
#define LENGTH(DX, DY)
Definition: buffer2.c:26
Definition: dgraph.h:6
#define LEFT_SIDE
Definition: buffer2.c:35
double y
Definition: dgraph.h:17
Vector map info.
Definition: dig_structs.h:1259
double * y
Array of Y coordinates.
Definition: dig_structs.h:1684
void Vect_line_parallel2(struct line_pnts *InPoints, double da, double db, double dalpha, int side, int round, double tol, struct line_pnts *OutPoints)
Definition: buffer2.c:1188
#define NON_LOOPED_LINE
Definition: buffer2.c:37
void G_warning(const char *,...) __attribute__((format(printf
int v2
Definition: dgraph.h:8
#define G_realloc(p, n)
Definition: defs/gis.h:114
double * z
Array of Z coordinates.
Definition: dig_structs.h:1688
int v1
Definition: dgraph.h:7
#define _(str)
Definition: glocale.h:10
void Vect_line_buffer2(const struct line_pnts *Points, double da, double db, double dalpha, int round, int caps, double tol, struct line_pnts **oPoints, struct line_pnts ***iPoints, int *inner_count)
Creates buffer around line.
Definition: buffer2.c:1005
int ecount
Definition: dgraph.h:18
void Vect_destroy_line_struct(struct line_pnts *)
Frees all memory associated with a line_pnts structure, including the structure itself.
Definition: line.c:77
int Vect_get_area_points(const struct Map_info *, int, struct line_pnts *)
Returns polygon array of points (outer ring) of given area.
int G_debug(int, const char *,...) __attribute__((format(printf
#define RIGHT_SIDE
Definition: buffer2.c:34
void Vect_area_buffer2(const struct Map_info *Map, int area, double da, double db, double dalpha, int round, int caps, double tol, struct line_pnts **oPoints, struct line_pnts ***iPoints, int *inner_count)
Creates buffer around area.
Definition: buffer2.c:1070
void Vect_point_buffer2(double px, double py, double da, double db, double dalpha, int round, double tol, struct line_pnts **oPoints)
Creates buffer around the point (px, py).
Definition: buffer2.c:1135
void Vect_reset_line(struct line_pnts *)
Reset line.
Definition: line.c:130
int Vect_line_delete_point(struct line_pnts *, int)
Delete point at given index and move all points above down.
Definition: line.c:211
int dig_find_area_poly(struct line_pnts *, double *)
Definition: diglib/poly.c:100