GRASS Programmer's Manual 6.4.1(2011)
|
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 }