GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-535c39c9fc
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  struct Slink *next;
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 *);
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 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;
211  struct Slink *Head;
212  static struct link_head *Token;
213  struct Slink *tmp;
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 */
244  Token = (struct link_head *)link_init(sizeof(struct Slink));
245  first_time = 0;
246  }
247 
248  Head = (struct Slink *)link_new(Token);
249  tmp = (struct Slink *)link_new(Token);
250 
251  Head->next = tmp;
252  tmp->next = NULL;
253 
254  Head->x = x_min;
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 
260  destroy_links(Token, Head);
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 
286  \param Head
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  */
295 static int Vect__divide_and_conquer(struct Slink *Head,
296  const struct line_pnts *Points,
297  struct link_head *Token, double *X,
298  double *Y, int levels)
299 {
300  struct Slink *A, *B, *C;
301 
302  G_debug(3, "Vect__divide_and_conquer(): LEVEL %d", levels);
303  A = Head;
304  B = Head->next;
305 
306  do {
307  C = (struct Slink *)link_new(Token);
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 
333 static void destroy_links(struct link_head *Token, struct Slink *Head)
334 {
335  struct Slink *p, *tmp;
336 
337  p = Head;
338 
339  while (p != NULL) {
340  tmp = p->next;
341  link_dispose(Token, (VOID_T *)p);
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;
861  const struct Plus_head *Plus;
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 
885  Vect_read_line(Map, Points, NULL, line);
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;
930  const struct Plus_head *Plus;
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 
953  Vect_read_line(Map, Points, NULL, line);
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
struct link_head * link_init(int)
Definition: linkm/init.c:42
void link_dispose(struct link_head *, VOID_T *)
Definition: dispose.c:10
void link_exit_on_error(int)
Definition: linkm/init.c:37
VOID_T * link_new(struct link_head *)
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
Definition: linkm.h:8
#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
struct Plus_head plus
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