GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
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>
20 #include <grass/linkm.h>
21 #include <grass/glocale.h>
22 
23 struct Slink
24 {
25  struct Slink *next;
26  double x;
27 };
28 
29 
30 /* function prototypes */
31 static int comp_double(double *, double *);
32 static int V__within(double, double, double);
35 static void destroy_links(struct link_head *, struct Slink *);
36 static int Vect__divide_and_conquer(struct Slink *, const struct line_pnts *,
37  struct link_head *, double *, double *,
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;
214  struct Slink *Head;
215  static struct link_head *Token;
216  struct Slink *tmp;
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 */
248  Token = (struct link_head *)link_init(sizeof(struct Slink));
249  first_time = 0;
250  }
251 
252  Head = (struct Slink *)link_new(Token);
253  tmp = (struct Slink *)link_new(Token);
254 
255  Head->next = tmp;
256  tmp->next = NULL;
257 
258  Head->x = x_min;
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 
264  destroy_links(Token, Head);
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 
290  \param Head
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
300 Vect__divide_and_conquer(struct Slink *Head,
301  const struct line_pnts *Points,
302  struct link_head *Token,
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);
308  A = Head;
309  B = Head->next;
310 
311  do {
312  C = (struct Slink *)link_new(Token);
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 
339 static void destroy_links(struct link_head *Token, struct Slink *Head)
340 {
341  struct Slink *p, *tmp;
342 
343  p = Head;
344 
345  while (p != NULL) {
346  tmp = p->next;
347  link_dispose(Token, (VOID_T *) p);
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;
863  const struct Plus_head *Plus;
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 
887  Vect_read_line(Map, Points, NULL, line);
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;
932  const struct Plus_head *Plus;
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 
956  Vect_read_line(Map, Points, NULL, line);
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
Definition: linkm.h:8
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
void link_dispose(struct link_head *, VOID_T *)
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
VOID_T * link_new(struct link_head *)
Definition: new.c:12
plus_t * lines
List of boundary lines.
Definition: dig_structs.h:1621
void link_exit_on_error(int)
Definition: linkm/init.c:35
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
struct Plus_head plus
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
struct link_head * link_init(int)
Definition: linkm/init.c:40
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