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