00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include <stdlib.h>
00019 #include <math.h>
00020 #include <grass/gis.h>
00021 #include <grass/Vect.h>
00022
00023 struct line_pnts *Vect__new_line_struct (void);
00024
00034 struct line_pnts *
00035 Vect_new_line_struct ()
00036 {
00037 struct line_pnts *p;
00038
00039 if (NULL == (p = Vect__new_line_struct ()))
00040 G_fatal_error ("New_line: Out of memory");
00041
00042 return p;
00043 }
00044
00045 struct line_pnts *
00046 Vect__new_line_struct ()
00047 {
00048 struct line_pnts *p;
00049
00050 p = (struct line_pnts *) malloc (sizeof (struct line_pnts));
00051
00052
00053 if (p)
00054 p->alloc_points = p->n_points = 0;
00055
00056 if (p)
00057 p->x = p->y = p->z = NULL;
00058
00059 return p;
00060 }
00061
00068 int
00069 Vect_destroy_line_struct (struct line_pnts *p)
00070 {
00071 if (p)
00072 {
00073 if (p->alloc_points)
00074 {
00075 free ((char *) p->x);
00076 free ((char *) p->y);
00077 free ((char *) p->z);
00078 }
00079 free ((char *) p);
00080 }
00081
00082 return 0;
00083 }
00084
00092 int
00093 Vect_copy_xyz_to_pnts (
00094 struct line_pnts *Points, double *x, double *y, double *z, int n)
00095 {
00096 register int i;
00097
00098 if (0 > dig_alloc_points (Points, n))
00099 return (-1);
00100
00101 for (i = 0; i < n; i++)
00102 {
00103 Points->x[i] = x[i];
00104 Points->y[i] = y[i];
00105 if ( z != NULL )
00106 Points->z[i] = z[i];
00107 else
00108 Points->z[i] = 0;
00109 Points->n_points = n;
00110 }
00111
00112 return (0);
00113 }
00114
00115
00124
00125 int
00126 Vect_reset_line (struct line_pnts *Points)
00127 {
00128 Points->n_points = 0;
00129
00130 return 0;
00131 }
00132
00142 int
00143 Vect_append_point (struct line_pnts *Points, double x, double y, double z)
00144 {
00145 register int n;
00146
00147 if (0 > dig_alloc_points (Points, Points->n_points + 1))
00148 return (-1);
00149
00150 n = Points->n_points;
00151 Points->x[n] = x;
00152 Points->y[n] = y;
00153 Points->z[n] = z;
00154 return ++(Points->n_points);
00155 }
00156
00167 int
00168 Vect_line_insert_point (struct line_pnts *Points, int index, double x, double y, double z)
00169 {
00170 register int n;
00171
00172 if ( index < 0 || index > Points->n_points-1 )
00173 G_fatal_error ("Index out of range in Vect_line_insert_point()");
00174
00175 if (0 > dig_alloc_points (Points, Points->n_points + 1))
00176 return (-1);
00177
00178
00179 for ( n = Points->n_points; n > index; n-- ) {
00180 Points->x[n] = Points->x[n-1];
00181 Points->y[n] = Points->y[n-1];
00182 Points->z[n] = Points->z[n-1];
00183 }
00184
00185 Points->x[index] = x;
00186 Points->y[index] = y;
00187 Points->z[index] = z;
00188 return ++(Points->n_points);
00189 }
00190
00198 int
00199 Vect_line_delete_point (struct line_pnts *Points, int index)
00200 {
00201 register int n;
00202
00203 if ( index < 0 || index > Points->n_points-1 )
00204 G_fatal_error ("Index out of range in Vect_line_delete_point()");
00205
00206 if ( Points->n_points == 0 )
00207 return 0;
00208
00209
00210 for ( n = index; n < Points->n_points - 1; n++ ) {
00211 Points->x[n] = Points->x[n+1];
00212 Points->y[n] = Points->y[n+1];
00213 Points->z[n] = Points->z[n+1];
00214 }
00215
00216 return --(Points->n_points);
00217 }
00218
00225 int
00226 Vect_line_prune (struct line_pnts *Points)
00227 {
00228 int i, j;
00229
00230 if ( Points->n_points > 0 ) {
00231 j = 1;
00232 for (i = 1; i < Points->n_points; i++) {
00233 if ( Points->x[i] != Points->x[j-1] || Points->y[i] != Points->y[j-1]
00234 || Points->z[i] != Points->z[j-1])
00235 {
00236 Points->x[j] = Points->x[i];
00237 Points->y[j] = Points->y[i];
00238 Points->z[j] = Points->z[i];
00239 j++;
00240 }
00241 }
00242 Points->n_points = j;
00243 }
00244
00245 return (Points->n_points);
00246 }
00247
00254 int
00255 Vect_line_prune_thresh (struct line_pnts *Points, double threshold)
00256 {
00257 int ret;
00258
00259 ret = dig_prune(Points,threshold);
00260
00261 if ( ret < Points->n_points );
00262 Points->n_points = ret;
00263
00264 return ( Points->n_points );
00265 }
00266
00277 int
00278 Vect_append_points (struct line_pnts *Points, struct line_pnts *APoints,
00279 int direction)
00280 {
00281 int i, n, on, an;
00282
00283 on = Points->n_points;
00284 an = APoints->n_points;
00285 n = on + an;
00286
00287
00288 if (0 > dig_alloc_points (Points, n))
00289 return (-1);
00290
00291 if ( direction == GV_FORWARD ) {
00292 for (i = 0; i < an; i++) {
00293 Points->x[on + i] = APoints->x[i];
00294 Points->y[on + i] = APoints->y[i];
00295 Points->z[on + i] = APoints->z[i];
00296 }
00297 } else {
00298 for (i = 0; i < an; i++) {
00299 Points->x[on + i] = APoints->x[an - i - 1];
00300 Points->y[on + i] = APoints->y[an - i - 1];
00301 Points->z[on + i] = APoints->z[an - i - 1];
00302 }
00303 }
00304
00305 Points->n_points = n;
00306 return n;
00307 }
00308
00309
00317 int
00318 Vect_copy_pnts_to_xyz (
00319 struct line_pnts *Points, double *x, double *y, double *z, int *n)
00320 {
00321 register int i;
00322
00323 for (i = 0; i < *n; i++)
00324 {
00325 x[i] = Points->x[i];
00326 y[i] = Points->y[i];
00327 if ( z != NULL )
00328 z[i] = Points->z[i];
00329 *n = Points->n_points;
00330 }
00331
00332 return (Points->n_points);
00333 }
00334
00354 int
00355 Vect_point_on_line ( struct line_pnts *Points, double distance,
00356 double *x, double *y, double *z, double *angle, double *slope )
00357 {
00358 int j, np, seg=0;
00359 double dist = 0, length;
00360 double xp=0, yp=0, zp=0, dx=0, dy=0, dz=0, dxy=0, dxyz, k, rest;
00361
00362 G_debug ( 3, "Vect_point_on_line(): distance = %f", distance);
00363 if ( (distance < 0) || (Points->n_points < 2) ) return 0;
00364
00365
00366 length = Vect_line_length ( Points );
00367 G_debug ( 3, " length = %f", length);
00368 if ( distance < 0 || distance > length ) {
00369 G_debug ( 3, " -> outside line");
00370 return 0;
00371 }
00372
00373 np = Points->n_points;
00374 if ( distance == 0 ) {
00375 G_debug ( 3, " -> first point");
00376 xp = Points->x[0];
00377 yp = Points->y[0];
00378 zp = Points->z[0];
00379 dx = Points->x[1] - Points->x[0];
00380 dy = Points->y[1] - Points->y[0];
00381 dz = Points->z[1] - Points->z[0];
00382 dxy = hypot (dx, dy);
00383 seg = 1;
00384 } else if ( distance == length ) {
00385 G_debug ( 3, " -> last point");
00386 xp = Points->x[np-1];
00387 yp = Points->y[np-1];
00388 zp = Points->z[np-1];
00389 dx = Points->x[np-1] - Points->x[np-2];
00390 dy = Points->y[np-1] - Points->y[np-2];
00391 dz = Points->z[np-1] - Points->z[np-2];
00392 dxy = hypot (dx, dy);
00393 seg = np - 1;
00394 } else {
00395 for ( j = 0; j < Points->n_points - 1; j++) {
00396
00397
00398 dx = Points->x[j+1] - Points->x[j];
00399 dy = Points->y[j+1] - Points->y[j];
00400 dz = Points->z[j+1] - Points->z[j];
00401 dxy = hypot (dx, dy);
00402 dxyz = hypot (dxy, dz);
00403
00404 dist += dxyz;
00405 if ( dist >= distance ){
00406 rest = distance - dist + dxyz;
00407 k = rest / dxyz;
00408
00409 xp = Points->x[j] + k * dx;
00410 yp = Points->y[j] + k * dy;
00411 zp = Points->z[j] + k * dz;
00412 seg = j + 1;
00413 break;
00414 }
00415 }
00416 }
00417
00418 if ( x != NULL ) *x = xp;
00419 if ( y != NULL ) *y = yp;
00420 if ( z != NULL ) *z = zp;
00421
00422
00423 if ( angle != NULL ) *angle = atan2(dy, dx);
00424
00425
00426 if ( slope != NULL ) *slope = atan2(dz, dxy);
00427
00428 return seg;
00429 }
00430
00444 int
00445 Vect_line_segment ( struct line_pnts *InPoints, double start, double end,
00446 struct line_pnts *OutPoints )
00447 {
00448 int i, seg1, seg2;
00449 double length, tmp;
00450 double x1, y1, z1, x2, y2, z2;
00451
00452 G_debug ( 3, "Vect_line_segment(): start = %f, end = %f, n_points = %d", start, end, InPoints->n_points);
00453
00454 Vect_reset_line (OutPoints);
00455
00456 if ( start > end ) {
00457 tmp = start;
00458 start = end;
00459 end = tmp;
00460 }
00461
00462
00463 if ( end < 0 ) return 0;
00464 length = Vect_line_length ( InPoints );
00465 if ( start > length ) return 0;
00466
00467
00468 seg1 = Vect_point_on_line ( InPoints, start, &x1, &y1, &z1, NULL, NULL);
00469 seg2 = Vect_point_on_line ( InPoints, end, &x2, &y2, &z2, NULL, NULL);
00470
00471 G_debug ( 3, " -> seg1 = %d seg2 = %d", seg1, seg2);
00472
00473 if ( seg1 == 0 || seg2 == 0 ) {
00474 G_warning ("Segment outside line, no segment created");
00475 return 0;
00476 }
00477
00478 Vect_append_point ( OutPoints, x1, y1, z1 );
00479
00480 for ( i = seg1; i < seg2; i++ ) {
00481 Vect_append_point ( OutPoints, InPoints->x[i], InPoints->y[i], InPoints->z[i] );
00482 };
00483
00484 Vect_append_point ( OutPoints, x2, y2, z2 );
00485
00486 return 1;
00487 }
00488
00495 double
00496 Vect_line_length ( struct line_pnts *Points )
00497 {
00498 int j;
00499 double dx, dy, dz, len = 0;
00500
00501 if ( Points->n_points < 2 ) return 0;
00502
00503 for ( j = 0; j < Points->n_points - 1; j++) {
00504 dx = Points->x[j+1] - Points->x[j];
00505 dy = Points->y[j+1] - Points->y[j];
00506 dz = Points->z[j+1] - Points->z[j];
00507 len += hypot ( hypot (dx, dy), dz );
00508 }
00509
00510 return len;
00511 }
00512
00513
00520 double
00521 Vect_line_geodesic_length ( struct line_pnts *Points )
00522 {
00523 int j, dc;
00524 double dx, dy, dz, dxy, len = 0;
00525
00526 dc = G_begin_distance_calculations();
00527
00528 if ( Points->n_points < 2 ) return 0;
00529
00530 for ( j = 0; j < Points->n_points - 1; j++) {
00531 if ( dc == 2 )
00532 dxy = G_geodesic_distance ( Points->x[j], Points->y[j], Points->x[j+1], Points->y[j+1] );
00533 else {
00534 dx = Points->x[j+1] - Points->x[j];
00535 dy = Points->y[j+1] - Points->y[j];
00536 dxy = hypot (dx, dy);
00537 }
00538
00539 dz = Points->z[j+1] - Points->z[j];
00540 len += hypot ( dxy, dz );
00541 }
00542
00543 return len;
00544 }
00545
00565
00566 int
00567 Vect_line_distance (
00568 struct line_pnts *points,
00569 double ux, double uy, double uz,
00570 int with_z,
00571 double *px, double *py, double *pz,
00572 double *dist,
00573 double *spdist,
00574 double *lpdist)
00575 {
00576 register int i;
00577 register double distance;
00578 register double new_dist;
00579 double tpx, tpy, tpz, tdist, tspdist, tlpdist=0;
00580 double dx, dy, dz;
00581 register int n_points;
00582 int segment;
00583
00584 n_points = points->n_points;
00585
00586 if ( n_points == 1 ) {
00587 distance = dig_distance2_point_to_line (ux, uy, uz, points->x[0], points->y[0], points->z[0],
00588 points->x[0], points->y[0], points->z[0],
00589 with_z,
00590 NULL, NULL, NULL, NULL, NULL);
00591 tpx = points->x[0];
00592 tpy = points->y[0];
00593 tpz = points->z[0];
00594 tdist = sqrt(distance);
00595 tspdist = 0;
00596 tlpdist = 0;
00597 segment = 0;
00598
00599 } else {
00600
00601 distance = dig_distance2_point_to_line (ux, uy, uz, points->x[0], points->y[0], points->z[0],
00602 points->x[1], points->y[1], points->z[1],
00603 with_z,
00604 NULL, NULL, NULL, NULL, NULL);
00605 segment = 1;
00606
00607 for ( i = 1 ; i < n_points - 1; i++)
00608 {
00609 new_dist = dig_distance2_point_to_line (ux, uy, uz,
00610 points->x[i], points->y[i], points->z[i],
00611 points->x[i + 1], points->y[i + 1], points->z[i + 1],
00612 with_z,
00613 NULL, NULL, NULL, NULL, NULL);
00614 if (new_dist < distance)
00615 {
00616 distance = new_dist;
00617 segment = i + 1;
00618 }
00619 }
00620
00621
00622 new_dist = dig_distance2_point_to_line (ux, uy, uz,
00623 points->x[segment - 1], points->y[segment - 1], points->z[segment - 1],
00624 points->x[segment], points->y[segment], points->z[segment],
00625 with_z,
00626 &tpx, &tpy, &tpz, &tspdist, NULL);
00627
00628
00629 if ( lpdist ) {
00630 tlpdist = 0;
00631 for ( i = 0; i < segment - 1; i++) {
00632 dx = points->x[i+1] - points->x[i];
00633 dy = points->y[i+1] - points->y[i];
00634 if ( with_z )
00635 dz = points->z[i+1] - points->z[i];
00636 else
00637 dz = 0;
00638
00639 tlpdist += hypot ( hypot (dx, dy), dz );
00640 }
00641 tlpdist += tspdist;
00642 }
00643 tdist = sqrt(distance);
00644 }
00645
00646 if (px) *px = tpx;
00647 if (py) *py = tpy;
00648 if (pz && with_z) *pz = tpz;
00649 if (dist) *dist = tdist;
00650 if (spdist) *spdist = tspdist;
00651 if (lpdist) *lpdist = tlpdist;
00652
00653 return (segment);
00654 }
00655
00656
00665 double
00666 Vect_points_distance (
00667 double x1, double y1, double z1,
00668 double x2, double y2, double z2,
00669 int with_z)
00670 {
00671 double dx, dy, dz;
00672
00673
00674 dx = x2 - x1;
00675 dy = y2 - y1;
00676 dz = z2 - z1;
00677
00678 if (with_z)
00679 return hypot ( hypot (dx, dy), dz );
00680 else
00681 return hypot (dx, dy);
00682
00683 }
00684
00691 int
00692 Vect_line_box ( struct line_pnts *Points, BOUND_BOX *Box )
00693 {
00694 dig_line_box ( Points, Box );
00695 return 0;
00696 }
00697
00703
00704 void Vect_line_reverse ( struct line_pnts *Points )
00705 {
00706 int i, j, np;
00707 double x, y, z;
00708
00709 np = (int) Points->n_points/2 ;
00710
00711 for ( i=0; i<np; i++) {
00712 j = Points->n_points - i - 1;
00713 x = Points->x[i];
00714 y = Points->y[i];
00715 z = Points->z[i];
00716 Points->x[i] = Points->x[j];
00717 Points->y[i] = Points->y[j];
00718 Points->z[i] = Points->z[j];
00719 Points->x[j] = x;
00720 Points->y[j] = y;
00721 Points->z[j] = z;
00722 }
00723 }
00724
00725
00734 int
00735 Vect_get_line_cat ( struct Map_info *Map, int line, int field ) {
00736
00737 static struct line_cats *cats = NULL;
00738 int cat, ltype;
00739
00740 if ( cats == NULL )
00741 cats = Vect_new_cats_struct ();
00742
00743 ltype=Vect_read_line (Map, NULL, cats, line );
00744 Vect_cat_get(cats, field, &cat);
00745 G_debug (3, "Vect_get_line_cat: display line %d, ltype %d, cat %d", line, ltype, cat);
00746
00747 return cat;
00748 }
00749