00001
00021 #include <stdlib.h>
00022 #include <grass/Vect.h>
00023 #include <grass/gis.h>
00024 #include <grass/glocale.h>
00025
00026 #include "dgraph.h"
00027
00028 #define LENGTH(DX, DY) (sqrt((DX*DX)+(DY*DY)))
00029 #ifndef MIN
00030 #define MIN(X,Y) ((X<Y)?X:Y)
00031 #endif
00032 #ifndef MAX
00033 #define MAX(X,Y) ((X>Y)?X:Y)
00034 #endif
00035 #define PI M_PI
00036 #define RIGHT_SIDE 1
00037 #define LEFT_SIDE -1
00038 #define LOOPED_LINE 1
00039 #define NON_LOOPED_LINE 0
00040
00041
00042 static void norm_vector(double x1, double y1, double x2, double y2, double *x,
00043 double *y)
00044 {
00045 double dx, dy, l;
00046
00047 dx = x2 - x1;
00048 dy = y2 - y1;
00049 if ((dx == 0) && (dy == 0)) {
00050
00051
00052 *x = 0;
00053 *y = 0;
00054 return;
00055 }
00056 l = LENGTH(dx, dy);
00057 *x = dx / l;
00058 *y = dy / l;
00059 return;
00060 }
00061
00062 static void rotate_vector(double x, double y, double cosa, double sina,
00063 double *nx, double *ny)
00064 {
00065 *nx = x * cosa - y * sina;
00066 *ny = x * sina + y * cosa;
00067 return;
00068 }
00069
00070
00071
00072
00073
00074 static void elliptic_transform(double x, double y, double da, double db,
00075 double dalpha, double *nx, double *ny)
00076 {
00077 double cosa = cos(dalpha);
00078
00079 double sina = sin(dalpha);
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089 double va, vb;
00090
00091 va = (x * cosa + y * sina) * da;
00092 vb = (x * (-sina) + y * cosa) * db;
00093 *nx = va * cosa + vb * (-sina);
00094 *ny = va * sina + vb * cosa;
00095 return;
00096 }
00097
00098
00099
00100
00101
00102
00103
00104 static void elliptic_tangent(double x, double y, double da, double db,
00105 double dalpha, double *px, double *py)
00106 {
00107 double cosa = cos(dalpha);
00108
00109 double sina = sin(dalpha);
00110
00111 double u, v, len;
00112
00113
00114 rotate_vector(x, y, cosa, -sina, &x, &y);
00115
00116
00117 u = da * da * y;
00118 v = -db * db * x;
00119 len = da * db / sqrt(da * da * v * v + db * db * u * u);
00120 u *= len;
00121 v *= len;
00122 rotate_vector(u, v, cosa, sina, px, py);
00123 return;
00124 }
00125
00126
00127
00128
00129
00130 static void line_coefficients(double x1, double y1, double x2, double y2,
00131 double *a, double *b, double *c)
00132 {
00133 *a = y2 - y1;
00134 *b = x1 - x2;
00135 *c = x2 * y1 - x1 * y2;
00136 return;
00137 }
00138
00139
00140
00141
00142
00143
00144 static int line_intersection(double a1, double b1, double c1, double a2,
00145 double b2, double c2, double *x, double *y)
00146 {
00147 double d;
00148
00149 if (fabs(a2 * b1 - a1 * b2) == 0) {
00150 if (fabs(a2 * c1 - a1 * c2) == 0)
00151 return 2;
00152 else
00153 return 0;
00154 }
00155 else {
00156 d = a1 * b2 - a2 * b1;
00157 *x = (b1 * c2 - b2 * c1) / d;
00158 *y = (c1 * a2 - c2 * a1) / d;
00159 return 1;
00160 }
00161 }
00162
00163 static double angular_tolerance(double tol, double da, double db)
00164 {
00165 double a = MAX(da, db);
00166
00167 if (tol > a)
00168 tol = a;
00169 return 2 * acos(1 - tol / a);
00170 }
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184 static void parallel_line(struct line_pnts *Points, double da, double db,
00185 double dalpha, int side, int round, int caps, int looped,
00186 double tol, struct line_pnts *nPoints)
00187 {
00188 int i, j, res, np;
00189
00190 double *x, *y;
00191
00192 double tx, ty, vx, vy, wx, wy, nx, ny, mx, my, rx, ry;
00193
00194 double vx1, vy1, wx1, wy1;
00195
00196 double a0, b0, c0, a1, b1, c1;
00197
00198 double phi1, phi2, delta_phi;
00199
00200 double nsegments, angular_tol, angular_step;
00201
00202 int inner_corner, turns360;
00203
00204 G_debug(3, "parallel_line()");
00205
00206 if (looped && 0) {
00207
00208 return;
00209 }
00210
00211 Vect_reset_line(nPoints);
00212
00213 if (looped) {
00214 Vect_append_point(Points, Points->x[1], Points->y[1], Points->z[1]);
00215 }
00216 np = Points->n_points;
00217 x = Points->x;
00218 y = Points->y;
00219
00220 if ((np == 0) || (np == 1))
00221 return;
00222
00223 if ((da == 0) || (db == 0)) {
00224 Vect_copy_xyz_to_pnts(nPoints, x, y, NULL, np);
00225 return;
00226 }
00227
00228 side = (side >= 0) ? (1) : (-1);
00229 dalpha *= PI / 180;
00230 angular_tol = angular_tolerance(tol, da, db);
00231
00232 for (i = 0; i < np - 1; i++) {
00233
00234 a0 = a1;
00235 b0 = b1;
00236 c0 = c1;
00237 wx = vx;
00238 wy = vy;
00239
00240
00241 norm_vector(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);
00242 if ((tx == 0) && (ty == 0))
00243 continue;
00244
00245 elliptic_tangent(side * tx, side * ty, da, db, dalpha, &vx, &vy);
00246
00247 nx = x[i] + vx;
00248 ny = y[i] + vy;
00249
00250 mx = x[i + 1] + vx;
00251 my = y[i + 1] + vy;
00252
00253 line_coefficients(nx, ny, mx, my, &a1, &b1, &c1);
00254
00255 if (i == 0) {
00256 if (!looped)
00257 Vect_append_point(nPoints, nx, ny, 0);
00258 continue;
00259 }
00260
00261 delta_phi = atan2(ty, tx) - atan2(y[i] - y[i - 1], x[i] - x[i - 1]);
00262 if (delta_phi > PI)
00263 delta_phi -= 2 * PI;
00264 else if (delta_phi <= -PI)
00265 delta_phi += 2 * PI;
00266
00267 turns360 = (fabs(fabs(delta_phi) - PI) < 1e-15);
00268 inner_corner = (side * delta_phi <= 0) && (!turns360);
00269
00270 if ((turns360) && (!(caps && round))) {
00271 if (caps) {
00272 norm_vector(0, 0, vx, vy, &tx, &ty);
00273 elliptic_tangent(side * tx, side * ty, da, db, dalpha, &tx,
00274 &ty);
00275 }
00276 else {
00277 tx = 0;
00278 ty = 0;
00279 }
00280 Vect_append_point(nPoints, x[i] + wx + tx, y[i] + wy + ty, 0);
00281 Vect_append_point(nPoints, nx + tx, ny + ty, 0);
00282 }
00283 else if ((!round) || inner_corner) {
00284 res = line_intersection(a0, b0, c0, a1, b1, c1, &rx, &ry);
00285
00286
00287
00288
00289
00290 if (res == 1) {
00291 if (!round)
00292 Vect_append_point(nPoints, rx, ry, 0);
00293 else {
00294
00295
00296
00297 Vect_append_point(nPoints, rx, ry, 0);
00298 }
00299 }
00300 }
00301 else {
00302
00303
00304
00305 elliptic_transform(wx, wy, 1 / da, 1 / db, dalpha, &wx1, &wy1);
00306 elliptic_transform(vx, vy, 1 / da, 1 / db, dalpha, &vx1, &vy1);
00307
00308 phi1 = atan2(wy1, wx1);
00309 phi2 = atan2(vy1, vx1);
00310 delta_phi = side * (phi2 - phi1);
00311
00312
00313 if (delta_phi < 0)
00314 delta_phi += 2 * PI;
00315
00316 nsegments = (int)(delta_phi / angular_tol) + 1;
00317 angular_step = side * (delta_phi / nsegments);
00318
00319 for (j = 0; j <= nsegments; j++) {
00320 elliptic_transform(cos(phi1), sin(phi1), da, db, dalpha, &tx,
00321 &ty);
00322 Vect_append_point(nPoints, x[i] + tx, y[i] + ty, 0);
00323 phi1 += angular_step;
00324 }
00325 }
00326
00327 if ((!looped) && (i == np - 2)) {
00328 Vect_append_point(nPoints, mx, my, 0);
00329 }
00330 }
00331
00332 if (looped) {
00333 Vect_append_point(nPoints, nPoints->x[0], nPoints->y[0],
00334 nPoints->z[0]);
00335 }
00336
00337 Vect_line_prune(nPoints);
00338
00339 if (looped) {
00340 Vect_line_delete_point(Points, Points->n_points - 1);
00341 }
00342 }
00343
00344
00345 static void convolution_line(struct line_pnts *Points, double da, double db,
00346 double dalpha, int side, int round, int caps,
00347 double tol, struct line_pnts *nPoints)
00348 {
00349 int i, j, res, np;
00350
00351 double *x, *y;
00352
00353 double tx, ty, vx, vy, wx, wy, nx, ny, mx, my, rx, ry;
00354
00355 double vx1, vy1, wx1, wy1;
00356
00357 double a0, b0, c0, a1, b1, c1;
00358
00359 double phi1, phi2, delta_phi;
00360
00361 double nsegments, angular_tol, angular_step;
00362
00363 double angle0, angle1;
00364
00365 int inner_corner, turns360;
00366
00367 G_debug(3, "convolution_line() side = %d", side);
00368
00369 np = Points->n_points;
00370 x = Points->x;
00371 y = Points->y;
00372 if ((np == 0) || (np == 1))
00373 return;
00374 if ((x[0] != x[np - 1]) || (y[0] != y[np - 1])) {
00375 G_fatal_error(_("Line is not looped"));
00376 return;
00377 }
00378
00379 Vect_reset_line(nPoints);
00380
00381 if ((da == 0) || (db == 0)) {
00382 Vect_copy_xyz_to_pnts(nPoints, x, y, NULL, np);
00383 return;
00384 }
00385
00386 side = (side >= 0) ? (1) : (-1);
00387 dalpha *= PI / 180;
00388 angular_tol = angular_tolerance(tol, da, db);
00389
00390 i = np - 2;
00391 norm_vector(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);
00392 elliptic_tangent(side * tx, side * ty, da, db, dalpha, &vx, &vy);
00393 angle1 = atan2(ty, tx);
00394 nx = x[i] + vx;
00395 ny = y[i] + vy;
00396 mx = x[i + 1] + vx;
00397 my = y[i + 1] + vy;
00398 if (!round)
00399 line_coefficients(nx, ny, mx, my, &a1, &b1, &c1);
00400
00401 for (i = 0; i <= np - 2; i++) {
00402 G_debug(4, "point %d, segment %d-%d", i, i, i + 1);
00403
00404 if (!round) {
00405 a0 = a1;
00406 b0 = b1;
00407 c0 = c1;
00408 }
00409 wx = vx;
00410 wy = vy;
00411 angle0 = angle1;
00412
00413 norm_vector(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);
00414 if ((tx == 0) && (ty == 0))
00415 continue;
00416 elliptic_tangent(side * tx, side * ty, da, db, dalpha, &vx, &vy);
00417 angle1 = atan2(ty, tx);
00418 nx = x[i] + vx;
00419 ny = y[i] + vy;
00420 mx = x[i + 1] + vx;
00421 my = y[i + 1] + vy;
00422 if (!round)
00423 line_coefficients(nx, ny, mx, my, &a1, &b1, &c1);
00424
00425
00426 delta_phi = angle1 - angle0;
00427 if (delta_phi > PI)
00428 delta_phi -= 2 * PI;
00429 else if (delta_phi <= -PI)
00430 delta_phi += 2 * PI;
00431
00432 turns360 = (fabs(fabs(delta_phi) - PI) < 1e-15);
00433 inner_corner = (side * delta_phi <= 0) && (!turns360);
00434
00435
00436
00437 if (turns360 && caps && (!round)) {
00438 norm_vector(0, 0, vx, vy, &tx, &ty);
00439 elliptic_tangent(side * tx, side * ty, da, db, dalpha, &tx, &ty);
00440 Vect_append_point(nPoints, x[i] + wx + tx, y[i] + wy + ty, 0);
00441 G_debug(4, " append point (c) x=%.16f y=%.16f", x[i] + wx + tx,
00442 y[i] + wy + ty);
00443 Vect_append_point(nPoints, nx + tx, ny + ty, 0);
00444 G_debug(4, " append point (c) x=%.16f y=%.16f", nx + tx, ny + ty);
00445 }
00446
00447 if ((!turns360) && (!round) && (!inner_corner)) {
00448 res = line_intersection(a0, b0, c0, a1, b1, c1, &rx, &ry);
00449 if (res == 1) {
00450 Vect_append_point(nPoints, rx, ry, 0);
00451 G_debug(4, " append point (o) x=%.16f y=%.16f", rx, ry);
00452 }
00453 else if (res == 2) {
00454
00455 }
00456 else
00457 G_fatal_error
00458 ("unexpected result of line_intersection() res = %d",
00459 res);
00460 }
00461
00462 if (round && (!inner_corner) && (!turns360 || caps)) {
00463
00464
00465
00466 elliptic_transform(wx, wy, 1 / da, 1 / db, dalpha, &wx1, &wy1);
00467 elliptic_transform(vx, vy, 1 / da, 1 / db, dalpha, &vx1, &vy1);
00468
00469 phi1 = atan2(wy1, wx1);
00470 phi2 = atan2(vy1, vx1);
00471 delta_phi = side * (phi2 - phi1);
00472
00473
00474 if (delta_phi < 0)
00475 delta_phi += 2 * PI;
00476
00477 nsegments = (int)(delta_phi / angular_tol) + 1;
00478 angular_step = side * (delta_phi / nsegments);
00479
00480 phi1 += angular_step;
00481 for (j = 1; j <= nsegments - 1; j++) {
00482 elliptic_transform(cos(phi1), sin(phi1), da, db, dalpha, &tx,
00483 &ty);
00484 Vect_append_point(nPoints, x[i] + tx, y[i] + ty, 0);
00485 G_debug(4, " append point (r) x=%.16f y=%.16f", x[i] + tx,
00486 y[i] + ty);
00487 phi1 += angular_step;
00488 }
00489 }
00490
00491 Vect_append_point(nPoints, nx, ny, 0);
00492 G_debug(4, " append point (s) x=%.16f y=%.16f", nx, ny);
00493 Vect_append_point(nPoints, mx, my, 0);
00494 G_debug(4, " append point (s) x=%.16f y=%.16f", mx, my);
00495 }
00496
00497
00498 Vect_append_point(nPoints, nPoints->x[0], nPoints->y[0], nPoints->z[0]);
00499
00500 }
00501
00502
00503
00504
00505
00506
00507 static void extract_contour(struct planar_graph *pg, struct pg_edge *first,
00508 int side, int winding, int stop_at_line_end,
00509 struct line_pnts *nPoints)
00510 {
00511 int j;
00512
00513 int v;
00514
00515 int v0;
00516
00517 int eside;
00518
00519 double eangle;
00520
00521 struct pg_vertex *vert;
00522
00523 struct pg_vertex *vert0;
00524
00525 struct pg_edge *edge;
00526
00527
00528
00529 double opt_angle, tangle;
00530
00531 int opt_j, opt_side, opt_flag;
00532
00533 G_debug(3,
00534 "extract_contour(): v1=%d, v2=%d, side=%d, stop_at_line_end=%d",
00535 first->v1, first->v2, side, stop_at_line_end);
00536
00537 Vect_reset_line(nPoints);
00538
00539 edge = first;
00540 if (side >= 0) {
00541 eside = 1;
00542 v0 = edge->v1;
00543 v = edge->v2;
00544 }
00545 else {
00546 eside = -1;
00547 v0 = edge->v2;
00548 v = edge->v1;
00549 }
00550 vert0 = &(pg->v[v0]);
00551 vert = &(pg->v[v]);
00552 eangle = atan2(vert->y - vert0->y, vert->x - vert0->x);
00553
00554 while (1) {
00555 Vect_append_point(nPoints, vert0->x, vert0->y, 0);
00556 G_debug(4, "ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d", v0,
00557 v, eside, edge->v1, edge->v2);
00558 G_debug(4, "ec: append point x=%.18f y=%.18f", vert0->x, vert0->y);
00559
00560
00561 if (eside == 1) {
00562 edge->visited_right = 1;
00563 edge->winding_right = winding;
00564 }
00565 else {
00566 edge->visited_left = 1;
00567 edge->winding_left = winding;
00568 }
00569
00570 opt_flag = 1;
00571 for (j = 0; j < vert->ecount; j++) {
00572
00573 if (vert->edges[j] != edge) {
00574 tangle = vert->angles[j] - eangle;
00575 if (tangle < -PI)
00576 tangle += 2 * PI;
00577 else if (tangle > PI)
00578 tangle -= 2 * PI;
00579
00580
00581 if (opt_flag || (tangle < opt_angle)) {
00582 opt_j = j;
00583 opt_side = (vert->edges[j]->v1 == v) ? (1) : (-1);
00584 opt_angle = tangle;
00585 opt_flag = 0;
00586 }
00587 }
00588 }
00589
00590
00591
00592
00593 if (opt_flag) {
00594 if (stop_at_line_end) {
00595 G_debug(3, " end has been reached, will stop here");
00596 break;
00597 }
00598 else {
00599 opt_j = 0;
00600 opt_side = -eside;
00601 G_debug(3, " end has been reached, turning around");
00602 }
00603 }
00604
00605
00606 if ((vert->edges[opt_j] == first) && (opt_side == side))
00607 break;
00608 if (opt_side == 1) {
00609 if (vert->edges[opt_j]->visited_right) {
00610 G_warning(_("Next edge was visited but it is not the first one !!! breaking loop"));
00611 G_debug(4,
00612 "ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d",
00613 v, (edge->v1 == v) ? (edge->v2) : (edge->v1),
00614 opt_side, vert->edges[opt_j]->v1,
00615 vert->edges[opt_j]->v2);
00616 break;
00617 }
00618 }
00619 else {
00620 if (vert->edges[opt_j]->visited_left) {
00621 G_warning(_("Next edge was visited but it is not the first one !!! breaking loop"));
00622 G_debug(4,
00623 "ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d",
00624 v, (edge->v1 == v) ? (edge->v2) : (edge->v1),
00625 opt_side, vert->edges[opt_j]->v1,
00626 vert->edges[opt_j]->v2);
00627 break;
00628 }
00629 }
00630
00631 edge = vert->edges[opt_j];
00632 eside = opt_side;
00633 v0 = v;
00634 v = (edge->v1 == v) ? (edge->v2) : (edge->v1);
00635 vert0 = vert;
00636 vert = &(pg->v[v]);
00637 eangle = vert0->angles[opt_j];
00638 }
00639 Vect_append_point(nPoints, vert->x, vert->y, 0);
00640 G_debug(4, "ec: append point x=%.18f y=%.18f", vert->x, vert->y);
00641
00642 return;
00643 }
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655 static void extract_outer_contour(struct planar_graph *pg, int side,
00656 struct line_pnts *nPoints)
00657 {
00658 int i;
00659
00660 int flag;
00661
00662 int v;
00663
00664 struct pg_vertex *vert;
00665
00666 struct pg_edge *edge;
00667
00668 double min_x, min_angle;
00669
00670 G_debug(3, "extract_outer_contour()");
00671
00672 if (side != 0) {
00673 G_fatal_error(_("side != 0 feature not implemented"));
00674 return;
00675 }
00676
00677
00678 flag = 1;
00679 for (i = 0; i < pg->vcount; i++) {
00680 if (flag || (pg->v[i].x < min_x)) {
00681 v = i;
00682 min_x = pg->v[i].x;
00683 flag = 0;
00684 }
00685 }
00686 vert = &(pg->v[v]);
00687
00688 flag = 1;
00689 for (i = 0; i < vert->ecount; i++) {
00690 if (flag || (vert->angles[i] < min_angle)) {
00691 edge = vert->edges[i];
00692 min_angle = vert->angles[i];
00693 flag = 0;
00694 }
00695 }
00696
00697
00698 extract_contour(pg, edge, (edge->v1 == v) ? RIGHT_SIDE : LEFT_SIDE, 0, 0,
00699 nPoints);
00700
00701 return;
00702 }
00703
00704
00705
00706
00707
00708
00709
00710
00711 static int extract_inner_contour(struct planar_graph *pg, int *winding,
00712 struct line_pnts *nPoints)
00713 {
00714 int i, w;
00715
00716 struct pg_edge *edge;
00717
00718 G_debug(3, "extract_inner_contour()");
00719
00720 for (i = 0; i < pg->ecount; i++) {
00721 edge = &(pg->e[i]);
00722 if (edge->visited_left) {
00723 if (!(pg->e[i].visited_right)) {
00724 w = edge->winding_left - 1;
00725 extract_contour(pg, &(pg->e[i]), RIGHT_SIDE, w, 0, nPoints);
00726 *winding = w;
00727 return 1;
00728 }
00729 }
00730 else {
00731 if (pg->e[i].visited_right) {
00732 w = edge->winding_right + 1;
00733 extract_contour(pg, &(pg->e[i]), LEFT_SIDE, w, 0, nPoints);
00734 *winding = w;
00735 return 1;
00736 }
00737 }
00738 }
00739
00740 return 0;
00741 }
00742
00743
00744
00745
00746
00747
00748 static int point_in_buf(struct line_pnts *Points, double px, double py, double da,
00749 double db, double dalpha)
00750 {
00751 int i, np;
00752
00753 double cx, cy;
00754
00755 double delta, delta_k, k;
00756
00757 double vx, vy, wx, wy, mx, my, nx, ny;
00758
00759 double len, tx, ty, d, da2;
00760
00761 G_debug(3, "point_in_buf()");
00762
00763 dalpha *= PI / 180;
00764
00765 np = Points->n_points;
00766 da2 = da * da;
00767 for (i = 0; i < np - 1; i++) {
00768 vx = Points->x[i];
00769 vy = Points->y[i];
00770 wx = Points->x[i + 1];
00771 wy = Points->y[i + 1];
00772
00773 if (da != db) {
00774 mx = wx - vx;
00775 my = wy - vy;
00776 len = LENGTH(mx, my);
00777 elliptic_tangent(mx / len, my / len, da, db, dalpha, &cx, &cy);
00778
00779 delta = mx * cy - my * cx;
00780 delta_k = (px - vx) * cy - (py - vy) * cx;
00781 k = delta_k / delta;
00782
00783 if (k <= 0) {
00784 nx = vx;
00785 ny = vy;
00786 }
00787 else if (k >= 1) {
00788 nx = wx;
00789 ny = wy;
00790 }
00791 else {
00792 nx = vx + k * mx;
00793 ny = vy + k * my;
00794 }
00795
00796
00797 elliptic_transform(px - nx, py - ny, 1 / da, 1 / db, dalpha, &tx,
00798 &ty);
00799
00800 d = dig_distance2_point_to_line(nx + tx, ny + ty, 0, vx, vy, 0,
00801 wx, wy, 0, 0, NULL, NULL, NULL,
00802 NULL, NULL);
00803
00804
00805 if (d <= 1) {
00806
00807 return 1;
00808 }
00809 }
00810 else {
00811 d = dig_distance2_point_to_line(px, py, 0, vx, vy, 0, wx, wy, 0,
00812 0, NULL, NULL, NULL, NULL, NULL);
00813
00814 if (d <= da2) {
00815 return 1;
00816 }
00817 }
00818 }
00819 return 0;
00820 }
00821
00822
00823
00824 static int get_polygon_orientation(const double *x, const double *y, int n)
00825 {
00826 double x1, y1, x2, y2;
00827
00828 double area;
00829
00830 x2 = x[n - 1];
00831 y2 = y[n - 1];
00832
00833 area = 0;
00834 while (--n >= 0) {
00835 x1 = x2;
00836 y1 = y2;
00837
00838 x2 = *x++;
00839 y2 = *y++;
00840
00841 area += (y2 + y1) * (x2 - x1);
00842 }
00843 return (area > 0);
00844 }
00845
00846
00847 static void add_line_to_array(struct line_pnts *Points,
00848 struct line_pnts ***arrPoints, int *count,
00849 int *allocated, int more)
00850 {
00851 if (*allocated == *count) {
00852 *allocated += more;
00853 *arrPoints =
00854 G_realloc(*arrPoints, (*allocated) * sizeof(struct line_pnts *));
00855 }
00856 (*arrPoints)[*count] = Points;
00857 (*count)++;
00858 return;
00859 }
00860
00861 static void destroy_lines_array(struct line_pnts **arr, int count)
00862 {
00863 int i;
00864
00865 for (i = 0; i < count; i++)
00866 Vect_destroy_line_struct(arr[i]);
00867 G_free(arr);
00868 }
00869
00870
00871
00872
00873 static void buffer_lines(struct line_pnts *area_outer, struct line_pnts **area_isles,
00874 int isles_count, int side, double da, double db,
00875 double dalpha, int round, int caps, double tol,
00876 struct line_pnts **oPoints, struct line_pnts ***iPoints,
00877 int *inner_count)
00878 {
00879 struct planar_graph *pg2;
00880
00881 struct line_pnts *sPoints, *cPoints;
00882
00883 struct line_pnts **arrPoints;
00884
00885 int i, count = 0;
00886
00887 int res, winding;
00888
00889 int auto_side;
00890
00891 int more = 8;
00892
00893 int allocated = 0;
00894
00895 double px, py;
00896
00897 G_debug(3, "buffer_lines()");
00898
00899 auto_side = (side == 0);
00900
00901
00902 sPoints = Vect_new_line_struct();
00903 cPoints = Vect_new_line_struct();
00904 arrPoints = NULL;
00905
00906
00907 G_debug(3, " processing outer contour");
00908 *oPoints = Vect_new_line_struct();
00909 if (auto_side)
00910 side =
00911 get_polygon_orientation(area_outer->x, area_outer->y,
00912 area_outer->n_points -
00913 1) ? LEFT_SIDE : RIGHT_SIDE;
00914 convolution_line(area_outer, da, db, dalpha, side, round, caps, tol,
00915 sPoints);
00916 pg2 = pg_create(sPoints);
00917 extract_outer_contour(pg2, 0, *oPoints);
00918 res = extract_inner_contour(pg2, &winding, cPoints);
00919 while (res != 0) {
00920 if (winding == 0) {
00921 if (!Vect_point_in_poly(cPoints->x[0], cPoints->y[0], area_outer)) {
00922 if (Vect_get_point_in_poly(cPoints, &px, &py) != 0)
00923 G_fatal_error(_("Vect_get_point_in_poly() failed"));
00924 if (!point_in_buf(area_outer, px, py, da, db, dalpha)) {
00925 add_line_to_array(cPoints, &arrPoints, &count, &allocated,
00926 more);
00927 cPoints = Vect_new_line_struct();
00928 }
00929 }
00930 }
00931 res = extract_inner_contour(pg2, &winding, cPoints);
00932 }
00933 pg_destroy_struct(pg2);
00934
00935
00936 G_debug(3, " processing inner contours");
00937 for (i = 0; i < isles_count; i++) {
00938 if (auto_side)
00939 side =
00940 get_polygon_orientation(area_isles[i]->x, area_isles[i]->y,
00941 area_isles[i]->n_points -
00942 1) ? RIGHT_SIDE : LEFT_SIDE;
00943 convolution_line(area_isles[i], da, db, dalpha, side, round, caps,
00944 tol, sPoints);
00945 pg2 = pg_create(sPoints);
00946 extract_outer_contour(pg2, 0, cPoints);
00947 res = extract_inner_contour(pg2, &winding, cPoints);
00948 while (res != 0) {
00949 if (winding == -1) {
00950
00951
00952
00953 if (Vect_point_in_poly
00954 (cPoints->x[0], cPoints->y[0], area_isles[i])) {
00955 if (Vect_get_point_in_poly(cPoints, &px, &py) != 0)
00956 G_fatal_error(_("Vect_get_point_in_poly() failed"));
00957 if (!point_in_buf(area_isles[i], px, py, da, db, dalpha)) {
00958 add_line_to_array(cPoints, &arrPoints, &count,
00959 &allocated, more);
00960 cPoints = Vect_new_line_struct();
00961 }
00962 }
00963 }
00964 res = extract_inner_contour(pg2, &winding, cPoints);
00965 }
00966 pg_destroy_struct(pg2);
00967 }
00968
00969 arrPoints = G_realloc(arrPoints, count * sizeof(struct line_pnts *));
00970 *inner_count = count;
00971 *iPoints = arrPoints;
00972
00973 Vect_destroy_line_struct(sPoints);
00974 Vect_destroy_line_struct(cPoints);
00975
00976 G_debug(3, "buffer_lines() ... done");
00977
00978 return;
00979 }
00980
00981
00998 void Vect_line_buffer2(struct line_pnts *Points, double da, double db,
00999 double dalpha, int round, int caps, double tol,
01000 struct line_pnts **oPoints,
01001 struct line_pnts ***iPoints, int *inner_count)
01002 {
01003 struct planar_graph *pg;
01004
01005 struct line_pnts *tPoints, *outer;
01006
01007 struct line_pnts **isles;
01008
01009 int isles_count = 0;
01010
01011 int res, winding;
01012
01013 int more = 8;
01014
01015 int isles_allocated = 0;
01016
01017 G_debug(2, "Vect_line_buffer()");
01018
01019
01020 tPoints = Vect_new_line_struct();
01021 isles = NULL;
01022 pg = pg_create(Points);
01023
01024
01025 outer = Vect_new_line_struct();
01026 extract_outer_contour(pg, 0, outer);
01027
01028
01029 res = extract_inner_contour(pg, &winding, tPoints);
01030 while (res != 0) {
01031 add_line_to_array(tPoints, &isles, &isles_count, &isles_allocated,
01032 more);
01033 tPoints = Vect_new_line_struct();
01034 res = extract_inner_contour(pg, &winding, tPoints);
01035 }
01036
01037 buffer_lines(outer, isles, isles_count, RIGHT_SIDE, da, db, dalpha, round,
01038 caps, tol, oPoints, iPoints, inner_count);
01039
01040 Vect_destroy_line_struct(tPoints);
01041 Vect_destroy_line_struct(outer);
01042 destroy_lines_array(isles, isles_count);
01043 pg_destroy_struct(pg);
01044
01045 return;
01046 }
01047
01063 void Vect_area_buffer2(struct Map_info *Map, int area, double da, double db,
01064 double dalpha, int round, int caps, double tol,
01065 struct line_pnts **oPoints,
01066 struct line_pnts ***iPoints, int *inner_count)
01067 {
01068 struct line_pnts *tPoints, *outer;
01069
01070 struct line_pnts **isles;
01071
01072 int isles_count = 0;
01073
01074 int i, isle;
01075
01076 int more = 8;
01077
01078 int isles_allocated = 0;
01079
01080 G_debug(2, "Vect_area_buffer()");
01081
01082
01083 tPoints = Vect_new_line_struct();
01084 isles_count = Vect_get_area_num_isles(Map, area);
01085 isles_allocated = isles_count;
01086 isles = G_malloc(isles_allocated * sizeof(struct line_pnts *));
01087
01088
01089 outer = Vect_new_line_struct();
01090 Vect_get_area_points(Map, area, outer);
01091 Vect_append_point(outer, outer->x[0], outer->y[0], outer->z[0]);
01092
01093
01094 for (i = 0; i < isles_count; i++) {
01095 isle = Vect_get_area_isle(Map, area, i);
01096 Vect_get_isle_points(Map, isle, tPoints);
01097
01098
01099
01100
01101
01102
01103 Vect_append_point(tPoints, tPoints->x[0], tPoints->y[0],
01104 tPoints->z[0]);
01105 add_line_to_array(tPoints, &isles, &isles_count, &isles_allocated,
01106 more);
01107 tPoints = Vect_new_line_struct();
01108 }
01109
01110 buffer_lines(outer, isles, isles_count, 0, da, db, dalpha, round, caps,
01111 tol, oPoints, iPoints, inner_count);
01112
01113 Vect_destroy_line_struct(tPoints);
01114 Vect_destroy_line_struct(outer);
01115 destroy_lines_array(isles, isles_count);
01116
01117 return;
01118 }
01119
01132 void Vect_point_buffer2(double px, double py, double da, double db,
01133 double dalpha, int round, double tol,
01134 struct line_pnts **oPoints)
01135 {
01136 double tx, ty;
01137
01138 double angular_tol, angular_step, phi1;
01139
01140 int j, nsegments;
01141
01142 G_debug(2, "Vect_point_buffer()");
01143
01144 *oPoints = Vect_new_line_struct();
01145
01146 dalpha *= PI / 180;
01147
01148 if (round || (!round)) {
01149 angular_tol = angular_tolerance(tol, da, db);
01150
01151 nsegments = (int)(2 * PI / angular_tol) + 1;
01152 angular_step = 2 * PI / nsegments;
01153
01154 phi1 = 0;
01155 for (j = 0; j < nsegments; j++) {
01156 elliptic_transform(cos(phi1), sin(phi1), da, db, dalpha, &tx,
01157 &ty);
01158 Vect_append_point(*oPoints, px + tx, py + ty, 0);
01159 phi1 += angular_step;
01160 }
01161 }
01162 else {
01163
01164 }
01165
01166
01167 Vect_append_point(*oPoints, (*oPoints)->x[0], (*oPoints)->y[0],
01168 (*oPoints)->z[0]);
01169
01170 return;
01171 }
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187 void Vect_line_parallel2(struct line_pnts *InPoints, double da, double db,
01188 double dalpha, int side, int round, double tol,
01189 struct line_pnts *OutPoints)
01190 {
01191 G_debug(2, "Vect_line_parallel(): npoints = %d, da = %f, "
01192 "db = %f, dalpha = %f, side = %d, round_corners = %d, tol = %f",
01193 InPoints->n_points, da, db, dalpha, side, round, tol);
01194
01195 parallel_line(InPoints, da, db, dalpha, side, round, 1, NON_LOOPED_LINE,
01196 tol, OutPoints);
01197
01198
01199
01200
01201
01202 return;
01203 }