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