GRASS Programmer's Manual
6.4.2(2012)
|
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 if (nPoints->n_points > 0 && (nPoints->x[nPoints->n_points - 1] != vert0->x || 00537 nPoints->y[nPoints->n_points - 1] != vert0->y)) 00538 Vect_append_point(nPoints, vert0->x, vert0->y, 0); 00539 else 00540 Vect_append_point(nPoints, vert0->x, vert0->y, 0); 00541 G_debug(4, "ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d", v0, 00542 v, eside, edge->v1, edge->v2); 00543 G_debug(4, "ec: append point x=%.18f y=%.18f", vert0->x, vert0->y); 00544 00545 /* mark current edge as visited on the appropriate side */ 00546 if (eside == 1) { 00547 edge->visited_right = 1; 00548 edge->winding_right = winding; 00549 } 00550 else { 00551 edge->visited_left = 1; 00552 edge->winding_left = winding; 00553 } 00554 00555 opt_flag = 1; 00556 for (j = 0; j < vert->ecount; j++) { 00557 /* exclude current edge */ 00558 if (vert->edges[j] != edge) { 00559 tangle = vert->angles[j] - eangle; 00560 if (tangle < -PI) 00561 tangle += 2 * PI; 00562 else if (tangle > PI) 00563 tangle -= 2 * PI; 00564 /* now tangle is in (-PI, PI) */ 00565 00566 if (opt_flag || (tangle < opt_angle)) { 00567 opt_j = j; 00568 opt_side = (vert->edges[j]->v1 == v) ? (1) : (-1); 00569 opt_angle = tangle; 00570 opt_flag = 0; 00571 } 00572 } 00573 } 00574 00575 /* 00576 G_debug(4, "ec: opt: side=%d opt_flag=%d opt_angle=%.18f opt_j=%d opt_step=%d", 00577 side, opt_flag, opt_angle, opt_j, opt_step); 00578 */ 00579 00580 /* if line end is reached (no other edges at curr vertex) */ 00581 if (opt_flag) { 00582 if (stop_at_line_end) { 00583 G_debug(3, " end has been reached, will stop here"); 00584 break; 00585 } 00586 else { 00587 opt_j = 0; /* the only edge of vert is vert->edges[0] */ 00588 opt_side = -eside; /* go to the other side of the current edge */ 00589 G_debug(3, " end has been reached, turning around"); 00590 } 00591 } 00592 00593 /* break condition */ 00594 if ((vert->edges[opt_j] == first) && (opt_side == side)) 00595 break; 00596 if (opt_side == 1) { 00597 if (vert->edges[opt_j]->visited_right) { 00598 G_warning(_("Next edge was visited but it is not the first one !!! breaking loop")); 00599 G_debug(4, 00600 "ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d", 00601 v, (edge->v1 == v) ? (edge->v2) : (edge->v1), 00602 opt_side, vert->edges[opt_j]->v1, 00603 vert->edges[opt_j]->v2); 00604 break; 00605 } 00606 } 00607 else { 00608 if (vert->edges[opt_j]->visited_left) { 00609 G_warning(_("Next edge was visited but it is not the first one !!! breaking loop")); 00610 G_debug(4, 00611 "ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d", 00612 v, (edge->v1 == v) ? (edge->v2) : (edge->v1), 00613 opt_side, vert->edges[opt_j]->v1, 00614 vert->edges[opt_j]->v2); 00615 break; 00616 } 00617 } 00618 00619 edge = vert->edges[opt_j]; 00620 eside = opt_side; 00621 v0 = v; 00622 v = (edge->v1 == v) ? (edge->v2) : (edge->v1); 00623 vert0 = vert; 00624 vert = &(pg->v[v]); 00625 eangle = vert0->angles[opt_j]; 00626 } 00627 Vect_append_point(nPoints, vert->x, vert->y, 0); 00628 Vect_line_prune(nPoints); 00629 G_debug(4, "ec: append point x=%.18f y=%.18f", vert->x, vert->y); 00630 00631 return; 00632 } 00633 00634 /* 00635 * This function extracts the outer contour of a (self crossing) line. 00636 * It can generate left/right contour if none of the line ends are in a loop. 00637 * If one or both of them is in a loop, then there's only one contour 00638 * 00639 * side: side > 0 - right contour, side < 0 - left contour, side = 0 - outer contour 00640 * if side != 0 and there's only one contour, the function returns it 00641 * 00642 * TODO: Implement side != 0 feature; 00643 */ 00644 static void extract_outer_contour(struct planar_graph *pg, int side, 00645 struct line_pnts *nPoints) 00646 { 00647 int i; 00648 int flag; 00649 int v; 00650 struct pg_vertex *vert; 00651 struct pg_edge *edge; 00652 double min_x, min_angle; 00653 00654 G_debug(3, "extract_outer_contour()"); 00655 00656 if (side != 0) { 00657 G_fatal_error(_("side != 0 feature not implemented")); 00658 return; 00659 } 00660 00661 /* find a line segment which is on the outer contour */ 00662 flag = 1; 00663 for (i = 0; i < pg->vcount; i++) { 00664 if (flag || (pg->v[i].x < min_x)) { 00665 v = i; 00666 min_x = pg->v[i].x; 00667 flag = 0; 00668 } 00669 } 00670 vert = &(pg->v[v]); 00671 00672 flag = 1; 00673 for (i = 0; i < vert->ecount; i++) { 00674 if (flag || (vert->angles[i] < min_angle)) { 00675 edge = vert->edges[i]; 00676 min_angle = vert->angles[i]; 00677 flag = 0; 00678 } 00679 } 00680 00681 /* the winding on the outer contour is 0 */ 00682 extract_contour(pg, edge, (edge->v1 == v) ? RIGHT_SIDE : LEFT_SIDE, 0, 0, 00683 nPoints); 00684 00685 return; 00686 } 00687 00688 /* 00689 * Extracts contours which are not visited. 00690 * IMPORTANT: the outer contour must be visited (you should call extract_outer_contour() to do that), 00691 * so that extract_inner_contour() doesn't return it 00692 * 00693 * returns: 0 when there are no more inner contours; otherwise, 1 00694 */ 00695 static int extract_inner_contour(struct planar_graph *pg, int *winding, 00696 struct line_pnts *nPoints) 00697 { 00698 int i, w; 00699 struct pg_edge *edge; 00700 00701 G_debug(3, "extract_inner_contour()"); 00702 00703 for (i = 0; i < pg->ecount; i++) { 00704 edge = &(pg->e[i]); 00705 if (edge->visited_left) { 00706 if (!(pg->e[i].visited_right)) { 00707 w = edge->winding_left - 1; 00708 extract_contour(pg, &(pg->e[i]), RIGHT_SIDE, w, 0, nPoints); 00709 *winding = w; 00710 return 1; 00711 } 00712 } 00713 else { 00714 if (pg->e[i].visited_right) { 00715 w = edge->winding_right + 1; 00716 extract_contour(pg, &(pg->e[i]), LEFT_SIDE, w, 0, nPoints); 00717 *winding = w; 00718 return 1; 00719 } 00720 } 00721 } 00722 00723 return 0; 00724 } 00725 00726 /* point_in_buf - test if point px,py is in d buffer of Points 00727 ** dalpha is in degrees 00728 ** returns: 1 in buffer 00729 ** 0 not in buffer 00730 */ 00731 static int point_in_buf(struct line_pnts *Points, double px, double py, double da, 00732 double db, double dalpha) 00733 { 00734 int i, np; 00735 double cx, cy; 00736 double delta, delta_k, k; 00737 double vx, vy, wx, wy, mx, my, nx, ny; 00738 double len, tx, ty, d, da2; 00739 00740 G_debug(3, "point_in_buf()"); 00741 00742 dalpha *= PI / 180; /* convert dalpha from degrees to radians */ 00743 00744 np = Points->n_points; 00745 da2 = da * da; 00746 for (i = 0; i < np - 1; i++) { 00747 vx = Points->x[i]; 00748 vy = Points->y[i]; 00749 wx = Points->x[i + 1]; 00750 wy = Points->y[i + 1]; 00751 00752 if (da != db) { 00753 mx = wx - vx; 00754 my = wy - vy; 00755 len = LENGTH(mx, my); 00756 elliptic_tangent(mx / len, my / len, da, db, dalpha, &cx, &cy); 00757 00758 delta = mx * cy - my * cx; 00759 delta_k = (px - vx) * cy - (py - vy) * cx; 00760 k = delta_k / delta; 00761 /* G_debug(4, "k = %g, k1 = %g", k, (mx * (px - vx) + my * (py - vy)) / (mx * mx + my * my)); */ 00762 if (k <= 0) { 00763 nx = vx; 00764 ny = vy; 00765 } 00766 else if (k >= 1) { 00767 nx = wx; 00768 ny = wy; 00769 } 00770 else { 00771 nx = vx + k * mx; 00772 ny = vy + k * my; 00773 } 00774 00775 /* inverse transform */ 00776 elliptic_transform(px - nx, py - ny, 1 / da, 1 / db, dalpha, &tx, 00777 &ty); 00778 00779 d = dig_distance2_point_to_line(nx + tx, ny + ty, 0, vx, vy, 0, 00780 wx, wy, 0, 0, NULL, NULL, NULL, 00781 NULL, NULL); 00782 00783 /* G_debug(4, "sqrt(d)*da = %g, len' = %g, olen = %g", sqrt(d)*da, da*LENGTH(tx,ty), LENGTH((px-nx),(py-ny))); */ 00784 if (d <= 1) { 00785 /* G_debug(1, "d=%g", d); */ 00786 return 1; 00787 } 00788 } 00789 else { 00790 d = dig_distance2_point_to_line(px, py, 0, vx, vy, 0, wx, wy, 0, 00791 0, NULL, NULL, NULL, NULL, NULL); 00792 /* G_debug(4, "sqrt(d) = %g", sqrt(d)); */ 00793 if (d <= da2) { 00794 return 1; 00795 } 00796 } 00797 } 00798 00799 return 0; 00800 } 00801 00802 /* returns 0 for ccw, non-zero for cw 00803 */ 00804 static int get_polygon_orientation(const double *x, const double *y, int n) 00805 { 00806 double x1, y1, x2, y2; 00807 double area; 00808 00809 x2 = x[n - 1]; 00810 y2 = y[n - 1]; 00811 00812 area = 0; 00813 while (--n >= 0) { 00814 x1 = x2; 00815 y1 = y2; 00816 00817 x2 = *x++; 00818 y2 = *y++; 00819 00820 area += (y2 + y1) * (x2 - x1); 00821 } 00822 00823 return (area > 0); 00824 } 00825 00826 /* internal */ 00827 static void add_line_to_array(struct line_pnts *Points, 00828 struct line_pnts ***arrPoints, int *count, 00829 int *allocated, int more) 00830 { 00831 if (*allocated == *count) { 00832 *allocated += more; 00833 *arrPoints = 00834 G_realloc(*arrPoints, (*allocated) * sizeof(struct line_pnts *)); 00835 } 00836 (*arrPoints)[*count] = Points; 00837 (*count)++; 00838 00839 return; 00840 } 00841 00842 static void destroy_lines_array(struct line_pnts **arr, int count) 00843 { 00844 int i; 00845 00846 for (i = 0; i < count; i++) 00847 Vect_destroy_line_struct(arr[i]); 00848 G_free(arr); 00849 } 00850 00851 /* area_outer and area_isles[i] must be closed non self-intersecting lines 00852 side: 0 - auto, 1 - right, -1 left 00853 */ 00854 static void buffer_lines(struct line_pnts *area_outer, struct line_pnts **area_isles, 00855 int isles_count, int side, double da, double db, 00856 double dalpha, int round, int caps, double tol, 00857 struct line_pnts **oPoints, struct line_pnts ***iPoints, 00858 int *inner_count) 00859 { 00860 struct planar_graph *pg2; 00861 struct line_pnts *sPoints, *cPoints; 00862 struct line_pnts **arrPoints; 00863 int i, count = 0; 00864 int res, winding; 00865 int auto_side; 00866 int more = 8; 00867 int allocated = 0; 00868 double px, py; 00869 00870 G_debug(3, "buffer_lines()"); 00871 00872 auto_side = (side == 0); 00873 00874 /* initializations */ 00875 sPoints = Vect_new_line_struct(); 00876 cPoints = Vect_new_line_struct(); 00877 arrPoints = NULL; 00878 00879 /* outer contour */ 00880 G_debug(3, " processing outer contour"); 00881 *oPoints = Vect_new_line_struct(); 00882 if (auto_side) 00883 side = 00884 get_polygon_orientation(area_outer->x, area_outer->y, 00885 area_outer->n_points - 00886 1) ? LEFT_SIDE : RIGHT_SIDE; 00887 convolution_line(area_outer, da, db, dalpha, side, round, caps, tol, 00888 sPoints); 00889 pg2 = pg_create(sPoints); 00890 extract_outer_contour(pg2, 0, *oPoints); 00891 res = extract_inner_contour(pg2, &winding, cPoints); 00892 while (res != 0) { 00893 if (winding == 0) { 00894 int check_poly = 1; 00895 double area_size; 00896 00897 dig_find_area_poly(cPoints, &area_size); 00898 if (area_size == 0) { 00899 G_warning(_("zero area size")); 00900 check_poly = 0; 00901 } 00902 if (cPoints->x[0] != cPoints->x[cPoints->n_points - 1] || 00903 cPoints->y[0] != cPoints->y[cPoints->n_points - 1]) { 00904 00905 G_warning(_("Line was not closed")); 00906 check_poly = 0; 00907 } 00908 00909 if (check_poly && !Vect_point_in_poly(cPoints->x[0], cPoints->y[0], area_outer)) { 00910 if (Vect_get_point_in_poly(cPoints, &px, &py) == 0) { 00911 if (!point_in_buf(area_outer, px, py, da, db, dalpha)) { 00912 add_line_to_array(cPoints, &arrPoints, &count, &allocated, 00913 more); 00914 cPoints = Vect_new_line_struct(); 00915 } 00916 } 00917 else { 00918 G_warning(_("Vect_get_point_in_poly() failed")); 00919 } 00920 } 00921 } 00922 res = extract_inner_contour(pg2, &winding, cPoints); 00923 } 00924 pg_destroy_struct(pg2); 00925 00926 /* inner contours */ 00927 G_debug(3, " processing inner contours"); 00928 for (i = 0; i < isles_count; i++) { 00929 if (auto_side) 00930 side = 00931 get_polygon_orientation(area_isles[i]->x, area_isles[i]->y, 00932 area_isles[i]->n_points - 00933 1) ? RIGHT_SIDE : LEFT_SIDE; 00934 convolution_line(area_isles[i], da, db, dalpha, side, round, caps, 00935 tol, sPoints); 00936 pg2 = pg_create(sPoints); 00937 extract_outer_contour(pg2, 0, cPoints); 00938 res = extract_inner_contour(pg2, &winding, cPoints); 00939 while (res != 0) { 00940 if (winding == -1) { 00941 int check_poly = 1; 00942 double area_size; 00943 00944 dig_find_area_poly(cPoints, &area_size); 00945 if (area_size == 0) { 00946 G_warning(_("zero area size")); 00947 check_poly = 0; 00948 } 00949 if (cPoints->x[0] != cPoints->x[cPoints->n_points - 1] || 00950 cPoints->y[0] != cPoints->y[cPoints->n_points - 1]) { 00951 00952 G_warning(_("Line was not closed")); 00953 check_poly = 0; 00954 } 00955 00956 /* we need to check if the area is in the buffer. 00957 I've simplfied convolution_line(), so that it runs faster, 00958 however that leads to ocasional problems */ 00959 if (check_poly && Vect_point_in_poly 00960 (cPoints->x[0], cPoints->y[0], area_isles[i])) { 00961 if (Vect_get_point_in_poly(cPoints, &px, &py) == 0) { 00962 if (!point_in_buf(area_isles[i], px, py, da, db, dalpha)) { 00963 add_line_to_array(cPoints, &arrPoints, &count, 00964 &allocated, more); 00965 cPoints = Vect_new_line_struct(); 00966 } 00967 } 00968 else { 00969 G_warning(_("Vect_get_point_in_poly() failed")); 00970 } 00971 } 00972 } 00973 res = extract_inner_contour(pg2, &winding, cPoints); 00974 } 00975 pg_destroy_struct(pg2); 00976 } 00977 00978 arrPoints = G_realloc(arrPoints, count * sizeof(struct line_pnts *)); 00979 *inner_count = count; 00980 *iPoints = arrPoints; 00981 00982 Vect_destroy_line_struct(sPoints); 00983 Vect_destroy_line_struct(cPoints); 00984 00985 G_debug(3, "buffer_lines() ... done"); 00986 00987 return; 00988 } 00989 00990 01007 void Vect_line_buffer2(struct line_pnts *Points, double da, double db, 01008 double dalpha, int round, int caps, double tol, 01009 struct line_pnts **oPoints, 01010 struct line_pnts ***iPoints, int *inner_count) 01011 { 01012 struct planar_graph *pg; 01013 struct line_pnts *tPoints, *outer; 01014 struct line_pnts **isles; 01015 int isles_count = 0; 01016 int res, winding; 01017 int more = 8; 01018 int isles_allocated = 0; 01019 01020 G_debug(2, "Vect_line_buffer()"); 01021 01022 Vect_line_prune(Points); 01023 01024 if (Points->n_points == 1) 01025 return Vect_point_buffer2(Points->x[0], Points->y[0], da, db, 01026 dalpha, round, tol, oPoints); 01027 01028 /* initializations */ 01029 tPoints = Vect_new_line_struct(); 01030 isles = NULL; 01031 pg = pg_create(Points); 01032 01033 /* outer contour */ 01034 outer = Vect_new_line_struct(); 01035 extract_outer_contour(pg, 0, outer); 01036 01037 /* inner contours */ 01038 res = extract_inner_contour(pg, &winding, tPoints); 01039 while (res != 0) { 01040 add_line_to_array(tPoints, &isles, &isles_count, &isles_allocated, 01041 more); 01042 tPoints = Vect_new_line_struct(); 01043 res = extract_inner_contour(pg, &winding, tPoints); 01044 } 01045 01046 buffer_lines(outer, isles, isles_count, RIGHT_SIDE, da, db, dalpha, round, 01047 caps, tol, oPoints, iPoints, inner_count); 01048 01049 Vect_destroy_line_struct(tPoints); 01050 Vect_destroy_line_struct(outer); 01051 destroy_lines_array(isles, isles_count); 01052 pg_destroy_struct(pg); 01053 01054 return; 01055 } 01056 01072 void Vect_area_buffer2(struct Map_info *Map, int area, double da, double db, 01073 double dalpha, int round, int caps, double tol, 01074 struct line_pnts **oPoints, 01075 struct line_pnts ***iPoints, int *inner_count) 01076 { 01077 struct line_pnts *tPoints, *outer; 01078 struct line_pnts **isles; 01079 int isles_count = 0, n_isles; 01080 int i, isle; 01081 int more = 8; 01082 int isles_allocated = 0; 01083 01084 G_debug(2, "Vect_area_buffer()"); 01085 01086 /* initializations */ 01087 tPoints = Vect_new_line_struct(); 01088 n_isles = Vect_get_area_num_isles(Map, area); 01089 isles_allocated = n_isles; 01090 isles = G_malloc(isles_allocated * sizeof(struct line_pnts *)); 01091 01092 /* outer contour */ 01093 outer = Vect_new_line_struct(); 01094 Vect_get_area_points(Map, area, outer); 01095 /* does not work with zero length line segments */ 01096 Vect_line_prune(outer); 01097 01098 /* inner contours */ 01099 for (i = 0; i < n_isles; i++) { 01100 isle = Vect_get_area_isle(Map, area, i); 01101 Vect_get_isle_points(Map, isle, tPoints); 01102 01103 /* Check if the isle is big enough */ 01104 /* 01105 if (Vect_line_length(tPoints) < 2*PI*max) 01106 continue; 01107 */ 01108 /* does not work with zero length line segments */ 01109 Vect_line_prune(tPoints); 01110 add_line_to_array(tPoints, &isles, &isles_count, &isles_allocated, 01111 more); 01112 tPoints = Vect_new_line_struct(); 01113 } 01114 01115 buffer_lines(outer, isles, isles_count, 0, da, db, dalpha, round, caps, 01116 tol, oPoints, iPoints, inner_count); 01117 01118 Vect_destroy_line_struct(tPoints); 01119 Vect_destroy_line_struct(outer); 01120 destroy_lines_array(isles, isles_count); 01121 01122 return; 01123 } 01124 01137 void Vect_point_buffer2(double px, double py, double da, double db, 01138 double dalpha, int round, double tol, 01139 struct line_pnts **oPoints) 01140 { 01141 double tx, ty; 01142 double angular_tol, angular_step, phi1; 01143 int j, nsegments; 01144 01145 G_debug(2, "Vect_point_buffer()"); 01146 01147 *oPoints = Vect_new_line_struct(); 01148 01149 dalpha *= PI / 180; /* convert dalpha from degrees to radians */ 01150 01151 if (round || (!round)) { 01152 angular_tol = angular_tolerance(tol, da, db); 01153 01154 nsegments = (int)(2 * PI / angular_tol) + 1; 01155 angular_step = 2 * PI / nsegments; 01156 01157 phi1 = 0; 01158 for (j = 0; j < nsegments; j++) { 01159 elliptic_transform(cos(phi1), sin(phi1), da, db, dalpha, &tx, 01160 &ty); 01161 Vect_append_point(*oPoints, px + tx, py + ty, 0); 01162 phi1 += angular_step; 01163 } 01164 } 01165 else { 01166 01167 } 01168 01169 /* close the output line */ 01170 Vect_append_point(*oPoints, (*oPoints)->x[0], (*oPoints)->y[0], 01171 (*oPoints)->z[0]); 01172 01173 return; 01174 } 01175 01176 01177 /* 01178 \brief Create parrallel line 01179 01180 See also Vect_line_parallel(). 01181 01182 \param InPoints input line geometry 01183 \param da distance along major axis 01184 \param da distance along minor axis 01185 \param dalpha angle between 0x and major axis 01186 \param round make corners round 01187 \param tol maximum distance between theoretical arc and output segments 01188 \param[out] OutPoints output line 01189 */ 01190 void Vect_line_parallel2(struct line_pnts *InPoints, double da, double db, 01191 double dalpha, int side, int round, double tol, 01192 struct line_pnts *OutPoints) 01193 { 01194 G_debug(2, "Vect_line_parallel(): npoints = %d, da = %f, " 01195 "db = %f, dalpha = %f, side = %d, round_corners = %d, tol = %f", 01196 InPoints->n_points, da, db, dalpha, side, round, tol); 01197 01198 parallel_line(InPoints, da, db, dalpha, side, round, 1, NON_LOOPED_LINE, 01199 tol, OutPoints); 01200 01201 /* if (!loops) 01202 clean_parallel(OutPoints, InPoints, distance, rm_end); 01203 */ 01204 01205 return; 01206 }