00001
00002
00003
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
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
00024
00025 dx = dy = 0;
00026 }
00027 *x = dx/l;
00028 *y = dy/l;
00029 }
00030
00031
00032
00033
00034
00035
00036
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
00079
00080
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
00103
00104
00105
00106
00107
00108
00109
00110
00111
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
00135 while( first < np-2 ){
00136
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; }
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; }
00148
00149
00150 if (sa > sa_max ) sa_max = sa;
00151 if (sa < sa_max ) break;
00152
00153
00154 if ( (sb-sa) == 1 ){
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++ ) {
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 ) ){
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 {
00173 first=sb;
00174 continue;
00175 }
00176 }
00177
00178 for (i=j;i<Points->n_points;i++)
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
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
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
00231
00232
00233
00234
00235
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 );
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 ) {
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
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) );
00354
00355 if ( Points == NULL )
00356 Points = Vect_new_line_struct();
00357
00358 if ( PPoints == NULL )
00359 PPoints = Vect_new_line_struct();
00360
00361
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 ) {
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
00380 Vect_append_point ( OutPoints, OutPoints->x[0], OutPoints->y[0], 0 );
00381 } else {
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
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
00397
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
00411 vect ( lx1, ly1, lx2, ly2, &nx, &ny);
00412
00413
00414 sangle = atan2 ( -nx, ny );
00415 sx = lx2 + ny * distance;
00416 sy = ly2 - nx * distance;
00417
00418
00419 ex = lx2 - ny * distance;
00420 ey = ly2 + nx * distance;
00421
00422 Vect_append_point ( OutPoints, sx, sy, 0 );
00423
00424
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
00435 Vect_append_point ( OutPoints, OutPoints->x[0], OutPoints->y[0], 0 );
00436 }
00437 Vect_line_prune ( OutPoints );
00438 }
00439
00440
00441