00001
00020 #include <stdlib.h>
00021 #include <math.h>
00022 #include <grass/Vect.h>
00023 #include <grass/gis.h>
00024
00025 #define LENGTH(DX, DY) ( sqrt( (DX*DX)+(DY*DY) ) )
00026 #define PI M_PI
00027
00028
00029 static void vect(double x1, double y1, double x2, double y2, double *x,
00030 double *y)
00031 {
00032 double dx, dy, l;
00033
00034 dx = x2 - x1;
00035 dy = y2 - y1;
00036 l = LENGTH(dx, dy);
00037 if (l == 0) {
00038
00039
00040 dx = dy = 0;
00041 }
00042 *x = dx / l;
00043 *y = dy / l;
00044 }
00045
00046
00047
00048
00049
00050
00051
00052
00053 static int find_cross(struct line_pnts *Points, int s1, int s2, int s3,
00054 int s4, int *s5, int *s6)
00055 {
00056 int i, j, np, ret;
00057 double *x, *y;
00058
00059 G_debug(5,
00060 "find_cross(): npoints = %d, s1 = %d, s2 = %d, s3 = %d, s4 = %d",
00061 Points->n_points, s1, s2, s3, s4);
00062
00063 x = Points->x;
00064 y = Points->y;
00065 np = Points->n_points;
00066
00067 for (i = s1; i <= s2; i++) {
00068 for (j = s3; j <= s4; j++) {
00069 if (j == i) {
00070 continue;
00071 }
00072 ret =
00073 dig_test_for_intersection(x[i], y[i], x[i + 1], y[i + 1],
00074 x[j], y[j], x[j + 1], y[j + 1]);
00075 if (ret == 1 && ((i - j) > 1 || (i - j) < -1)) {
00076 *s5 = i;
00077 *s6 = j;
00078 G_debug(5, " intersection: s5 = %d, s6 = %d", *s5, *s6);
00079 return 1;
00080 }
00081 if (ret == -1) {
00082 *s5 = i;
00083 *s6 = j;
00084 G_debug(5, " overlap: s5 = %d, s6 = %d", *s5, *s6);
00085 return -1;
00086 }
00087 }
00088 }
00089 G_debug(5, " no intersection");
00090 return 0;
00091 }
00092
00093
00094
00095
00096
00097 static int point_in_buf(struct line_pnts *Points, double px, double py,
00098 double d)
00099 {
00100 int i, np;
00101 double sd;
00102
00103 np = Points->n_points;
00104 d *= d;
00105 for (i = 0; i < np - 1; i++) {
00106 sd = dig_distance2_point_to_line(px, py, 0,
00107 Points->x[i], Points->y[i], 0,
00108 Points->x[i + 1], Points->y[i + 1],
00109 0, 0, NULL, NULL, NULL, NULL, NULL);
00110 if (sd <= d) {
00111 return 1;
00112 }
00113 }
00114 return 0;
00115 }
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132 static void clean_parallel(struct line_pnts *Points,
00133 struct line_pnts *origPoints, double d, int rm_end)
00134 {
00135 int i, j, np, npn, sa, sb;
00136 int sa_max = 0;
00137 int first = 0, current, last, lcount;
00138 double *x, *y, px, py, ix, iy;
00139 static struct line_pnts *sPoints = NULL;
00140
00141 G_debug(4, "clean_parallel(): npoints = %d, d = %f, rm_end = %d",
00142 Points->n_points, d, rm_end);
00143
00144 x = Points->x;
00145 y = Points->y;
00146 np = Points->n_points;
00147
00148 if (sPoints == NULL)
00149 sPoints = Vect_new_line_struct();
00150
00151 Vect_reset_line(sPoints);
00152
00153 npn = 1;
00154
00155
00156 while (first < np - 2) {
00157
00158 current = first;
00159 last = Points->n_points - 2;
00160 lcount = 0;
00161 while (find_cross
00162 (Points, current, last - 1, current + 1, last, &sa,
00163 &sb) != 0) {
00164 if (lcount == 0) {
00165 first = sa;
00166 }
00167
00168 current = sa + 1;
00169 last = sb;
00170 lcount++;
00171 G_debug(5, " current = %d, last = %d, lcount = %d", current,
00172 last, lcount);
00173 }
00174 if (lcount == 0) {
00175 break;
00176 }
00177
00178
00179 if (sa > sa_max)
00180 sa_max = sa;
00181 if (sa < sa_max)
00182 break;
00183
00184
00185 if ((sb - sa) == 1) {
00186 j = sb + 1;
00187 npn = sa + 1;
00188 }
00189 else {
00190 Vect_reset_line(sPoints);
00191 dig_find_intersection(x[sa], y[sa], x[sa + 1], y[sa + 1], x[sb],
00192 y[sb], x[sb + 1], y[sb + 1], &ix, &iy);
00193 Vect_append_point(sPoints, ix, iy, 0);
00194 for (i = sa + 1; i < sb + 1; i++) {
00195 Vect_append_point(sPoints, x[i], y[i], 0);
00196 }
00197 Vect_find_poly_centroid(sPoints, &px, &py);
00198 if (point_in_buf(origPoints, px, py, d)) {
00199 npn = sa + 1;
00200 x[npn] = ix;
00201 y[npn] = iy;
00202 j = sb + 1;
00203 npn++;
00204 if (lcount == 0) {
00205 first = sb;
00206 }
00207 }
00208 else {
00209 first = sb;
00210 continue;
00211 }
00212 }
00213
00214 for (i = j; i < Points->n_points; i++) {
00215 x[npn] = x[i];
00216 y[npn] = y[i];
00217 npn++;
00218 }
00219 Points->n_points = npn;
00220 }
00221
00222 if (rm_end) {
00223
00224 j = 0;
00225 for (i = 0; i < Points->n_points - 1; i++) {
00226 px = (x[i] + x[i + 1]) / 2;
00227 py = (y[i] + y[i + 1]) / 2;
00228 if (point_in_buf(origPoints, x[i], y[i], d * 0.9999)
00229 && point_in_buf(origPoints, px, py, d * 0.9999)) {
00230 j++;
00231 }
00232 else {
00233 break;
00234 }
00235 }
00236 if (j > 0) {
00237 npn = 0;
00238 for (i = j; i < Points->n_points; i++) {
00239 x[npn] = x[i];
00240 y[npn] = y[i];
00241 npn++;
00242 }
00243 Points->n_points = npn;
00244 }
00245
00246 j = 0;
00247 for (i = Points->n_points - 1; i >= 1; i--) {
00248 px = (x[i] + x[i - 1]) / 2;
00249 py = (y[i] + y[i - 1]) / 2;
00250 if (point_in_buf(origPoints, x[i], y[i], d * 0.9999)
00251 && point_in_buf(origPoints, px, py, d * 0.9999)) {
00252 j++;
00253 }
00254 else {
00255 break;
00256 }
00257 }
00258 if (j > 0) {
00259 Points->n_points -= j;
00260 }
00261 }
00262 }
00263
00264
00265
00266
00267
00268
00269
00270
00271 static void parallel_line(struct line_pnts *Points, double d, double tol,
00272 struct line_pnts *nPoints)
00273 {
00274 int i, j, np, na, side;
00275 double *x, *y, nx, ny, tx, ty, vx, vy, ux, uy, wx, wy;
00276 double atol, atol2, a, av, aw;
00277
00278 G_debug(4, "parallel_line()");
00279
00280 Vect_reset_line(nPoints);
00281
00282 Vect_line_prune(Points);
00283 np = Points->n_points;
00284 x = Points->x;
00285 y = Points->y;
00286
00287 if (np == 0)
00288 return;
00289
00290 if (np == 1) {
00291 Vect_append_point(nPoints, x[0], y[0], 0);
00292 return;
00293 }
00294
00295 if (d == 0) {
00296 Vect_copy_xyz_to_pnts(nPoints, x, y, NULL, np);
00297 return;
00298 }
00299
00300 side = (int)(d / fabs(d));
00301 atol = 2 * acos(1 - tol / fabs(d));
00302
00303 for (i = 0; i < np - 1; i++) {
00304 vect(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);
00305 vx = ty * d;
00306 vy = -tx * d;
00307
00308 nx = x[i] + vx;
00309 ny = y[i] + vy;
00310 Vect_append_point(nPoints, nx, ny, 0);
00311
00312 nx = x[i + 1] + vx;
00313 ny = y[i + 1] + vy;
00314 Vect_append_point(nPoints, nx, ny, 0);
00315
00316 if (i < np - 2) {
00317 vect(x[i + 1], y[i + 1], x[i + 2], y[i + 2], &ux, &uy);
00318 wx = uy * d;
00319 wy = -ux * d;
00320 av = atan2(vy, vx);
00321 aw = atan2(wy, wx);
00322 a = (aw - av) * side;
00323 if (a < 0)
00324 a += 2 * PI;
00325
00326
00327 if (a <= PI && a > atol) {
00328 na = (int)(a / atol);
00329 atol2 = a / (na + 1) * side;
00330 for (j = 0; j < na; j++) {
00331 av += atol2;
00332 nx = x[i + 1] + fabs(d) * cos(av);
00333 ny = y[i + 1] + fabs(d) * sin(av);
00334 Vect_append_point(nPoints, nx, ny, 0);
00335 }
00336 }
00337 }
00338 }
00339 Vect_line_prune(nPoints);
00340 }
00341
00353 void
00354 Vect_line_parallel(struct line_pnts *InPoints, double distance,
00355 double tolerance, int rm_end, struct line_pnts *OutPoints)
00356 {
00357 G_debug(4,
00358 "Vect_line_parallel(): npoints = %d, distance = %f, tolerance = %f",
00359 InPoints->n_points, distance, tolerance);
00360
00361 parallel_line(InPoints, distance, tolerance, OutPoints);
00362
00363 clean_parallel(OutPoints, InPoints, distance, rm_end);
00364
00365 return;
00366 }
00367
00379 void
00380 Vect_line_buffer(struct line_pnts *InPoints, double distance,
00381 double tolerance, struct line_pnts *OutPoints)
00382 {
00383 double dangle;
00384 int side, npoints;
00385 static struct line_pnts *Points = NULL;
00386 static struct line_pnts *PPoints = NULL;
00387
00388 distance = fabs(distance);
00389
00390 dangle = 2 * acos(1 - tolerance / fabs(distance));
00391
00392 if (Points == NULL)
00393 Points = Vect_new_line_struct();
00394
00395 if (PPoints == NULL)
00396 PPoints = Vect_new_line_struct();
00397
00398
00399 Vect_reset_line(Points);
00400 Vect_append_points(Points, InPoints, GV_FORWARD);
00401 Vect_line_prune(Points);
00402
00403 Vect_reset_line(OutPoints);
00404
00405 npoints = Points->n_points;
00406 if (npoints <= 0) {
00407 return;
00408 }
00409 else if (npoints == 1) {
00410 double angle, x, y;
00411
00412 for (angle = 0; angle < 2 * PI; angle += dangle) {
00413 x = Points->x[0] + distance * cos(angle);
00414 y = Points->y[0] + distance * sin(angle);
00415 Vect_append_point(OutPoints, x, y, 0);
00416 }
00417
00418 Vect_append_point(OutPoints, OutPoints->x[0], OutPoints->y[0], 0);
00419 }
00420 else {
00421 for (side = 0; side < 2; side++) {
00422 double angle, sangle;
00423 double lx1, ly1, lx2, ly2;
00424 double x, y, nx, ny, sx, sy, ex, ey;
00425
00426
00427 if (side == 0) {
00428 Vect_line_parallel(Points, distance, tolerance, 0, PPoints);
00429 Vect_append_points(OutPoints, PPoints, GV_FORWARD);
00430 }
00431 else {
00432 Vect_line_parallel(Points, -distance, tolerance, 0, PPoints);
00433 Vect_append_points(OutPoints, PPoints, GV_BACKWARD);
00434 }
00435
00436
00437
00438 if (side == 0) {
00439 lx1 = Points->x[npoints - 2];
00440 ly1 = Points->y[npoints - 2];
00441 lx2 = Points->x[npoints - 1];
00442 ly2 = Points->y[npoints - 1];
00443 }
00444 else {
00445 lx1 = Points->x[1];
00446 ly1 = Points->y[1];
00447 lx2 = Points->x[0];
00448 ly2 = Points->y[0];
00449 }
00450
00451
00452 vect(lx1, ly1, lx2, ly2, &nx, &ny);
00453
00454
00455 sangle = atan2(-nx, ny);
00456 sx = lx2 + ny * distance;
00457 sy = ly2 - nx * distance;
00458
00459
00460 ex = lx2 - ny * distance;
00461 ey = ly2 + nx * distance;
00462
00463 Vect_append_point(OutPoints, sx, sy, 0);
00464
00465
00466 for (angle = dangle; angle < PI; angle += dangle) {
00467 x = lx2 + distance * cos(sangle + angle);
00468 y = ly2 + distance * sin(sangle + angle);
00469 Vect_append_point(OutPoints, x, y, 0);
00470 }
00471
00472 Vect_append_point(OutPoints, ex, ey, 0);
00473 }
00474
00475
00476 Vect_append_point(OutPoints, OutPoints->x[0], OutPoints->y[0], 0);
00477 }
00478 Vect_line_prune(OutPoints);
00479
00480 return;
00481 }