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