buffer2.c

Go to the documentation of this file.
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     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  * (x,y) shoud be normalized vector for common transforms; This func transforms (x,y) to a vector corresponding to da, db, dalpha params
00073  * dalpha is in radians
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     /*    double cc = cosa*cosa;
00083        double ss = sina*sina;
00084        double t = (da-db)*sina*cosa;
00085 
00086        *nx = (da*cc + db*ss)*x + t*y;
00087        *ny = (da*ss + db*cc)*y + t*x;
00088        return; */
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  * vect(x,y) must be normalized
00101  * gives the tangent point of the tangent to ellpise(da,db,dalpha) parallel to vect(x,y)
00102  * dalpha is in radians
00103  * ellipse center is in (0,0)
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     /* 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     return;
00125 }
00126 
00127 
00128 /*
00129  * !!! This is not line in GRASS' sense. See http://en.wikipedia.org/wiki/Line_%28mathematics%29
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  * Finds intersection of two straight lines. Returns 0 if the lines are parallel, 1 if they cross,
00142  * 2 if they are the same line.
00143  * !!!!!!!!!!!!!!!! FIX THIS TOLLERANCE CONSTANTS BAD (and UGLY) CODE !!!!!!!!!
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  * This function generates parallel line (with loops, but not like the old ones).
00175  * It is not to be used directly for creating buffers.
00176  * + added elliptical buffers/par.lines support
00177  *
00178  * dalpha - direction of elliptical buffer major axis in degrees
00179  * da - distance along major axis
00180  * db: distance along minor (perp.) axis
00181  * side: side >= 0 - right side, side < 0 - left side
00182  * when (da == db) we have plain distances (old case)
00183  * round - 1 for round corners, 0 for sharp corners. (tol is used only if round == 1)
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         /* start point != end point */
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);    /* normalize variable */
00230     dalpha *= PI / 180;         /* convert dalpha from degrees to radians */
00231     angular_tol = angular_tolerance(tol, da, db);
00232 
00233     for (i = 0; i < np - 1; i++) {
00234         /* save the old values */
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         /* now delta_phi is in [-pi;pi] */
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);    /* nx == x[i] + vx, ny == y[i] + vy */
00283         }
00284         else if ((!round) || inner_corner) {
00285             res = line_intersection(a0, b0, c0, a1, b1, c1, &rx, &ry);
00286             /*                if (res == 0) {
00287                G_debug(4, "a0=%.18f, b0=%.18f, c0=%.18f, a1=%.18f, b1=%.18f, c1=%.18f", a0, b0, c0, a1, b1, c1);
00288                G_fatal_error("Two consequtive line segments are parallel, but not on one straight line! This should never happen.");
00289                return;
00290                }  */
00291             if (res == 1) {
00292                 if (!round)
00293                     Vect_append_point(nPoints, rx, ry, 0);
00294                 else {
00295                     /*                    d = dig_distance2_point_to_line(rx, ry, 0, x[i-1], y[i-1], 0, x[i], y[i], 0,
00296                        0, NULL, NULL, NULL, NULL, NULL);
00297                        if ( */
00298                     Vect_append_point(nPoints, rx, ry, 0);
00299                 }
00300             }
00301         }
00302         else {
00303             /* we should draw elliptical arc for outside corner */
00304 
00305             /* inverse transforms */
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             /* make delta_phi in [0, 2pi] */
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 /* input line must be looped */
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);    /* normalize variable */
00388     dalpha *= PI / 180;         /* convert dalpha from degrees to radians */
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         /* save the old values */
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         /* now delta_phi is in [-pi;pi] */
00433         turns360 = (fabs(fabs(delta_phi) - PI) < 1e-15);
00434         inner_corner = (side * delta_phi <= 0) && (!turns360);
00435 
00436 
00437         /* if <line turns 360> and (<caps> and <not round>) */
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);    /* nx == x[i] + vx, ny == y[i] + vy */
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                 /* no need to append point in this case */
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             /* we should draw elliptical arc for outside corner */
00465 
00466             /* inverse transforms */
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             /* make delta_phi in [0, 2pi] */
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     /* close the output line */
00499     Vect_append_point(nPoints, nPoints->x[0], nPoints->y[0], nPoints->z[0]);
00500     /*    Vect_line_prune ( nPoints ); */
00501 }
00502 
00503 /*
00504  * side: side >= 0 - extracts contour on right side of edge, side < 0 - extracts contour on left side of edge
00505  * if the extracted contour is the outer contour, it is returned in ccw order
00506  * else if it is inner contour, it is returned in cw order
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;                      /* current vertex number */
00515 
00516     int v0;
00517 
00518     int eside;                  /* side of the current edge */
00519 
00520     double eangle;              /* current edge angle with Ox (according to the current direction) */
00521 
00522     struct pg_vertex *vert;     /* current vertex */
00523 
00524     struct pg_vertex *vert0;    /* last vertex */
00525 
00526     struct pg_edge *edge;       /* current edge; must be edge of vert */
00527 
00528     /*    int cs; *//* on which side are we turning along the contour */
00529     /* we will always turn right and dont need that one */
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         /* mark current edge as visited on the appropriate side */
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             /* exclude current edge */
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                 /* now tangle is in (-PI, PI) */
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         //        G_debug(4, "ec: opt: side=%d opt_flag=%d opt_angle=%.18f opt_j=%d opt_step=%d", side, opt_flag, opt_angle, opt_j, opt_step);
00592 
00593         /* if line end is reached (no other edges at curr vertex) */
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;      /* the only edge of vert is vert->edges[0] */
00601                 opt_side = -eside;      /* go to the other side of the current edge */
00602                 G_debug(3, "    end has been reached, turning around");
00603             }
00604         }
00605 
00606         /* break condition */
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  * This function extracts the outer contour of a (self crossing) line.
00648  * It can generate left/right contour if none of the line ends are in a loop.
00649  * If one or both of them is in a loop, then there's only one contour
00650  * 
00651  * side: side > 0 - right contour, side < 0 - left contour, side = 0 - outer contour
00652  *       if side != 0 and there's only one contour, the function returns it
00653  * 
00654  * TODO: Implement side != 0 feature;
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     /* find a line segment which is on the outer contour */
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     /* the winding on the outer contour is 0 */
00699     extract_contour(pg, edge, (edge->v1 == v) ? RIGHT_SIDE : LEFT_SIDE, 0, 0,
00700                     nPoints);
00701 
00702     return;
00703 }
00704 
00705 /*
00706  * Extracts contours which are not visited.
00707  * IMPORTANT: the outer contour must be visited (you should call extract_outer_contour() to do that),
00708  *            so that extract_inner_contour() doesn't return it
00709  *
00710  * returns: 0 when there are no more inner contours; otherwise, 1
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 /* point_in_buf - test if point px,py is in d buffer of Points
00745  ** dalpha is in degrees
00746  ** returns:  1 in buffer
00747  **           0 not in buffer
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;         /* convert dalpha from degrees to radians */
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             /*            G_debug(4, "k = %g, k1 = %g", k, (mx * (px - vx) + my * (py - vy)) / (mx * mx + my * my)); */
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             /* inverse transform */
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             /*            G_debug(4, "sqrt(d)*da = %g, len' = %g, olen = %g", sqrt(d)*da, da*LENGTH(tx,ty), LENGTH((px-nx),(py-ny))); */
00806             if (d <= 1) {
00807                 //G_debug(1, "d=%g", d);
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             /*            G_debug(4, "sqrt(d)     = %g", sqrt(d)); */
00815             if (d <= da2) {
00816                 return 1;
00817             }
00818         }
00819     }
00820     return 0;
00821 }
00822 
00823 /* returns 0 for ccw, non-zero for cw
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 /* internal */
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 /* area_outer and area_isles[i] must be closed non self-intersecting lines
00872    side: 0 - auto, 1 - right, -1 left
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     /* initializations */
00903     sPoints = Vect_new_line_struct();
00904     cPoints = Vect_new_line_struct();
00905     arrPoints = NULL;
00906 
00907     /* outer contour */
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     /* inner contours */
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                 /* we need to check if the area is in the buffer.
00952                    I've simplfied convolution_line(), so that it runs faster,
00953                    however that leads to ocasional problems */
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     /* initializations */
01021     tPoints = Vect_new_line_struct();
01022     isles = NULL;
01023     pg = pg_create(Points);
01024 
01025     /* outer contour */
01026     outer = Vect_new_line_struct();
01027     extract_outer_contour(pg, 0, outer);
01028 
01029     /* inner contours */
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     /* initializations */
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     /* outer contour */
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     /* inner contours */
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         /* Check if the isle is big enough */
01100         /*
01101            if (Vect_line_length(tPoints) < 2*PI*max)
01102            continue;
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;         /* convert dalpha from degrees to radians */
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     /* close the output line */
01168     Vect_append_point(*oPoints, (*oPoints)->x[0], (*oPoints)->y[0],
01169                       (*oPoints)->z[0]);
01170 
01171     return;
01172 }
01173 
01174 
01175 /*
01176    \brief Create parrallel line
01177 
01178    See also Vect_line_parallel().
01179    
01180    \param InPoints input line geometry
01181    \param da distance along major axis
01182    \param da distance along minor axis
01183    \param dalpha angle between 0x and major axis
01184    \param round make corners round
01185    \param tol maximum distance between theoretical arc and output segments
01186    \param[out] OutPoints output line
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     /*    if (!loops)
01200        clean_parallel(OutPoints, InPoints, distance, rm_end);
01201      */
01202     
01203     return;
01204 }
Generated on Tue Apr 6 13:28:09 2010 for GRASS Programmer's Manual by  doxygen 1.6.3