GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-24365e1d47
Vlib/poly.c
Go to the documentation of this file.
1 /*!
2  \file lib/vector/Vlib/poly.c
3
4  \brief Vector library - polygon related fns
5
6  Higher level functions for reading/writing/manipulating vectors.
7
8  (C) 2001-2009 by the GRASS Development Team
9
10  This program is free software under the GNU General Public License
11  (>=v2). Read the file COPYING that comes with GRASS for details.
12
13  \author Original author CERL, probably Dave Gerdes or Mike Higgins.
14  \author Update to GRASS 5.7 Radim Blazek and David D. Gray.
15  */
16
17 #include <math.h>
18 #include <stdlib.h>
19 #include <grass/vector.h>
21 #include <grass/glocale.h>
22
25  double x;
26 };
27
28 /* function prototypes */
29 static int comp_double(const void *, const void *);
30 static int V__within(double, double, double);
31 int Vect__intersect_y_line_with_poly(const struct line_pnts *, double,
32  struct line_pnts *);
33 int Vect__intersect_x_line_with_poly(const struct line_pnts *, double,
34  struct line_pnts *);
36 static int Vect__divide_and_conquer(struct Slink *, const struct line_pnts *,
38  int);
39
40 /*!
41  \brief Get point inside area and outside all islands.
42
43  Take a line and intersect it with the polygon and any islands.
44  sort the list of X values from these intersections. This will
45  be a list of segments alternating IN/OUT/IN/OUT of the polygon.
46  Pick the largest IN segment and take the midpoint.
47
48  \param Map vector map
49  \param area area id
50  \param[out] X,Y point coordinateds
51
52  \return 0 on success
53  \return -1 on error
54  */
55 int Vect_get_point_in_area(struct Map_info *Map, int area, double *X, double *Y)
56 {
57  static struct line_pnts *Points;
58  static struct line_pnts **IPoints;
59  static int first_time = 1;
60  static int isl_allocated = 0;
61  int i, n_isles;
62
63  G_debug(3, "Vect_get_point_in_area()");
64
65  if (first_time) {
66  Points = Vect_new_line_struct();
67  IPoints = NULL;
68  first_time = 0;
69  }
70  n_isles = Vect_get_area_num_isles(Map, area);
71  if (n_isles > isl_allocated) {
72  IPoints = (struct line_pnts **)G_realloc(
73  IPoints, (1 + n_isles) * sizeof(struct line_pnts *));
74  for (i = isl_allocated; i < n_isles; i++)
75  IPoints[i] = Vect_new_line_struct();
76  isl_allocated = n_isles;
77  }
78
79  if (0 > Vect_get_area_points(Map, area, Points))
80  return -1;
81
82  for (i = 0; i < n_isles; i++) {
83  IPoints[i]->alloc_points = 0;
84  if (0 > Vect_get_isle_points(Map, Vect_get_area_isle(Map, area, i),
85  IPoints[i]))
86  return -1;
87  }
88  return (Vect_get_point_in_poly_isl((const struct line_pnts *)Points,
89  (const struct line_pnts **)IPoints,
90  n_isles, X, Y));
91
92  return -1;
93 }
94
95 static int comp_double(const void *i, const void *j)
96 {
97  if (*(const double *)i < *(const double *)j)
98  return -1;
99
100  return (*(const double *)i > *(const double *)j);
101 }
102
103 static int V__within(double a, double x, double b)
104 {
105  if (a < b)
106  return (x >= a && x < b);
107
108  return (x > b && x <= a);
109 }
110
111 /*
112  \brief Intersects line with polygon
113
114  For each intersection of a polygon w/ a line, stuff the X value in
115  the Inter Points array. I used line_pnts, just cuz the memory
116  management was already there. I am getting real tired of managing
117  realloc stuff.
118
119  \param Points line
120  \param y y coordinate of horizontal line
121  \param Inter intersections of horizontal line with points line
122
123  \return 0 on success
124  \return -1 on error
125  */
126 int Vect__intersect_y_line_with_poly(const struct line_pnts *Points, double y,
127  struct line_pnts *Inter)
128 {
129  int i;
130  double a, b, c, d, x;
131  double perc;
132
133  for (i = 1; i < Points->n_points; i++) {
134  a = Points->y[i - 1];
135  b = Points->y[i];
136
137  c = Points->x[i - 1];
138  d = Points->x[i];
139
140  if (V__within(a, y, b)) {
141  if (a == b)
142  continue;
143
144  perc = (y - a) / (b - a);
145  x = perc * (d - c) + c; /* interp X */
146
147  if (0 > Vect_append_point(Inter, x, y, 0))
148  return -1;
149  }
150  }
151  return 0;
152 }
153
154 /*
155  \brief Intersects line with polygon
156
157  For each intersection of a polygon w/ a line, stuff the Y value in
158  the Inter Points array. I used line_pnts, just cuz the memory
159  management was already there. I am getting real tired of managing
160  realloc stuff.
161
162  \param Points line
163  \param x x coordinate of vertical line
164  \param Inter intersections of horizontal line with points line
165
166  \return 0 on success
167  \return -1 on error
168  */
169 int Vect__intersect_x_line_with_poly(const struct line_pnts *Points, double x,
170  struct line_pnts *Inter)
171 {
172  int i;
173  double a, b, c, d, y;
174  double perc;
175
176  for (i = 1; i < Points->n_points; i++) {
177  a = Points->x[i - 1];
178  b = Points->x[i];
179
180  c = Points->y[i - 1];
181  d = Points->y[i];
182
183  if (V__within(a, x, b)) {
184  if (a == b)
185  continue;
186
187  perc = (x - a) / (b - a);
188  y = perc * (d - c) + c; /* interp Y */
189
190  if (0 > Vect_append_point(Inter, x, y, 0))
191  return -1;
192  }
193  }
194  return 0;
195 }
196
197 /*!
198  \brief Get point inside polygon.
199
200  This does NOT consider ISLANDS!
201
202  \param Points polygon
203  \param[out] X,Y point coordinates
204
205  \return 0 on success
206  \return -1 on error
207  */
208 int Vect_get_point_in_poly(const struct line_pnts *Points, double *X, double *Y)
209 {
210  double cent_x, cent_y;
214  static int first_time = 1;
215  register int i;
216  double x_max, x_min;
217  int ret;
218
219  /* get centroid */
220  Vect_find_poly_centroid(Points, &cent_x, &cent_y);
221  /* is it w/in poly? */
222  if (Vect_point_in_poly(cent_x, cent_y, Points) == 1) {
223  *X = cent_x;
224  *Y = cent_y;
225  return 0;
226  }
227
228  /* guess we have to do it the hard way... */
229  G_debug(3, "Vect_get_point_in_poly(): divide and conquer");
230
231  /* get min and max x values */
232  x_max = x_min = Points->x[0];
233  for (i = 0; i < Points->n_points; i++) {
234  if (x_min > Points->x[i])
235  x_min = Points->x[i];
236  if (x_max < Points->x[i])
237  x_max = Points->x[i];
238  }
239
240  /* init the linked list */
241  if (first_time) {
242  /* will never call link_cleanup () */
243  link_exit_on_error(1); /* kill program if out of memory */
245  first_time = 0;
246  }
247
250
252  tmp->next = NULL;
253
255  tmp->x = x_max;
256
257  *Y = cent_y; /* pick line segment (x_min, cent_y) - (x_max, cent_y) */
258  ret = Vect__divide_and_conquer(Head, Points, Token, X, Y, 10);
259
261
262  if (ret < 0) {
263  G_warning("Vect_get_point_in_poly(): %s",
264  _("Unable to find point in polygon"));
265  return -1;
266  }
267
268  G_debug(3, "Found point in %d iterations", 10 - ret);
269
270  return 0;
271 }
272
273 /*
274  \brief Provide a breadth first binary division of real space along line
275  segment.
276
277  Looking for a point w/in the polygon.
278
279  This routine walks along the list of points on line segment
280  and divides each pair in half. It sticks that new point right into
281  the list, and then checks to see if it is inside the poly.
282
283  After going through the whole list, it calls itself. The list
284  now has a whole extra set of points to divide again.
285
287  \param Points
288  \param Token
289  \param X,Y
290  \param levels
291
292  \return # levels it took
293  \return -1 if exceeded # of levels
294  */
296  const struct line_pnts *Points,
298  double *Y, int levels)
299 {
300  struct Slink *A, *B, *C;
301
302  G_debug(3, "Vect__divide_and_conquer(): LEVEL %d", levels);
305
306  do {
308  A->next = C;
309  C->next = B;
310
311  C->x = (A->x + B->x) / 2.;
312
313  if (Vect_point_in_poly(C->x, *Y, Points) == 1) {
314  *X = C->x;
315  return levels;
316  }
317
318  A = B;
319  B = B->next;
320  } while (B != NULL);
321
322  /*
323  ** If it got through the entire loop and still no hits,
324  ** then lets go a level deeper and divide again.
325  */
326
327  if (levels <= 0)
328  return -1;
329
330  return Vect__divide_and_conquer(Head, Points, Token, X, Y, --levels);
331 }
332
334 {
336
338
339  while (p != NULL) {
340  tmp = p->next;
342  p = tmp;
343  }
344 }
345
346 /*!
347  \brief Get centroid of polygon
348
349  \param points polygon
350  \param[out] cent_x,cent_y centroid coordinates
351
352  \return 0 on success
353  \return -1 on error
354  */
355 int Vect_find_poly_centroid(const struct line_pnts *points, double *cent_x,
356  double *cent_y)
357 {
358  int i;
359  double *xptr1, *yptr1;
360  double *xptr2, *yptr2;
361  double cent_weight_x, cent_weight_y;
362  double len, tot_len;
363
364  tot_len = 0.0;
365  cent_weight_x = 0.0;
366  cent_weight_y = 0.0;
367
368  xptr1 = points->x;
369  yptr1 = points->y;
370  xptr2 = points->x + 1;
371  yptr2 = points->y + 1;
372
373  /* center of gravity of the polygon line, not area */
374  for (i = 1; i < points->n_points; i++) {
375  len = hypot(*xptr1 - *xptr2, *yptr1 - *yptr2);
376  cent_weight_x += len * ((*xptr1 + *xptr2) / 2.);
377  cent_weight_y += len * ((*yptr1 + *yptr2) / 2.);
378  tot_len += len;
379  xptr1++;
380  xptr2++;
381  yptr1++;
382  yptr2++;
383  }
384
385  if (tot_len == 0.0)
386  return -1;
387
388  *cent_x = cent_weight_x / tot_len;
389  *cent_y = cent_weight_y / tot_len;
390
391  return 0;
392 }
393
394 /*
395  ** returns true if point is in any of islands /w in area
396  ** returns 0 if not
397  ** returns -1 on error
398  */
399 /*
400  int
401  Vect_point_in_islands (
402  struct Map_info *Map,
403  int area,
404  double cent_x, double cent_y)
405  {
406  struct P_area *Area;
407  static struct line_pnts *TPoints;
408  static int first_time = 1;
409  int isle;
410
411  if (first_time == 1)
412  {
413  TPoints = Vect_new_line_struct ();
414  first_time = 0;
415  }
416
417  Area = &(Map->plus.Area[area]);
418
419  for (isle = 0; isle < Area->n_isles; isle++)
420  {
421  if (0 > Vect_get_isle_points (Map, Area->isles[isle], TPoints))
422  return -1;
423
424  if ( Vect_point_in_poly (cent_x, cent_y, TPoints) == 1 )
425  return 1;
426  }
427
428  return 0;
429  }
430  */
431
432 /*!
433  \brief Get point inside polygon but outside the islands specifiled in
434  IPoints.
435
436  Take a line and intersect it with the polygon and any islands.
437  sort the list of X values from these intersections. This will be a
438  list of segments alternating IN/OUT/IN/OUT of the polygon. Pick the
439  largest IN segment and take the midpoint.
440
441  \param Points polygon (boundary)
442  \param IPoints isles (list of isle boundaries)
443  \param n_isles number of isles
444  \param[out] att_x,att_y point coordinates
445
446  \return 0 on success
447  \return -1 on error
448  */
449 int Vect_get_point_in_poly_isl(const struct line_pnts *Points,
450  const struct line_pnts **IPoints, int n_isles,
451  double *att_x, double *att_y)
452 {
453  static struct line_pnts *Intersects;
454  static int first_time = 1;
455  double cent_x, cent_y;
456  register int i, j;
457  double max, hi_x, lo_x, hi_y, lo_y;
458  double fa, fb, dmax;
459  int exp;
460  int maxpos;
461  int point_in_sles = 0;
462  double diff;
463  int ret;
464
465  G_debug(3, "Vect_get_point_in_poly_isl(): n_isles = %d", n_isles);
466
467  if (first_time) {
468  Intersects = Vect_new_line_struct();
469  first_time = 0;
470  }
471
472  if (Points->n_points < 3) { /* test */
473  if (Points->n_points > 0) {
474  *att_x = Points->x[0];
475  *att_y = Points->y[0];
476  return 0;
477  }
478  return -1;
479  }
480
481  /* get centroid */
482  Vect_find_poly_centroid(Points, &cent_x, &cent_y);
483  /* is it w/in poly? */
484  if (Vect_point_in_poly(cent_x, cent_y, Points) == 1) {
485  /* if the point is inside the polygon */
486  for (i = 0; i < n_isles; i++) {
487  if (Vect_point_in_poly(cent_x, cent_y, IPoints[i]) >= 1) {
488  point_in_sles = 1;
489  break;
490  }
491  }
492  if (!point_in_sles) {
493  *att_x = cent_x;
494  *att_y = cent_y;
495  return 0;
496  }
497  }
498  /* guess we have to do it the hard way... */
499
500  /* first find att_y close to cent_y so that no points lie on the line */
501  /* find the point closest to line from below, and point close to line
502  from above and take average of their y-coordinates */
503  /* same for x */
504
505  /* first initializing lo_x,hi_x and lo_y,hi_y
506  * to be any 2 pnts on either side of cent_x and cent_y */
507  hi_y = cent_y - 1;
508  lo_y = cent_y + 1;
509
510  hi_x = cent_x - 1;
511  lo_x = cent_x + 1;
512  for (i = 0; i < Points->n_points; i++) {
513  if ((lo_y < cent_y) && (hi_y >= cent_y) && (lo_x < cent_x) &&
514  (hi_x >= cent_x))
515  break; /* already initialized */
516  if (Points->y[i] < cent_y)
517  lo_y = Points->y[i];
518  if (Points->y[i] >= cent_y)
519  hi_y = Points->y[i];
520
521  if (Points->x[i] < cent_x)
522  lo_x = Points->x[i];
523  if (Points->x[i] >= cent_x)
524  hi_x = Points->x[i];
525  }
526  /* first going through boundary points */
527  for (i = 0; i < Points->n_points; i++) {
528  if ((Points->y[i] < cent_y) &&
529  ((cent_y - Points->y[i]) < (cent_y - lo_y)))
530  lo_y = Points->y[i];
531  if ((Points->y[i] >= cent_y) &&
532  ((Points->y[i] - cent_y) < (hi_y - cent_y)))
533  hi_y = Points->y[i];
534
535  if ((Points->x[i] < cent_x) &&
536  ((cent_x - Points->x[i]) < (cent_x - lo_x)))
537  lo_x = Points->x[i];
538  if ((Points->x[i] >= cent_x) &&
539  ((Points->x[i] - cent_x) < (hi_x - cent_x)))
540  hi_x = Points->x[i];
541  }
542  for (i = 0; i < n_isles; i++) {
543  for (j = 0; j < IPoints[i]->n_points; j++) {
544  if ((IPoints[i]->y[j] < cent_y) &&
545  ((cent_y - IPoints[i]->y[j]) < (cent_y - lo_y)))
546  lo_y = IPoints[i]->y[j];
547  if ((IPoints[i]->y[j] >= cent_y) &&
548  ((IPoints[i]->y[j] - cent_y) < (hi_y - cent_y)))
549  hi_y = IPoints[i]->y[j];
550
551  if ((IPoints[i]->x[j] < cent_x) &&
552  ((cent_x - IPoints[i]->x[j]) < (cent_x - lo_x)))
553  lo_x = IPoints[i]->x[j];
554  if ((IPoints[i]->x[j] >= cent_x) &&
555  ((IPoints[i]->x[j] - cent_x) < (hi_x - cent_x)))
556  hi_x = IPoints[i]->x[j];
557  }
558  }
559
560  if (lo_y == hi_y)
561  return (-1); /* area is empty */
562
563  *att_y = (hi_y + lo_y) / 2.0;
564
565  Intersects->n_points = 0;
566  Vect__intersect_y_line_with_poly(Points, *att_y, Intersects);
567
568  /* add in intersections w/ holes */
569  for (i = 0; i < n_isles; i++) {
570  if (0 >
571  Vect__intersect_y_line_with_poly(IPoints[i], *att_y, Intersects))
572  return -1;
573  }
574
575  if (Intersects->n_points < 2) /* test */
576  return -1;
577
578  qsort(Intersects->x, (size_t)Intersects->n_points, sizeof(double),
579  comp_double);
580
581  max = 0;
582  maxpos = 0;
583
584  /* find area of MAX distance */
585  for (i = 0; i < Intersects->n_points; i += 2) {
586  diff = Intersects->x[i + 1] - Intersects->x[i];
587
588  if (diff > max) {
589  max = diff;
590  maxpos = i;
591  }
592  }
593  /* ULP single precision 23, double 52 bits, here 42 */
594  /* if the difference is too small, the point will be on a line
595  * ULP double is too small, ULP single too large */
596  fa = fabs(Intersects->x[maxpos]);
597  fb = fabs(Intersects->x[maxpos + 1]);
598  if (fa > fb)
599  dmax = frexp(fa, &exp);
600  else
601  dmax = frexp(fb, &exp);
602  exp -= 42;
603  dmax = ldexp(dmax, exp);
604
605  if (max > dmax) {
606  *att_x = (Intersects->x[maxpos] + Intersects->x[maxpos + 1]) / 2.;
607  }
608  else {
609  /* try x intersect */
610  G_debug(3, "Vect_get_point_in_poly_isl(): trying x intersect");
611
612  if (lo_x == hi_x)
613  return (-1); /* area is empty */
614
615  *att_x = (hi_x + lo_x) / 2.0;
616
617  Intersects->n_points = 0;
618  Vect__intersect_x_line_with_poly(Points, *att_x, Intersects);
619
620  /* add in intersections w/ holes */
621  for (i = 0; i < n_isles; i++) {
622  if (0 > Vect__intersect_x_line_with_poly(IPoints[i], *att_x,
623  Intersects))
624  return -1;
625  }
626
627  if (Intersects->n_points < 2) /* test */
628  return -1;
629
630  qsort(Intersects->y, (size_t)Intersects->n_points, sizeof(double),
631  comp_double);
632
633  max = 0;
634  maxpos = 0;
635
636  /* find area of MAX distance */
637  for (i = 0; i < Intersects->n_points; i += 2) {
638  diff = Intersects->y[i + 1] - Intersects->y[i];
639
640  if (diff > max) {
641  max = diff;
642  maxpos = i;
643  }
644  }
645  /* ULP single precision 23, double 52 bits, here 42 */
646  fa = fabs(Intersects->y[maxpos]);
647  fb = fabs(Intersects->y[maxpos + 1]);
648  if (fa > fb)
649  dmax = frexp(fa, &exp);
650  else
651  dmax = frexp(fb, &exp);
652  exp -= 42;
653  dmax = ldexp(dmax, exp);
654  if (max > dmax) {
655  *att_y = (Intersects->y[maxpos] + Intersects->y[maxpos + 1]) / 2.;
656  }
657  else {
658  /* area was (nearly) empty: example ((x1,y1), (x2,y2), (x1,y1)) */
659  G_warning("Vect_get_point_in_poly_isl(): collapsed area");
660  return -1;
661  }
662  }
663
664  /* is it now w/in poly? */
665  cent_x = *att_x;
666  cent_y = *att_y;
667  point_in_sles = 0;
668
669  ret = Vect_point_in_poly(cent_x, cent_y, Points);
670  if (ret == 2) {
671  /* point on outer ring, should not happen because of ULP test above */
672  G_warning("Vect_get_point_in_poly_isl(), the hard way: centroid is on "
673  "outer ring, max dist is %g",
674  max);
675  return -1;
676  }
677  if (ret == 1) {
678  /* if the point is inside the polygon, should not happen because of ULP
679  * test above */
680  for (i = 0; i < n_isles; i++) {
681  if (Vect_point_in_poly(cent_x, cent_y, IPoints[i]) >= 1) {
682  point_in_sles = 1;
683  G_warning("Vect_get_point_in_poly_isl(), the hard way: "
684  "centroid is in isle, max dist is %g",
685  max);
686  break;
687  }
688  }
689  if (!point_in_sles) {
690  return 0;
691  }
692  }
693
694  return -1;
695 }
696
697 /* Intersect segments of Points with ray from point X,Y to the right.
698  * Returns: -1 point exactly on segment
699  * number of intersections
700  */
701 static int segments_x_ray(double X, double Y, const struct line_pnts *Points)
702 {
703  double x1, x2, y1, y2;
704  double x_inter;
705  int n_intersects;
706  int n;
707
708  G_debug(3, "segments_x_ray(): x = %f y = %f n_points = %d", X, Y,
709  Points->n_points);
710
711  /* Follow the ray from X,Y along positive x and find number of
712  * intersections. Coordinates exactly on ray are considered to be slightly
713  * above. */
714
715  n_intersects = 0;
716  for (n = 1; n < Points->n_points; n++) {
717  x1 = Points->x[n - 1];
718  y1 = Points->y[n - 1];
719  x2 = Points->x[n];
720  y2 = Points->y[n];
721
722  /* G_debug() is slow, avoid it in loops over points,
723  * activate when needed */
724  /*
725  G_debug(3, "X = %f Y = %f x1 = %f y1 = %f x2 = %f y2 = %f", X, Y, x1,
726  y1, x2, y2);
727  */
728
729  /* I know, it should be possible to do that with less conditions,
730  * but it should be enough readable also! */
731
732  /* first, skip segments that obviously do not intersect with test ray */
733
734  /* segment above (X is not important) */
735  if (y1 > Y && y2 > Y)
736  continue;
737
738  /* segment below (X is not important) */
739  if (y1 < Y && y2 < Y)
740  continue;
741
742  /* segment left from X -> no intersection */
743  if (x1 < X && x2 < X)
744  continue;
745
746  /* point on vertex */
747  if ((x1 == X && y1 == Y) || (x2 == X && y2 == Y))
748  return -1;
749
750  /* on vertical boundary */
751  if (x1 == x2 && x1 == X) {
752  if ((y1 <= Y && y2 >= Y) || (y1 >= Y && y2 <= Y))
753  return -1;
754  }
755
756  /* on horizontal boundary */
757  if (y1 == y2 && y1 == Y) {
758  if ((x1 <= X && x2 >= X) || (x1 >= X && x2 <= X))
759  return -1;
760  else
761  continue; /* segment on ray (X is not important) */
762  }
763
764  /* segment on ray (X is not important) */
765  /* if (y1 == Y && y2 == Y)
766  continue; */
767
768  /* one end on Y second above (X is not important) */
769  if ((y1 == Y && y2 > Y) || (y2 == Y && y1 > Y))
770  continue;
771
772  /* For following cases we know that at least one of x1 and x2 is >= X */
773
774  /* one end of segment on Y second below Y */
775  if (y1 == Y && y2 < Y) {
776  if (x1 >= X) /* x of the end on the ray is >= X */
777  n_intersects++;
778  continue;
779  }
780  if (y2 == Y && y1 < Y) {
781  if (x2 >= X)
782  n_intersects++;
783  continue;
784  }
785
786  /* one end of segment above Y second below Y */
787  if ((y1 < Y && y2 > Y) || (y1 > Y && y2 < Y)) {
788  if (x1 >= X && x2 >= X) {
789  n_intersects++;
790  continue;
791  }
792
793  /* now either x1 < X && x2 > X or x1 > X && x2 < X -> calculate
794  * intersection */
795  x_inter = dig_x_intersect(x1, x2, y1, y2, Y);
796  G_debug(3, "x_inter = %f", x_inter);
797  if (x_inter == X)
798  return -1; /* point on segment, do not assume inside/outside */
799  else if (x_inter > X)
800  n_intersects++;
801
802  continue; /* would not be necessary, just to check, see below */
803  }
804  /* should not be reached (one condition is not necessary, but it is
805  * maybe better readable and it is a check) */
806  G_warning("segments_x_ray() %s: X = %f Y = %f x1 = %f y1 = %f x2 = %f "
807  "y2 = %f",
808  _("conditions failed"), X, Y, x1, y1, x2, y2);
809  }
810
811  return n_intersects;
812 }
813
814 /*!
815  \brief Determines if a point (X,Y) is inside a polygon.
816
817  \param X,Y point coordinates
818  \param Points polygon
819
820  \return 0 - outside
821  \return 1 - inside
822  \return 2 - on the boundary
823  */
824 int Vect_point_in_poly(double X, double Y, const struct line_pnts *Points)
825 {
826  int n_intersects;
827
828  G_debug(3, "Vect_point_in_poly(): x = %f y = %f n_points = %d", X, Y,
829  Points->n_points);
830
831  n_intersects = segments_x_ray(X, Y, Points);
832
833  if (n_intersects == -1)
834  return 2;
835
836  /* odd number of intersections: inside, return 1
837  * even number of intersections: outside, return 0 */
838  return (n_intersects & 1);
839 }
840
841 /*!
842  \brief Determines if a point (X,Y) is inside an area outer ring. Islands are
843  not considered.
844
845  \param X,Y point coordinates
846  \param Map vector map
847  \param area area id
848  \param box area bounding box
849
850  \return 0 - outside
851  \return 1 - inside
852  \return 2 - on the boundary
853  */
854 int Vect_point_in_area_outer_ring(double X, double Y, struct Map_info *Map,
855  int area, struct bound_box *box)
856 {
857  static int first = 1;
858  int n_intersects, inter;
859  int i, line;
860  static struct line_pnts *Points;
862  struct P_area *Area;
863
864  /* keep in sync with Vect_point_in_island() */
865
866  G_debug(3, "Vect_point_in_area_outer_ring(): x = %f y = %f area = %d", X, Y,
867  area);
868
869  if (first == 1) {
870  Points = Vect_new_line_struct();
871  first = 0;
872  }
873
874  Plus = &(Map->plus);
875  Area = Plus->Area[area];
876
877  /* First it must be in box */
878  if (X < box->W || X > box->E || Y > box->N || Y < box->S)
879  return 0;
880
881  n_intersects = 0;
882  for (i = 0; i < Area->n_lines; i++) {
883  line = abs(Area->lines[i]);
884
886
887  /* if the bbox of the line would be available,
888  * the bbox could be used for a first check: */
889
890  /* Vect_line_box(Points, &lbox);
891  * do not check lines that obviously do not intersect with test ray */
892  /* if ((lbox.N < Y) || (lbox.S > Y) || (lbox.E < X))
893  continue;
894  */
895
896  /* retrieving the bbox from the spatial index or
897  * calculating the box from the vertices is slower than
898  * just feeding the line to segments_x_ray() */
899
900  inter = segments_x_ray(X, Y, Points);
901  if (inter == -1)
902  return 2;
903  n_intersects += inter;
904  }
905
906  /* odd number of intersections: inside, return 1
907  * even number of intersections: outside, return 0 */
908  return (n_intersects & 1);
909 }
910
911 /*!
912  \brief Determines if a point (X,Y) is inside an island.
913
914  \param X,Y point coordinates
915  \param Map vector map
916  \param isle isle id
917  \param box isle bounding box
918
919  \return 0 - outside
920  \return 1 - inside
921  \return 2 - on the boundary
922  */
923 int Vect_point_in_island(double X, double Y, struct Map_info *Map, int isle,
924  struct bound_box *box)
925 {
926  static int first = 1;
927  int n_intersects, inter;
928  int i, line;
929  static struct line_pnts *Points;
931  struct P_isle *Isle;
932
933  /* keep in sync with Vect_point_in_area_outer_ring() */
934
935  G_debug(3, "Vect_point_in_island(): x = %f y = %f isle = %d", X, Y, isle);
936
937  if (first == 1) {
938  Points = Vect_new_line_struct();
939  first = 0;
940  }
941
942  Plus = &(Map->plus);
943  Isle = Plus->Isle[isle];
944
945  /* First it must be in box */
946  if (X < box->W || X > box->E || Y > box->N || Y < box->S)
947  return 0;
948
949  n_intersects = 0;
950  for (i = 0; i < Isle->n_lines; i++) {
951  line = abs(Isle->lines[i]);
952
954
955  /* if the bbox of the line would be available,
956  * the bbox could be used for a first check: */
957
958  /* Vect_line_box(Points, &lbox);
959  * don't check lines that obviously do not intersect with test ray */
960  /* if ((lbox.N < Y) || (lbox.S > Y) || (lbox.E < X))
961  continue;
962  */
963
964  /* retrieving the bbox from the spatial index or
965  * calculating the box from the vertices is slower than
966  * just feeding the line to segments_x_ray() */
967
968  inter = segments_x_ray(X, Y, Points);
969  if (inter == -1)
970  return 2;
971  n_intersects += inter;
972  }
973
974  /* odd number of intersections: inside, return 1
975  * even number of intersections: outside, return 0 */
976  return (n_intersects & 1);
977 }
int Vect_get_point_in_poly_isl(const struct line_pnts *Points, const struct line_pnts **IPoints, int n_isles, double *att_x, double *att_y)
Get point inside polygon but outside the islands specifiled in IPoints.
Definition: Vlib/poly.c:449
int Vect_point_in_poly(double X, double Y, const struct line_pnts *Points)
Determines if a point (X,Y) is inside a polygon.
Definition: Vlib/poly.c:824
int Vect_get_point_in_poly(const struct line_pnts *Points, double *X, double *Y)
Get point inside polygon.
Definition: Vlib/poly.c:208
int Vect_find_poly_centroid(const struct line_pnts *points, double *cent_x, double *cent_y)
Get centroid of polygon.
Definition: Vlib/poly.c:355
int Vect__intersect_x_line_with_poly(const struct line_pnts *, double, struct line_pnts *)
Definition: Vlib/poly.c:169
int Vect__intersect_y_line_with_poly(const struct line_pnts *, double, struct line_pnts *)
Definition: Vlib/poly.c:126
int Vect_point_in_island(double X, double Y, struct Map_info *Map, int isle, struct bound_box *box)
Determines if a point (X,Y) is inside an island.
Definition: Vlib/poly.c:923
int Vect_point_in_area_outer_ring(double X, double Y, struct Map_info *Map, int area, struct bound_box *box)
Determines if a point (X,Y) is inside an area outer ring. Islands are not considered.
Definition: Vlib/poly.c:854
int Vect_get_point_in_area(struct Map_info *Map, int area, double *X, double *Y)
Get point inside area and outside all islands.
Definition: Vlib/poly.c:55
#define NULL
Definition: ccmath.h:32
#define G_realloc(p, n)
Definition: defs/gis.h:96
void G_warning(const char *,...) __attribute__((format(printf
int G_debug(int, const char *,...) __attribute__((format(printf
Definition: dispose.c:10
Definition: new.c:12
int Vect_get_isle_points(struct Map_info *, int, struct line_pnts *)
Returns polygon array of points for given isle.
int Vect_get_area_points(struct Map_info *, int, struct line_pnts *)
Returns polygon array of points (outer ring) of given area.
int Vect_get_area_isle(struct Map_info *, int, int)
Returns isle id for area.
int Vect_read_line(struct Map_info *, struct line_pnts *, struct line_cats *, int)
Read vector feature (topological level required)
int Vect_get_area_num_isles(struct Map_info *, int)
Returns number of isles for given area.
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
int Vect_append_point(struct line_pnts *, double, double, double)
Appends one point to the end of a line.
Definition: line.c:148
double dig_x_intersect(double, double, double, double, double)
Definition: inside.c:20
#define max(x, y)
Definition: draw2.c:30
#define _(str)
Definition: glocale.h:10
#define VOID_T
#define X
Definition: ogsf.h:140
#define W
Definition: ogsf.h:143
#define Y
Definition: ogsf.h:141
double b
Definition: r_raster.c:39
Vector map info.
Definition: dig_structs.h:1243
Plus info (topology, version, ...)
Definition: dig_structs.h:1270
Area (topology) info.
Definition: dig_structs.h:1583
plus_t n_lines
Number of boundary lines.
Definition: dig_structs.h:1587
plus_t * lines
List of boundary lines.
Definition: dig_structs.h:1598
Isle (topology) info.
Definition: dig_structs.h:1623
plus_t * lines
List of boundary lines.
Definition: dig_structs.h:1638
plus_t n_lines
Number of boundary lines.
Definition: dig_structs.h:1627
Basic topology-related info.
Definition: dig_structs.h:769
struct P_area ** Area
Array of areas.
Definition: dig_structs.h:875
struct P_isle ** Isle
Array of isles.
Definition: dig_structs.h:879
Bounding box.
Definition: dig_structs.h:64
double N
North.
Definition: dig_structs.h:68
double E
East.
Definition: dig_structs.h:76
Feature geometry info - coordinates.
Definition: dig_structs.h:1651
double * y
Array of Y coordinates.
Definition: dig_structs.h:1659
int alloc_points
Allocated space for points.
Definition: dig_structs.h:1671
double * x
Array of X coordinates.
Definition: dig_structs.h:1655
int n_points
Number of points.
Definition: dig_structs.h:1667
#define x