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