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