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