GRASS Programmer's Manual 6.4.1(2011)
buffer.c
Go to the documentation of this file.
00001 
00020 #include <stdlib.h>
00021 #include <math.h>
00022 #include <grass/Vect.h>
00023 #include <grass/gis.h>
00024 
00025 #define LENGTH(DX, DY)  (  sqrt( (DX*DX)+(DY*DY) )  )
00026 #define PI M_PI
00027 
00028 /* vector() calculates normalized vector form two points */
00029 static void vect(double x1, double y1, double x2, double y2, double *x,
00030                  double *y)
00031 {
00032     double dx, dy, l;
00033 
00034     dx = x2 - x1;
00035     dy = y2 - y1;
00036     l = LENGTH(dx, dy);
00037     if (l == 0) {
00038         /* assume that dx == dy == 0, which should give (NaN,NaN) */
00039         /* without this, very small dx or dy could result in Infinity */
00040         dx = dy = 0;
00041     }
00042     *x = dx / l;
00043     *y = dy / l;
00044 }
00045 
00046 /* find_cross find first crossing between segments from s1 to s2 and from s3 to s4
00047  ** s5 is set to first segment and s6 to second
00048  ** neighbours are taken as crossing each other only if overlap
00049  ** returns: 1 found
00050  **         -1 found overlap
00051  **          0 not found
00052  */
00053 static int find_cross(struct line_pnts *Points, int s1, int s2, int s3,
00054                       int s4, int *s5, int *s6)
00055 {
00056     int i, j, np, ret;
00057     double *x, *y;
00058 
00059     G_debug(5,
00060             "find_cross(): npoints = %d, s1 = %d, s2 = %d, s3 = %d, s4 = %d",
00061             Points->n_points, s1, s2, s3, s4);
00062 
00063     x = Points->x;
00064     y = Points->y;
00065     np = Points->n_points;
00066 
00067     for (i = s1; i <= s2; i++) {
00068         for (j = s3; j <= s4; j++) {
00069             if (j == i) {
00070                 continue;
00071             }
00072             ret =
00073                 dig_test_for_intersection(x[i], y[i], x[i + 1], y[i + 1],
00074                                           x[j], y[j], x[j + 1], y[j + 1]);
00075             if (ret == 1 && ((i - j) > 1 || (i - j) < -1)) {
00076                 *s5 = i;
00077                 *s6 = j;
00078                 G_debug(5, "  intersection: s5 = %d, s6 = %d", *s5, *s6);
00079                 return 1;
00080             }
00081             if (ret == -1) {
00082                 *s5 = i;
00083                 *s6 = j;
00084                 G_debug(5, "  overlap: s5 = %d, s6 = %d", *s5, *s6);
00085                 return -1;
00086             }
00087         }
00088     }
00089     G_debug(5, "  no intersection");
00090     return 0;
00091 }
00092 
00093 /* point_in_buf - test if point px,py is in d buffer of Points
00094  ** returns:  1 in buffer
00095  **           0 not  in buffer
00096  */
00097 static int point_in_buf(struct line_pnts *Points, double px, double py,
00098                         double d)
00099 {
00100     int i, np;
00101     double sd;
00102 
00103     np = Points->n_points;
00104     d *= d;
00105     for (i = 0; i < np - 1; i++) {
00106         sd = dig_distance2_point_to_line(px, py, 0,
00107                                          Points->x[i], Points->y[i], 0,
00108                                          Points->x[i + 1], Points->y[i + 1],
00109                                          0, 0, NULL, NULL, NULL, NULL, NULL);
00110         if (sd <= d) {
00111             return 1;
00112         }
00113     }
00114     return 0;
00115 }
00116 
00117 /* clean_parallel - clean parallel line created by parallel_line:
00118  ** - looking for loops and if loop doesn't contain any other loop
00119  **   and centroid of loop is in buffer removes this loop (repeated)
00120  ** - optionally removes all end points in buffer
00121  *    parameters:
00122  *      Points - parallel line
00123  *      origPoints - original line
00124  *      d - offset
00125  *      rm_end - remove end points in buffer
00126  ** note1: on some lines (multiply selfcrossing; lines with end points
00127  **        in buffer of line other; some shapes of ends ) may create nosense
00128  ** note2: this function is stupid and slow, somebody more clever
00129  **        than I am should write paralle_line + clean_parallel
00130  **        better;    RB March 2000
00131  */
00132 static void clean_parallel(struct line_pnts *Points,
00133                            struct line_pnts *origPoints, double d, int rm_end)
00134 {
00135     int i, j, np, npn, sa, sb;
00136     int sa_max = 0;
00137     int first = 0, current, last, lcount;
00138     double *x, *y, px, py, ix, iy;
00139     static struct line_pnts *sPoints = NULL;
00140 
00141     G_debug(4, "clean_parallel(): npoints = %d, d = %f, rm_end = %d",
00142             Points->n_points, d, rm_end);
00143 
00144     x = Points->x;
00145     y = Points->y;
00146     np = Points->n_points;
00147 
00148     if (sPoints == NULL)
00149         sPoints = Vect_new_line_struct();
00150 
00151     Vect_reset_line(sPoints);
00152 
00153     npn = 1;
00154 
00155     /* remove loops */
00156     while (first < np - 2) {
00157         /* find first loop which doesn't contain any other loop */
00158         current = first;
00159         last = Points->n_points - 2;
00160         lcount = 0;
00161         while (find_cross
00162                (Points, current, last - 1, current + 1, last, &sa,
00163                 &sb) != 0) {
00164             if (lcount == 0) {
00165                 first = sa;
00166             }                   /* move first forward */
00167 
00168             current = sa + 1;
00169             last = sb;
00170             lcount++;
00171             G_debug(5, "  current = %d, last = %d, lcount = %d", current,
00172                     last, lcount);
00173         }
00174         if (lcount == 0) {
00175             break;
00176         }                       /* loop not found */
00177 
00178         /* ensure sa is monotonically increasing, so npn doesn't reset low */
00179         if (sa > sa_max)
00180             sa_max = sa;
00181         if (sa < sa_max)
00182             break;
00183 
00184         /* remove loop if in buffer */
00185         if ((sb - sa) == 1) {   /* neighbouring lines overlap */
00186             j = sb + 1;
00187             npn = sa + 1;
00188         }
00189         else {
00190             Vect_reset_line(sPoints);
00191             dig_find_intersection(x[sa], y[sa], x[sa + 1], y[sa + 1], x[sb],
00192                                   y[sb], x[sb + 1], y[sb + 1], &ix, &iy);
00193             Vect_append_point(sPoints, ix, iy, 0);
00194             for (i = sa + 1; i < sb + 1; i++) { /* create loop polygon */
00195                 Vect_append_point(sPoints, x[i], y[i], 0);
00196             }
00197             Vect_find_poly_centroid(sPoints, &px, &py);
00198             if (point_in_buf(origPoints, px, py, d)) {  /* is loop in buffer ? */
00199                 npn = sa + 1;
00200                 x[npn] = ix;
00201                 y[npn] = iy;
00202                 j = sb + 1;
00203                 npn++;
00204                 if (lcount == 0) {
00205                     first = sb;
00206                 }
00207             }
00208             else {              /* loop is not in buffer */
00209                 first = sb;
00210                 continue;
00211             }
00212         }
00213 
00214         for (i = j; i < Points->n_points; i++) {        /* move points down */
00215             x[npn] = x[i];
00216             y[npn] = y[i];
00217             npn++;
00218         }
00219         Points->n_points = npn;
00220     }
00221 
00222     if (rm_end) {
00223         /* remove points from start in buffer */
00224         j = 0;
00225         for (i = 0; i < Points->n_points - 1; i++) {
00226             px = (x[i] + x[i + 1]) / 2;
00227             py = (y[i] + y[i + 1]) / 2;
00228             if (point_in_buf(origPoints, x[i], y[i], d * 0.9999)
00229                 && point_in_buf(origPoints, px, py, d * 0.9999)) {
00230                 j++;
00231             }
00232             else {
00233                 break;
00234             }
00235         }
00236         if (j > 0) {
00237             npn = 0;
00238             for (i = j; i < Points->n_points; i++) {
00239                 x[npn] = x[i];
00240                 y[npn] = y[i];
00241                 npn++;
00242             }
00243             Points->n_points = npn;
00244         }
00245         /* remove points from end in buffer */
00246         j = 0;
00247         for (i = Points->n_points - 1; i >= 1; i--) {
00248             px = (x[i] + x[i - 1]) / 2;
00249             py = (y[i] + y[i - 1]) / 2;
00250             if (point_in_buf(origPoints, x[i], y[i], d * 0.9999)
00251                 && point_in_buf(origPoints, px, py, d * 0.9999)) {
00252                 j++;
00253             }
00254             else {
00255                 break;
00256             }
00257         }
00258         if (j > 0) {
00259             Points->n_points -= j;
00260         }
00261     }
00262 }
00263 
00264 /* parallel_line - remove duplicate points from input line and
00265  *  creates new parallel line in 'd' offset distance;
00266  *  'tol' is tolerance between arc and polyline;
00267  *  this function doesn't care about created loops;
00268  *
00269  *  New line is written to existing nPoints structure.
00270  */
00271 static void parallel_line(struct line_pnts *Points, double d, double tol,
00272                           struct line_pnts *nPoints)
00273 {
00274     int i, j, np, na, side;
00275     double *x, *y, nx, ny, tx, ty, vx, vy, ux, uy, wx, wy;
00276     double atol, atol2, a, av, aw;
00277 
00278     G_debug(4, "parallel_line()");
00279 
00280     Vect_reset_line(nPoints);
00281 
00282     Vect_line_prune(Points);
00283     np = Points->n_points;
00284     x = Points->x;
00285     y = Points->y;
00286 
00287     if (np == 0)
00288         return;
00289 
00290     if (np == 1) {
00291         Vect_append_point(nPoints, x[0], y[0], 0);      /* ? OK, should make circle for points ? */
00292         return;
00293     }
00294 
00295     if (d == 0) {
00296         Vect_copy_xyz_to_pnts(nPoints, x, y, NULL, np);
00297         return;
00298     }
00299 
00300     side = (int)(d / fabs(d));
00301     atol = 2 * acos(1 - tol / fabs(d));
00302 
00303     for (i = 0; i < np - 1; i++) {
00304         vect(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);
00305         vx = ty * d;
00306         vy = -tx * d;
00307 
00308         nx = x[i] + vx;
00309         ny = y[i] + vy;
00310         Vect_append_point(nPoints, nx, ny, 0);
00311 
00312         nx = x[i + 1] + vx;
00313         ny = y[i + 1] + vy;
00314         Vect_append_point(nPoints, nx, ny, 0);
00315 
00316         if (i < np - 2) {       /* use polyline instead of arc between line segments */
00317             vect(x[i + 1], y[i + 1], x[i + 2], y[i + 2], &ux, &uy);
00318             wx = uy * d;
00319             wy = -ux * d;
00320             av = atan2(vy, vx);
00321             aw = atan2(wy, wx);
00322             a = (aw - av) * side;
00323             if (a < 0)
00324                 a += 2 * PI;
00325 
00326             /* TODO: a <= PI can probably fail because of representation error */
00327             if (a <= PI && a > atol) {
00328                 na = (int)(a / atol);
00329                 atol2 = a / (na + 1) * side;
00330                 for (j = 0; j < na; j++) {
00331                     av += atol2;
00332                     nx = x[i + 1] + fabs(d) * cos(av);
00333                     ny = y[i + 1] + fabs(d) * sin(av);
00334                     Vect_append_point(nPoints, nx, ny, 0);
00335                 }
00336             }
00337         }
00338     }
00339     Vect_line_prune(nPoints);
00340 }
00341 
00353 void
00354 Vect_line_parallel(struct line_pnts *InPoints, double distance,
00355                    double tolerance, int rm_end, struct line_pnts *OutPoints)
00356 {
00357     G_debug(4,
00358             "Vect_line_parallel(): npoints = %d, distance = %f, tolerance = %f",
00359             InPoints->n_points, distance, tolerance);
00360 
00361     parallel_line(InPoints, distance, tolerance, OutPoints);
00362 
00363     clean_parallel(OutPoints, InPoints, distance, rm_end);
00364 
00365     return;
00366 }
00367 
00379 void
00380 Vect_line_buffer(struct line_pnts *InPoints, double distance,
00381                  double tolerance, struct line_pnts *OutPoints)
00382 {
00383     double dangle;
00384     int side, npoints;
00385     static struct line_pnts *Points = NULL;
00386     static struct line_pnts *PPoints = NULL;
00387 
00388     distance = fabs(distance);
00389 
00390     dangle = 2 * acos(1 - tolerance / fabs(distance));  /* angle step */
00391 
00392     if (Points == NULL)
00393         Points = Vect_new_line_struct();
00394 
00395     if (PPoints == NULL)
00396         PPoints = Vect_new_line_struct();
00397 
00398     /* Copy and prune input */
00399     Vect_reset_line(Points);
00400     Vect_append_points(Points, InPoints, GV_FORWARD);
00401     Vect_line_prune(Points);
00402 
00403     Vect_reset_line(OutPoints);
00404 
00405     npoints = Points->n_points;
00406     if (npoints <= 0) {
00407         return;
00408     }
00409     else if (npoints == 1) {    /* make a circle */
00410         double angle, x, y;
00411 
00412         for (angle = 0; angle < 2 * PI; angle += dangle) {
00413             x = Points->x[0] + distance * cos(angle);
00414             y = Points->y[0] + distance * sin(angle);
00415             Vect_append_point(OutPoints, x, y, 0);
00416         }
00417         /* Close polygon */
00418         Vect_append_point(OutPoints, OutPoints->x[0], OutPoints->y[0], 0);
00419     }
00420     else {                      /* 2 and more points */
00421         for (side = 0; side < 2; side++) {
00422             double angle, sangle;
00423             double lx1, ly1, lx2, ly2;
00424             double x, y, nx, ny, sx, sy, ex, ey;
00425 
00426             /* Parallel on one side */
00427             if (side == 0) {
00428                 Vect_line_parallel(Points, distance, tolerance, 0, PPoints);
00429                 Vect_append_points(OutPoints, PPoints, GV_FORWARD);
00430             }
00431             else {
00432                 Vect_line_parallel(Points, -distance, tolerance, 0, PPoints);
00433                 Vect_append_points(OutPoints, PPoints, GV_BACKWARD);
00434             }
00435 
00436             /* Arc at the end */
00437             /* 2 points at theend of original line */
00438             if (side == 0) {
00439                 lx1 = Points->x[npoints - 2];
00440                 ly1 = Points->y[npoints - 2];
00441                 lx2 = Points->x[npoints - 1];
00442                 ly2 = Points->y[npoints - 1];
00443             }
00444             else {
00445                 lx1 = Points->x[1];
00446                 ly1 = Points->y[1];
00447                 lx2 = Points->x[0];
00448                 ly2 = Points->y[0];
00449             }
00450 
00451             /* normalized vector */
00452             vect(lx1, ly1, lx2, ly2, &nx, &ny);
00453 
00454             /* starting point */
00455             sangle = atan2(-nx, ny);    /* starting angle */
00456             sx = lx2 + ny * distance;
00457             sy = ly2 - nx * distance;
00458 
00459             /* end point */
00460             ex = lx2 - ny * distance;
00461             ey = ly2 + nx * distance;
00462 
00463             Vect_append_point(OutPoints, sx, sy, 0);
00464 
00465             /* arc */
00466             for (angle = dangle; angle < PI; angle += dangle) {
00467                 x = lx2 + distance * cos(sangle + angle);
00468                 y = ly2 + distance * sin(sangle + angle);
00469                 Vect_append_point(OutPoints, x, y, 0);
00470             }
00471 
00472             Vect_append_point(OutPoints, ex, ey, 0);
00473         }
00474 
00475         /* Close polygon */
00476         Vect_append_point(OutPoints, OutPoints->x[0], OutPoints->y[0], 0);
00477     }
00478     Vect_line_prune(OutPoints);
00479 
00480     return;
00481 }
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Defines