buffer.c

Go to the documentation of this file.
00001 /* Functions: nearest, adjust_line, parallel_line 
00002 **
00003 ** Author: Radim Blazek Feb 2000
00004 ** 
00005 **
00006 */
00007 #include <stdlib.h>
00008 #include <math.h>
00009 #include <grass/Vect.h>
00010 #include <grass/gis.h>
00011 
00012 #define LENGTH(DX, DY)  (  sqrt( (DX*DX)+(DY*DY) )  )
00013 #define PI M_PI
00014 
00015 /* vector() calculates normalized vector form two points */
00016 static void vect(double x1, double y1, double x2, double y2, double *x, double *y )
00017 {
00018     double dx, dy, l;
00019     dx  = x2 - x1;
00020     dy  = y2 - y1;      
00021     l   = LENGTH ( dx, dy );
00022     if (l == 0) {
00023         /* assume that dx == dy == 0, which should give (NaN,NaN) */
00024         /* without this, very small dx or dy could result in Infinity */
00025         dx = dy = 0;
00026     }
00027     *x  = dx/l;
00028     *y  = dy/l;
00029 }
00030 
00031 /* find_cross find first crossing between segments from s1 to s2 and from s3 to s4 
00032 ** s5 is set to first segment and s6 to second
00033 ** neighbours are taken as crossing each other only if overlap
00034 ** returns: 1 found
00035 **         -1 found overlap 
00036 **          0 not found    
00037 */
00038 int find_cross ( struct line_pnts *Points, int s1, int s2, int s3, int s4,  int *s5, int *s6 )  
00039 {
00040     int i, j, np, ret;
00041     double *x, *y;
00042 
00043     G_debug (5, "find_cross(): npoints = %d, s1 = %d, s2 = %d, s3 = %d, s4 = %d", 
00044                          Points->n_points, s1, s2, s3, s4 ); 
00045     
00046     x = Points->x;
00047     y = Points->y;
00048     np = Points->n_points;
00049 
00050     for ( i=s1; i<=s2; i++) 
00051     {
00052         for ( j=s3; j<=s4; j++) 
00053         {
00054             if ( j==i ){
00055                 continue;           
00056             }
00057             ret = dig_test_for_intersection ( x[i], y[i], x[i+1], y[i+1], x[j], y[j], x[j+1], y[j+1] );     
00058             if ( ret == 1 &&  ( (i-j) > 1 || (i-j) < -1 )  )
00059             {
00060                 *s5 = i;
00061                 *s6 = j;
00062                 G_debug (5, "  intersection: s5 = %d, s6 = %d", *s5, *s6 ); 
00063                 return 1;
00064             }       
00065             if (  ret == -1  )
00066             {
00067                 *s5 = i;
00068                 *s6 = j;
00069                 G_debug (5, "  overlap: s5 = %d, s6 = %d", *s5, *s6 ); 
00070                 return -1;
00071             }
00072         }    
00073     }
00074     G_debug (5, "  no intersection" ); 
00075     return 0;
00076 }
00077 
00078 /* point_in_buf - test if point px,py is in d buffer of Points
00079 ** returns:  1 in buffer  
00080 **           0 not  in buffer   
00081 */
00082 int point_in_buf ( struct line_pnts *Points, double px, double py, double d )
00083 {
00084     int i, np;
00085     double sd;
00086 
00087     np = Points->n_points;
00088     d *= d;
00089     for ( i=0; i < np-1; i++) 
00090     {
00091         sd = dig_distance2_point_to_line ( px, py, 0, 
00092                 Points->x[i], Points->y[i], 0, Points->x[i+1], Points->y[i+1], 0,
00093                 0, NULL, NULL, NULL, NULL, NULL );     
00094         if ( sd <= d )
00095         {
00096             return 1;
00097         }           
00098     }
00099     return 0;
00100 }
00101 
00102 /* clean_parallel - clean parallel line created by parallel_line:
00103 ** - looking for loopes and if loop doesn't contain any other loop
00104 **   and centroid of loop is in buffer removes this loop (repeated)
00105 ** - optionaly removes all end points in buffer 
00106 ** Points - parallel line, oPoints - original line, d - offset, rmend - remove end points in buffer 
00107 ** note1: on some lines (multiply selfcrossing; lines with end points 
00108 **        in buffer of line other; some shapes of ends ) may create nosense
00109 ** note2: this function is stupid and slow, somebody more clever
00110 **        than I am should write paralle_line + clean_parallel 
00111 **        better;    RB March 2000
00112 */
00113 void clean_parallel ( struct line_pnts *Points, struct line_pnts *oPoints, double d , int rm_end )  
00114 {
00115     int i, j, np, npn, sa, sb;
00116     int sa_max = 0;
00117     int first=0, current, last, lcount;
00118     double *x, *y, px, py, ix, iy;
00119     static struct line_pnts *sPoints = NULL;    
00120 
00121     G_debug (4, "clean_parallel(): npoints = %d, d = %f, rm_end = %d", Points->n_points, d, rm_end ); 
00122 
00123     x = Points->x;
00124     y = Points->y;
00125     np = Points->n_points;    
00126 
00127     if ( sPoints == NULL )
00128         sPoints = Vect_new_line_struct(); 
00129 
00130     Vect_reset_line ( sPoints );
00131 
00132     npn=1;
00133 
00134     /* remove loopes */
00135     while( first < np-2 ){
00136         /* find first loop which doesn't contain any other loop */
00137         current=first;  last=Points->n_points-2;  lcount=0;
00138         while( find_cross ( Points, current, last-1, current+1, last,  &sa, &sb ) != 0 ) 
00139         {  
00140             if ( lcount == 0 ){ first=sa; } /* move first forward */    
00141 
00142             current=sa+1;
00143             last=sb;    
00144             lcount++;
00145             G_debug (5, "  current = %d, last = %d, lcount = %d", current, last, lcount); 
00146         }
00147         if ( lcount == 0 ) { break; }   /* loop not found */
00148 
00149         /* ensure sa is monotonically increasing, so npn doesn't reset low */
00150         if (sa > sa_max ) sa_max = sa;
00151         if (sa < sa_max ) break;
00152 
00153         /* remove loop if in buffer */          
00154         if ( (sb-sa) == 1 ){ /* neighbouring lines overlap */  
00155             j=sb+1;
00156             npn=sa+1;
00157         } else {
00158             Vect_reset_line ( sPoints );
00159             dig_find_intersection ( x[sa],y[sa],x[sa+1],y[sa+1],x[sb],y[sb],x[sb+1],y[sb+1], &ix,&iy);
00160             Vect_append_point ( sPoints, ix, iy, 0 );
00161             for ( i=sa+1 ; i < sb+1; i++ ) { /* create loop polygon */
00162                 Vect_append_point ( sPoints, x[i], y[i], 0 );
00163             }
00164             Vect_find_poly_centroid  ( sPoints, &px, &py);   
00165             if ( point_in_buf( oPoints, px, py, d )  ){ /* is loop in buffer ? */
00166                 npn=sa+1;
00167                 x[npn] = ix;
00168                 y[npn] = iy;
00169                 j=sb+1;
00170                 npn++;
00171                 if ( lcount == 0 ){ first=sb; } 
00172             } else {  /* loop is not in buffer */
00173                 first=sb;
00174                 continue;
00175             }
00176         }           
00177 
00178         for (i=j;i<Points->n_points;i++) /* move points down */ 
00179         {
00180             x[npn] = x[i];
00181             y[npn] = y[i];
00182             npn++;
00183         }    
00184         Points->n_points=npn;
00185     }
00186     
00187     if ( rm_end ) {
00188         /* remove points from start in buffer */
00189         j=0;
00190         for (i=0;i<Points->n_points-1;i++) 
00191         {
00192             px=(x[i]+x[i+1])/2;
00193             py=(y[i]+y[i+1])/2;
00194             if ( point_in_buf ( oPoints, x[i], y[i], d*0.9999) 
00195                  && point_in_buf ( oPoints, px, py, d*0.9999) ){
00196                 j++;
00197             } else {
00198                 break;
00199             }
00200         }
00201         if (j>0){
00202             npn=0;    
00203             for (i=j;i<Points->n_points;i++) 
00204             {
00205                 x[npn] = x[i];
00206                 y[npn] = y[i];
00207                 npn++;
00208             }
00209             Points->n_points = npn;
00210         }    
00211         /* remove points from end in buffer */
00212         j=0;
00213         for (i=Points->n_points-1 ;i>=1; i--) 
00214         {
00215             px=(x[i]+x[i-1])/2;
00216             py=(y[i]+y[i-1])/2;
00217             if ( point_in_buf ( oPoints, x[i], y[i], d*0.9999) 
00218                  && point_in_buf ( oPoints, px, py, d*0.9999) ){
00219                 j++;
00220             } else {
00221                 break;
00222             }
00223         }
00224         if (j>0){
00225             Points->n_points -= j;
00226         }    
00227     }
00228 }
00229 
00230 /* parallel_line - remove duplicate points from input line and 
00231 *  creates new parallel line in 'd' offset distance; 
00232 *  'tol' is tolerance between arc and polyline;
00233 *  this function doesn't care about created loopes;
00234 *  
00235 *  New line is written to existing nPoints structure.
00236 */
00237 void parallel_line (struct line_pnts *Points, double d, double tol, struct line_pnts *nPoints)  
00238 {
00239     int i, j, np, na, side;
00240     double *x, *y, nx, ny, tx, ty, vx, vy, ux, uy, wx, wy;
00241     double atol, atol2, a, av, aw;
00242 
00243     G_debug (4, "parallel_line()" ); 
00244 
00245     Vect_reset_line ( nPoints );
00246     
00247     Vect_line_prune ( Points );  
00248     np = Points->n_points;
00249     x = Points->x;
00250     y = Points->y;
00251     
00252     if ( np == 0 )
00253         return;
00254 
00255     if ( np == 1 ) {
00256         Vect_append_point ( nPoints, x[0], y[0], 0 ); /* ? OK, should make circle for points ? */
00257         return;
00258     }
00259 
00260     if ( d == 0 ) {
00261         Vect_copy_xyz_to_pnts ( nPoints, x, y, NULL, np );
00262         return;
00263     }
00264 
00265     side = (int)(d/fabs(d));
00266     atol = 2 * acos( 1-tol/fabs(d) );
00267 
00268     for (i = 0; i < np-1; i++)
00269     {
00270         vect ( x[i], y[i], x[i+1], y[i+1], &tx, &ty);
00271         vx  = ty * d;
00272         vy  = -tx * d;
00273 
00274         nx = x[i] + vx; 
00275         ny = y[i] + vy;
00276         Vect_append_point ( nPoints, nx, ny, 0 );       
00277 
00278         nx = x[i+1] + vx; 
00279         ny = y[i+1] + vy;
00280         Vect_append_point ( nPoints, nx, ny, 0 );
00281         
00282         if ( i < np-2 ) {  /* use polyline instead of arc between line segments */
00283             vect ( x[i+1], y[i+1], x[i+2], y[i+2], &ux, &uy);
00284             wx  =  uy * d;
00285             wy  = -ux * d;                  
00286             av = atan2 ( vy, vx );  
00287             aw = atan2 ( wy, wx );
00288             a = (aw - av) * side;
00289             if ( a < 0 )  a+=2*PI;
00290 
00291             /* TODO: a <= PI can probably fail because of representation error */
00292             if ( a <= PI && a > atol)
00293             {
00294                 na = (int) (a/atol);
00295                 atol2 = a/(na+1) * side;
00296                 for (j = 0; j < na; j++)
00297                 {
00298                     av+=atol2;
00299                     nx = x[i+1] + fabs(d) * cos(av);
00300                     ny = y[i+1] + fabs(d) * sin(av); 
00301                     Vect_append_point ( nPoints, nx, ny, 0 );
00302                 }
00303             }
00304         }   
00305     }
00306     Vect_line_prune ( nPoints );  
00307 }
00308 
00319 void
00320 Vect_line_parallel ( struct line_pnts *InPoints, double distance, double tolerance, int rm_end,
00321                      struct line_pnts *OutPoints )
00322 {
00323     G_debug (4, "Vect_line_parallel(): npoints = %d, distance = %f, tolerance = %f", 
00324                   InPoints->n_points, distance, tolerance);
00325 
00326     parallel_line ( InPoints, distance, tolerance, OutPoints);
00327             
00328     clean_parallel ( OutPoints, InPoints, distance, rm_end );
00329 }
00330 
00342 void
00343 Vect_line_buffer ( struct line_pnts *InPoints, double distance, double tolerance,
00344                      struct line_pnts *OutPoints )
00345 {
00346     double dangle;
00347     int    side, npoints;
00348     static struct line_pnts *Points = NULL;    
00349     static struct line_pnts *PPoints = NULL;    
00350 
00351     distance = fabs (distance );
00352 
00353     dangle = 2 * acos( 1-tolerance/fabs(distance) ); /* angle step */
00354     
00355     if ( Points == NULL )
00356         Points = Vect_new_line_struct(); 
00357     
00358     if ( PPoints == NULL )
00359         PPoints = Vect_new_line_struct(); 
00360 
00361     /* Copy and prune input */
00362     Vect_reset_line ( Points );
00363     Vect_append_points ( Points, InPoints, GV_FORWARD );
00364     Vect_line_prune ( Points );
00365 
00366     Vect_reset_line ( OutPoints );
00367 
00368     npoints = Points->n_points;
00369     if ( npoints <= 0 ) {
00370         return;
00371     } else if ( npoints == 1 ) { /* make a circle */
00372         double angle, x, y;
00373 
00374         for ( angle = 0; angle < 2*PI; angle += dangle ) {
00375             x = Points->x[0] + distance * cos( angle );
00376             y = Points->y[0] + distance * sin( angle ); 
00377             Vect_append_point ( OutPoints, x, y, 0 );
00378         }
00379         /* Close polygon */    
00380         Vect_append_point ( OutPoints, OutPoints->x[0], OutPoints->y[0], 0 );
00381     } else { /* 2 and more points */
00382         for ( side = 0; side < 2; side++ ) {
00383             double angle, sangle;
00384             double lx1, ly1, lx2, ly2;
00385             double x, y, nx, ny, sx, sy, ex, ey;
00386             
00387             /* Parallel on one side */
00388             if ( side == 0 ) {
00389                 Vect_line_parallel ( Points, distance, tolerance, 0, PPoints );
00390                 Vect_append_points ( OutPoints, PPoints, GV_FORWARD );
00391             } else {
00392                 Vect_line_parallel ( Points, -distance, tolerance, 0, PPoints );
00393                 Vect_append_points ( OutPoints, PPoints, GV_BACKWARD );
00394             }
00395 
00396             /* Arc at the end */
00397             /* 2 points at theend of original line */
00398             if ( side == 0 ) {
00399                 lx1 = Points->x[npoints-2];
00400                 ly1 = Points->y[npoints-2];
00401                 lx2 = Points->x[npoints-1];
00402                 ly2 = Points->y[npoints-1];
00403             } else {
00404                 lx1 = Points->x[1];
00405                 ly1 = Points->y[1];
00406                 lx2 = Points->x[0];
00407                 ly2 = Points->y[0];
00408             }
00409 
00410             /* normalized vector */
00411             vect ( lx1, ly1, lx2, ly2, &nx, &ny);
00412 
00413             /* starting point */
00414             sangle = atan2 ( -nx, ny ); /* starting angle */
00415             sx  = lx2 + ny * distance;
00416             sy  = ly2 - nx * distance;              
00417 
00418             /* end point */
00419             ex  = lx2 - ny * distance;
00420             ey  = ly2 + nx * distance;              
00421 
00422             Vect_append_point ( OutPoints, sx, sy, 0 );
00423          
00424             /* arc */
00425             for ( angle = dangle; angle < PI; angle += dangle ) {
00426                 x = lx2 + distance * cos( sangle + angle );
00427                 y = ly2 + distance * sin( sangle + angle ); 
00428                 Vect_append_point ( OutPoints, x, y, 0 );
00429             }
00430             
00431             Vect_append_point ( OutPoints, ex, ey, 0 );
00432         }
00433         
00434         /* Close polygon */    
00435         Vect_append_point ( OutPoints, OutPoints->x[0], OutPoints->y[0], 0 );
00436     }
00437     Vect_line_prune ( OutPoints );
00438 }
00439 
00440 
00441 

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