line.c

Go to the documentation of this file.
00001 /*
00002 ****************************************************************************
00003 *
00004 * MODULE:       Vector library 
00005 *               
00006 * AUTHOR(S):    Original author CERL, probably Dave Gerdes or Mike Higgins.
00007 *               Update to GRASS 5.7 Radim Blazek and David D. Gray.
00008 *
00009 * PURPOSE:      Higher level functions for reading/writing/manipulating vectors.
00010 *
00011 * COPYRIGHT:    (C) 2001 by the GRASS Development Team
00012 *
00013 *               This program is free software under the GNU General Public
00014 *               License (>=v2). Read the file COPYING that comes with GRASS
00015 *               for details.
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   /* alloc_points MUST be initialized to zero */
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)                        /* probably a moot test */
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 /* NEW after 4.0 alpha */
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   /* move up */
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   /* move down */
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) /* GV_FORWARD, GV_BACKWARD */
00280 {
00281   int i, n, on, an;
00282 
00283   on = Points->n_points;
00284   an = APoints->n_points;
00285   n = on + an;
00286   
00287   /* Should be OK, dig_alloc_points calls realloc */
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     /* Check if first or last */
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             /* dxyz = G_distance(Points->x[j], Points->y[j], 
00397                                  Points->x[j+1], Points->y[j+1]); */
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 ){ /* point is on the current line part */
00406                 rest = distance - dist + dxyz; /* from first point of segment to point */
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     /* calculate angle */
00423     if ( angle != NULL ) *angle = atan2(dy, dx);
00424 
00425     /* calculate slope */
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     /* Check start/end */
00463     if ( end < 0 ) return 0; 
00464     length = Vect_line_length ( InPoints );
00465     if ( start > length ) return 0;
00466     
00467     /* Find coordinates and segments of start/end */
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 /* original dig__check_dist () in grass50 */
00566 int 
00567 Vect_line_distance (
00568                   struct line_pnts *points, /* line */
00569                   double ux, double uy, double uz,    /* point */
00570                   int    with_z,   /* use z coordinate, (3D calculation) */
00571                   double *px, double *py, double *pz, /* point on line */
00572                   double *dist,             /* distance to line */
00573                   double *spdist,           /* distance of point from segment beginning */
00574                   double *lpdist)           /* distance of point from line beginning */
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       /* we have segment and now we can recalculate other values (speed) */
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       /* calculate distance from beginning of line */
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,    /* point 1 */
00668                   double x2, double y2, double z2,    /* point 2 */
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 /* reverse_line - reverse order of Points ( direction of line ) */
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 

Generated on Sun Apr 6 17:32:44 2008 for GRASS by  doxygen 1.5.5