GRASS Programmer's Manual
6.4.1(2011)
|
00001 00021 #include <math.h> 00022 #include <stdlib.h> 00023 #include <grass/Vect.h> 00024 #include <grass/gis.h> 00025 #include <grass/linkm.h> 00026 #include <grass/glocale.h> 00027 00028 struct Slink 00029 { 00030 double x; 00031 struct Slink *next; 00032 }; 00033 00034 00035 /* function prototypes */ 00036 static int comp_double(double *, double *); 00037 static int V__within(double, double, double); 00038 int Vect__intersect_line_with_poly(); 00039 static void destroy_links(struct Slink *); 00040 static int Vect__divide_and_conquer(struct Slink *, struct line_pnts *, 00041 struct link_head *, double *, double *, 00042 int); 00043 00044 00060 int 00061 Vect_get_point_in_area(struct Map_info *Map, int area, double *X, double *Y) 00062 { 00063 static struct line_pnts *Points; 00064 static struct line_pnts **IPoints; 00065 static int first_time = 1; 00066 static int isl_allocated = 0; 00067 int i, n_isles; 00068 00069 G_debug(3, "Vect_get_point_in_area()"); 00070 00071 if (first_time) { 00072 Points = Vect_new_line_struct(); 00073 IPoints = NULL; 00074 first_time = 0; 00075 } 00076 n_isles = Vect_get_area_num_isles(Map, area); 00077 if (n_isles > isl_allocated) { 00078 IPoints = (struct line_pnts **) 00079 G_realloc(IPoints, (1 + n_isles) * sizeof(struct line_pnts *)); 00080 for (i = isl_allocated; i < n_isles; i++) 00081 IPoints[i] = Vect_new_line_struct(); 00082 isl_allocated = n_isles; 00083 } 00084 00085 if (0 > Vect_get_area_points(Map, area, Points)) 00086 return -1; 00087 00088 for (i = 0; i < n_isles; i++) { 00089 IPoints[i]->alloc_points = 0; 00090 if (0 > 00091 Vect_get_isle_points(Map, Vect_get_area_isle(Map, area, i), 00092 IPoints[i])) 00093 return -1; 00094 } 00095 return (Vect_get_point_in_poly_isl(Points, IPoints, n_isles, X, Y)); 00096 00097 return -1; 00098 } 00099 00100 static int comp_double(double *i, double *j) 00101 { 00102 if (*i < *j) 00103 return -1; 00104 00105 if (*i > *j) 00106 return 1; 00107 00108 return 0; 00109 } 00110 00111 static int V__within(double a, double x, double b) 00112 { 00113 double tmp; 00114 00115 if (a > b) { 00116 tmp = a; 00117 a = b; 00118 b = tmp; 00119 } 00120 00121 return (x >= a && x <= b); 00122 } 00123 00124 /* 00125 \brief Intersects line with polygon 00126 00127 For each intersection of a polygon w/ a line, stuff the X value in 00128 the Inter Points array. I used line_pnts, just cuz the memory 00129 management was already there. I am getting real tired of managing 00130 realloc stuff. Assumes that no vertex of polygon lies on Y This is 00131 taken care of by functions calling this function 00132 00133 \param Points line 00134 \param y ? 00135 \param Iter intersection ? 00136 00137 \return 0 on success 00138 \return -1 on error 00139 */ 00140 int 00141 Vect__intersect_line_with_poly(struct line_pnts *Points, 00142 double y, struct line_pnts *Inter) 00143 { 00144 int i; 00145 double a, b, c, d, x; 00146 double perc; 00147 00148 for (i = 1; i < Points->n_points; i++) { 00149 a = Points->y[i - 1]; 00150 b = Points->y[i]; 00151 00152 c = Points->x[i - 1]; 00153 d = Points->x[i]; 00154 00155 if (V__within(a, y, b)) { 00156 if (a == b) 00157 continue; 00158 00159 perc = (y - a) / (b - a); 00160 x = perc * (d - c) + c; /* interp X */ 00161 00162 if (0 > Vect_append_point(Inter, x, y, 0)) 00163 return -1; 00164 } 00165 } 00166 return 0; 00167 } 00168 00180 int Vect_get_point_in_poly(struct line_pnts *Points, double *X, double *Y) 00181 { 00182 double cent_x, cent_y; 00183 struct Slink *Head; 00184 static struct link_head *Token; 00185 struct Slink *tmp; 00186 static int first_time = 1; 00187 register int i; 00188 double x_max, x_min; 00189 int ret; 00190 00191 /* get centroid */ 00192 Vect_find_poly_centroid(Points, ¢_x, ¢_y); 00193 /* is it w/in poly? */ 00194 if (Vect_point_in_poly(cent_x, cent_y, Points) == 1) { 00195 *X = cent_x; 00196 *Y = cent_y; 00197 return 0; 00198 } 00199 00200 /* guess we have to do it the hard way... */ 00201 /* get min and max x values */ 00202 x_max = x_min = Points->x[0]; 00203 for (i = 0; i < Points->n_points; i++) { 00204 if (x_min > Points->x[i]) 00205 x_min = Points->x[i]; 00206 if (x_max < Points->x[i]) 00207 x_max = Points->x[i]; 00208 } 00209 00210 00211 /* init the linked list */ 00212 if (first_time) { 00213 /* will never call link_cleanup () */ 00214 link_exit_on_error(1); /* kill program if out of memory */ 00215 Token = (struct link_head *)link_init(sizeof(struct Slink)); 00216 first_time = 0; 00217 } 00218 00219 Head = (struct Slink *)link_new(Token); 00220 tmp = (struct Slink *)link_new(Token); 00221 00222 Head->next = tmp; 00223 tmp->next = NULL; 00224 00225 Head->x = x_min; 00226 tmp->x = x_max; 00227 00228 *Y = cent_y; /* pick line segment (x_min, cent_y) - (x_max, cent_y) */ 00229 ret = Vect__divide_and_conquer(Head, Points, Token, X, Y, 10); 00230 00231 destroy_links(Head); 00232 00233 if (ret < 0) { 00234 G_warning("Vect_get_point_in_poly(): %s", 00235 _("Unable to find point in polygon")); 00236 return -1; 00237 } 00238 00239 G_debug(3, "Found point in %d iterations", 10 - ret); 00240 00241 return 0; 00242 } 00243 00244 00245 /* 00246 \brief Provide a breadth first binary division of real space along line segment. 00247 00248 Looking for a point w/in the polygon. 00249 00250 This routine walks along the list of points on line segment 00251 and divides each pair in half. It sticks that new point right into 00252 the list, and then checks to see if it is inside the poly. 00253 00254 After going through the whole list, it calls itself. The list 00255 now has a whole extra set of points to divide again. 00256 00257 \param Head 00258 \param Points 00259 \param Token 00260 \param X,Y 00261 \param levels 00262 00263 \return # levels it took 00264 \return -1 if exceeded # of levels 00265 */ 00266 static int 00267 Vect__divide_and_conquer(struct Slink *Head, 00268 struct line_pnts *Points, 00269 struct link_head *Token, 00270 double *X, double *Y, int levels) 00271 { 00272 struct Slink *A, *B, *C; 00273 00274 G_debug(3, "Vect__divide_and_conquer(): LEVEL %d", levels); 00275 A = Head; 00276 B = Head->next; 00277 00278 do { 00279 C = (struct Slink *)link_new(Token); 00280 A->next = C; 00281 C->next = B; 00282 00283 C->x = (A->x + B->x) / 2.; 00284 00285 if (Vect_point_in_poly(C->x, *Y, Points) == 1) { 00286 *X = C->x; 00287 return levels; 00288 } 00289 00290 A = B; 00291 B = B->next; 00292 } 00293 while (B != NULL); 00294 00295 /* 00296 ** If it got through the entire loop and still no hits, 00297 ** then lets go a level deeper and divide again. 00298 */ 00299 00300 if (levels <= 0) 00301 return -1; 00302 00303 return Vect__divide_and_conquer(Head, Points, Token, X, Y, --levels); 00304 } 00305 00306 static void destroy_links(struct Slink *Head) 00307 { 00308 struct Slink *p, *tmp; 00309 00310 p = Head; 00311 00312 while (p != NULL) { 00313 tmp = p->next; 00314 link_dispose((struct link_head *)Head, (VOID_T *) p); 00315 p = tmp; 00316 } 00317 } 00318 00319 00329 int 00330 Vect_find_poly_centroid(struct line_pnts *points, 00331 double *cent_x, double *cent_y) 00332 { 00333 int i; 00334 double *xptr1, *yptr1; 00335 double *xptr2, *yptr2; 00336 double cent_weight_x, cent_weight_y; 00337 double len, tot_len; 00338 00339 tot_len = 0.0; 00340 cent_weight_x = 0.0; 00341 cent_weight_y = 0.0; 00342 00343 xptr1 = points->x; 00344 yptr1 = points->y; 00345 xptr2 = points->x + 1; 00346 yptr2 = points->y + 1; 00347 00348 for (i = 1; i < points->n_points; i++) { 00349 len = hypot(*xptr1 - *xptr2, *yptr1 - *yptr2); 00350 cent_weight_x += len * ((*xptr1 + *xptr2) / 2.); 00351 cent_weight_y += len * ((*yptr1 + *yptr2) / 2.); 00352 tot_len += len; 00353 xptr1++; 00354 xptr2++; 00355 yptr1++; 00356 yptr2++; 00357 } 00358 00359 if (tot_len == 0.0) 00360 return -1; 00361 00362 *cent_x = cent_weight_x / tot_len; 00363 *cent_y = cent_weight_y / tot_len; 00364 00365 return 0; 00366 } 00367 00368 00369 /* 00370 ** returns true if point is in any of islands /w in area 00371 ** returns 0 if not 00372 ** returns -1 on error 00373 */ 00374 /* 00375 int 00376 Vect_point_in_islands ( 00377 struct Map_info *Map, 00378 int area, 00379 double cent_x, double cent_y) 00380 { 00381 P_AREA *Area; 00382 static struct line_pnts *TPoints; 00383 static int first_time = 1; 00384 int isle; 00385 00386 if (first_time == 1) 00387 { 00388 TPoints = Vect_new_line_struct (); 00389 first_time = 0; 00390 } 00391 00392 Area = &(Map->plus.Area[area]); 00393 00394 for (isle = 0; isle < Area->n_isles; isle++) 00395 { 00396 if (0 > Vect_get_isle_points (Map, Area->isles[isle], TPoints)) 00397 return -1; 00398 00399 if ( Vect_point_in_poly (cent_x, cent_y, TPoints) == 1 ) 00400 return 1; 00401 } 00402 00403 return 0; 00404 } 00405 */ 00406 00423 int 00424 Vect_get_point_in_poly_isl(struct line_pnts *Points, 00425 struct line_pnts **IPoints, int n_isles, 00426 double *att_x, double *att_y) 00427 { 00428 static struct line_pnts *Intersects; 00429 static int first_time = 1; 00430 double cent_x, cent_y; 00431 register int i, j; 00432 double max, hi_y, lo_y; 00433 int maxpos; 00434 int point_in_sles = 0; 00435 double diff; 00436 00437 G_debug(3, "Vect_get_point_in_poly_isl(): n_isles = %d", n_isles); 00438 00439 if (first_time) { 00440 Intersects = Vect_new_line_struct(); 00441 first_time = 0; 00442 } 00443 00444 if (Points->n_points < 3) { /* test */ 00445 if (Points->n_points > 0) { 00446 *att_x = Points->x[0]; 00447 *att_y = Points->y[0]; 00448 return 0; 00449 } 00450 return -1; 00451 } 00452 00453 /* get centroid */ 00454 Vect_find_poly_centroid(Points, ¢_x, ¢_y); 00455 /* is it w/in poly? */ 00456 if (Vect_point_in_poly(cent_x, cent_y, Points) == 1) 00457 /* if the point is iside the polygon */ 00458 { 00459 for (i = 0; i < n_isles; i++) { 00460 if (Vect_point_in_poly(cent_x, cent_y, IPoints[i]) >= 1) { 00461 point_in_sles = 1; 00462 break; 00463 } 00464 } 00465 if (!point_in_sles) { 00466 *att_x = cent_x; 00467 *att_y = cent_y; 00468 return 0; 00469 } 00470 } 00471 /* guess we have to do it the hard way... */ 00472 00473 /* first find att_y close to cent_y so that no points lie on the line */ 00474 /* find the point closest to line from below, and point close to line 00475 from above and take average of their y-coordinates */ 00476 00477 /* first initializing lo_y,hi_y to be any 2 pnts on either side of cent_y */ 00478 hi_y = cent_y - 1; 00479 lo_y = cent_y + 1; 00480 for (i = 0; i < Points->n_points; i++) { 00481 if ((lo_y < cent_y) && (hi_y >= cent_y)) 00482 break; /* already initialized */ 00483 if (Points->y[i] < cent_y) 00484 lo_y = Points->y[i]; 00485 if (Points->y[i] >= cent_y) 00486 hi_y = Points->y[i]; 00487 } 00488 /* first going throught boundary points */ 00489 for (i = 0; i < Points->n_points; i++) { 00490 if ((Points->y[i] < cent_y) && 00491 ((cent_y - Points->y[i]) < (cent_y - lo_y))) 00492 lo_y = Points->y[i]; 00493 if ((Points->y[i] >= cent_y) && 00494 ((Points->y[i] - cent_y) < (hi_y - cent_y))) 00495 hi_y = Points->y[i]; 00496 } 00497 for (i = 0; i < n_isles; i++) 00498 for (j = 0; j < IPoints[i]->n_points; j++) { 00499 if ((IPoints[i]->y[j] < cent_y) && 00500 ((cent_y - IPoints[i]->y[j]) < (cent_y - lo_y))) 00501 lo_y = IPoints[i]->y[j]; 00502 00503 if ((IPoints[i]->y[j] >= cent_y) && 00504 ((IPoints[i]->y[j] - cent_y) < (hi_y - cent_y))) 00505 hi_y = IPoints[i]->y[j]; 00506 } 00507 00508 if (lo_y == hi_y) 00509 return (-1); /* area is empty */ 00510 else 00511 *att_y = (hi_y + lo_y) / 2.0; 00512 00513 Intersects->n_points = 0; 00514 Vect__intersect_line_with_poly(Points, *att_y, Intersects); 00515 00516 /* add in intersections w/ holes */ 00517 for (i = 0; i < n_isles; i++) { 00518 if (0 > 00519 Vect__intersect_line_with_poly(IPoints[i], *att_y, Intersects)) 00520 return -1; 00521 } 00522 00523 if (Intersects->n_points < 2) /* test */ 00524 return -1; 00525 00526 qsort(Intersects->x, (size_t) Intersects->n_points, sizeof(double), 00527 (void *)comp_double); 00528 00529 max = 0; 00530 maxpos = 0; 00531 00532 /* find area of MAX distance */ 00533 for (i = 0; i < Intersects->n_points; i += 2) { 00534 diff = Intersects->x[i + 1] - Intersects->x[i]; 00535 00536 if (diff > max) { 00537 max = diff; 00538 maxpos = i; 00539 } 00540 } 00541 if (max == 0.0) /* area was empty: example ((x1,y1), (x2,y2), (x1,y1)) */ 00542 return -1; 00543 00544 *att_x = (Intersects->x[maxpos] + Intersects->x[maxpos + 1]) / 2.; 00545 00546 return 0; 00547 } 00548 00549 00550 /* Intersect segments of Points with ray from point X,Y to the right. 00551 * Returns: -1 point exactly on segment 00552 * number of intersections 00553 */ 00554 static int segments_x_ray(double X, double Y, struct line_pnts *Points) 00555 { 00556 double x1, x2, y1, y2; 00557 double x_inter; 00558 int n_intersects; 00559 int n; 00560 00561 G_debug(3, "segments_x_ray(): x = %f y = %f n_points = %d", X, Y, 00562 Points->n_points); 00563 00564 /* Follow the ray from X,Y along positive x and find number of intersections. 00565 * Coordinates exactly on ray are considered to be slightly above. */ 00566 00567 n_intersects = 0; 00568 for (n = 0; n < Points->n_points - 1; n++) { 00569 x1 = Points->x[n]; 00570 y1 = Points->y[n]; 00571 x2 = Points->x[n + 1]; 00572 y2 = Points->y[n + 1]; 00573 00574 G_debug(3, "X = %f Y = %f x1 = %f y1 = %f x2 = %f y2 = %f", X, Y, x1, 00575 y1, x2, y2); 00576 00577 /* I know, it should be possible to do that with less conditions, but it should be 00578 * enough readable also! */ 00579 00580 /* segment left from X -> no intersection */ 00581 if (x1 < X && x2 < X) 00582 continue; 00583 00584 /* point on vertex */ 00585 if ((x1 == X && y1 == Y) || (x2 == X && y2 == Y)) 00586 return -1; 00587 00588 /* on vertical boundary */ 00589 if ((x1 == x2 && x1 == X) && 00590 ((y1 <= Y && y2 >= Y) || (y1 >= Y && y2 <= Y))) 00591 return -1; 00592 00593 /* on horizontal boundary */ 00594 if ((y1 == y2 && y1 == Y) && 00595 ((x1 <= X && x2 >= X) || (x1 >= X && x2 <= X))) 00596 return -1; 00597 00598 /* segment on ray (X is not important) */ 00599 if (y1 == Y && y2 == Y) 00600 continue; 00601 00602 /* segment above (X is not important) */ 00603 if (y1 > Y && y2 > Y) 00604 continue; 00605 00606 /* segment below (X is not important) */ 00607 if (y1 < Y && y2 < Y) 00608 continue; 00609 00610 /* one end on Y second above (X is not important) */ 00611 if ((y1 == Y && y2 > Y) || (y2 == Y && y1 > Y)) 00612 continue; 00613 00614 /* For following cases we know that at least one of x1 and x2 is >= X */ 00615 00616 /* one end of segment on Y second below Y */ 00617 if (y1 == Y && y2 < Y) { 00618 if (x1 >= X) /* x of the end on the ray is >= X */ 00619 n_intersects++; 00620 continue; 00621 } 00622 if (y2 == Y && y1 < Y) { 00623 if (x2 >= X) 00624 n_intersects++; 00625 continue; 00626 } 00627 00628 /* one end of segment above Y second below Y */ 00629 if ((y1 < Y && y2 > Y) || (y1 > Y && y2 < Y)) { 00630 if (x1 >= X && x2 >= X) { 00631 n_intersects++; 00632 continue; 00633 } 00634 00635 /* now either x1 < X && x2 > X or x1 > X && x2 < X -> calculate intersection */ 00636 x_inter = dig_x_intersect(x1, x2, y1, y2, Y); 00637 G_debug(3, "x_inter = %f", x_inter); 00638 if (x_inter == X) 00639 return 1; 00640 else if (x_inter > X) 00641 n_intersects++; 00642 00643 continue; /* would not be necessary, just to check, see below */ 00644 } 00645 /* should not be reached (one condition is not necessary, but it is may be better readable 00646 * and it is a check) */ 00647 G_warning 00648 ("segments_x_ray() %s: X = %f Y = %f x1 = %f y1 = %f x2 = %f y2 = %f", 00649 _("conditions failed"), X, Y, x1, y1, x2, y2); 00650 } 00651 00652 return n_intersects; 00653 } 00654 00665 int Vect_point_in_poly(double X, double Y, struct line_pnts *Points) 00666 { 00667 int n_intersects; 00668 00669 G_debug(3, "Vect_point_in_poly(): x = %f y = %f n_points = %d", X, Y, 00670 Points->n_points); 00671 00672 n_intersects = segments_x_ray(X, Y, Points); 00673 00674 if (n_intersects == -1) 00675 return 2; 00676 00677 if (n_intersects % 2) 00678 return 1; 00679 else 00680 return 0; 00681 } 00682 00694 int 00695 Vect_point_in_area_outer_ring(double X, double Y, struct Map_info *Map, 00696 int area) 00697 { 00698 static int first = 1; 00699 int n_intersects, inter; 00700 int i, line; 00701 static struct line_pnts *Points; 00702 struct Plus_head *Plus; 00703 P_LINE *Line; 00704 P_AREA *Area; 00705 00706 G_debug(3, "Vect_point_in_area_outer_ring(): x = %f y = %f area = %d", X, 00707 Y, area); 00708 00709 if (first == 1) { 00710 Points = Vect_new_line_struct(); 00711 first = 0; 00712 } 00713 00714 Plus = &(Map->plus); 00715 Area = Plus->Area[area]; 00716 00717 /* First it must be in box */ 00718 if (X < Area->W || X > Area->E || Y > Area->N || Y < Area->S) 00719 return 0; 00720 00721 n_intersects = 0; 00722 for (i = 0; i < Area->n_lines; i++) { 00723 line = abs(Area->lines[i]); 00724 G_debug(3, " line[%d] = %d", i, line); 00725 00726 Line = Plus->Line[line]; 00727 00728 /* dont check lines that obviously do not intersect with test ray */ 00729 if ((Line->N < Y) || (Line->S > Y) || (Line->E < X)) 00730 continue; 00731 00732 Vect_read_line(Map, Points, NULL, line); 00733 00734 inter = segments_x_ray(X, Y, Points); 00735 G_debug(3, " inter = %d", inter); 00736 00737 if (inter == -1) 00738 return 2; 00739 n_intersects += inter; 00740 G_debug(3, " n_intersects = %d", n_intersects); 00741 } 00742 00743 if (n_intersects % 2) 00744 return 1; 00745 else 00746 return 0; 00747 } 00748 00760 int Vect_point_in_island(double X, double Y, struct Map_info *Map, int isle) 00761 { 00762 static int first = 1; 00763 int n_intersects, inter; 00764 int i, line; 00765 static struct line_pnts *Points; 00766 struct Plus_head *Plus; 00767 P_LINE *Line; 00768 P_ISLE *Isle; 00769 00770 G_debug(3, "Vect_point_in_island(): x = %f y = %f isle = %d", X, Y, isle); 00771 00772 if (first == 1) { 00773 Points = Vect_new_line_struct(); 00774 first = 0; 00775 } 00776 00777 Plus = &(Map->plus); 00778 Isle = Plus->Isle[isle]; 00779 00780 if (X < Isle->W || X > Isle->E || Y > Isle->N || Y < Isle->S) 00781 return 0; 00782 00783 n_intersects = 0; 00784 for (i = 0; i < Isle->n_lines; i++) { 00785 line = abs(Isle->lines[i]); 00786 00787 Line = Plus->Line[line]; 00788 00789 /* dont check lines that obviously do not intersect with test ray */ 00790 if ((Line->N < Y) || (Line->S > Y) || (Line->E < X)) 00791 continue; 00792 00793 Vect_read_line(Map, Points, NULL, line); 00794 00795 inter = segments_x_ray(X, Y, Points); 00796 if (inter == -1) 00797 return 2; 00798 n_intersects += inter; 00799 } 00800 00801 if (n_intersects % 2) 00802 return 1; 00803 else 00804 return 0; 00805 }