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
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
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
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;
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
00192 Vect_find_poly_centroid(Points, ¢_x, ¢_y);
00193
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
00201
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
00212 if (first_time) {
00213
00214 link_exit_on_error(1);
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;
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
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
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
00297
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
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
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) {
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
00454 Vect_find_poly_centroid(Points, ¢_x, ¢_y);
00455
00456 if (Vect_point_in_poly(cent_x, cent_y, Points) == 1)
00457
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
00472
00473
00474
00475
00476
00477
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;
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
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);
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
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)
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
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)
00542 return -1;
00543
00544 *att_x = (Intersects->x[maxpos] + Intersects->x[maxpos + 1]) / 2.;
00545
00546 return 0;
00547 }
00548
00549
00550
00551
00552
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
00565
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
00578
00579
00580
00581 if (x1 < X && x2 < X)
00582 continue;
00583
00584
00585 if ((x1 == X && y1 == Y) || (x2 == X && y2 == Y))
00586 return -1;
00587
00588
00589 if ((x1 == x2 && x1 == X) &&
00590 ((y1 <= Y && y2 >= Y) || (y1 >= Y && y2 <= Y)))
00591 return -1;
00592
00593
00594 if ((y1 == y2 && y1 == Y) &&
00595 ((x1 <= X && x2 >= X) || (x1 >= X && x2 <= X)))
00596 return -1;
00597
00598
00599 if (y1 == Y && y2 == Y)
00600 continue;
00601
00602
00603 if (y1 > Y && y2 > Y)
00604 continue;
00605
00606
00607 if (y1 < Y && y2 < Y)
00608 continue;
00609
00610
00611 if ((y1 == Y && y2 > Y) || (y2 == Y && y1 > Y))
00612 continue;
00613
00614
00615
00616
00617 if (y1 == Y && y2 < Y) {
00618 if (x1 >= 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
00629 if ((y1 < Y && y2 > Y) || (y1 > Y && y2 < Y)) {
00630 if (x1 >= X && x2 >= X) {
00631 n_intersects++;
00632 continue;
00633 }
00634
00635
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;
00644 }
00645
00646
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
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
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
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 }