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 "Vect.h"
00010 #include "gis.h"
00011 
00012 #define LENGTH(DX, DY)  (  sqrt( (DX*DX)+(DY*DY) )  )
00013 #define PI 3.141592653589793116
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 first=0, current, last, lcount;
00117     double *x, *y, px, py, ix, iy;
00118     static struct line_pnts *sPoints = NULL;    
00119 
00120     G_debug (4, "clean_parallel(): npoints = %d, d = %f, rm_end = %d", Points->n_points, d, rm_end ); 
00121 
00122     x = Points->x;
00123     y = Points->y;
00124     np = Points->n_points;    
00125 
00126     if ( sPoints == NULL )
00127         sPoints = Vect_new_line_struct(); 
00128 
00129     Vect_reset_line ( sPoints );
00130 
00131     npn=1;
00132 
00133     /* remove loopes */
00134     while( first < np-2 ){
00135         /* find first loop which doesn't contain any other loop */
00136         current=first;  last=Points->n_points-2;  lcount=0;
00137         while( find_cross ( Points, current, last-1, current+1, last,  &sa, &sb ) != 0 ) 
00138         {  
00139             if ( lcount == 0 ){ first=sa; } /* move first forward */    
00140 
00141             current=sa+1;
00142             last=sb;    
00143             lcount++;
00144             G_debug (5, "  current = %d, last = %d, lcount = %d", current, last, lcount); 
00145         }
00146         if ( lcount == 0 ) { break; }   /* loop not found */
00147 
00148         /* remove loop if in buffer */          
00149         if ( (sb-sa) == 1 ){ /* neighbouring lines overlap */  
00150             j=sb+1;
00151             npn=sa+1;
00152         } else {
00153             Vect_reset_line ( sPoints );
00154             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);
00155             Vect_append_point ( sPoints, ix, iy, 0 );
00156             for ( i=sa+1 ; i < sb+1; i++ ) { /* create loop polygon */
00157                 Vect_append_point ( sPoints, x[i], y[i], 0 );
00158             }
00159             Vect_find_poly_centroid  ( sPoints, &px, &py);   
00160             if ( point_in_buf( oPoints, px, py, d )  ){ /* is loop in buffer ? */
00161                 npn=sa+1;
00162                 x[npn] = ix;
00163                 y[npn] = iy;
00164                 j=sb+1;
00165                 npn++;
00166                 if ( lcount == 0 ){ first=sb; } 
00167             } else {  /* loop is not in buffer */
00168                 first=sb;
00169                 continue;
00170             }
00171         }           
00172 
00173         for (i=j;i<Points->n_points;i++) /* move points down */ 
00174         {
00175             x[npn] = x[i];
00176             y[npn] = y[i];
00177             npn++;
00178         }    
00179         Points->n_points=npn;
00180     }
00181     
00182     if ( rm_end ) {
00183         /* remove points from start in buffer */
00184         j=0;
00185         for (i=0;i<Points->n_points-1;i++) 
00186         {
00187             px=(x[i]+x[i+1])/2;
00188             py=(y[i]+y[i+1])/2;
00189             if ( point_in_buf ( oPoints, x[i], y[i], d*0.9999) 
00190                  && point_in_buf ( oPoints, px, py, d*0.9999) ){
00191                 j++;
00192             } else {
00193                 break;
00194             }
00195         }
00196         if (j>0){
00197             npn=0;    
00198             for (i=j;i<Points->n_points;i++) 
00199             {
00200                 x[npn] = x[i];
00201                 y[npn] = y[i];
00202                 npn++;
00203             }
00204             Points->n_points = npn;
00205         }    
00206         /* remove points from end in buffer */
00207         j=0;
00208         for (i=Points->n_points-1 ;i>=1; i--) 
00209         {
00210             px=(x[i]+x[i-1])/2;
00211             py=(y[i]+y[i-1])/2;
00212             if ( point_in_buf ( oPoints, x[i], y[i], d*0.9999) 
00213                  && point_in_buf ( oPoints, px, py, d*0.9999) ){
00214                 j++;
00215             } else {
00216                 break;
00217             }
00218         }
00219         if (j>0){
00220             Points->n_points -= j;
00221         }    
00222     }
00223 }
00224 
00225 /* parallel_line - remove duplicate points from input line and 
00226 *  creates new parallel line in 'd' offset distance; 
00227 *  'tol' is tolerance between arc and polyline;
00228 *  this function doesn't care about created loopes;
00229 *  
00230 *  New line is written to existing nPoints structure.
00231 */
00232 void parallel_line (struct line_pnts *Points, double d, double tol, struct line_pnts *nPoints)  
00233 {
00234     int i, j, np, na, side;
00235     double *x, *y, nx, ny, tx, ty, vx, vy, ux, uy, wx, wy;
00236     double atol, atol2, a, av, aw;
00237 
00238     G_debug (4, "parallel_line()" ); 
00239 
00240     Vect_reset_line ( nPoints );
00241     
00242     Vect_line_prune ( Points );  
00243     np = Points->n_points;
00244     x = Points->x;
00245     y = Points->y;
00246     
00247     if ( np == 0 )
00248         return;
00249 
00250     if ( np == 1 ) {
00251         Vect_append_point ( nPoints, x[0], y[0], 0 ); /* ? OK, should make circle for points ? */
00252         return;
00253     }
00254 
00255     if ( d == 0 ) {
00256         Vect_copy_xyz_to_pnts ( nPoints, x, y, NULL, np );
00257         return;
00258     }
00259 
00260     side = d/abs(d);
00261     atol = 2 * acos( 1-tol/fabs(d) );
00262 
00263     for (i = 0; i < np-1; i++)
00264     {
00265         vect ( x[i], y[i], x[i+1], y[i+1], &tx, &ty);
00266         vx  = ty * d;
00267         vy  = -tx * d;
00268 
00269         nx = x[i] + vx; 
00270         ny = y[i] + vy;
00271         Vect_append_point ( nPoints, nx, ny, 0 );       
00272 
00273         nx = x[i+1] + vx; 
00274         ny = y[i+1] + vy;
00275         Vect_append_point ( nPoints, nx, ny, 0 );
00276         
00277         if ( i < np-2 ) {  /* use polyline instead of arc between line segments */
00278             vect ( x[i+1], y[i+1], x[i+2], y[i+2], &ux, &uy);
00279             wx  =  uy * d;
00280             wy  = -ux * d;                  
00281             av = atan2 ( vy, vx );  
00282             aw = atan2 ( wy, wx );
00283             a = (aw - av) * side;
00284             if ( a < 0 )  a+=2*PI;
00285             if ( a < PI && a > atol)
00286             {
00287                 na = (int) (a/atol);
00288                 atol2 = a/(na+1) * side;
00289                 for (j = 0; j < na; j++)
00290                 {
00291                     av+=atol2;
00292                     nx = x[i+1] + fabs(d) * cos(av);
00293                     ny = y[i+1] + fabs(d) * sin(av); 
00294                     Vect_append_point ( nPoints, nx, ny, 0 );
00295                 }
00296             }
00297         }   
00298     }
00299     Vect_line_prune ( nPoints );  
00300 }
00301 
00312 void
00313 Vect_line_parallel ( struct line_pnts *InPoints, double distance, double tolerance, int rm_end,
00314                      struct line_pnts *OutPoints )
00315 {
00316     G_debug (4, "Vect_line_parallel(): npoints = %d, distance = %f, tolerance = %f", 
00317                   InPoints->n_points, distance, tolerance);
00318 
00319     parallel_line ( InPoints, distance, tolerance, OutPoints);
00320             
00321     clean_parallel ( OutPoints, InPoints, distance, rm_end );
00322 }
00323 
00335 void
00336 Vect_line_buffer ( struct line_pnts *InPoints, double distance, double tolerance,
00337                      struct line_pnts *OutPoints )
00338 {
00339     double dangle;
00340     int    side, npoints;
00341     static struct line_pnts *Points = NULL;    
00342     static struct line_pnts *PPoints = NULL;    
00343 
00344     distance = fabs (distance );
00345 
00346     dangle = 2 * acos( 1-tolerance/fabs(distance) ); /* angle step */
00347     
00348     if ( Points == NULL )
00349         Points = Vect_new_line_struct(); 
00350     
00351     if ( PPoints == NULL )
00352         PPoints = Vect_new_line_struct(); 
00353 
00354     /* Copy and prune input */
00355     Vect_reset_line ( Points );
00356     Vect_append_points ( Points, InPoints, GV_FORWARD );
00357     Vect_line_prune ( Points );
00358 
00359     Vect_reset_line ( OutPoints );
00360 
00361     npoints = Points->n_points;
00362     if ( npoints <= 0 ) {
00363         return;
00364     } else if ( npoints == 1 ) { /* make a circle */
00365         double angle, x, y;
00366 
00367         for ( angle = 0; angle < 2*PI; angle += dangle ) {
00368             x = Points->x[0] + distance * cos( angle );
00369             y = Points->y[0] + distance * sin( angle ); 
00370             Vect_append_point ( OutPoints, x, y, 0 );
00371         }
00372         /* Close polygon */    
00373         Vect_append_point ( OutPoints, OutPoints->x[0], OutPoints->y[0], 0 );
00374     } else { /* 2 and more points */
00375         for ( side = 0; side < 2; side++ ) {
00376             double angle, sangle;
00377             double lx1, ly1, lx2, ly2;
00378             double x, y, nx, ny, sx, sy, ex, ey;
00379             
00380             /* Parallel on one side */
00381             if ( side == 0 ) {
00382                 Vect_line_parallel ( Points, distance, tolerance, 0, PPoints );
00383                 Vect_append_points ( OutPoints, PPoints, GV_FORWARD );
00384             } else {
00385                 Vect_line_parallel ( Points, -distance, tolerance, 0, PPoints );
00386                 Vect_append_points ( OutPoints, PPoints, GV_BACKWARD );
00387             }
00388 
00389             /* Arc at the end */
00390             /* 2 points at theend of original line */
00391             if ( side == 0 ) {
00392                 lx1 = Points->x[npoints-2];
00393                 ly1 = Points->y[npoints-2];
00394                 lx2 = Points->x[npoints-1];
00395                 ly2 = Points->y[npoints-1];
00396             } else {
00397                 lx1 = Points->x[1];
00398                 ly1 = Points->y[1];
00399                 lx2 = Points->x[0];
00400                 ly2 = Points->y[0];
00401             }
00402 
00403             /* normalized vector */
00404             vect ( lx1, ly1, lx2, ly2, &nx, &ny);
00405 
00406             /* starting point */
00407             sangle = atan2 ( -nx, ny ); /* starting angle */
00408             sx  = lx2 + ny * distance;
00409             sy  = ly2 - nx * distance;              
00410 
00411             /* end point */
00412             ex  = lx2 - ny * distance;
00413             ey  = ly2 + nx * distance;              
00414 
00415             Vect_append_point ( OutPoints, sx, sy, 0 );
00416          
00417             /* arc */
00418             for ( angle = dangle; angle < PI; angle += dangle ) {
00419                 x = lx2 + distance * cos( sangle + angle );
00420                 y = ly2 + distance * sin( sangle + angle ); 
00421                 Vect_append_point ( OutPoints, x, y, 0 );
00422             }
00423             
00424             Vect_append_point ( OutPoints, ex, ey, 0 );
00425         }
00426         
00427         /* Close polygon */    
00428         Vect_append_point ( OutPoints, OutPoints->x[0], OutPoints->y[0], 0 );
00429     }
00430     Vect_line_prune ( OutPoints );
00431 }
00432 
00433 
00434 

Generated on Sat Jul 22 22:05:56 2006 for GRASS by  doxygen 1.4.7