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