27 #define LENGTH(DX, DY) (sqrt((DX*DX)+(DY*DY))) 32 #define NON_LOOPED_LINE 0 35 static void norm_vector(
double x1,
double y1,
double x2,
double y2,
double *
x,
42 if ((dx == 0) && (dy == 0)) {
56 static void rotate_vector(
double x,
double y,
double cosa,
double sina,
57 double *nx,
double *ny)
59 *nx =
x * cosa -
y * sina;
60 *ny =
x * sina +
y * cosa;
69 static void elliptic_transform(
double x,
double y,
double da,
double db,
70 double dalpha,
double *nx,
double *ny)
72 double cosa = cos(dalpha);
73 double sina = sin(dalpha);
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;
99 static void elliptic_tangent(
double x,
double y,
double da,
double db,
100 double dalpha,
double *px,
double *py)
102 double cosa = cos(dalpha);
103 double sina = sin(dalpha);
107 rotate_vector(
x,
y, cosa, -sina, &
x, &
y);
112 len = da * db / sqrt(da * da * v * v + db * db * u * u);
115 rotate_vector(u, v, cosa, sina, px, py);
124 static void line_coefficients(
double x1,
double y1,
double x2,
double y2,
125 double *a,
double *
b,
double *c)
129 *c = x2 * y1 - x1 * y2;
139 static int line_intersection(
double a1,
double b1,
double c1,
double a2,
140 double b2,
double c2,
double *
x,
double *
y)
144 if (fabs(a2 * b1 - a1 * b2) == 0) {
145 if (fabs(a2 * c1 - a1 * c2) == 0)
151 d = a1 * b2 - a2 * b1;
152 *
x = (b1 * c2 - b2 * c1) / d;
153 *y = (c1 * a2 - c2 * a1) / d;
158 static double angular_tolerance(
double tol,
double da,
double db)
160 double a =
MAX(da, db);
165 return 2 * acos(1 - tol / a);
180 static void parallel_line(
struct line_pnts *Points,
double da,
double db,
181 double dalpha,
int side,
int round,
int caps,
int looped,
186 double tx, ty, vx, vy, wx, wy, nx, ny, mx, my, rx, ry;
187 double vx1, vy1, wx1, wy1;
188 double a0, b0, c0, a1, b1, c1;
189 double phi1, phi2, delta_phi;
190 double nsegments, angular_tol, angular_step;
191 int inner_corner, turns360;
209 if ((np == 0) || (np == 1))
212 if ((da == 0) || (db == 0)) {
217 side = (side >= 0) ? (1) : (-1);
219 angular_tol = angular_tolerance(tol, da, db);
221 for (i = 0; i < np - 1; i++) {
230 norm_vector(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);
231 if ((tx == 0) && (ty == 0))
234 elliptic_tangent(side * tx, side * ty, da, db, dalpha, &vx, &vy);
242 line_coefficients(nx, ny, mx, my, &a1, &b1, &c1);
250 delta_phi = atan2(ty, tx) - atan2(y[i] - y[i - 1], x[i] - x[i - 1]);
253 else if (delta_phi <= -
PI)
256 turns360 = (fabs(fabs(delta_phi) -
PI) < 1e-15);
257 inner_corner = (side * delta_phi <= 0) && (!turns360);
259 if ((turns360) && (!(caps && round))) {
261 norm_vector(0, 0, vx, vy, &tx, &ty);
262 elliptic_tangent(side * tx, side * ty, da, db, dalpha, &tx,
272 else if ((!round) || inner_corner) {
273 res = line_intersection(a0, b0, c0, a1, b1, c1, &rx, &ry);
294 elliptic_transform(wx, wy, 1 / da, 1 / db, dalpha, &wx1, &wy1);
295 elliptic_transform(vx, vy, 1 / da, 1 / db, dalpha, &vx1, &vy1);
297 phi1 = atan2(wy1, wx1);
298 phi2 = atan2(vy1, vx1);
299 delta_phi = side * (phi2 - phi1);
305 nsegments = (int)(delta_phi / angular_tol) + 1;
306 angular_step = side * (delta_phi / nsegments);
308 for (j = 0; j <= nsegments; j++) {
309 elliptic_transform(cos(phi1), sin(phi1), da, db, dalpha, &tx,
312 phi1 += angular_step;
316 if ((!looped) && (i == np - 2)) {
334 static void convolution_line(
struct line_pnts *Points,
double da,
double db,
335 double dalpha,
int side,
int round,
int caps,
340 double tx, ty, vx, vy, wx, wy, nx, ny, mx, my, rx, ry;
341 double vx1, vy1, wx1, wy1;
342 double a0, b0, c0, a1, b1, c1;
343 double phi1, phi2, delta_phi;
344 double nsegments, angular_tol, angular_step;
345 double angle0, angle1;
346 int inner_corner, turns360;
348 G_debug(3,
"convolution_line() side = %d", side);
353 if ((np == 0) || (np == 1))
355 if ((x[0] != x[np - 1]) || (y[0] != y[np - 1])) {
362 if ((da == 0) || (db == 0)) {
367 side = (side >= 0) ? (1) : (-1);
369 angular_tol = angular_tolerance(tol, da, db);
372 norm_vector(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);
373 elliptic_tangent(side * tx, side * ty, da, db, dalpha, &vx, &vy);
374 angle1 = atan2(ty, tx);
380 line_coefficients(nx, ny, mx, my, &a1, &b1, &c1);
382 for (i = 0; i <= np - 2; i++) {
383 G_debug(4,
"point %d, segment %d-%d", i, i, i + 1);
394 norm_vector(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);
395 if ((tx == 0) && (ty == 0))
397 elliptic_tangent(side * tx, side * ty, da, db, dalpha, &vx, &vy);
398 angle1 = atan2(ty, tx);
404 line_coefficients(nx, ny, mx, my, &a1, &b1, &c1);
407 delta_phi = angle1 - angle0;
410 else if (delta_phi <= -
PI)
413 turns360 = (fabs(fabs(delta_phi) -
PI) < 1e-15);
414 inner_corner = (side * delta_phi <= 0) && (!turns360);
418 if (turns360 && caps && (!round)) {
419 norm_vector(0, 0, vx, vy, &tx, &ty);
420 elliptic_tangent(side * tx, side * ty, da, db, dalpha, &tx, &ty);
422 G_debug(4,
" append point (c) x=%.16f y=%.16f", x[i] + wx + tx,
425 G_debug(4,
" append point (c) x=%.16f y=%.16f", nx + tx, ny + ty);
428 if ((!turns360) && (!round) && (!inner_corner)) {
429 res = line_intersection(a0, b0, c0, a1, b1, c1, &rx, &ry);
432 G_debug(4,
" append point (o) x=%.16f y=%.16f", rx, ry);
438 G_fatal_error(
_(
"Unexpected result of line_intersection() res = %d"),
442 if (round && (!inner_corner) && (!turns360 || caps)) {
446 elliptic_transform(wx, wy, 1 / da, 1 / db, dalpha, &wx1, &wy1);
447 elliptic_transform(vx, vy, 1 / da, 1 / db, dalpha, &vx1, &vy1);
449 phi1 = atan2(wy1, wx1);
450 phi2 = atan2(vy1, vx1);
451 delta_phi = side * (phi2 - phi1);
457 nsegments = (int)(delta_phi / angular_tol) + 1;
458 angular_step = side * (delta_phi / nsegments);
460 phi1 += angular_step;
461 for (j = 1; j <= nsegments - 1; j++) {
462 elliptic_transform(cos(phi1), sin(phi1), da, db, dalpha, &tx,
465 G_debug(4,
" append point (r) x=%.16f y=%.16f", x[i] + tx,
467 phi1 += angular_step;
472 G_debug(4,
" append point (s) x=%.16f y=%.16f", nx, ny);
474 G_debug(4,
" append point (s) x=%.16f y=%.16f", mx, my);
488 int side,
int winding,
int stop_at_line_end,
502 double opt_angle, tangle;
503 int opt_j, opt_side, opt_flag;
505 G_debug(3,
"extract_contour(): v1=%d, v2=%d, side=%d, stop_at_line_end=%d",
506 first->
v1, first->
v2, side, stop_at_line_end);
521 vert0 = &(pg->
v[v0]);
523 eangle = atan2(vert->
y - vert0->
y, vert->
x - vert0->
x);
527 G_debug(4,
"ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d", v0,
528 v, eside, edge->
v1, edge->
v2);
529 G_debug(4,
"ec: append point x=%.18f y=%.18f", vert0->
x, vert0->
y);
542 for (j = 0; j < vert->
ecount; j++) {
544 if (vert->
edges[j] != edge) {
545 tangle = vert->
angles[j] - eangle;
548 else if (tangle >
PI)
552 if (opt_flag || (tangle < opt_angle)) {
554 opt_side = (vert->
edges[j]->
v1 == v) ? (1) : (-1);
568 if (stop_at_line_end) {
569 G_debug(3,
" end has been reached, will stop here");
575 G_debug(3,
" end has been reached, turning around");
580 if ((vert->
edges[opt_j] == first) && (opt_side == side))
584 G_warning(
_(
"Next edge was visited (right) but it is not the first one !!! breaking loop"));
586 "ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d",
587 v, (edge->
v1 == v) ? (edge->
v2) : (edge->
v1),
588 opt_side, vert->
edges[opt_j]->
v1,
595 G_warning(
_(
"Next edge was visited (left) but it is not the first one !!! breaking loop"));
597 "ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d",
598 v, (edge->
v1 == v) ? (edge->
v2) : (edge->
v1),
599 opt_side, vert->
edges[opt_j]->
v1,
605 edge = vert->
edges[opt_j];
608 v = (edge->
v1 == v) ? (edge->
v2) : (edge->
v1);
611 eangle = vert0->
angles[opt_j];
615 G_debug(4,
"ec: append point x=%.18f y=%.18f", vert->
x, vert->
y);
630 static void extract_outer_contour(
struct planar_graph *pg,
int side,
638 double min_x, min_angle;
640 G_debug(3,
"extract_outer_contour()");
649 for (i = 0; i < pg->
vcount; i++) {
650 if (flag || (pg->
v[i].
x < min_x)) {
659 for (i = 0; i < vert->
ecount; i++) {
660 if (flag || (vert->
angles[i] < min_angle)) {
661 edge = vert->
edges[i];
662 min_angle = vert->
angles[i];
681 static int extract_inner_contour(
struct planar_graph *pg,
int *winding,
687 G_debug(3,
"extract_inner_contour()");
689 for (i = 0; i < pg->
ecount; i++) {
694 extract_contour(pg, &(pg->
e[i]),
RIGHT_SIDE, w, 0, nPoints);
702 extract_contour(pg, &(pg->
e[i]),
LEFT_SIDE, w, 0, nPoints);
717 static int point_in_buf(
struct line_pnts *Points,
double px,
double py,
double da,
718 double db,
double dalpha)
722 double delta, delta_k, k;
723 double vx, vy, wx, wy, mx, my, nx, ny;
724 double len, tx, ty, d, da2;
732 for (i = 0; i < np - 1; i++) {
735 wx = Points->
x[i + 1];
736 wy = Points->
y[i + 1];
742 elliptic_tangent(mx / len, my / len, da, db, dalpha, &cx, &cy);
744 delta = mx * cy - my * cx;
745 delta_k = (px - vx) * cy - (py - vy) * cx;
762 elliptic_transform(px - nx, py - ny, 1 / da, 1 / db, dalpha, &tx,
790 static int get_polygon_orientation(
const double *
x,
const double *y,
int n)
792 double x1, y1, x2, y2;
806 area += (y2 + y1) * (x2 - x1);
813 static void add_line_to_array(
struct line_pnts *Points,
815 int *allocated,
int more)
817 if (*allocated == *
count) {
822 (*arrPoints)[*
count] = Points;
828 static void destroy_lines_array(
struct line_pnts **arr,
int count)
832 for (i = 0; i <
count; i++)
841 int isles_count,
int side,
double da,
double db,
842 double dalpha,
int round,
int caps,
double tol,
858 auto_side = (side == 0);
866 G_debug(3,
" processing outer contour");
870 get_polygon_orientation(area_outer->
x, area_outer->
y,
873 convolution_line(area_outer, da, db, dalpha, side, round, caps, tol,
876 extract_outer_contour(pg2, 0, *oPoints);
877 res = extract_inner_contour(pg2, &winding, cPoints);
884 if (area_size == 0) {
888 if (cPoints->
x[0] != cPoints->
x[cPoints->
n_points - 1] ||
889 cPoints->
y[0] != cPoints->
y[cPoints->
n_points - 1]) {
897 if (!point_in_buf(area_outer, px, py, da, db, dalpha)) {
898 add_line_to_array(cPoints, &arrPoints, &count, &allocated,
904 G_warning(
_(
"Vect_get_point_in_poly() failed"));
908 res = extract_inner_contour(pg2, &winding, cPoints);
913 G_debug(3,
" processing inner contours");
914 for (i = 0; i < isles_count; i++) {
917 get_polygon_orientation(area_isles[i]->
x, area_isles[i]->
y,
920 convolution_line(area_isles[i], da, db, dalpha, side, round, caps,
923 extract_outer_contour(pg2, 0, cPoints);
924 res = extract_inner_contour(pg2, &winding, cPoints);
931 if (area_size == 0) {
935 if (cPoints->
x[0] != cPoints->
x[cPoints->
n_points - 1] ||
936 cPoints->
y[0] != cPoints->
y[cPoints->
n_points - 1]) {
946 (cPoints->
x[0], cPoints->
y[0], area_isles[i])) {
948 if (!point_in_buf(area_isles[i], px, py, da, db, dalpha)) {
949 add_line_to_array(cPoints, &arrPoints, &count,
955 G_warning(
_(
"Vect_get_point_in_poly() failed"));
959 res = extract_inner_contour(pg2, &winding, cPoints);
965 *inner_count =
count;
966 *iPoints = arrPoints;
971 G_debug(3,
"buffer_lines() ... done");
1001 double dalpha,
int round,
int caps,
double tol,
1003 struct line_pnts ***iPoints,
int *inner_count)
1008 int isles_count = 0;
1011 int isles_allocated = 0;
1013 G_debug(2,
"Vect_line_buffer()");
1019 dalpha, round, tol, oPoints);
1028 extract_outer_contour(pg, 0, outer);
1031 res = extract_inner_contour(pg, &winding, tPoints);
1033 add_line_to_array(tPoints, &isles, &isles_count, &isles_allocated,
1036 res = extract_inner_contour(pg, &winding, tPoints);
1039 buffer_lines(outer, isles, isles_count,
RIGHT_SIDE, da, db, dalpha, round,
1040 caps, tol, oPoints, iPoints, inner_count);
1044 destroy_lines_array(isles, isles_count);
1066 double dalpha,
int round,
int caps,
double tol,
1068 struct line_pnts ***iPoints,
int *inner_count)
1072 int isles_count = 0, n_isles;
1075 int isles_allocated = 0;
1077 G_debug(2,
"Vect_area_buffer()");
1082 isles_allocated = n_isles;
1092 for (i = 0; i < n_isles; i++) {
1103 add_line_to_array(tPoints, &isles, &isles_count, &isles_allocated,
1108 buffer_lines(outer, isles, isles_count, 0, da, db, dalpha, round, caps,
1109 tol, oPoints, iPoints, inner_count);
1113 destroy_lines_array(isles, isles_count);
1131 double dalpha,
int round,
double tol,
1135 double angular_tol, angular_step, phi1;
1138 G_debug(2,
"Vect_point_buffer()");
1144 if (round || (!round)) {
1145 angular_tol = angular_tolerance(tol, da, db);
1147 nsegments = (int)(2 *
PI / angular_tol) + 1;
1148 angular_step = 2 *
PI / nsegments;
1151 for (j = 0; j < nsegments; j++) {
1152 elliptic_transform(cos(phi1), sin(phi1), da, db, dalpha, &tx,
1155 phi1 += angular_step;
1184 double dalpha,
int side,
int round,
double tol,
1187 G_debug(2,
"Vect_line_parallel(): npoints = %d, da = %f, " 1188 "db = %f, dalpha = %f, side = %d, round_corners = %d, tol = %f",
1189 InPoints->
n_points, da, db, dalpha, side, round, tol);
1191 parallel_line(InPoints, da, db, dalpha, side, round, 1,
NON_LOOPED_LINE,
int Vect_point_in_poly(double, double, const struct line_pnts *)
Determines if a point (X,Y) is inside a polygon.
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
void pg_destroy_struct(struct planar_graph *pg)
int Vect_get_area_isle(const struct Map_info *, int, int)
Returns isle id for area.
int n_points
Number of points.
int Vect_copy_xyz_to_pnts(struct line_pnts *, const double *, const double *, const double *, int)
Copy points from array to line_pnts structure.
void G_free(void *)
Free allocated memory.
int Vect_get_area_num_isles(const struct Map_info *, int)
Returns number of isles for given area.
int Vect_line_prune(struct line_pnts *)
Remove duplicate points, i.e. zero length segments.
double * x
Array of X coordinates.
Feature geometry info - coordinates.
struct planar_graph * pg_create(const struct line_pnts *Points)
int Vect_get_point_in_poly(const struct line_pnts *, double *, double *)
Get point inside polygon.
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
double dig_distance2_point_to_line(double, double, double, double, double, double, double, double, double, int, double *, double *, double *, double *, int *)
int Vect_get_isle_points(const struct Map_info *, int, struct line_pnts *)
Returns polygon array of points for given isle.
int Vect_append_point(struct line_pnts *, double, double, double)
Appends one point to the end of a line.
double * y
Array of Y coordinates.
void Vect_line_parallel2(struct line_pnts *InPoints, double da, double db, double dalpha, int side, int round, double tol, struct line_pnts *OutPoints)
void G_warning(const char *,...) __attribute__((format(printf
double * z
Array of Z coordinates.
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.
void Vect_destroy_line_struct(struct line_pnts *)
Frees all memory associated with a line_pnts structure, including the structure itself.
int Vect_get_area_points(const struct Map_info *, int, struct line_pnts *)
Returns polygon array of points (outer ring) of given area.
int G_debug(int, const char *,...) __attribute__((format(printf
void Vect_area_buffer2(const struct Map_info *Map, int area, double da, double db, double dalpha, int round, int caps, double tol, struct line_pnts **oPoints, struct line_pnts ***iPoints, int *inner_count)
Creates buffer around area.
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).
void Vect_reset_line(struct line_pnts *)
Reset line.
int Vect_line_delete_point(struct line_pnts *, int)
Delete point at given index and move all points above down.
int dig_find_area_poly(struct line_pnts *, double *)