line.c

Go to the documentation of this file.
00001 
00019 #include <stdlib.h>
00020 #include <math.h>
00021 #include <grass/gis.h>
00022 #include <grass/Vect.h>
00023 #include <grass/glocale.h>
00024 
00040 struct line_pnts *Vect__new_line_struct(void);
00041 
00055 struct line_pnts *Vect_new_line_struct()
00056 {
00057     struct line_pnts *p;
00058 
00059     if (NULL == (p = Vect__new_line_struct()))
00060         G_fatal_error("Vect_new_line_struct(): %s", _("Out of memory"));
00061 
00062     return p;
00063 }
00064 
00065 struct line_pnts *Vect__new_line_struct()
00066 {
00067     struct line_pnts *p;
00068 
00069     p = (struct line_pnts *)malloc(sizeof(struct line_pnts));
00070 
00071     /* alloc_points MUST be initialized to zero */
00072     if (p)
00073         p->alloc_points = p->n_points = 0;
00074 
00075     if (p)
00076         p->x = p->y = p->z = NULL;
00077 
00078     return p;
00079 }
00080 
00088 int Vect_destroy_line_struct(struct line_pnts *p)
00089 {
00090     if (p) {                    /* probably a moot test */
00091         if (p->alloc_points) {
00092             G_free((char *)p->x);
00093             G_free((char *)p->y);
00094             G_free((char *)p->z);
00095         }
00096         G_free((char *)p);
00097     }
00098 
00099     return 0;
00100 }
00101 
00112 int
00113 Vect_copy_xyz_to_pnts(struct line_pnts *Points, double *x, double *y,
00114                       double *z, int n)
00115 {
00116     register int i;
00117 
00118     if (0 > dig_alloc_points(Points, n))
00119         return (-1);
00120 
00121     for (i = 0; i < n; i++) {
00122         Points->x[i] = x[i];
00123         Points->y[i] = y[i];
00124         if (z != NULL)
00125             Points->z[i] = z[i];
00126         else
00127             Points->z[i] = 0;
00128         Points->n_points = n;
00129     }
00130 
00131     return (0);
00132 }
00133 
00134 
00146 int Vect_reset_line(struct line_pnts *Points)
00147 {
00148     Points->n_points = 0;
00149 
00150     return 0;
00151 }
00152 
00166 int Vect_append_point(struct line_pnts *Points, double x, double y, double z)
00167 {
00168     register int n;
00169 
00170     if (0 > dig_alloc_points(Points, Points->n_points + 1))
00171         return (-1);
00172 
00173     n = Points->n_points;
00174     Points->x[n] = x;
00175     Points->y[n] = y;
00176     Points->z[n] = z;
00177 
00178     return ++(Points->n_points);
00179 }
00180 
00191 int
00192 Vect_line_insert_point(struct line_pnts *Points, int index, double x,
00193                        double y, double z)
00194 {
00195     register int n;
00196 
00197     if (index < 0 || index > Points->n_points - 1)
00198         G_fatal_error("%s Vect_line_insert_point()",
00199                       _("Index out of range in"));
00200 
00201     if (0 > dig_alloc_points(Points, Points->n_points + 1))
00202         return (-1);
00203 
00204     /* move up */
00205     for (n = Points->n_points; n > index; n--) {
00206         Points->x[n] = Points->x[n - 1];
00207         Points->y[n] = Points->y[n - 1];
00208         Points->z[n] = Points->z[n - 1];
00209     }
00210 
00211     Points->x[index] = x;
00212     Points->y[index] = y;
00213     Points->z[index] = z;
00214     return ++(Points->n_points);
00215 }
00216 
00225 int Vect_line_delete_point(struct line_pnts *Points, int index)
00226 {
00227     register int n;
00228 
00229     if (index < 0 || index > Points->n_points - 1)
00230         G_fatal_error("%s Vect_line_insert_point()",
00231                       _("Index out of range in"));
00232 
00233     if (Points->n_points == 0)
00234         return 0;
00235 
00236     /* move down */
00237     for (n = index; n < Points->n_points - 1; n++) {
00238         Points->x[n] = Points->x[n + 1];
00239         Points->y[n] = Points->y[n + 1];
00240         Points->z[n] = Points->z[n + 1];
00241     }
00242 
00243     return --(Points->n_points);
00244 }
00245 
00253 int Vect_line_prune(struct line_pnts *Points)
00254 {
00255     int i, j;
00256 
00257     if (Points->n_points > 0) {
00258         j = 1;
00259         for (i = 1; i < Points->n_points; i++) {
00260             if (Points->x[i] != Points->x[j - 1] ||
00261                 Points->y[i] != Points->y[j - 1]
00262                 || Points->z[i] != Points->z[j - 1]) {
00263                 Points->x[j] = Points->x[i];
00264                 Points->y[j] = Points->y[i];
00265                 Points->z[j] = Points->z[i];
00266                 j++;
00267             }
00268         }
00269         Points->n_points = j;
00270     }
00271 
00272     return (Points->n_points);
00273 }
00274 
00283 int Vect_line_prune_thresh(struct line_pnts *Points, double threshold)
00284 {
00285     int ret;
00286 
00287     ret = dig_prune(Points, threshold);
00288 
00289     if (ret < Points->n_points)
00290         Points->n_points = ret;
00291 
00292     return (Points->n_points);
00293 }
00294 
00309 int
00310 Vect_append_points(struct line_pnts *Points, struct line_pnts *APoints,
00311                    int direction)
00312 {
00313     int i, n, on, an;
00314 
00315     on = Points->n_points;
00316     an = APoints->n_points;
00317     n = on + an;
00318 
00319     /* Should be OK, dig_alloc_points calls realloc */
00320     if (0 > dig_alloc_points(Points, n))
00321         return (-1);
00322 
00323     if (direction == GV_FORWARD) {
00324         for (i = 0; i < an; i++) {
00325             Points->x[on + i] = APoints->x[i];
00326             Points->y[on + i] = APoints->y[i];
00327             Points->z[on + i] = APoints->z[i];
00328         }
00329     }
00330     else {
00331         for (i = 0; i < an; i++) {
00332             Points->x[on + i] = APoints->x[an - i - 1];
00333             Points->y[on + i] = APoints->y[an - i - 1];
00334             Points->z[on + i] = APoints->z[an - i - 1];
00335         }
00336     }
00337 
00338     Points->n_points = n;
00339     return n;
00340 }
00341 
00342 
00356 int
00357 Vect_copy_pnts_to_xyz(struct line_pnts *Points, double *x, double *y,
00358                       double *z, int *n)
00359 {
00360     register int i;
00361 
00362     for (i = 0; i < *n; i++) {
00363         x[i] = Points->x[i];
00364         y[i] = Points->y[i];
00365         if (z != NULL)
00366             z[i] = Points->z[i];
00367         *n = Points->n_points;
00368     }
00369 
00370     return (Points->n_points);
00371 }
00372 
00390 int
00391 Vect_point_on_line(struct line_pnts *Points, double distance,
00392                    double *x, double *y, double *z, double *angle,
00393                    double *slope)
00394 {
00395     int j, np, seg = 0;
00396     double dist = 0, length;
00397     double xp = 0, yp = 0, zp = 0, dx = 0, dy = 0, dz = 0, dxy =
00398         0, dxyz, k, rest;
00399 
00400     G_debug(3, "Vect_point_on_line(): distance = %f", distance);
00401     if ((distance < 0) || (Points->n_points < 2))
00402         return 0;
00403 
00404     /* Check if first or last */
00405     length = Vect_line_length(Points);
00406     G_debug(3, "  length = %f", length);
00407     if (distance < 0 || distance > length) {
00408         G_debug(3, "  -> outside line");
00409         return 0;
00410     }
00411 
00412     np = Points->n_points;
00413     if (distance == 0) {
00414         G_debug(3, "  -> first point");
00415         xp = Points->x[0];
00416         yp = Points->y[0];
00417         zp = Points->z[0];
00418         dx = Points->x[1] - Points->x[0];
00419         dy = Points->y[1] - Points->y[0];
00420         dz = Points->z[1] - Points->z[0];
00421         dxy = hypot(dx, dy);
00422         seg = 1;
00423     }
00424     else if (distance == length) {
00425         G_debug(3, "  -> last point");
00426         xp = Points->x[np - 1];
00427         yp = Points->y[np - 1];
00428         zp = Points->z[np - 1];
00429         dx = Points->x[np - 1] - Points->x[np - 2];
00430         dy = Points->y[np - 1] - Points->y[np - 2];
00431         dz = Points->z[np - 1] - Points->z[np - 2];
00432         dxy = hypot(dx, dy);
00433         seg = np - 1;
00434     }
00435     else {
00436         for (j = 0; j < Points->n_points - 1; j++) {
00437             /* dxyz = G_distance(Points->x[j], Points->y[j], 
00438                Points->x[j+1], Points->y[j+1]); */
00439             dx = Points->x[j + 1] - Points->x[j];
00440             dy = Points->y[j + 1] - Points->y[j];
00441             dz = Points->z[j + 1] - Points->z[j];
00442             dxy = hypot(dx, dy);
00443             dxyz = hypot(dxy, dz);
00444 
00445             dist += dxyz;
00446             if (dist >= distance) {     /* point is on the current line part */
00447                 rest = distance - dist + dxyz;  /* from first point of segment to point */
00448                 k = rest / dxyz;
00449 
00450                 xp = Points->x[j] + k * dx;
00451                 yp = Points->y[j] + k * dy;
00452                 zp = Points->z[j] + k * dz;
00453                 seg = j + 1;
00454                 break;
00455             }
00456         }
00457     }
00458 
00459     if (x != NULL)
00460         *x = xp;
00461     if (y != NULL)
00462         *y = yp;
00463     if (z != NULL)
00464         *z = zp;
00465 
00466     /* calculate angle */
00467     if (angle != NULL)
00468         *angle = atan2(dy, dx);
00469 
00470     /* calculate slope */
00471     if (slope != NULL)
00472         *slope = atan2(dz, dxy);
00473 
00474     return seg;
00475 }
00476 
00494 int
00495 Vect_line_segment(struct line_pnts *InPoints, double start, double end,
00496                   struct line_pnts *OutPoints)
00497 {
00498     int i, seg1, seg2;
00499     double length, tmp;
00500     double x1, y1, z1, x2, y2, z2;
00501 
00502     G_debug(3, "Vect_line_segment(): start = %f, end = %f, n_points = %d",
00503             start, end, InPoints->n_points);
00504 
00505     Vect_reset_line(OutPoints);
00506 
00507     if (start > end) {
00508         tmp = start;
00509         start = end;
00510         end = tmp;
00511     }
00512 
00513     /* Check start/end */
00514     if (end < 0)
00515         return 0;
00516     length = Vect_line_length(InPoints);
00517     if (start > length)
00518         return 0;
00519 
00520     /* Find coordinates and segments of start/end */
00521     seg1 = Vect_point_on_line(InPoints, start, &x1, &y1, &z1, NULL, NULL);
00522     seg2 = Vect_point_on_line(InPoints, end, &x2, &y2, &z2, NULL, NULL);
00523 
00524     G_debug(3, "  -> seg1 = %d seg2 = %d", seg1, seg2);
00525 
00526     if (seg1 == 0 || seg2 == 0) {
00527         G_warning(_("Segment outside line, no segment created"));
00528         return 0;
00529     }
00530 
00531     Vect_append_point(OutPoints, x1, y1, z1);
00532 
00533     for (i = seg1; i < seg2; i++) {
00534         Vect_append_point(OutPoints, InPoints->x[i], InPoints->y[i],
00535                           InPoints->z[i]);
00536     };
00537 
00538     Vect_append_point(OutPoints, x2, y2, z2);
00539 
00540     return 1;
00541 }
00542 
00552 double Vect_line_length(struct line_pnts *Points)
00553 {
00554     int j;
00555     double dx, dy, dz, len = 0;
00556 
00557     if (Points->n_points < 2)
00558         return 0;
00559 
00560     for (j = 0; j < Points->n_points - 1; j++) {
00561         dx = Points->x[j + 1] - Points->x[j];
00562         dy = Points->y[j + 1] - Points->y[j];
00563         dz = Points->z[j + 1] - Points->z[j];
00564         len += hypot(hypot(dx, dy), dz);
00565     }
00566 
00567     return len;
00568 }
00569 
00570 
00580 double Vect_line_geodesic_length(struct line_pnts *Points)
00581 {
00582     int j, dc;
00583     double dx, dy, dz, dxy, len = 0;
00584 
00585     dc = G_begin_distance_calculations();
00586 
00587     if (Points->n_points < 2)
00588         return 0;
00589 
00590     for (j = 0; j < Points->n_points - 1; j++) {
00591         if (dc == 2)
00592             dxy =
00593                 G_geodesic_distance(Points->x[j], Points->y[j],
00594                                     Points->x[j + 1], Points->y[j + 1]);
00595         else {
00596             dx = Points->x[j + 1] - Points->x[j];
00597             dy = Points->y[j + 1] - Points->y[j];
00598             dxy = hypot(dx, dy);
00599         }
00600 
00601         dz = Points->z[j + 1] - Points->z[j];
00602         len += hypot(dxy, dz);
00603     }
00604 
00605     return len;
00606 }
00607 
00627 int
00628 Vect_line_distance(struct line_pnts *points,
00629                    double ux, double uy, double uz,
00630                    int with_z,
00631                    double *px, double *py, double *pz,
00632                    double *dist, double *spdist, double *lpdist)
00633 {
00634     register int i;
00635     register double distance;
00636     register double new_dist;
00637     double tpx, tpy, tpz, tdist, tspdist, tlpdist = 0;
00638     double dx, dy, dz;
00639     register int n_points;
00640     int segment;
00641 
00642     n_points = points->n_points;
00643 
00644     if (n_points == 1) {
00645         distance =
00646             dig_distance2_point_to_line(ux, uy, uz, points->x[0],
00647                                         points->y[0], points->z[0],
00648                                         points->x[0], points->y[0],
00649                                         points->z[0], with_z, NULL, NULL,
00650                                         NULL, NULL, NULL);
00651         tpx = points->x[0];
00652         tpy = points->y[0];
00653         tpz = points->z[0];
00654         tdist = sqrt(distance);
00655         tspdist = 0;
00656         tlpdist = 0;
00657         segment = 0;
00658 
00659     }
00660     else {
00661 
00662         distance =
00663             dig_distance2_point_to_line(ux, uy, uz, points->x[0],
00664                                         points->y[0], points->z[0],
00665                                         points->x[1], points->y[1],
00666                                         points->z[1], with_z, NULL, NULL,
00667                                         NULL, NULL, NULL);
00668         segment = 1;
00669 
00670         for (i = 1; i < n_points - 1; i++) {
00671             new_dist = dig_distance2_point_to_line(ux, uy, uz,
00672                                                    points->x[i], points->y[i],
00673                                                    points->z[i],
00674                                                    points->x[i + 1],
00675                                                    points->y[i + 1],
00676                                                    points->z[i + 1], with_z,
00677                                                    NULL, NULL, NULL, NULL,
00678                                                    NULL);
00679             if (new_dist < distance) {
00680                 distance = new_dist;
00681                 segment = i + 1;
00682             }
00683         }
00684 
00685         /* we have segment and now we can recalculate other values (speed) */
00686         new_dist = dig_distance2_point_to_line(ux, uy, uz,
00687                                                points->x[segment - 1],
00688                                                points->y[segment - 1],
00689                                                points->z[segment - 1],
00690                                                points->x[segment],
00691                                                points->y[segment],
00692                                                points->z[segment], with_z,
00693                                                &tpx, &tpy, &tpz, &tspdist,
00694                                                NULL);
00695 
00696         /* calculate distance from beginning of line */
00697         if (lpdist) {
00698             tlpdist = 0;
00699             for (i = 0; i < segment - 1; i++) {
00700                 dx = points->x[i + 1] - points->x[i];
00701                 dy = points->y[i + 1] - points->y[i];
00702                 if (with_z)
00703                     dz = points->z[i + 1] - points->z[i];
00704                 else
00705                     dz = 0;
00706 
00707                 tlpdist += hypot(hypot(dx, dy), dz);
00708             }
00709             tlpdist += tspdist;
00710         }
00711         tdist = sqrt(distance);
00712     }
00713 
00714     if (px)
00715         *px = tpx;
00716     if (py)
00717         *py = tpy;
00718     if (pz && with_z)
00719         *pz = tpz;
00720     if (dist)
00721         *dist = tdist;
00722     if (spdist)
00723         *spdist = tspdist;
00724     if (lpdist)
00725         *lpdist = tlpdist;
00726 
00727     return (segment);
00728 }
00729 
00730 
00742 double Vect_points_distance(double x1, double y1, double z1,    /* point 1 */
00743                             double x2, double y2, double z2,    /* point 2 */
00744                             int with_z)
00745 {
00746     double dx, dy, dz;
00747 
00748 
00749     dx = x2 - x1;
00750     dy = y2 - y1;
00751     dz = z2 - z1;
00752 
00753     if (with_z)
00754         return hypot(hypot(dx, dy), dz);
00755     else
00756         return hypot(dx, dy);
00757 
00758 }
00759 
00768 int Vect_line_box(struct line_pnts *Points, BOUND_BOX * Box)
00769 {
00770     dig_line_box(Points, Box);
00771     return 0;
00772 }
00773 
00781 void Vect_line_reverse(struct line_pnts *Points)
00782 {
00783     int i, j, np;
00784     double x, y, z;
00785 
00786     np = (int)Points->n_points / 2;
00787 
00788     for (i = 0; i < np; i++) {
00789         j = Points->n_points - i - 1;
00790         x = Points->x[i];
00791         y = Points->y[i];
00792         z = Points->z[i];
00793         Points->x[i] = Points->x[j];
00794         Points->y[i] = Points->y[j];
00795         Points->z[i] = Points->z[j];
00796         Points->x[j] = x;
00797         Points->y[j] = y;
00798         Points->z[j] = z;
00799     }
00800 }
00801 
00812 int Vect_get_line_cat(struct Map_info *Map, int line, int field)
00813 {
00814 
00815     static struct line_cats *cats = NULL;
00816     int cat, ltype;
00817 
00818     if (cats == NULL)
00819         cats = Vect_new_cats_struct();
00820 
00821     ltype = Vect_read_line(Map, NULL, cats, line);
00822     Vect_cat_get(cats, field, &cat);
00823     G_debug(3, "Vect_get_line_cat: display line %d, ltype %d, cat %d", line,
00824             ltype, cat);
00825 
00826     return cat;
00827 }