00001
00002
00003
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
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 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
00134 while( first < np-2 ){
00135
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; }
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; }
00147
00148
00149 if ( (sb-sa) == 1 ){
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++ ) {
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 ) ){
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 {
00168 first=sb;
00169 continue;
00170 }
00171 }
00172
00173 for (i=j;i<Points->n_points;i++)
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
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
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
00226
00227
00228
00229
00230
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 );
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 ) {
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) );
00347
00348 if ( Points == NULL )
00349 Points = Vect_new_line_struct();
00350
00351 if ( PPoints == NULL )
00352 PPoints = Vect_new_line_struct();
00353
00354
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 ) {
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
00373 Vect_append_point ( OutPoints, OutPoints->x[0], OutPoints->y[0], 0 );
00374 } else {
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
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
00390
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
00404 vect ( lx1, ly1, lx2, ly2, &nx, &ny);
00405
00406
00407 sangle = atan2 ( -nx, ny );
00408 sx = lx2 + ny * distance;
00409 sy = ly2 - nx * distance;
00410
00411
00412 ex = lx2 - ny * distance;
00413 ey = ly2 + nx * distance;
00414
00415 Vect_append_point ( OutPoints, sx, sy, 0 );
00416
00417
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
00428 Vect_append_point ( OutPoints, OutPoints->x[0], OutPoints->y[0], 0 );
00429 }
00430 Vect_line_prune ( OutPoints );
00431 }
00432
00433
00434