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