buffer2.c

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

Generated on Thu Jul 16 13:21:17 2009 for GRASS Programmer's Manual by  doxygen 1.5.6