GRASS Programmer's Manual  6.5.svn(2014)-r66266
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)) {
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)) {
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) {
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);
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
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
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