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