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
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) {
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
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
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
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
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
00438
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) {
00447 rest = distance - dist + dxyz;
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
00467 if (angle != NULL)
00468 *angle = atan2(dy, dx);
00469
00470
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
00514 if (end < 0)
00515 return 0;
00516 length = Vect_line_length(InPoints);
00517 if (start > length)
00518 return 0;
00519
00520
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
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
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,
00743 double x2, double y2, double z2,
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 }