PostGIS  2.2.8dev-r@@SVN_REVISION@@
lwgeom_functions_analytic.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  *
6  * Copyright (C) 2001-2005 Refractions Research Inc.
7  *
8  * This is free software; you can redistribute and/or modify it under
9  * the terms of the GNU General Public Licence. See the COPYING file.
10  *
11  **********************************************************************/
12 
13 #include "postgres.h"
14 #include "fmgr.h"
15 #include "liblwgeom.h"
16 #include "liblwgeom_internal.h" /* For FP comparators. */
17 #include "lwgeom_pg.h"
18 #include "math.h"
19 #include "lwgeom_rtree.h"
21 
22 
23 /***********************************************************************
24  * Simple Douglas-Peucker line simplification.
25  * No checks are done to avoid introduction of self-intersections.
26  * No topology relations are considered.
27  *
28  * --strk@keybit.net;
29  ***********************************************************************/
30 
31 
32 /* Prototypes */
33 Datum LWGEOM_simplify2d(PG_FUNCTION_ARGS);
34 Datum LWGEOM_SetEffectiveArea(PG_FUNCTION_ARGS);
35 Datum ST_LineCrossingDirection(PG_FUNCTION_ARGS);
36 
37 
38 double determineSide(POINT2D *seg1, POINT2D *seg2, POINT2D *point);
39 int isOnSegment(POINT2D *seg1, POINT2D *seg2, POINT2D *point);
40 int point_in_ring(POINTARRAY *pts, POINT2D *point);
41 int point_in_ring_rtree(RTREE_NODE *root, POINT2D *point);
42 
43 
45 Datum LWGEOM_simplify2d(PG_FUNCTION_ARGS)
46 {
47  GSERIALIZED *geom = PG_GETARG_GSERIALIZED_P(0);
48  double dist = PG_GETARG_FLOAT8(1);
49  GSERIALIZED *result;
50  int type = gserialized_get_type(geom);
51  LWGEOM *in, *out;
52  bool preserve_collapsed = false;
53 
54  /* Handle optional argument to preserve collapsed features */
55  if ( PG_NARGS() > 2 && ! PG_ARGISNULL(2) )
56  preserve_collapsed = true;
57 
58  /* Can't simplify points! */
59  if ( type == POINTTYPE || type == MULTIPOINTTYPE )
60  PG_RETURN_POINTER(geom);
61 
62  in = lwgeom_from_gserialized(geom);
63 
64  out = lwgeom_simplify(in, dist, preserve_collapsed);
65  if ( ! out ) PG_RETURN_NULL();
66 
67  /* COMPUTE_BBOX TAINTING */
68  if ( in->bbox ) lwgeom_add_bbox(out);
69 
70  result = geometry_serialize(out);
71  lwgeom_free(out);
72  PG_FREE_IF_COPY(geom, 0);
73  PG_RETURN_POINTER(result);
74 }
75 
77 Datum LWGEOM_SetEffectiveArea(PG_FUNCTION_ARGS)
78 {
79  GSERIALIZED *geom = PG_GETARG_GSERIALIZED_P(0);
80  GSERIALIZED *result;
81  int type = gserialized_get_type(geom);
82  LWGEOM *in;
83  LWGEOM *out;
84  double area=0;
85  int set_area=0;
86 
87  if ( type == POINTTYPE || type == MULTIPOINTTYPE )
88  PG_RETURN_POINTER(geom);
89 
90  if ( (PG_NARGS()>1) && (!PG_ARGISNULL(1)) )
91  area = PG_GETARG_FLOAT8(1);
92 
93  if ( (PG_NARGS()>2) && (!PG_ARGISNULL(2)) )
94  set_area = PG_GETARG_INT32(2);
95 
96  in = lwgeom_from_gserialized(geom);
97 
98  out = lwgeom_set_effective_area(in,set_area, area);
99  if ( ! out ) PG_RETURN_NULL();
100 
101  /* COMPUTE_BBOX TAINTING */
102  if ( in->bbox ) lwgeom_add_bbox(out);
103 
104  result = geometry_serialize(out);
105  lwgeom_free(out);
106  PG_FREE_IF_COPY(geom, 0);
107  PG_RETURN_POINTER(result);
108 }
109 
110 
111 /***********************************************************************
112  * --strk@keybit.net;
113  ***********************************************************************/
114 
115 /***********************************************************************
116  * Interpolate a point along a line, useful for Geocoding applications
117  * SELECT line_interpolate_point( 'LINESTRING( 0 0, 2 2'::geometry, .5 )
118  * returns POINT( 1 1 ).
119  * Works in 2d space only.
120  *
121  * Initially written by: jsunday@rochgrp.com
122  * Ported to LWGEOM by: strk@refractions.net
123  ***********************************************************************/
124 
125 Datum LWGEOM_line_interpolate_point(PG_FUNCTION_ARGS);
126 
128 Datum LWGEOM_line_interpolate_point(PG_FUNCTION_ARGS)
129 {
130  GSERIALIZED *gser = PG_GETARG_GSERIALIZED_P(0);
131  GSERIALIZED *result;
132  double distance = PG_GETARG_FLOAT8(1);
133  LWLINE *line;
134  LWGEOM *geom;
135  LWPOINT *point;
136  POINTARRAY *ipa, *opa;
137  POINT4D pt;
138  int nsegs, i;
139  double length, slength, tlength;
140 
141  if ( distance < 0 || distance > 1 )
142  {
143  elog(ERROR,"line_interpolate_point: 2nd arg isn't within [0,1]");
144  PG_RETURN_NULL();
145  }
146 
147  if ( gserialized_get_type(gser) != LINETYPE )
148  {
149  elog(ERROR,"line_interpolate_point: 1st arg isn't a line");
150  PG_RETURN_NULL();
151  }
152 
153  /* Empty.InterpolatePoint == Point Empty */
154  if ( gserialized_is_empty(gser) )
155  {
157  result = geometry_serialize(lwpoint_as_lwgeom(point));
158  lwpoint_free(point);
159  PG_RETURN_POINTER(result);
160  }
161 
162  geom = lwgeom_from_gserialized(gser);
163  line = lwgeom_as_lwline(geom);
164  ipa = line->points;
165 
166  /* If distance is one of the two extremes, return the point on that
167  * end rather than doing any expensive computations
168  */
169  if ( distance == 0.0 || distance == 1.0 )
170  {
171  if ( distance == 0.0 )
172  getPoint4d_p(ipa, 0, &pt);
173  else
174  getPoint4d_p(ipa, ipa->npoints-1, &pt);
175 
176  opa = ptarray_construct(lwgeom_has_z(geom), lwgeom_has_m(geom), 1);
177  ptarray_set_point4d(opa, 0, &pt);
178 
179  point = lwpoint_construct(line->srid, NULL, opa);
180  PG_RETURN_POINTER(geometry_serialize(lwpoint_as_lwgeom(point)));
181  }
182 
183  /* Interpolate a point on the line */
184  nsegs = ipa->npoints - 1;
185  length = ptarray_length_2d(ipa);
186  tlength = 0;
187  for ( i = 0; i < nsegs; i++ )
188  {
189  POINT4D p1, p2;
190  POINT4D *p1ptr=&p1, *p2ptr=&p2; /* don't break
191  * strict-aliasing rules
192  */
193 
194  getPoint4d_p(ipa, i, &p1);
195  getPoint4d_p(ipa, i+1, &p2);
196 
197  /* Find the relative length of this segment */
198  slength = distance2d_pt_pt((POINT2D*)p1ptr, (POINT2D*)p2ptr)/length;
199 
200  /* If our target distance is before the total length we've seen
201  * so far. create a new point some distance down the current
202  * segment.
203  */
204  if ( distance < tlength + slength )
205  {
206  double dseg = (distance - tlength) / slength;
207  interpolate_point4d(&p1, &p2, &pt, dseg);
208  opa = ptarray_construct(lwgeom_has_z(geom), lwgeom_has_m(geom), 1);
209  ptarray_set_point4d(opa, 0, &pt);
210  point = lwpoint_construct(line->srid, NULL, opa);
211  PG_RETURN_POINTER(geometry_serialize(lwpoint_as_lwgeom(point)));
212  }
213  tlength += slength;
214  }
215 
216  /* Return the last point on the line. This shouldn't happen, but
217  * could if there's some floating point rounding errors. */
218  getPoint4d_p(ipa, ipa->npoints-1, &pt);
219  opa = ptarray_construct(lwgeom_has_z(geom), lwgeom_has_m(geom), 1);
220  ptarray_set_point4d(opa, 0, &pt);
221  point = lwpoint_construct(line->srid, NULL, opa);
222  PG_FREE_IF_COPY(gser, 0);
223  PG_RETURN_POINTER(geometry_serialize(lwpoint_as_lwgeom(point)));
224 }
225 /***********************************************************************
226  * --jsunday@rochgrp.com;
227  ***********************************************************************/
228 
229 /***********************************************************************
230  *
231  * Grid application function for postgis.
232  *
233  * WHAT IS
234  * -------
235  *
236  * This is a grid application function for postgis.
237  * You use it to stick all of a geometry points to
238  * a custom grid defined by its origin and cell size
239  * in geometry units.
240  *
241  * Points reduction is obtained by collapsing all
242  * consecutive points falling on the same grid node and
243  * removing all consecutive segments S1,S2 having
244  * S2.startpoint = S1.endpoint and S2.endpoint = S1.startpoint.
245  *
246  * ISSUES
247  * ------
248  *
249  * Only 2D is contemplated in grid application.
250  *
251  * Consecutive coincident segments removal does not work
252  * on first/last segments (They are not considered consecutive).
253  *
254  * Grid application occurs on a geometry subobject basis, thus no
255  * points reduction occurs for multipoint geometries.
256  *
257  * USAGE TIPS
258  * ----------
259  *
260  * Run like this:
261  *
262  * SELECT SnapToGrid(<geometry>, <originX>, <originY>, <sizeX>, <sizeY>);
263  *
264  * Grid cells are not visible on a screen as long as the screen
265  * point size is equal or bigger then the grid cell size.
266  * This assumption may be used to reduce the number of points in
267  * a map for a given display scale.
268  *
269  * Keeping multiple resolution versions of the same data may be used
270  * in conjunction with MINSCALE/MAXSCALE keywords of mapserv to speed
271  * up rendering.
272  *
273  * Check also the use of DP simplification to reduce grid visibility.
274  * I'm still researching about the relationship between grid size and
275  * DP epsilon values - please tell me if you know more about this.
276  *
277  *
278  * --strk@keybit.net;
279  *
280  ***********************************************************************/
281 
282 #define CHECK_RING_IS_CLOSE
283 #define SAMEPOINT(a,b) ((a)->x==(b)->x&&(a)->y==(b)->y)
284 
285 
286 
287 /* Forward declarations */
288 Datum LWGEOM_snaptogrid(PG_FUNCTION_ARGS);
289 Datum LWGEOM_snaptogrid_pointoff(PG_FUNCTION_ARGS);
290 
291 
292 
294 Datum LWGEOM_snaptogrid(PG_FUNCTION_ARGS)
295 {
296  LWGEOM *in_lwgeom;
297  GSERIALIZED *out_geom = NULL;
298  LWGEOM *out_lwgeom;
299  gridspec grid;
300 
301  GSERIALIZED *in_geom = PG_GETARG_GSERIALIZED_P(0);
302 
303  /* Set grid values to zero to start */
304  memset(&grid, 0, sizeof(gridspec));
305 
306  grid.ipx = PG_GETARG_FLOAT8(1);
307  grid.ipy = PG_GETARG_FLOAT8(2);
308  grid.xsize = PG_GETARG_FLOAT8(3);
309  grid.ysize = PG_GETARG_FLOAT8(4);
310 
311  /* Return input geometry if input geometry is empty */
312  if ( gserialized_is_empty(in_geom) )
313  {
314  PG_RETURN_POINTER(in_geom);
315  }
316 
317  /* Return input geometry if input grid is meaningless */
318  if ( grid.xsize==0 && grid.ysize==0 && grid.zsize==0 && grid.msize==0 )
319  {
320  PG_RETURN_POINTER(in_geom);
321  }
322 
323  in_lwgeom = lwgeom_from_gserialized(in_geom);
324 
325  POSTGIS_DEBUGF(3, "SnapToGrid got a %s", lwtype_name(in_lwgeom->type));
326 
327  out_lwgeom = lwgeom_grid(in_lwgeom, &grid);
328  if ( out_lwgeom == NULL ) PG_RETURN_NULL();
329 
330  /* COMPUTE_BBOX TAINTING */
331  if ( in_lwgeom->bbox )
332  lwgeom_add_bbox(out_lwgeom);
333 
334 
335  POSTGIS_DEBUGF(3, "SnapToGrid made a %s", lwtype_name(out_lwgeom->type));
336 
337  out_geom = geometry_serialize(out_lwgeom);
338 
339  PG_RETURN_POINTER(out_geom);
340 }
341 
342 
343 #if POSTGIS_DEBUG_LEVEL >= 4
344 /* Print grid using given reporter */
345 static void
346 grid_print(const gridspec *grid)
347 {
348  lwpgnotice("GRID(%g %g %g %g, %g %g %g %g)",
349  grid->ipx, grid->ipy, grid->ipz, grid->ipm,
350  grid->xsize, grid->ysize, grid->zsize, grid->msize);
351 }
352 #endif
353 
354 
356 Datum LWGEOM_snaptogrid_pointoff(PG_FUNCTION_ARGS)
357 {
358  GSERIALIZED *in_geom, *in_point;
359  LWGEOM *in_lwgeom;
360  LWPOINT *in_lwpoint;
361  GSERIALIZED *out_geom = NULL;
362  LWGEOM *out_lwgeom;
363  gridspec grid;
364  /* BOX3D box3d; */
365  POINT4D offsetpoint;
366 
367  in_geom = PG_GETARG_GSERIALIZED_P(0);
368 
369  /* Return input geometry if input geometry is empty */
370  if ( gserialized_is_empty(in_geom) )
371  {
372  PG_RETURN_POINTER(in_geom);
373  }
374 
375  in_point = PG_GETARG_GSERIALIZED_P(1);
376  in_lwpoint = lwgeom_as_lwpoint(lwgeom_from_gserialized(in_point));
377  if ( in_lwpoint == NULL )
378  {
379  lwpgerror("Offset geometry must be a point");
380  }
381 
382  grid.xsize = PG_GETARG_FLOAT8(2);
383  grid.ysize = PG_GETARG_FLOAT8(3);
384  grid.zsize = PG_GETARG_FLOAT8(4);
385  grid.msize = PG_GETARG_FLOAT8(5);
386 
387  /* Take offsets from point geometry */
388  getPoint4d_p(in_lwpoint->point, 0, &offsetpoint);
389  grid.ipx = offsetpoint.x;
390  grid.ipy = offsetpoint.y;
391  if (FLAGS_GET_Z(in_lwpoint->flags) ) grid.ipz = offsetpoint.z;
392  else grid.ipz=0;
393  if (FLAGS_GET_M(in_lwpoint->flags) ) grid.ipm = offsetpoint.m;
394  else grid.ipm=0;
395 
396 #if POSTGIS_DEBUG_LEVEL >= 4
397  grid_print(&grid);
398 #endif
399 
400  /* Return input geometry if input grid is meaningless */
401  if ( grid.xsize==0 && grid.ysize==0 && grid.zsize==0 && grid.msize==0 )
402  {
403  PG_RETURN_POINTER(in_geom);
404  }
405 
406  in_lwgeom = lwgeom_from_gserialized(in_geom);
407 
408  POSTGIS_DEBUGF(3, "SnapToGrid got a %s", lwtype_name(in_lwgeom->type));
409 
410  out_lwgeom = lwgeom_grid(in_lwgeom, &grid);
411  if ( out_lwgeom == NULL ) PG_RETURN_NULL();
412 
413  /* COMPUTE_BBOX TAINTING */
414  if ( in_lwgeom->bbox ) lwgeom_add_bbox(out_lwgeom);
415 
416  POSTGIS_DEBUGF(3, "SnapToGrid made a %s", lwtype_name(out_lwgeom->type));
417 
418  out_geom = geometry_serialize(out_lwgeom);
419 
420  PG_RETURN_POINTER(out_geom);
421 }
422 
423 
424 /*
425 ** crossingDirection(line1, line2)
426 **
427 ** Determines crossing direction of line2 relative to line1.
428 ** Only accepts LINESTRING as parameters!
429 */
431 Datum ST_LineCrossingDirection(PG_FUNCTION_ARGS)
432 {
433  int type1, type2, rv;
434  LWLINE *l1 = NULL;
435  LWLINE *l2 = NULL;
436  GSERIALIZED *geom1 = PG_GETARG_GSERIALIZED_P(0);
437  GSERIALIZED *geom2 = PG_GETARG_GSERIALIZED_P(1);
438 
440 
441  type1 = gserialized_get_type(geom1);
442  type2 = gserialized_get_type(geom2);
443 
444  if ( type1 != LINETYPE || type2 != LINETYPE )
445  {
446  elog(ERROR,"This function only accepts LINESTRING as arguments.");
447  PG_RETURN_NULL();
448  }
449 
452 
453  rv = lwline_crossing_direction(l1, l2);
454 
455  PG_FREE_IF_COPY(geom1, 0);
456  PG_FREE_IF_COPY(geom2, 1);
457 
458  PG_RETURN_INT32(rv);
459 
460 }
461 
462 
463 
464 /***********************************************************************
465  * --strk@keybit.net
466  ***********************************************************************/
467 
468 Datum LWGEOM_line_substring(PG_FUNCTION_ARGS);
469 
471 Datum LWGEOM_line_substring(PG_FUNCTION_ARGS)
472 {
473  GSERIALIZED *geom = PG_GETARG_GSERIALIZED_P(0);
474  double from = PG_GETARG_FLOAT8(1);
475  double to = PG_GETARG_FLOAT8(2);
476  LWGEOM *olwgeom;
477  POINTARRAY *ipa, *opa;
478  GSERIALIZED *ret;
479  int type = gserialized_get_type(geom);
480 
481  if ( from < 0 || from > 1 )
482  {
483  elog(ERROR,"line_interpolate_point: 2nd arg isn't within [0,1]");
484  PG_RETURN_NULL();
485  }
486 
487  if ( to < 0 || to > 1 )
488  {
489  elog(ERROR,"line_interpolate_point: 3rd arg isn't within [0,1]");
490  PG_RETURN_NULL();
491  }
492 
493  if ( from > to )
494  {
495  elog(ERROR, "2nd arg must be smaller then 3rd arg");
496  PG_RETURN_NULL();
497  }
498 
499  if ( type == LINETYPE )
500  {
502 
503  if ( lwgeom_is_empty((LWGEOM*)iline) )
504  {
505  /* TODO return empty line */
506  lwline_release(iline);
507  PG_FREE_IF_COPY(geom, 0);
508  PG_RETURN_NULL();
509  }
510 
511  ipa = iline->points;
512 
513  opa = ptarray_substring(ipa, from, to, 0);
514 
515  if ( opa->npoints == 1 ) /* Point returned */
516  olwgeom = (LWGEOM *)lwpoint_construct(iline->srid, NULL, opa);
517  else
518  olwgeom = (LWGEOM *)lwline_construct(iline->srid, NULL, opa);
519 
520  }
521  else if ( type == MULTILINETYPE )
522  {
523  LWMLINE *iline;
524  int i = 0, g = 0;
525  int homogeneous = LW_TRUE;
526  LWGEOM **geoms = NULL;
527  double length = 0.0, sublength = 0.0, minprop = 0.0, maxprop = 0.0;
528 
530 
531  if ( lwgeom_is_empty((LWGEOM*)iline) )
532  {
533  /* TODO return empty collection */
534  lwmline_release(iline);
535  PG_FREE_IF_COPY(geom, 0);
536  PG_RETURN_NULL();
537  }
538 
539  /* Calculate the total length of the mline */
540  for ( i = 0; i < iline->ngeoms; i++ )
541  {
542  LWLINE *subline = (LWLINE*)iline->geoms[i];
543  if ( subline->points && subline->points->npoints > 1 )
544  length += ptarray_length_2d(subline->points);
545  }
546 
547  geoms = lwalloc(sizeof(LWGEOM*) * iline->ngeoms);
548 
549  /* Slice each sub-geometry of the multiline */
550  for ( i = 0; i < iline->ngeoms; i++ )
551  {
552  LWLINE *subline = (LWLINE*)iline->geoms[i];
553  double subfrom = 0.0, subto = 0.0;
554 
555  if ( subline->points && subline->points->npoints > 1 )
556  sublength += ptarray_length_2d(subline->points);
557 
558  /* Calculate proportions for this subline */
559  minprop = maxprop;
560  maxprop = sublength / length;
561 
562  /* This subline doesn't reach the lowest proportion requested
563  or is beyond the highest proporton */
564  if ( from > maxprop || to < minprop )
565  continue;
566 
567  if ( from <= minprop )
568  subfrom = 0.0;
569  if ( to >= maxprop )
570  subto = 1.0;
571 
572  if ( from > minprop && from <= maxprop )
573  subfrom = (from - minprop) / (maxprop - minprop);
574 
575  if ( to < maxprop && to >= minprop )
576  subto = (to - minprop) / (maxprop - minprop);
577 
578 
579  opa = ptarray_substring(subline->points, subfrom, subto, 0);
580  if ( opa && opa->npoints > 0 )
581  {
582  if ( opa->npoints == 1 ) /* Point returned */
583  {
584  geoms[g] = (LWGEOM *)lwpoint_construct(SRID_UNKNOWN, NULL, opa);
585  homogeneous = LW_FALSE;
586  }
587  else
588  {
589  geoms[g] = (LWGEOM *)lwline_construct(SRID_UNKNOWN, NULL, opa);
590  }
591  g++;
592  }
593 
594 
595 
596  }
597  /* If we got any points, we need to return a GEOMETRYCOLLECTION */
598  if ( ! homogeneous )
599  type = COLLECTIONTYPE;
600 
601  olwgeom = (LWGEOM*)lwcollection_construct(type, iline->srid, NULL, g, geoms);
602  }
603  else
604  {
605  elog(ERROR,"line_substring: 1st arg isn't a line");
606  PG_RETURN_NULL();
607  }
608 
609  ret = geometry_serialize(olwgeom);
610  lwgeom_free(olwgeom);
611  PG_FREE_IF_COPY(geom, 0);
612  PG_RETURN_POINTER(ret);
613 
614 }
615 
616 /*******************************************************************************
617  * The following is based on the "Fast Winding Number Inclusion of a Point
618  * in a Polygon" algorithm by Dan Sunday.
619  * http://softsurfer.com/Archive/algorithm_0103/algorithm_0103.htm#Winding%20Number
620  ******************************************************************************/
621 
622 /*
623  * returns: >0 for a point to the left of the segment,
624  * <0 for a point to the right of the segment,
625  * 0 for a point on the segment
626  */
627 double determineSide(POINT2D *seg1, POINT2D *seg2, POINT2D *point)
628 {
629  return ((seg2->x-seg1->x)*(point->y-seg1->y)-(point->x-seg1->x)*(seg2->y-seg1->y));
630 }
631 
632 /*
633  * This function doesn't test that the point falls on the line defined by
634  * the two points. It assumes that that has already been determined
635  * by having determineSide return within the tolerance. It simply checks
636  * that if the point is on the line, it is within the endpoints.
637  *
638  * returns: 1 if the point is not outside the bounds of the segment
639  * 0 if it is
640  */
641 int isOnSegment(POINT2D *seg1, POINT2D *seg2, POINT2D *point)
642 {
643  double maxX;
644  double maxY;
645  double minX;
646  double minY;
647 
648  if (seg1->x > seg2->x)
649  {
650  maxX = seg1->x;
651  minX = seg2->x;
652  }
653  else
654  {
655  maxX = seg2->x;
656  minX = seg1->x;
657  }
658  if (seg1->y > seg2->y)
659  {
660  maxY = seg1->y;
661  minY = seg2->y;
662  }
663  else
664  {
665  maxY = seg2->y;
666  minY = seg1->y;
667  }
668 
669  POSTGIS_DEBUGF(3, "maxX minX/maxY minY: %.8f %.8f/%.8f %.8f", maxX, minX, maxY, minY);
670 
671  if (maxX < point->x || minX > point->x)
672  {
673  POSTGIS_DEBUGF(3, "X value %.8f falls outside the range %.8f-%.8f", point->x, minX, maxX);
674 
675  return 0;
676  }
677  else if (maxY < point->y || minY > point->y)
678  {
679  POSTGIS_DEBUGF(3, "Y value %.8f falls outside the range %.8f-%.8f", point->y, minY, maxY);
680 
681  return 0;
682  }
683  return 1;
684 }
685 
686 /*
687  * return -1 iff point is outside ring pts
688  * return 1 iff point is inside ring pts
689  * return 0 iff point is on ring pts
690  */
692 {
693  int wn = 0;
694  int i;
695  double side;
696  POINT2D seg1;
697  POINT2D seg2;
698  LWMLINE *lines;
699 
700  POSTGIS_DEBUG(2, "point_in_ring called.");
701 
702  lines = RTreeFindLineSegments(root, point->y);
703  if (!lines)
704  return -1;
705 
706  for (i=0; i<lines->ngeoms; i++)
707  {
708  getPoint2d_p(lines->geoms[i]->points, 0, &seg1);
709  getPoint2d_p(lines->geoms[i]->points, 1, &seg2);
710 
711 
712  side = determineSide(&seg1, &seg2, point);
713 
714  POSTGIS_DEBUGF(3, "segment: (%.8f, %.8f),(%.8f, %.8f)", seg1.x, seg1.y, seg2.x, seg2.y);
715  POSTGIS_DEBUGF(3, "side result: %.8f", side);
716  POSTGIS_DEBUGF(3, "counterclockwise wrap %d, clockwise wrap %d", FP_CONTAINS_BOTTOM(seg1.y,point->y,seg2.y), FP_CONTAINS_BOTTOM(seg2.y,point->y,seg1.y));
717 
718  /* zero length segments are ignored. */
719  if (((seg2.x-seg1.x)*(seg2.x-seg1.x)+(seg2.y-seg1.y)*(seg2.y-seg1.y)) < 1e-12*1e-12)
720  {
721  POSTGIS_DEBUG(3, "segment is zero length... ignoring.");
722 
723  continue;
724  }
725 
726  /* a point on the boundary of a ring is not contained. */
727  /* WAS: if (fabs(side) < 1e-12), see #852 */
728  if (side == 0.0)
729  {
730  if (isOnSegment(&seg1, &seg2, point) == 1)
731  {
732  POSTGIS_DEBUGF(3, "point on ring boundary between points %d, %d", i, i+1);
733 
734  return 0;
735  }
736  }
737 
738  /*
739  * If the point is to the left of the line, and it's rising,
740  * then the line is to the right of the point and
741  * circling counter-clockwise, so incremement.
742  */
743  if (FP_CONTAINS_BOTTOM(seg1.y,point->y,seg2.y) && side>0)
744  {
745  POSTGIS_DEBUG(3, "incrementing winding number.");
746 
747  ++wn;
748  }
749  /*
750  * If the point is to the right of the line, and it's falling,
751  * then the line is to the right of the point and circling
752  * clockwise, so decrement.
753  */
754  else if (FP_CONTAINS_BOTTOM(seg2.y,point->y,seg1.y) && side<0)
755  {
756  POSTGIS_DEBUG(3, "decrementing winding number.");
757 
758  --wn;
759  }
760  }
761 
762  POSTGIS_DEBUGF(3, "winding number %d", wn);
763 
764  if (wn == 0)
765  return -1;
766  return 1;
767 }
768 
769 
770 /*
771  * return -1 iff point is outside ring pts
772  * return 1 iff point is inside ring pts
773  * return 0 iff point is on ring pts
774  */
776 {
777  int wn = 0;
778  int i;
779  double side;
780  POINT2D seg1;
781  POINT2D seg2;
782 
783  POSTGIS_DEBUG(2, "point_in_ring called.");
784 
785 
786  for (i=0; i<pts->npoints-1; i++)
787  {
788  getPoint2d_p(pts, i, &seg1);
789  getPoint2d_p(pts, i+1, &seg2);
790 
791 
792  side = determineSide(&seg1, &seg2, point);
793 
794  POSTGIS_DEBUGF(3, "segment: (%.8f, %.8f),(%.8f, %.8f)", seg1.x, seg1.y, seg2.x, seg2.y);
795  POSTGIS_DEBUGF(3, "side result: %.8f", side);
796  POSTGIS_DEBUGF(3, "counterclockwise wrap %d, clockwise wrap %d", FP_CONTAINS_BOTTOM(seg1.y,point->y,seg2.y), FP_CONTAINS_BOTTOM(seg2.y,point->y,seg1.y));
797 
798  /* zero length segments are ignored. */
799  if (((seg2.x-seg1.x)*(seg2.x-seg1.x)+(seg2.y-seg1.y)*(seg2.y-seg1.y)) < 1e-12*1e-12)
800  {
801  POSTGIS_DEBUG(3, "segment is zero length... ignoring.");
802 
803  continue;
804  }
805 
806  /* a point on the boundary of a ring is not contained. */
807  /* WAS: if (fabs(side) < 1e-12), see #852 */
808  if (side == 0.0)
809  {
810  if (isOnSegment(&seg1, &seg2, point) == 1)
811  {
812  POSTGIS_DEBUGF(3, "point on ring boundary between points %d, %d", i, i+1);
813 
814  return 0;
815  }
816  }
817 
818  /*
819  * If the point is to the left of the line, and it's rising,
820  * then the line is to the right of the point and
821  * circling counter-clockwise, so incremement.
822  */
823  if (FP_CONTAINS_BOTTOM(seg1.y,point->y,seg2.y) && side>0)
824  {
825  POSTGIS_DEBUG(3, "incrementing winding number.");
826 
827  ++wn;
828  }
829  /*
830  * If the point is to the right of the line, and it's falling,
831  * then the line is to the right of the point and circling
832  * clockwise, so decrement.
833  */
834  else if (FP_CONTAINS_BOTTOM(seg2.y,point->y,seg1.y) && side<0)
835  {
836  POSTGIS_DEBUG(3, "decrementing winding number.");
837 
838  --wn;
839  }
840  }
841 
842  POSTGIS_DEBUGF(3, "winding number %d", wn);
843 
844  if (wn == 0)
845  return -1;
846  return 1;
847 }
848 
849 /*
850  * return 0 iff point outside polygon or on boundary
851  * return 1 iff point inside polygon
852  */
853 int point_in_polygon_rtree(RTREE_NODE **root, int ringCount, LWPOINT *point)
854 {
855  int i;
856  POINT2D pt;
857 
858  POSTGIS_DEBUGF(2, "point_in_polygon called for %p %d %p.", root, ringCount, point);
859 
860  getPoint2d_p(point->point, 0, &pt);
861  /* assume bbox short-circuit has already been attempted */
862 
863  if (point_in_ring_rtree(root[0], &pt) != 1)
864  {
865  POSTGIS_DEBUG(3, "point_in_polygon_rtree: outside exterior ring.");
866 
867  return 0;
868  }
869 
870  for (i=1; i<ringCount; i++)
871  {
872  if (point_in_ring_rtree(root[i], &pt) != -1)
873  {
874  POSTGIS_DEBUGF(3, "point_in_polygon_rtree: within hole %d.", i);
875 
876  return 0;
877  }
878  }
879  return 1;
880 }
881 
882 /*
883  * return -1 if point outside polygon
884  * return 0 if point on boundary
885  * return 1 if point inside polygon
886  *
887  * Expected **root order is each exterior ring followed by its holes, eg. EIIEIIEI
888  */
889 int point_in_multipolygon_rtree(RTREE_NODE **root, int polyCount, int *ringCounts, LWPOINT *point)
890 {
891  int i, p, r, in_ring;
892  POINT2D pt;
893  int result = -1;
894 
895  POSTGIS_DEBUGF(2, "point_in_multipolygon_rtree called for %p %d %p.", root, polyCount, point);
896 
897  getPoint2d_p(point->point, 0, &pt);
898  /* assume bbox short-circuit has already been attempted */
899 
900  i = 0; /* the current index into the root array */
901 
902  /* is the point inside any of the sub-polygons? */
903  for ( p = 0; p < polyCount; p++ )
904  {
905  in_ring = point_in_ring_rtree(root[i], &pt);
906  POSTGIS_DEBUGF(4, "point_in_multipolygon_rtree: exterior ring (%d), point_in_ring returned %d", p, in_ring);
907  if ( in_ring == -1 ) /* outside the exterior ring */
908  {
909  POSTGIS_DEBUG(3, "point_in_multipolygon_rtree: outside exterior ring.");
910  }
911  else if ( in_ring == 0 ) /* on the boundary */
912  {
913  POSTGIS_DEBUGF(3, "point_in_multipolygon_rtree: on edge of exterior ring %d", p);
914  return 0;
915  } else {
916  result = in_ring;
917 
918  for(r=1; r<ringCounts[p]; r++)
919  {
920  in_ring = point_in_ring_rtree(root[i+r], &pt);
921  POSTGIS_DEBUGF(4, "point_in_multipolygon_rtree: interior ring (%d), point_in_ring returned %d", r, in_ring);
922  if (in_ring == 1) /* inside a hole => outside the polygon */
923  {
924  POSTGIS_DEBUGF(3, "point_in_multipolygon_rtree: within hole %d of exterior ring %d", r, p);
925  result = -1;
926  break;
927  }
928  if (in_ring == 0) /* on the edge of a hole */
929  {
930  POSTGIS_DEBUGF(3, "point_in_multipolygon_rtree: on edge of hole %d of exterior ring %d", r, p);
931  return 0;
932  }
933  }
934  /* if we have a positive result, we can short-circuit and return it */
935  if ( result != -1)
936  {
937  return result;
938  }
939  }
940  /* increment the index by the total number of rings in the sub-poly */
941  /* we do this here in case we short-cutted out of the poly before looking at all the rings */
942  i += ringCounts[p];
943  }
944 
945  return result; /* -1 = outside, 0 = boundary, 1 = inside */
946 
947 }
948 
949 /*
950  * return -1 iff point outside polygon
951  * return 0 iff point on boundary
952  * return 1 iff point inside polygon
953  */
954 int point_in_polygon(LWPOLY *polygon, LWPOINT *point)
955 {
956  int i, result, in_ring;
957  POINT2D pt;
958 
959  POSTGIS_DEBUG(2, "point_in_polygon called.");
960 
961  getPoint2d_p(point->point, 0, &pt);
962  /* assume bbox short-circuit has already been attempted */
963 
964  /* everything is outside of an empty polygon */
965  if ( polygon->nrings == 0 ) return -1;
966 
967  in_ring = point_in_ring(polygon->rings[0], &pt);
968  if ( in_ring == -1) /* outside the exterior ring */
969  {
970  POSTGIS_DEBUG(3, "point_in_polygon: outside exterior ring.");
971  return -1;
972  }
973  result = in_ring;
974 
975  for (i=1; i<polygon->nrings; i++)
976  {
977  in_ring = point_in_ring(polygon->rings[i], &pt);
978  if (in_ring == 1) /* inside a hole => outside the polygon */
979  {
980  POSTGIS_DEBUGF(3, "point_in_polygon: within hole %d.", i);
981  return -1;
982  }
983  if (in_ring == 0) /* on the edge of a hole */
984  {
985  POSTGIS_DEBUGF(3, "point_in_polygon: on edge of hole %d.", i);
986  return 0;
987  }
988  }
989  return result; /* -1 = outside, 0 = boundary, 1 = inside */
990 }
991 
992 /*
993  * return -1 iff point outside multipolygon
994  * return 0 iff point on multipolygon boundary
995  * return 1 iff point inside multipolygon
996  */
997 int point_in_multipolygon(LWMPOLY *mpolygon, LWPOINT *point)
998 {
999  int i, j, result, in_ring;
1000  POINT2D pt;
1001 
1002  POSTGIS_DEBUG(2, "point_in_polygon called.");
1003 
1004  getPoint2d_p(point->point, 0, &pt);
1005  /* assume bbox short-circuit has already been attempted */
1006 
1007  result = -1;
1008 
1009  for (j = 0; j < mpolygon->ngeoms; j++ )
1010  {
1011 
1012  LWPOLY *polygon = mpolygon->geoms[j];
1013 
1014  /* everything is outside of an empty polygon */
1015  if ( polygon->nrings == 0 ) continue;
1016 
1017  in_ring = point_in_ring(polygon->rings[0], &pt);
1018  if ( in_ring == -1) /* outside the exterior ring */
1019  {
1020  POSTGIS_DEBUG(3, "point_in_polygon: outside exterior ring.");
1021  continue;
1022  }
1023  if ( in_ring == 0 )
1024  {
1025  return 0;
1026  }
1027 
1028  result = in_ring;
1029 
1030  for (i=1; i<polygon->nrings; i++)
1031  {
1032  in_ring = point_in_ring(polygon->rings[i], &pt);
1033  if (in_ring == 1) /* inside a hole => outside the polygon */
1034  {
1035  POSTGIS_DEBUGF(3, "point_in_polygon: within hole %d.", i);
1036  result = -1;
1037  break;
1038  }
1039  if (in_ring == 0) /* on the edge of a hole */
1040  {
1041  POSTGIS_DEBUGF(3, "point_in_polygon: on edge of hole %d.", i);
1042  return 0;
1043  }
1044  }
1045  if ( result != -1)
1046  {
1047  return result;
1048  }
1049  }
1050  return result;
1051 }
1052 
1053 
1054 /*******************************************************************************
1055  * End of "Fast Winding Number Inclusion of a Point in a Polygon" derivative.
1056  ******************************************************************************/
1057 
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:536
double x
Definition: liblwgeom.h:336
#define LINETYPE
Definition: liblwgeom.h:71
uint32_t gserialized_get_type(const GSERIALIZED *s)
Extract the geometry type from the serialized form (it hides in the anonymous data area...
Definition: g_serialized.c:55
Datum LWGEOM_simplify2d(PG_FUNCTION_ARGS)
POINTARRAY * ptarray_construct(char hasz, char hasm, uint32_t npoints)
Construct an empty pointarray, allocating storage and setting the npoints, but not filling in any inf...
Definition: ptarray.c:62
GBOX * bbox
Definition: liblwgeom.h:382
double m
Definition: liblwgeom.h:336
LWGEOM * lwgeom_grid(const LWGEOM *lwgeom, const gridspec *grid)
Definition: lwgeom.c:1868
LWCOLLECTION * lwcollection_construct(uint8_t type, int srid, GBOX *bbox, uint32_t ngeoms, LWGEOM **geoms)
Definition: lwcollection.c:30
char * r
Definition: cu_in_wkt.c:24
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
int npoints
Definition: liblwgeom.h:355
int gserialized_has_m(const GSERIALIZED *gser)
Check if a GSERIALIZED has an M ordinate.
Definition: g_serialized.c:29
int lwline_crossing_direction(const LWLINE *l1, const LWLINE *l2)
Given two lines, characterize how (and if) they cross each other.
Definition: lwalgorithm.c:446
The following struct and methods are used for a 1D RTree implementation, described at: http://lin-ear...
Definition: lwgeom_rtree.h:21
Datum area(PG_FUNCTION_ARGS)
void lwpoint_free(LWPOINT *pt)
Definition: lwpoint.c:182
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1050
POINTARRAY * ptarray_substring(POINTARRAY *pa, double d1, double d2, double tolerance)
start location (distance from start / total distance) end location (distance from start / total dist...
Definition: ptarray.c:1058
#define MULTIPOINTTYPE
Definition: liblwgeom.h:73
double ptarray_length_2d(const POINTARRAY *pts)
Find the 2d length of the given POINTARRAY (even if it&#39;s 3d)
Definition: ptarray.c:1645
Datum LWGEOM_snaptogrid(PG_FUNCTION_ARGS)
int point_in_multipolygon_rtree(RTREE_NODE **root, int polyCount, int *ringCounts, LWPOINT *point)
LWMLINE * RTreeFindLineSegments(RTREE_NODE *root, double value)
Retrieves a collection of line segments given the root and crossing value.
Definition: lwgeom_rtree.c:436
double distance2d_pt_pt(const POINT2D *p1, const POINT2D *p2)
The old function nessecary for ptarray_segmentize2d in ptarray.c.
Definition: measures.c:2301
void error_if_srid_mismatch(int srid1, int srid2)
Definition: lwutil.c:341
LWPOINT * lwpoint_construct_empty(int srid, char hasz, char hasm)
Definition: lwpoint.c:120
LWPOINT * lwgeom_as_lwpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:80
int32_t srid
Definition: liblwgeom.h:405
Datum LWGEOM_SetEffectiveArea(PG_FUNCTION_ARGS)
POINTARRAY * point
Definition: liblwgeom.h:395
int gserialized_has_z(const GSERIALIZED *gser)
Check if a GSERIALIZED has a Z ordinate.
Definition: g_serialized.c:24
int point_in_multipolygon(LWMPOLY *mpolygon, LWPOINT *point)
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:836
int32_t srid
Definition: liblwgeom.h:383
void interpolate_point4d(POINT4D *A, POINT4D *B, POINT4D *I, double F)
Find interpolation point I between point A and point B so that the len(AI) == len(AB)*F and I falls o...
Definition: lwgeom_api.c:806
int ngeoms
Definition: liblwgeom.h:465
LWGEOM * lwgeom_set_effective_area(const LWGEOM *igeom, int set_area, double trshld)
int gserialized_is_empty(const GSERIALIZED *g)
Check if a GSERIALIZED is empty without deserializing first.
Definition: g_serialized.c:139
Datum LWGEOM_line_substring(PG_FUNCTION_ARGS)
double x
Definition: liblwgeom.h:312
Datum LWGEOM_snaptogrid_pointoff(PG_FUNCTION_ARGS)
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:188
LWGEOM * geom
#define LW_FALSE
Definition: liblwgeom.h:62
LWPOLY ** geoms
Definition: liblwgeom.h:480
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:61
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:29
LWGEOM * lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed)
Definition: lwgeom.c:1558
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:172
LWMLINE * lwgeom_as_lwmline(const LWGEOM *lwgeom)
Definition: lwgeom.c:161
POINTARRAY ** rings
Definition: liblwgeom.h:441
int isOnSegment(POINT2D *seg1, POINT2D *seg2, POINT2D *point)
int nrings
Definition: liblwgeom.h:439
double y
Definition: liblwgeom.h:312
int getPoint2d_p(const POINTARRAY *pa, int n, POINT2D *point)
Definition: lwgeom_api.c:448
#define FLAGS_GET_Z(flags)
Macros for manipulating the &#39;flags&#39; byte.
Definition: liblwgeom.h:124
double z
Definition: liblwgeom.h:336
int point_in_ring_rtree(RTREE_NODE *root, POINT2D *point)
Datum distance(PG_FUNCTION_ARGS)
int ngeoms
Definition: liblwgeom.h:478
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:89
uint8_t flags
Definition: liblwgeom.h:392
LWLINE ** geoms
Definition: liblwgeom.h:467
int point_in_polygon(LWPOLY *polygon, LWPOINT *point)
GSERIALIZED * geometry_serialize(LWGEOM *lwgeom)
Datum LWGEOM_line_interpolate_point(PG_FUNCTION_ARGS)
Datum ST_LineCrossingDirection(PG_FUNCTION_ARGS)
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:70
void lwgeom_add_bbox(LWGEOM *lwgeom)
Compute a bbox if not already computed.
Definition: lwgeom.c:599
#define FLAGS_GET_M(flags)
Definition: liblwgeom.h:125
int point_in_ring(POINTARRAY *pts, POINT2D *point)
uint8_t type
Definition: liblwgeom.h:380
type
Definition: ovdump.py:41
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:254
int point_in_polygon_rtree(RTREE_NODE **root, int ringCount, LWPOINT *point)
LWPOINT * lwpoint_construct(int srid, GBOX *bbox, POINTARRAY *point)
Definition: lwpoint.c:98
void * lwalloc(size_t size)
Definition: lwutil.c:199
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members) ...
Definition: lwgeom.c:1297
void lwmline_release(LWMLINE *lwline)
Definition: lwmline.c:18
#define FP_CONTAINS_BOTTOM(A, X, B)
PG_FUNCTION_INFO_V1(LWGEOM_simplify2d)
double y
Definition: liblwgeom.h:336
#define MULTILINETYPE
Definition: liblwgeom.h:74
void lwline_release(LWLINE *lwline)
Definition: lwline.c:121
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:843
int32_t gserialized_get_srid(const GSERIALIZED *s)
Extract the SRID from the serialized form (it is packed into three bytes so this is a handy function)...
Definition: g_serialized.c:69
if(!(yy_init))
int getPoint4d_p(const POINTARRAY *pa, int n, POINT4D *point)
Definition: lwgeom_api.c:231
#define COLLECTIONTYPE
Definition: liblwgeom.h:76
This library is the generic geometry handling section of PostGIS.
double determineSide(POINT2D *seg1, POINT2D *seg2, POINT2D *point)
POINTARRAY * points
Definition: liblwgeom.h:406
Snap to grid.