GRASS Programmer's Manual
6.4.2(2012)
|
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 /* vector() calculates normalized vector form two points */ 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 /* assume that dx == dy == 0, which should give (NaN,NaN) */ 00039 /* without this, very small dx or dy could result in Infinity */ 00040 dx = dy = 0; 00041 } 00042 *x = dx / l; 00043 *y = dy / l; 00044 } 00045 00046 /* find_cross find first crossing between segments from s1 to s2 and from s3 to s4 00047 ** s5 is set to first segment and s6 to second 00048 ** neighbours are taken as crossing each other only if overlap 00049 ** returns: 1 found 00050 ** -1 found overlap 00051 ** 0 not found 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 /* point_in_buf - test if point px,py is in d buffer of Points 00094 ** returns: 1 in buffer 00095 ** 0 not in buffer 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 /* clean_parallel - clean parallel line created by parallel_line: 00118 ** - looking for loops and if loop doesn't contain any other loop 00119 ** and centroid of loop is in buffer removes this loop (repeated) 00120 ** - optionally removes all end points in buffer 00121 * parameters: 00122 * Points - parallel line 00123 * origPoints - original line 00124 * d - offset 00125 * rm_end - remove end points in buffer 00126 ** note1: on some lines (multiply selfcrossing; lines with end points 00127 ** in buffer of line other; some shapes of ends ) may create nosense 00128 ** note2: this function is stupid and slow, somebody more clever 00129 ** than I am should write paralle_line + clean_parallel 00130 ** better; RB March 2000 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 /* remove loops */ 00156 while (first < np - 2) { 00157 /* find first loop which doesn't contain any other loop */ 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 } /* move first forward */ 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 } /* loop not found */ 00177 00178 /* ensure sa is monotonically increasing, so npn doesn't reset low */ 00179 if (sa > sa_max) 00180 sa_max = sa; 00181 if (sa < sa_max) 00182 break; 00183 00184 /* remove loop if in buffer */ 00185 if ((sb - sa) == 1) { /* neighbouring lines overlap */ 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++) { /* create loop polygon */ 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)) { /* is loop in buffer ? */ 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 { /* loop is not in buffer */ 00209 first = sb; 00210 continue; 00211 } 00212 } 00213 00214 for (i = j; i < Points->n_points; i++) { /* move points down */ 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 /* remove points from start in buffer */ 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 /* remove points from end in buffer */ 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 /* parallel_line - remove duplicate points from input line and 00265 * creates new parallel line in 'd' offset distance; 00266 * 'tol' is tolerance between arc and polyline; 00267 * this function doesn't care about created loops; 00268 * 00269 * New line is written to existing nPoints structure. 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); /* ? OK, should make circle for points ? */ 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) { /* use polyline instead of arc between line segments */ 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 /* TODO: a <= PI can probably fail because of representation error */ 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)); /* angle step */ 00391 00392 if (Points == NULL) 00393 Points = Vect_new_line_struct(); 00394 00395 if (PPoints == NULL) 00396 PPoints = Vect_new_line_struct(); 00397 00398 /* Copy and prune input */ 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) { /* make a circle */ 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 /* Close polygon */ 00418 Vect_append_point(OutPoints, OutPoints->x[0], OutPoints->y[0], 0); 00419 } 00420 else { /* 2 and more points */ 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 /* Parallel on one side */ 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 /* Arc at the end */ 00437 /* 2 points at theend of original line */ 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 /* normalized vector */ 00452 vect(lx1, ly1, lx2, ly2, &nx, &ny); 00453 00454 /* starting point */ 00455 sangle = atan2(-nx, ny); /* starting angle */ 00456 sx = lx2 + ny * distance; 00457 sy = ly2 - nx * distance; 00458 00459 /* end point */ 00460 ex = lx2 - ny * distance; 00461 ey = ly2 + nx * distance; 00462 00463 Vect_append_point(OutPoints, sx, sy, 0); 00464 00465 /* arc */ 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 /* Close polygon */ 00476 Vect_append_point(OutPoints, OutPoints->x[0], OutPoints->y[0], 0); 00477 } 00478 Vect_line_prune(OutPoints); 00479 00480 return; 00481 }