00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include <math.h>
00019 #include <stdlib.h>
00020 #include <grass/Vect.h>
00021 #include <grass/gis.h>
00022 #include <grass/linkm.h>
00023
00024 struct Slink
00025 {
00026 double x;
00027 struct Slink *next;
00028 };
00029
00030
00031
00032 static int comp_double (double *, double *);
00033 static int V__within (double, double, double);
00034 int Vect__intersect_line_with_poly ();
00035 static void destroy_links (struct Slink *);
00036 static int Vect__divide_and_conquer (struct Slink *, struct line_pnts *,
00037 struct link_head *, double *, double *, int);
00038
00039
00050 int
00051 Vect_get_point_in_area ( struct Map_info *Map, int area, double *X, double *Y)
00052 {
00053 static struct line_pnts *Points;
00054 static struct line_pnts **IPoints;
00055 static int first_time = 1;
00056 static int isl_allocated = 0;
00057 int i, n_isles;
00058
00059 G_debug ( 3, "Vect_get_point_in_area()" );
00060
00061 if (first_time) {
00062 Points = Vect_new_line_struct ();
00063 IPoints = NULL;
00064 first_time = 0;
00065 }
00066 n_isles = Vect_get_area_num_isles ( Map, area);
00067 if ( n_isles > isl_allocated) {
00068 IPoints = (struct line_pnts **)
00069 G_realloc (IPoints, (1 + n_isles) * sizeof (struct line_pnts *));
00070 for (i = isl_allocated; i < n_isles; i++)
00071 IPoints[i] = Vect_new_line_struct ();
00072 isl_allocated = n_isles;
00073 }
00074
00075 if (0 > Vect_get_area_points (Map, area, Points))
00076 return -1;
00077
00078 for (i = 0; i < n_isles; i++) {
00079 IPoints[i]->alloc_points = 0;
00080 if (0 > Vect_get_isle_points (Map, Vect_get_area_isle(Map, area, i), IPoints[i]))
00081 return -1;
00082 }
00083 return (Vect_get_point_in_poly_isl (Points, IPoints, n_isles, X, Y));
00084
00085 return -1;
00086 }
00087
00088 static int
00089 comp_double (double *i, double *j)
00090 {
00091 if (*i < *j)
00092 return -1;
00093
00094 if (*i > *j)
00095 return 1;
00096
00097 return 0;
00098 }
00099
00100 static int
00101 V__within (double a, double x, double b)
00102 {
00103 double tmp;
00104
00105 if (a > b)
00106 {
00107 tmp = a;
00108 a = b;
00109 b = tmp;
00110 }
00111
00112 return (x >= a && x <= b);
00113 }
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126 int
00127 Vect__intersect_line_with_poly (
00128 struct line_pnts *Points,
00129 double y,
00130 struct line_pnts *Inter)
00131 {
00132 int i;
00133 double a, b, c, d, x;
00134 double perc;
00135
00136 for (i = 1; i < Points->n_points; i++)
00137 {
00138 a = Points->y[i - 1];
00139 b = Points->y[i];
00140
00141 c = Points->x[i - 1];
00142 d = Points->x[i];
00143
00144 if (V__within (a, y, b))
00145 {
00146 if (a == b)
00147 continue;
00148
00149 perc = (y - a) / (b - a);
00150 x = perc * (d - c) + c;
00151
00152 if (0 > Vect_append_point (Inter, x, y, 0))
00153 return -1;
00154 }
00155 }
00156 return 0;
00157 }
00158
00165 int
00166 Vect_get_point_in_poly (struct line_pnts *Points, double *X, double *Y)
00167 {
00168 double cent_x, cent_y;
00169 struct Slink *Head;
00170 static struct link_head *Token;
00171 struct Slink *tmp;
00172 static int first_time = 1;
00173 register int i;
00174 double x_max, x_min;
00175 int ret;
00176
00177
00178 Vect_find_poly_centroid (Points, ¢_x, ¢_y);
00179
00180 if ( Vect_point_in_poly (cent_x, cent_y, Points) == 1 )
00181 {
00182 *X = cent_x;
00183 *Y = cent_y;
00184 return 0;
00185 }
00186
00187
00188
00189 x_max = x_min = Points->x[0];
00190 for (i = 0; i < Points->n_points; i++)
00191 {
00192 if (x_min > Points->x[i])
00193 x_min = Points->x[i];
00194 if (x_max < Points->x[i])
00195 x_max = Points->x[i];
00196 }
00197
00198
00199
00200 if (first_time)
00201 {
00202
00203 link_exit_on_error (1);
00204 Token = (struct link_head *) link_init (sizeof (struct Slink));
00205 first_time = 0;
00206 }
00207
00208 Head = (struct Slink *) link_new (Token);
00209 tmp = (struct Slink *) link_new (Token);
00210
00211 Head->next = tmp;
00212 tmp->next = NULL;
00213
00214 Head->x = x_min;
00215 tmp->x = x_max;
00216
00217 *Y = cent_y;
00218 ret = Vect__divide_and_conquer (Head, Points, Token, X, Y, 10);
00219
00220 destroy_links (Head);
00221
00222 if (ret < 0)
00223 {
00224 fprintf (stderr, "Could not find point in polygon\n");
00225 return -1;
00226 }
00227
00228
00229
00230 return 0;
00231 }
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247 static int
00248 Vect__divide_and_conquer (
00249 struct Slink *Head,
00250 struct line_pnts *Points,
00251 struct link_head *Token,
00252 double *X, double *Y,
00253 int levels)
00254 {
00255 struct Slink *A, *B, *C;
00256
00257
00258 A = Head;
00259 B = Head->next;
00260
00261 do
00262 {
00263 C = (struct Slink *) link_new (Token);
00264 A->next = C;
00265 C->next = B;
00266
00267 C->x = (A->x + B->x) / 2.;
00268
00269 if ( Vect_point_in_poly (C->x, *Y, Points) == 1 )
00270 {
00271 *X = C->x;
00272 return levels;
00273 }
00274
00275 A = B;
00276 B = B->next;
00277 }
00278 while (B != NULL);
00279
00280
00281
00282
00283
00284
00285 if (levels <= 0)
00286 return -1;
00287
00288 return Vect__divide_and_conquer (Head, Points, Token, X, Y, --levels);
00289 }
00290
00291
00292 static void
00293 destroy_links (struct Slink *Head)
00294 {
00295 struct Slink *p, *tmp;
00296
00297 p = Head;
00298
00299 while (p != NULL)
00300 {
00301 tmp = p->next;
00302 link_dispose ((struct link_head *) Head, (VOID_T *) p);
00303 p = tmp;
00304 }
00305 }
00306
00307
00315 int
00316 Vect_find_poly_centroid (
00317 struct line_pnts *points,
00318 double *cent_x, double *cent_y)
00319 {
00320 int i;
00321 double *xptr1, *yptr1;
00322 double *xptr2, *yptr2;
00323 double cent_weight_x, cent_weight_y;
00324 double len, tot_len;
00325
00326 tot_len = 0.0;
00327 cent_weight_x = 0.0;
00328 cent_weight_y = 0.0;
00329
00330 xptr1 = points->x;
00331 yptr1 = points->y;
00332 xptr2 = points->x + 1;
00333 yptr2 = points->y + 1;
00334
00335 for (i = 1; i < points->n_points; i++)
00336 {
00337 len = hypot (*xptr1 - *xptr2, *yptr1 - *yptr2);
00338 cent_weight_x += len * ((*xptr1 + *xptr2) / 2.);
00339 cent_weight_y += len * ((*yptr1 + *yptr2) / 2.);
00340 tot_len += len;
00341 xptr1++;
00342 xptr2++;
00343 yptr1++;
00344 yptr2++;
00345 }
00346
00347 if (tot_len == 0.0)
00348 return -1;
00349
00350 *cent_x = cent_weight_x / tot_len;
00351 *cent_y = cent_weight_y / tot_len;
00352
00353 return 0;
00354 }
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
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
00407 int
00408 Vect_get_point_in_poly_isl (
00409 struct line_pnts *Points, struct line_pnts **IPoints,
00410 int n_isles,
00411 double *att_x, double *att_y)
00412 {
00413 static struct line_pnts *Intersects;
00414 static int first_time = 1;
00415 double cent_x, cent_y;
00416 register int i, j;
00417 double max, hi_y, lo_y;
00418 int maxpos;
00419 int point_in_sles = 0;
00420 double diff;
00421
00422 G_debug ( 3, "Vect_get_point_in_poly_isl(): n_isles = %d", n_isles );
00423
00424 if (first_time)
00425 {
00426 Intersects = Vect_new_line_struct ();
00427 first_time = 0;
00428 }
00429
00430 if (Points->n_points < 3)
00431 {
00432 if (Points->n_points > 0)
00433 {
00434 *att_x = Points->x[0];
00435 *att_y = Points->y[0];
00436 return 0;
00437 }
00438 return -1;
00439 }
00440
00441
00442 Vect_find_poly_centroid (Points, ¢_x, ¢_y);
00443
00444 if ( Vect_point_in_poly (cent_x, cent_y, Points) == 1)
00445
00446 {
00447 for (i = 0; i < n_isles; i++)
00448 {
00449 if (Vect_point_in_poly (cent_x, cent_y, IPoints[i]) >= 1) {
00450 point_in_sles = 1;
00451 break;
00452 }
00453 }
00454 if (!point_in_sles)
00455 {
00456 *att_x = cent_x;
00457 *att_y = cent_y;
00458 return 0;
00459 }
00460 }
00461
00462
00463
00464
00465
00466
00467
00468 hi_y = cent_y - 1;
00469 lo_y = cent_y + 1;
00470 for (i = 0; i < Points->n_points; i++)
00471 {
00472 if ((lo_y < cent_y) && (hi_y >= cent_y))
00473 break;
00474 if (Points->y[i] < cent_y)
00475 lo_y = Points->y[i];
00476 if (Points->y[i] >= cent_y)
00477 hi_y = Points->y[i];
00478 }
00479
00480 for (i = 0; i < Points->n_points; i++)
00481 {
00482 if ((Points->y[i] < cent_y) && ((cent_y - Points->y[i]) < (cent_y - lo_y)))
00483 lo_y = Points->y[i];
00484 if ((Points->y[i] >= cent_y) && ((Points->y[i] - cent_y) < (hi_y - cent_y)))
00485 hi_y = Points->y[i];
00486 }
00487 for (i = 0; i < n_isles; i++)
00488 for (j = 0; j < IPoints[i]->n_points; j++)
00489 {
00490 if ((IPoints[i]->y[j] < cent_y) &&
00491 ((cent_y - IPoints[i]->y[j]) < (cent_y - lo_y)))
00492 lo_y = IPoints[i]->y[j];
00493
00494 if ((IPoints[i]->y[j] >= cent_y) &&
00495 ((IPoints[i]->y[j] - cent_y) < (hi_y - cent_y)))
00496 hi_y = IPoints[i]->y[j];
00497 }
00498
00499 if (lo_y == hi_y)
00500 return (-1);
00501 else
00502 *att_y = (hi_y + lo_y) / 2.0;
00503
00504 Intersects->n_points = 0;
00505 Vect__intersect_line_with_poly (Points, *att_y, Intersects);
00506
00507
00508 for (i = 0; i < n_isles; i++)
00509 {
00510 if (0 > Vect__intersect_line_with_poly (IPoints[i], *att_y, Intersects))
00511 return -1;
00512 }
00513
00514 if (Intersects->n_points < 2)
00515 return -1;
00516
00517 qsort (Intersects->x, (size_t)Intersects->n_points, sizeof (double), (void *) comp_double);
00518
00519 max = 0;
00520 maxpos = 0;
00521
00522
00523 for (i = 0; i < Intersects->n_points; i += 2)
00524 {
00525 diff = Intersects->x[i + 1] - Intersects->x[i];
00526
00527 if (diff > max)
00528 {
00529 max = diff;
00530 maxpos = i;
00531 }
00532 }
00533 if (max == 0.0)
00534 return -1;
00535
00536 *att_x = (Intersects->x[maxpos] + Intersects->x[maxpos + 1]) / 2.;
00537
00538 return 0;
00539 }
00540
00541
00542
00543
00544
00545
00546 int
00547 segments_x_ray ( double X, double Y, struct line_pnts *Points)
00548 {
00549 double x1, x2, y1, y2;
00550 double x_inter;
00551 int n_intersects;
00552 int n;
00553
00554 G_debug ( 3, "segments_x_ray(): x = %f y = %f n_points = %d", X, Y, Points->n_points );
00555
00556
00557
00558
00559 n_intersects = 0;
00560 for ( n = 0; n < Points->n_points-1; n++) {
00561 x1 = Points->x[n];
00562 y1 = Points->y[n];
00563 x2 = Points->x[n+1];
00564 y2 = Points->y[n+1];
00565
00566 G_debug ( 3, "X = %f Y = %f x1 = %f y1 = %f x2 = %f y2 = %f", X, Y, x1, y1, x2, y2 );
00567
00568
00569
00570
00571
00572 if ( x1 < X && x2 < X ) continue;
00573
00574
00575 if ( (x1 == X && y1 == Y) || (x2 == X && y2 == Y) ) return -1;
00576
00577
00578 if ( (x1 == x2 && x1 == X) && ( (y1 <= Y && y2 >= Y) || (y1 >= Y && y2 <= Y) ) ) return -1;
00579
00580
00581 if ( (y1 == y2 && y1 == Y) && ( (x1 <= X && x2 >= X) || (x1 >= X && x2 <= X) ) ) return -1;
00582
00583
00584 if ( y1 == Y && y2 == Y ) continue;
00585
00586
00587 if ( y1 > Y && y2 > Y ) continue;
00588
00589
00590 if ( y1 < Y && y2 < Y ) continue;
00591
00592
00593 if ( (y1 == Y && y2 > Y) || (y2 == Y && y1 > Y) ) continue;
00594
00595
00596
00597
00598 if ( y1 == Y && y2 < Y) {
00599 if ( x1 >= X)
00600 n_intersects++;
00601 continue;
00602 }
00603 if ( y2 == Y && y1 < Y ) {
00604 if ( x2 >= X)
00605 n_intersects++;
00606 continue;
00607 }
00608
00609
00610 if ( (y1 < Y && y2 > Y) || (y1 > Y && y2 < Y) ) {
00611 if ( x1 >= X && x2 >= X ) {
00612 n_intersects++;
00613 continue;
00614 }
00615
00616
00617 x_inter = dig_x_intersect ( x1, x2, y1, y2, Y);
00618 G_debug ( 3, "x_inter = %f", x_inter );
00619 if ( x_inter == X )
00620 return 1;
00621 else if (x_inter > X)
00622 n_intersects++;
00623
00624 continue;
00625 }
00626
00627
00628 G_warning ( "segments_x_ray() conditions failed:" );
00629 G_warning ( "X = %f Y = %f x1 = %f y1 = %f x2 = %f y2 = %f", X, Y, x1, y1, x2, y2 );
00630 }
00631
00632 return n_intersects;
00633 }
00634
00635
00636
00637
00638
00639
00640
00641
00642 int
00643 Vect_point_in_poly ( double X, double Y, struct line_pnts *Points)
00644 {
00645 int n_intersects;
00646
00647 G_debug ( 3, "Vect_point_in_poly(): x = %f y = %f n_points = %d", X, Y, Points->n_points );
00648
00649 n_intersects = segments_x_ray ( X, Y, Points);
00650
00651 if ( n_intersects == -1 ) return 2;
00652
00653 if (n_intersects % 2)
00654 return 1;
00655 else
00656 return 0;
00657 }
00658
00659
00660
00661
00662
00663
00664
00665
00666 int
00667 Vect_point_in_area_outer_ring ( double X, double Y, struct Map_info *Map, int area)
00668 {
00669 static int first = 1;
00670 int n_intersects, inter;
00671 int i, line;
00672 static struct line_pnts *Points;
00673 struct Plus_head *Plus;
00674 P_LINE *Line;
00675 P_AREA *Area;
00676
00677 G_debug ( 3, "Vect_point_in_area_outer_ring(): x = %f y = %f area = %d", X, Y, area );
00678
00679 if (first == 1) {
00680 Points = Vect_new_line_struct();
00681 first = 0;
00682 }
00683
00684 Plus = &(Map->plus);
00685 Area = Plus->Area[area];
00686
00687
00688 if ( X < Area->W || X > Area->E || Y > Area->N || Y < Area->S ) return 0;
00689
00690 n_intersects = 0;
00691 for (i = 0; i < Area->n_lines; i++) {
00692 line = abs(Area->lines[i]);
00693 G_debug ( 3, " line[%d] = %d", i, line );
00694
00695 Line = Plus->Line[line];
00696
00697
00698 if ((Line->N < Y) || (Line->S > Y) || (Line->E < X)) continue;
00699
00700 Vect_read_line (Map, Points, NULL, line );
00701
00702 inter = segments_x_ray ( X, Y, Points);
00703 G_debug ( 3, " inter = %d", inter );
00704
00705 if ( inter == -1 ) return 2;
00706 n_intersects += inter;
00707 G_debug ( 3, " n_intersects = %d", n_intersects );
00708 }
00709
00710 if (n_intersects % 2)
00711 return 1;
00712 else
00713 return 0;
00714 }
00715
00716
00717
00718
00719
00720
00721
00722
00723 int
00724 Vect_point_in_island ( double X, double Y, struct Map_info *Map, int isle)
00725 {
00726 static int first = 1;
00727 int n_intersects, inter;
00728 int i, line;
00729 static struct line_pnts *Points;
00730 struct Plus_head *Plus;
00731 P_LINE *Line;
00732 P_ISLE *Isle;
00733
00734 G_debug ( 3, "Vect_point_in_island(): x = %f y = %f isle = %d", X, Y, isle );
00735
00736 if (first == 1) {
00737 Points = Vect_new_line_struct();
00738 first = 0;
00739 }
00740
00741 Plus = &(Map->plus);
00742 Isle = Plus->Isle[isle];
00743
00744 if ( X < Isle->W || X > Isle->E || Y > Isle->N || Y < Isle->S ) return 0;
00745
00746 n_intersects = 0;
00747 for (i = 0; i < Isle->n_lines; i++) {
00748 line = abs(Isle->lines[i]);
00749
00750 Line = Plus->Line[line];
00751
00752
00753 if ((Line->N < Y) || (Line->S > Y) || (Line->E < X)) continue;
00754
00755 Vect_read_line (Map, Points, NULL, line );
00756
00757 inter = segments_x_ray ( X, Y, Points);
00758 if ( inter == -1 ) return 2;
00759 n_intersects += inter;
00760 }
00761
00762 if (n_intersects % 2)
00763 return 1;
00764 else
00765 return 0;
00766 }
00767