GRASS GIS 7 Programmer's Manual  7.5.svn(2017)-r71933
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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)) {
903  add_line_to_array(cPoints, &arrPoints, &count, &allocated,
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)) {
954  add_line_to_array(cPoints, &arrPoints, &count,
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 
985  See also Vect_line_buffer().
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) {
1038  add_line_to_array(tPoints, &isles, &isles_count, &isles_allocated,
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);
1108  add_line_to_array(tPoints, &isles, &isles_count, &isles_allocated,
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 
1178  See also Vect_line_parallel().
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 }
int Vect_point_in_poly(double X, double Y, const struct line_pnts *Points)
Determines if a point (X,Y) is inside a polygon.
Definition: Vlib/poly.c:826
#define PI
Definition: buffer2.c:33
double x
Definition: dgraph.h:16
char winding_right
Definition: dgraph.h:12
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:149
void pg_destroy_struct(struct planar_graph *pg)
Definition: dgraph.c:368
#define MAX(X, Y)
Definition: buffer2.c:31
struct line_pnts * Vect_new_line_struct()
Creates and initializes a line_pnts structure.
Definition: line.c:45
int n_points
Number of points.
Definition: dig_structs.h:1692
char visited_left
Definition: dgraph.h:9
int count
void Vect_destroy_line_struct(struct line_pnts *p)
Frees all memory associated with a line_pnts structure, including the structure itself.
Definition: line.c:77
int Vect_copy_xyz_to_pnts(struct line_pnts *Points, const double *x, const double *y, const double *z, int n)
Copy points from array to line_pnts structure.
Definition: line.c:99
int ecount
Definition: dgraph.h:27
struct pg_edge * e
Definition: dgraph.h:29
#define NULL
Definition: ccmath.h:32
#define x
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:149
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:159
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
struct pg_vertex * v
Definition: dgraph.h:26
struct pg_edge ** edges
Definition: dgraph.h:20
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
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 *)
char winding_left
Definition: dgraph.h:11
#define LENGTH(DX, DY)
Definition: buffer2.c:26
int Vect_get_point_in_poly(const struct line_pnts *Points, double *X, double *Y)
Get point inside polygon.
Definition: Vlib/poly.c:211
Definition: dgraph.h:6
#define LEFT_SIDE
Definition: buffer2.c:35
double y
Definition: dgraph.h:17
int Vect_get_area_num_isles(const struct Map_info *Map, int area)
Returns number of isles for given area.
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
int Vect_line_prune(struct line_pnts *Points)
Remove duplicate points, i.e. zero length segments.
Definition: line.c:281
int v2
Definition: dgraph.h:8
double * z
Array of Z coordinates.
Definition: dig_structs.h:1688
int v1
Definition: dgraph.h:7
#define _(str)
Definition: glocale.h:13
int Vect_line_delete_point(struct line_pnts *Points, int index)
Delete point at given index and move all points above down.
Definition: line.c:211
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 Vect_get_area_points(const struct Map_info *Map, int area, struct line_pnts *BPoints)
Returns polygon array of points (outer ring) of given area.
void Vect_reset_line(struct line_pnts *Points)
Reset line.
Definition: line.c:130
int ecount
Definition: dgraph.h:18
int
Reads the categories file for map name in mapset and stores the categories in the pcats structure...
int Vect_get_isle_points(const struct Map_info *Map, int isle, struct line_pnts *BPoints)
Returns polygon array of points for given isle.
int Vect_get_area_isle(const struct Map_info *Map, int area, int isle)
Returns isle id for area.
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:203
#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
int dig_find_area_poly(struct line_pnts *, double *)
Definition: diglib/poly.c:100