GRASS Programmer's Manual 6.4.1(2011)
|
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 /* norm_vector() calculates normalized vector form two points */ 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 /* assume that dx == dy == 0, which should give (NaN,NaN) */ 00052 /* without this, very small dx or dy could result in Infinity */ 00053 *x = 0; 00054 *y = 0; 00055 return; 00056 } 00057 l = LENGTH(dx, dy); 00058 *x = dx / l; 00059 *y = dy / l; 00060 00061 return; 00062 } 00063 00064 static void rotate_vector(double x, double y, double cosa, double sina, 00065 double *nx, double *ny) 00066 { 00067 *nx = x * cosa - y * sina; 00068 *ny = x * sina + y * cosa; 00069 00070 return; 00071 } 00072 00073 /* 00074 * (x,y) shoud be normalized vector for common transforms; This func transforms (x,y) to a vector corresponding to da, db, dalpha params 00075 * dalpha is in radians 00076 */ 00077 static void elliptic_transform(double x, double y, double da, double db, 00078 double dalpha, double *nx, double *ny) 00079 { 00080 double cosa = cos(dalpha); 00081 double sina = sin(dalpha); 00082 00083 /* double cc = cosa*cosa; 00084 double ss = sina*sina; 00085 double t = (da-db)*sina*cosa; 00086 00087 *nx = (da*cc + db*ss)*x + t*y; 00088 *ny = (da*ss + db*cc)*y + t*x; 00089 return; */ 00090 00091 double va, vb; 00092 00093 va = (x * cosa + y * sina) * da; 00094 vb = (x * (-sina) + y * cosa) * db; 00095 *nx = va * cosa + vb * (-sina); 00096 *ny = va * sina + vb * cosa; 00097 00098 return; 00099 } 00100 00101 /* 00102 * vect(x,y) must be normalized 00103 * gives the tangent point of the tangent to ellpise(da,db,dalpha) parallel to vect(x,y) 00104 * dalpha is in radians 00105 * ellipse center is in (0,0) 00106 */ 00107 static void elliptic_tangent(double x, double y, double da, double db, 00108 double dalpha, double *px, double *py) 00109 { 00110 double cosa = cos(dalpha); 00111 double sina = sin(dalpha); 00112 double u, v, len; 00113 00114 /* rotate (x,y) -dalpha radians */ 00115 rotate_vector(x, y, cosa, -sina, &x, &y); 00116 /*u = (x + da*y/db)/2; 00117 v = (y - db*x/da)/2; */ 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 00125 return; 00126 } 00127 00128 00129 /* 00130 * !!! This is not line in GRASS' sense. See http://en.wikipedia.org/wiki/Line_%28mathematics%29 00131 */ 00132 static void line_coefficients(double x1, double y1, double x2, double y2, 00133 double *a, double *b, double *c) 00134 { 00135 *a = y2 - y1; 00136 *b = x1 - x2; 00137 *c = x2 * y1 - x1 * y2; 00138 00139 return; 00140 } 00141 00142 /* 00143 * Finds intersection of two straight lines. Returns 0 if the lines are parallel, 1 if they cross, 00144 * 2 if they are the same line. 00145 * !!!!!!!!!!!!!!!! FIX THIS TOLLERANCE CONSTANTS BAD (and UGLY) CODE !!!!!!!!! 00146 */ 00147 static int line_intersection(double a1, double b1, double c1, double a2, 00148 double b2, double c2, double *x, double *y) 00149 { 00150 double d; 00151 00152 if (fabs(a2 * b1 - a1 * b2) == 0) { 00153 if (fabs(a2 * c1 - a1 * c2) == 0) 00154 return 2; 00155 else 00156 return 0; 00157 } 00158 else { 00159 d = a1 * b2 - a2 * b1; 00160 *x = (b1 * c2 - b2 * c1) / d; 00161 *y = (c1 * a2 - c2 * a1) / d; 00162 return 1; 00163 } 00164 } 00165 00166 static double angular_tolerance(double tol, double da, double db) 00167 { 00168 double a = MAX(da, db); 00169 00170 if (tol > a) 00171 tol = a; 00172 00173 return 2 * acos(1 - tol / a); 00174 } 00175 00176 /* 00177 * This function generates parallel line (with loops, but not like the old ones). 00178 * It is not to be used directly for creating buffers. 00179 * + added elliptical buffers/par.lines support 00180 * 00181 * dalpha - direction of elliptical buffer major axis in degrees 00182 * da - distance along major axis 00183 * db: distance along minor (perp.) axis 00184 * side: side >= 0 - right side, side < 0 - left side 00185 * when (da == db) we have plain distances (old case) 00186 * round - 1 for round corners, 0 for sharp corners. (tol is used only if round == 1) 00187 */ 00188 static void parallel_line(struct line_pnts *Points, double da, double db, 00189 double dalpha, int side, int round, int caps, int looped, 00190 double tol, struct line_pnts *nPoints) 00191 { 00192 int i, j, res, np; 00193 double *x, *y; 00194 double tx, ty, vx, vy, wx, wy, nx, ny, mx, my, rx, ry; 00195 double vx1, vy1, wx1, wy1; 00196 double a0, b0, c0, a1, b1, c1; 00197 double phi1, phi2, delta_phi; 00198 double nsegments, angular_tol, angular_step; 00199 int inner_corner, turns360; 00200 00201 G_debug(3, "parallel_line()"); 00202 00203 if (looped && 0) { 00204 /* start point != end point */ 00205 return; 00206 } 00207 00208 Vect_reset_line(nPoints); 00209 00210 if (looped) { 00211 Vect_append_point(Points, Points->x[1], Points->y[1], Points->z[1]); 00212 } 00213 np = Points->n_points; 00214 x = Points->x; 00215 y = Points->y; 00216 00217 if ((np == 0) || (np == 1)) 00218 return; 00219 00220 if ((da == 0) || (db == 0)) { 00221 Vect_copy_xyz_to_pnts(nPoints, x, y, NULL, np); 00222 return; 00223 } 00224 00225 side = (side >= 0) ? (1) : (-1); /* normalize variable */ 00226 dalpha *= PI / 180; /* convert dalpha from degrees to radians */ 00227 angular_tol = angular_tolerance(tol, da, db); 00228 00229 for (i = 0; i < np - 1; i++) { 00230 /* save the old values */ 00231 a0 = a1; 00232 b0 = b1; 00233 c0 = c1; 00234 wx = vx; 00235 wy = vy; 00236 00237 00238 norm_vector(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty); 00239 if ((tx == 0) && (ty == 0)) 00240 continue; 00241 00242 elliptic_tangent(side * tx, side * ty, da, db, dalpha, &vx, &vy); 00243 00244 nx = x[i] + vx; 00245 ny = y[i] + vy; 00246 00247 mx = x[i + 1] + vx; 00248 my = y[i + 1] + vy; 00249 00250 line_coefficients(nx, ny, mx, my, &a1, &b1, &c1); 00251 00252 if (i == 0) { 00253 if (!looped) 00254 Vect_append_point(nPoints, nx, ny, 0); 00255 continue; 00256 } 00257 00258 delta_phi = atan2(ty, tx) - atan2(y[i] - y[i - 1], x[i] - x[i - 1]); 00259 if (delta_phi > PI) 00260 delta_phi -= 2 * PI; 00261 else if (delta_phi <= -PI) 00262 delta_phi += 2 * PI; 00263 /* now delta_phi is in [-pi;pi] */ 00264 turns360 = (fabs(fabs(delta_phi) - PI) < 1e-15); 00265 inner_corner = (side * delta_phi <= 0) && (!turns360); 00266 00267 if ((turns360) && (!(caps && round))) { 00268 if (caps) { 00269 norm_vector(0, 0, vx, vy, &tx, &ty); 00270 elliptic_tangent(side * tx, side * ty, da, db, dalpha, &tx, 00271 &ty); 00272 } 00273 else { 00274 tx = 0; 00275 ty = 0; 00276 } 00277 Vect_append_point(nPoints, x[i] + wx + tx, y[i] + wy + ty, 0); 00278 Vect_append_point(nPoints, nx + tx, ny + ty, 0); /* nx == x[i] + vx, ny == y[i] + vy */ 00279 } 00280 else if ((!round) || inner_corner) { 00281 res = line_intersection(a0, b0, c0, a1, b1, c1, &rx, &ry); 00282 /* if (res == 0) { 00283 G_debug(4, "a0=%.18f, b0=%.18f, c0=%.18f, a1=%.18f, b1=%.18f, c1=%.18f", a0, b0, c0, a1, b1, c1); 00284 G_fatal_error("Two consequtive line segments are parallel, but not on one straight line! This should never happen."); 00285 return; 00286 } */ 00287 if (res == 1) { 00288 if (!round) 00289 Vect_append_point(nPoints, rx, ry, 0); 00290 else { 00291 /* d = dig_distance2_point_to_line(rx, ry, 0, x[i-1], y[i-1], 0, x[i], y[i], 0, 00292 0, NULL, NULL, NULL, NULL, NULL); 00293 if ( */ 00294 Vect_append_point(nPoints, rx, ry, 0); 00295 } 00296 } 00297 } 00298 else { 00299 /* we should draw elliptical arc for outside corner */ 00300 00301 /* inverse transforms */ 00302 elliptic_transform(wx, wy, 1 / da, 1 / db, dalpha, &wx1, &wy1); 00303 elliptic_transform(vx, vy, 1 / da, 1 / db, dalpha, &vx1, &vy1); 00304 00305 phi1 = atan2(wy1, wx1); 00306 phi2 = atan2(vy1, vx1); 00307 delta_phi = side * (phi2 - phi1); 00308 00309 /* make delta_phi in [0, 2pi] */ 00310 if (delta_phi < 0) 00311 delta_phi += 2 * PI; 00312 00313 nsegments = (int)(delta_phi / angular_tol) + 1; 00314 angular_step = side * (delta_phi / nsegments); 00315 00316 for (j = 0; j <= nsegments; j++) { 00317 elliptic_transform(cos(phi1), sin(phi1), da, db, dalpha, &tx, 00318 &ty); 00319 Vect_append_point(nPoints, x[i] + tx, y[i] + ty, 0); 00320 phi1 += angular_step; 00321 } 00322 } 00323 00324 if ((!looped) && (i == np - 2)) { 00325 Vect_append_point(nPoints, mx, my, 0); 00326 } 00327 } 00328 00329 if (looped) { 00330 Vect_append_point(nPoints, nPoints->x[0], nPoints->y[0], 00331 nPoints->z[0]); 00332 } 00333 00334 Vect_line_prune(nPoints); 00335 00336 if (looped) { 00337 Vect_line_delete_point(Points, Points->n_points - 1); 00338 } 00339 } 00340 00341 /* input line must be looped */ 00342 static void convolution_line(struct line_pnts *Points, double da, double db, 00343 double dalpha, int side, int round, int caps, 00344 double tol, struct line_pnts *nPoints) 00345 { 00346 int i, j, res, np; 00347 double *x, *y; 00348 double tx, ty, vx, vy, wx, wy, nx, ny, mx, my, rx, ry; 00349 double vx1, vy1, wx1, wy1; 00350 double a0, b0, c0, a1, b1, c1; 00351 double phi1, phi2, delta_phi; 00352 double nsegments, angular_tol, angular_step; 00353 double angle0, angle1; 00354 int inner_corner, turns360; 00355 00356 G_debug(3, "convolution_line() side = %d", side); 00357 00358 np = Points->n_points; 00359 x = Points->x; 00360 y = Points->y; 00361 if ((np == 0) || (np == 1)) 00362 return; 00363 if ((x[0] != x[np - 1]) || (y[0] != y[np - 1])) { 00364 G_fatal_error(_("Line is not looped")); 00365 return; 00366 } 00367 00368 Vect_reset_line(nPoints); 00369 00370 if ((da == 0) || (db == 0)) { 00371 Vect_copy_xyz_to_pnts(nPoints, x, y, NULL, np); 00372 return; 00373 } 00374 00375 side = (side >= 0) ? (1) : (-1); /* normalize variable */ 00376 dalpha *= PI / 180; /* convert dalpha from degrees to radians */ 00377 angular_tol = angular_tolerance(tol, da, db); 00378 00379 i = np - 2; 00380 norm_vector(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty); 00381 elliptic_tangent(side * tx, side * ty, da, db, dalpha, &vx, &vy); 00382 angle1 = atan2(ty, tx); 00383 nx = x[i] + vx; 00384 ny = y[i] + vy; 00385 mx = x[i + 1] + vx; 00386 my = y[i + 1] + vy; 00387 if (!round) 00388 line_coefficients(nx, ny, mx, my, &a1, &b1, &c1); 00389 00390 for (i = 0; i <= np - 2; i++) { 00391 G_debug(4, "point %d, segment %d-%d", i, i, i + 1); 00392 /* save the old values */ 00393 if (!round) { 00394 a0 = a1; 00395 b0 = b1; 00396 c0 = c1; 00397 } 00398 wx = vx; 00399 wy = vy; 00400 angle0 = angle1; 00401 00402 norm_vector(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty); 00403 if ((tx == 0) && (ty == 0)) 00404 continue; 00405 elliptic_tangent(side * tx, side * ty, da, db, dalpha, &vx, &vy); 00406 angle1 = atan2(ty, tx); 00407 nx = x[i] + vx; 00408 ny = y[i] + vy; 00409 mx = x[i + 1] + vx; 00410 my = y[i + 1] + vy; 00411 if (!round) 00412 line_coefficients(nx, ny, mx, my, &a1, &b1, &c1); 00413 00414 00415 delta_phi = angle1 - angle0; 00416 if (delta_phi > PI) 00417 delta_phi -= 2 * PI; 00418 else if (delta_phi <= -PI) 00419 delta_phi += 2 * PI; 00420 /* now delta_phi is in [-pi;pi] */ 00421 turns360 = (fabs(fabs(delta_phi) - PI) < 1e-15); 00422 inner_corner = (side * delta_phi <= 0) && (!turns360); 00423 00424 00425 /* if <line turns 360> and (<caps> and <not round>) */ 00426 if (turns360 && caps && (!round)) { 00427 norm_vector(0, 0, vx, vy, &tx, &ty); 00428 elliptic_tangent(side * tx, side * ty, da, db, dalpha, &tx, &ty); 00429 Vect_append_point(nPoints, x[i] + wx + tx, y[i] + wy + ty, 0); 00430 G_debug(4, " append point (c) x=%.16f y=%.16f", x[i] + wx + tx, 00431 y[i] + wy + ty); 00432 Vect_append_point(nPoints, nx + tx, ny + ty, 0); /* nx == x[i] + vx, ny == y[i] + vy */ 00433 G_debug(4, " append point (c) x=%.16f y=%.16f", nx + tx, ny + ty); 00434 } 00435 00436 if ((!turns360) && (!round) && (!inner_corner)) { 00437 res = line_intersection(a0, b0, c0, a1, b1, c1, &rx, &ry); 00438 if (res == 1) { 00439 Vect_append_point(nPoints, rx, ry, 0); 00440 G_debug(4, " append point (o) x=%.16f y=%.16f", rx, ry); 00441 } 00442 else if (res == 2) { 00443 /* no need to append point in this case */ 00444 } 00445 else 00446 G_fatal_error 00447 ("unexpected result of line_intersection() res = %d", 00448 res); 00449 } 00450 00451 if (round && (!inner_corner) && (!turns360 || caps)) { 00452 /* we should draw elliptical arc for outside corner */ 00453 00454 /* inverse transforms */ 00455 elliptic_transform(wx, wy, 1 / da, 1 / db, dalpha, &wx1, &wy1); 00456 elliptic_transform(vx, vy, 1 / da, 1 / db, dalpha, &vx1, &vy1); 00457 00458 phi1 = atan2(wy1, wx1); 00459 phi2 = atan2(vy1, vx1); 00460 delta_phi = side * (phi2 - phi1); 00461 00462 /* make delta_phi in [0, 2pi] */ 00463 if (delta_phi < 0) 00464 delta_phi += 2 * PI; 00465 00466 nsegments = (int)(delta_phi / angular_tol) + 1; 00467 angular_step = side * (delta_phi / nsegments); 00468 00469 phi1 += angular_step; 00470 for (j = 1; j <= nsegments - 1; j++) { 00471 elliptic_transform(cos(phi1), sin(phi1), da, db, dalpha, &tx, 00472 &ty); 00473 Vect_append_point(nPoints, x[i] + tx, y[i] + ty, 0); 00474 G_debug(4, " append point (r) x=%.16f y=%.16f", x[i] + tx, 00475 y[i] + ty); 00476 phi1 += angular_step; 00477 } 00478 } 00479 00480 Vect_append_point(nPoints, nx, ny, 0); 00481 G_debug(4, " append point (s) x=%.16f y=%.16f", nx, ny); 00482 Vect_append_point(nPoints, mx, my, 0); 00483 G_debug(4, " append point (s) x=%.16f y=%.16f", mx, my); 00484 } 00485 00486 /* close the output line */ 00487 Vect_append_point(nPoints, nPoints->x[0], nPoints->y[0], nPoints->z[0]); 00488 /* Vect_line_prune ( nPoints ); */ 00489 } 00490 00491 /* 00492 * side: side >= 0 - extracts contour on right side of edge, side < 0 - extracts contour on left side of edge 00493 * if the extracted contour is the outer contour, it is returned in ccw order 00494 * else if it is inner contour, it is returned in cw order 00495 */ 00496 static void extract_contour(struct planar_graph *pg, struct pg_edge *first, 00497 int side, int winding, int stop_at_line_end, 00498 struct line_pnts *nPoints) 00499 { 00500 int j; 00501 int v; /* current vertex number */ 00502 int v0; 00503 int eside; /* side of the current edge */ 00504 double eangle; /* current edge angle with Ox (according to the current direction) */ 00505 struct pg_vertex *vert; /* current vertex */ 00506 struct pg_vertex *vert0; /* last vertex */ 00507 struct pg_edge *edge; /* current edge; must be edge of vert */ 00508 00509 /* int cs; *//* on which side are we turning along the contour */ 00510 /* we will always turn right and dont need that one */ 00511 double opt_angle, tangle; 00512 int opt_j, opt_side, opt_flag; 00513 00514 G_debug(3, 00515 "extract_contour(): v1=%d, v2=%d, side=%d, stop_at_line_end=%d", 00516 first->v1, first->v2, side, stop_at_line_end); 00517 00518 Vect_reset_line(nPoints); 00519 00520 edge = first; 00521 if (side >= 0) { 00522 eside = 1; 00523 v0 = edge->v1; 00524 v = edge->v2; 00525 } 00526 else { 00527 eside = -1; 00528 v0 = edge->v2; 00529 v = edge->v1; 00530 } 00531 vert0 = &(pg->v[v0]); 00532 vert = &(pg->v[v]); 00533 eangle = atan2(vert->y - vert0->y, vert->x - vert0->x); 00534 00535 while (1) { 00536 Vect_append_point(nPoints, vert0->x, vert0->y, 0); 00537 G_debug(4, "ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d", v0, 00538 v, eside, edge->v1, edge->v2); 00539 G_debug(4, "ec: append point x=%.18f y=%.18f", vert0->x, vert0->y); 00540 00541 /* mark current edge as visited on the appropriate side */ 00542 if (eside == 1) { 00543 edge->visited_right = 1; 00544 edge->winding_right = winding; 00545 } 00546 else { 00547 edge->visited_left = 1; 00548 edge->winding_left = winding; 00549 } 00550 00551 opt_flag = 1; 00552 for (j = 0; j < vert->ecount; j++) { 00553 /* exclude current edge */ 00554 if (vert->edges[j] != edge) { 00555 tangle = vert->angles[j] - eangle; 00556 if (tangle < -PI) 00557 tangle += 2 * PI; 00558 else if (tangle > PI) 00559 tangle -= 2 * PI; 00560 /* now tangle is in (-PI, PI) */ 00561 00562 if (opt_flag || (tangle < opt_angle)) { 00563 opt_j = j; 00564 opt_side = (vert->edges[j]->v1 == v) ? (1) : (-1); 00565 opt_angle = tangle; 00566 opt_flag = 0; 00567 } 00568 } 00569 } 00570 00571 /* 00572 G_debug(4, "ec: opt: side=%d opt_flag=%d opt_angle=%.18f opt_j=%d opt_step=%d", 00573 side, opt_flag, opt_angle, opt_j, opt_step); 00574 */ 00575 00576 /* if line end is reached (no other edges at curr vertex) */ 00577 if (opt_flag) { 00578 if (stop_at_line_end) { 00579 G_debug(3, " end has been reached, will stop here"); 00580 break; 00581 } 00582 else { 00583 opt_j = 0; /* the only edge of vert is vert->edges[0] */ 00584 opt_side = -eside; /* go to the other side of the current edge */ 00585 G_debug(3, " end has been reached, turning around"); 00586 } 00587 } 00588 00589 /* break condition */ 00590 if ((vert->edges[opt_j] == first) && (opt_side == side)) 00591 break; 00592 if (opt_side == 1) { 00593 if (vert->edges[opt_j]->visited_right) { 00594 G_warning(_("Next edge was visited but it is not the first one !!! breaking loop")); 00595 G_debug(4, 00596 "ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d", 00597 v, (edge->v1 == v) ? (edge->v2) : (edge->v1), 00598 opt_side, vert->edges[opt_j]->v1, 00599 vert->edges[opt_j]->v2); 00600 break; 00601 } 00602 } 00603 else { 00604 if (vert->edges[opt_j]->visited_left) { 00605 G_warning(_("Next edge was visited but it is not the first one !!! breaking loop")); 00606 G_debug(4, 00607 "ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d", 00608 v, (edge->v1 == v) ? (edge->v2) : (edge->v1), 00609 opt_side, vert->edges[opt_j]->v1, 00610 vert->edges[opt_j]->v2); 00611 break; 00612 } 00613 } 00614 00615 edge = vert->edges[opt_j]; 00616 eside = opt_side; 00617 v0 = v; 00618 v = (edge->v1 == v) ? (edge->v2) : (edge->v1); 00619 vert0 = vert; 00620 vert = &(pg->v[v]); 00621 eangle = vert0->angles[opt_j]; 00622 } 00623 Vect_append_point(nPoints, vert->x, vert->y, 0); 00624 G_debug(4, "ec: append point x=%.18f y=%.18f", vert->x, vert->y); 00625 00626 return; 00627 } 00628 00629 /* 00630 * This function extracts the outer contour of a (self crossing) line. 00631 * It can generate left/right contour if none of the line ends are in a loop. 00632 * If one or both of them is in a loop, then there's only one contour 00633 * 00634 * side: side > 0 - right contour, side < 0 - left contour, side = 0 - outer contour 00635 * if side != 0 and there's only one contour, the function returns it 00636 * 00637 * TODO: Implement side != 0 feature; 00638 */ 00639 static void extract_outer_contour(struct planar_graph *pg, int side, 00640 struct line_pnts *nPoints) 00641 { 00642 int i; 00643 int flag; 00644 int v; 00645 struct pg_vertex *vert; 00646 struct pg_edge *edge; 00647 double min_x, min_angle; 00648 00649 G_debug(3, "extract_outer_contour()"); 00650 00651 if (side != 0) { 00652 G_fatal_error(_("side != 0 feature not implemented")); 00653 return; 00654 } 00655 00656 /* find a line segment which is on the outer contour */ 00657 flag = 1; 00658 for (i = 0; i < pg->vcount; i++) { 00659 if (flag || (pg->v[i].x < min_x)) { 00660 v = i; 00661 min_x = pg->v[i].x; 00662 flag = 0; 00663 } 00664 } 00665 vert = &(pg->v[v]); 00666 00667 flag = 1; 00668 for (i = 0; i < vert->ecount; i++) { 00669 if (flag || (vert->angles[i] < min_angle)) { 00670 edge = vert->edges[i]; 00671 min_angle = vert->angles[i]; 00672 flag = 0; 00673 } 00674 } 00675 00676 /* the winding on the outer contour is 0 */ 00677 extract_contour(pg, edge, (edge->v1 == v) ? RIGHT_SIDE : LEFT_SIDE, 0, 0, 00678 nPoints); 00679 00680 return; 00681 } 00682 00683 /* 00684 * Extracts contours which are not visited. 00685 * IMPORTANT: the outer contour must be visited (you should call extract_outer_contour() to do that), 00686 * so that extract_inner_contour() doesn't return it 00687 * 00688 * returns: 0 when there are no more inner contours; otherwise, 1 00689 */ 00690 static int extract_inner_contour(struct planar_graph *pg, int *winding, 00691 struct line_pnts *nPoints) 00692 { 00693 int i, w; 00694 struct pg_edge *edge; 00695 00696 G_debug(3, "extract_inner_contour()"); 00697 00698 for (i = 0; i < pg->ecount; i++) { 00699 edge = &(pg->e[i]); 00700 if (edge->visited_left) { 00701 if (!(pg->e[i].visited_right)) { 00702 w = edge->winding_left - 1; 00703 extract_contour(pg, &(pg->e[i]), RIGHT_SIDE, w, 0, nPoints); 00704 *winding = w; 00705 return 1; 00706 } 00707 } 00708 else { 00709 if (pg->e[i].visited_right) { 00710 w = edge->winding_right + 1; 00711 extract_contour(pg, &(pg->e[i]), LEFT_SIDE, w, 0, nPoints); 00712 *winding = w; 00713 return 1; 00714 } 00715 } 00716 } 00717 00718 return 0; 00719 } 00720 00721 /* point_in_buf - test if point px,py is in d buffer of Points 00722 ** dalpha is in degrees 00723 ** returns: 1 in buffer 00724 ** 0 not in buffer 00725 */ 00726 static int point_in_buf(struct line_pnts *Points, double px, double py, double da, 00727 double db, double dalpha) 00728 { 00729 int i, np; 00730 double cx, cy; 00731 double delta, delta_k, k; 00732 double vx, vy, wx, wy, mx, my, nx, ny; 00733 double len, tx, ty, d, da2; 00734 00735 G_debug(3, "point_in_buf()"); 00736 00737 dalpha *= PI / 180; /* convert dalpha from degrees to radians */ 00738 00739 np = Points->n_points; 00740 da2 = da * da; 00741 for (i = 0; i < np - 1; i++) { 00742 vx = Points->x[i]; 00743 vy = Points->y[i]; 00744 wx = Points->x[i + 1]; 00745 wy = Points->y[i + 1]; 00746 00747 if (da != db) { 00748 mx = wx - vx; 00749 my = wy - vy; 00750 len = LENGTH(mx, my); 00751 elliptic_tangent(mx / len, my / len, da, db, dalpha, &cx, &cy); 00752 00753 delta = mx * cy - my * cx; 00754 delta_k = (px - vx) * cy - (py - vy) * cx; 00755 k = delta_k / delta; 00756 /* G_debug(4, "k = %g, k1 = %g", k, (mx * (px - vx) + my * (py - vy)) / (mx * mx + my * my)); */ 00757 if (k <= 0) { 00758 nx = vx; 00759 ny = vy; 00760 } 00761 else if (k >= 1) { 00762 nx = wx; 00763 ny = wy; 00764 } 00765 else { 00766 nx = vx + k * mx; 00767 ny = vy + k * my; 00768 } 00769 00770 /* inverse transform */ 00771 elliptic_transform(px - nx, py - ny, 1 / da, 1 / db, dalpha, &tx, 00772 &ty); 00773 00774 d = dig_distance2_point_to_line(nx + tx, ny + ty, 0, vx, vy, 0, 00775 wx, wy, 0, 0, NULL, NULL, NULL, 00776 NULL, NULL); 00777 00778 /* G_debug(4, "sqrt(d)*da = %g, len' = %g, olen = %g", sqrt(d)*da, da*LENGTH(tx,ty), LENGTH((px-nx),(py-ny))); */ 00779 if (d <= 1) { 00780 /* G_debug(1, "d=%g", d); */ 00781 return 1; 00782 } 00783 } 00784 else { 00785 d = dig_distance2_point_to_line(px, py, 0, vx, vy, 0, wx, wy, 0, 00786 0, NULL, NULL, NULL, NULL, NULL); 00787 /* G_debug(4, "sqrt(d) = %g", sqrt(d)); */ 00788 if (d <= da2) { 00789 return 1; 00790 } 00791 } 00792 } 00793 00794 return 0; 00795 } 00796 00797 /* returns 0 for ccw, non-zero for cw 00798 */ 00799 static int get_polygon_orientation(const double *x, const double *y, int n) 00800 { 00801 double x1, y1, x2, y2; 00802 double area; 00803 00804 x2 = x[n - 1]; 00805 y2 = y[n - 1]; 00806 00807 area = 0; 00808 while (--n >= 0) { 00809 x1 = x2; 00810 y1 = y2; 00811 00812 x2 = *x++; 00813 y2 = *y++; 00814 00815 area += (y2 + y1) * (x2 - x1); 00816 } 00817 00818 return (area > 0); 00819 } 00820 00821 /* internal */ 00822 static void add_line_to_array(struct line_pnts *Points, 00823 struct line_pnts ***arrPoints, int *count, 00824 int *allocated, int more) 00825 { 00826 if (*allocated == *count) { 00827 *allocated += more; 00828 *arrPoints = 00829 G_realloc(*arrPoints, (*allocated) * sizeof(struct line_pnts *)); 00830 } 00831 (*arrPoints)[*count] = Points; 00832 (*count)++; 00833 00834 return; 00835 } 00836 00837 static void destroy_lines_array(struct line_pnts **arr, int count) 00838 { 00839 int i; 00840 00841 for (i = 0; i < count; i++) 00842 Vect_destroy_line_struct(arr[i]); 00843 G_free(arr); 00844 } 00845 00846 /* area_outer and area_isles[i] must be closed non self-intersecting lines 00847 side: 0 - auto, 1 - right, -1 left 00848 */ 00849 static void buffer_lines(struct line_pnts *area_outer, struct line_pnts **area_isles, 00850 int isles_count, int side, double da, double db, 00851 double dalpha, int round, int caps, double tol, 00852 struct line_pnts **oPoints, struct line_pnts ***iPoints, 00853 int *inner_count) 00854 { 00855 struct planar_graph *pg2; 00856 struct line_pnts *sPoints, *cPoints; 00857 struct line_pnts **arrPoints; 00858 int i, count = 0; 00859 int res, winding; 00860 int auto_side; 00861 int more = 8; 00862 int allocated = 0; 00863 double px, py; 00864 00865 G_debug(3, "buffer_lines()"); 00866 00867 auto_side = (side == 0); 00868 00869 /* initializations */ 00870 sPoints = Vect_new_line_struct(); 00871 cPoints = Vect_new_line_struct(); 00872 arrPoints = NULL; 00873 00874 /* outer contour */ 00875 G_debug(3, " processing outer contour"); 00876 *oPoints = Vect_new_line_struct(); 00877 if (auto_side) 00878 side = 00879 get_polygon_orientation(area_outer->x, area_outer->y, 00880 area_outer->n_points - 00881 1) ? LEFT_SIDE : RIGHT_SIDE; 00882 convolution_line(area_outer, da, db, dalpha, side, round, caps, tol, 00883 sPoints); 00884 pg2 = pg_create(sPoints); 00885 extract_outer_contour(pg2, 0, *oPoints); 00886 res = extract_inner_contour(pg2, &winding, cPoints); 00887 while (res != 0) { 00888 if (winding == 0) { 00889 if (!Vect_point_in_poly(cPoints->x[0], cPoints->y[0], area_outer)) { 00890 if (Vect_get_point_in_poly(cPoints, &px, &py) != 0) 00891 G_fatal_error(_("Vect_get_point_in_poly() failed")); 00892 if (!point_in_buf(area_outer, px, py, da, db, dalpha)) { 00893 add_line_to_array(cPoints, &arrPoints, &count, &allocated, 00894 more); 00895 cPoints = Vect_new_line_struct(); 00896 } 00897 } 00898 } 00899 res = extract_inner_contour(pg2, &winding, cPoints); 00900 } 00901 pg_destroy_struct(pg2); 00902 00903 /* inner contours */ 00904 G_debug(3, " processing inner contours"); 00905 for (i = 0; i < isles_count; i++) { 00906 if (auto_side) 00907 side = 00908 get_polygon_orientation(area_isles[i]->x, area_isles[i]->y, 00909 area_isles[i]->n_points - 00910 1) ? RIGHT_SIDE : LEFT_SIDE; 00911 convolution_line(area_isles[i], da, db, dalpha, side, round, caps, 00912 tol, sPoints); 00913 pg2 = pg_create(sPoints); 00914 extract_outer_contour(pg2, 0, cPoints); 00915 res = extract_inner_contour(pg2, &winding, cPoints); 00916 while (res != 0) { 00917 if (winding == -1) { 00918 /* we need to check if the area is in the buffer. 00919 I've simplfied convolution_line(), so that it runs faster, 00920 however that leads to ocasional problems */ 00921 if (Vect_point_in_poly 00922 (cPoints->x[0], cPoints->y[0], area_isles[i])) { 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_isles[i], px, py, da, db, dalpha)) { 00926 add_line_to_array(cPoints, &arrPoints, &count, 00927 &allocated, 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 arrPoints = G_realloc(arrPoints, count * sizeof(struct line_pnts *)); 00938 *inner_count = count; 00939 *iPoints = arrPoints; 00940 00941 Vect_destroy_line_struct(sPoints); 00942 Vect_destroy_line_struct(cPoints); 00943 00944 G_debug(3, "buffer_lines() ... done"); 00945 00946 return; 00947 } 00948 00949 00966 void Vect_line_buffer2(struct line_pnts *Points, double da, double db, 00967 double dalpha, int round, int caps, double tol, 00968 struct line_pnts **oPoints, 00969 struct line_pnts ***iPoints, int *inner_count) 00970 { 00971 struct planar_graph *pg; 00972 struct line_pnts *tPoints, *outer; 00973 struct line_pnts **isles; 00974 int isles_count = 0; 00975 int res, winding; 00976 int more = 8; 00977 int isles_allocated = 0; 00978 00979 G_debug(2, "Vect_line_buffer()"); 00980 00981 Vect_line_prune(Points); 00982 00983 if (Points->n_points == 1) 00984 return Vect_point_buffer2(Points->x[0], Points->y[0], da, db, 00985 dalpha, round, tol, oPoints); 00986 00987 /* initializations */ 00988 tPoints = Vect_new_line_struct(); 00989 isles = NULL; 00990 pg = pg_create(Points); 00991 00992 /* outer contour */ 00993 outer = Vect_new_line_struct(); 00994 extract_outer_contour(pg, 0, outer); 00995 00996 /* inner contours */ 00997 res = extract_inner_contour(pg, &winding, tPoints); 00998 while (res != 0) { 00999 add_line_to_array(tPoints, &isles, &isles_count, &isles_allocated, 01000 more); 01001 tPoints = Vect_new_line_struct(); 01002 res = extract_inner_contour(pg, &winding, tPoints); 01003 } 01004 01005 buffer_lines(outer, isles, isles_count, RIGHT_SIDE, da, db, dalpha, round, 01006 caps, tol, oPoints, iPoints, inner_count); 01007 01008 Vect_destroy_line_struct(tPoints); 01009 Vect_destroy_line_struct(outer); 01010 destroy_lines_array(isles, isles_count); 01011 pg_destroy_struct(pg); 01012 01013 return; 01014 } 01015 01031 void Vect_area_buffer2(struct Map_info *Map, int area, double da, double db, 01032 double dalpha, int round, int caps, double tol, 01033 struct line_pnts **oPoints, 01034 struct line_pnts ***iPoints, int *inner_count) 01035 { 01036 struct line_pnts *tPoints, *outer; 01037 struct line_pnts **isles; 01038 int isles_count = 0, n_isles; 01039 int i, isle; 01040 int more = 8; 01041 int isles_allocated = 0; 01042 01043 G_debug(2, "Vect_area_buffer()"); 01044 01045 /* initializations */ 01046 tPoints = Vect_new_line_struct(); 01047 n_isles = Vect_get_area_num_isles(Map, area); 01048 isles_allocated = n_isles; 01049 isles = G_malloc(isles_allocated * sizeof(struct line_pnts *)); 01050 01051 /* outer contour */ 01052 outer = Vect_new_line_struct(); 01053 Vect_get_area_points(Map, area, outer); 01054 /* does not work with zero length line segments */ 01055 Vect_line_prune(outer); 01056 01057 /* inner contours */ 01058 for (i = 0; i < n_isles; i++) { 01059 isle = Vect_get_area_isle(Map, area, i); 01060 Vect_get_isle_points(Map, isle, tPoints); 01061 01062 /* Check if the isle is big enough */ 01063 /* 01064 if (Vect_line_length(tPoints) < 2*PI*max) 01065 continue; 01066 */ 01067 /* does not work with zero length line segments */ 01068 Vect_line_prune(tPoints); 01069 add_line_to_array(tPoints, &isles, &isles_count, &isles_allocated, 01070 more); 01071 tPoints = Vect_new_line_struct(); 01072 } 01073 01074 buffer_lines(outer, isles, isles_count, 0, da, db, dalpha, round, caps, 01075 tol, oPoints, iPoints, inner_count); 01076 01077 Vect_destroy_line_struct(tPoints); 01078 Vect_destroy_line_struct(outer); 01079 destroy_lines_array(isles, isles_count); 01080 01081 return; 01082 } 01083 01096 void Vect_point_buffer2(double px, double py, double da, double db, 01097 double dalpha, int round, double tol, 01098 struct line_pnts **oPoints) 01099 { 01100 double tx, ty; 01101 double angular_tol, angular_step, phi1; 01102 int j, nsegments; 01103 01104 G_debug(2, "Vect_point_buffer()"); 01105 01106 *oPoints = Vect_new_line_struct(); 01107 01108 dalpha *= PI / 180; /* convert dalpha from degrees to radians */ 01109 01110 if (round || (!round)) { 01111 angular_tol = angular_tolerance(tol, da, db); 01112 01113 nsegments = (int)(2 * PI / angular_tol) + 1; 01114 angular_step = 2 * PI / nsegments; 01115 01116 phi1 = 0; 01117 for (j = 0; j < nsegments; j++) { 01118 elliptic_transform(cos(phi1), sin(phi1), da, db, dalpha, &tx, 01119 &ty); 01120 Vect_append_point(*oPoints, px + tx, py + ty, 0); 01121 phi1 += angular_step; 01122 } 01123 } 01124 else { 01125 01126 } 01127 01128 /* close the output line */ 01129 Vect_append_point(*oPoints, (*oPoints)->x[0], (*oPoints)->y[0], 01130 (*oPoints)->z[0]); 01131 01132 return; 01133 } 01134 01135 01136 /* 01137 \brief Create parrallel line 01138 01139 See also Vect_line_parallel(). 01140 01141 \param InPoints input line geometry 01142 \param da distance along major axis 01143 \param da distance along minor axis 01144 \param dalpha angle between 0x and major axis 01145 \param round make corners round 01146 \param tol maximum distance between theoretical arc and output segments 01147 \param[out] OutPoints output line 01148 */ 01149 void Vect_line_parallel2(struct line_pnts *InPoints, double da, double db, 01150 double dalpha, int side, int round, double tol, 01151 struct line_pnts *OutPoints) 01152 { 01153 G_debug(2, "Vect_line_parallel(): npoints = %d, da = %f, " 01154 "db = %f, dalpha = %f, side = %d, round_corners = %d, tol = %f", 01155 InPoints->n_points, da, db, dalpha, side, round, tol); 01156 01157 parallel_line(InPoints, da, db, dalpha, side, round, 1, NON_LOOPED_LINE, 01158 tol, OutPoints); 01159 01160 /* if (!loops) 01161 clean_parallel(OutPoints, InPoints, distance, rm_end); 01162 */ 01163 01164 return; 01165 }