PostGIS  2.5.0beta1dev-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  * PostGIS is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * PostGIS is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18  *
19  **********************************************************************
20  *
21  * Copyright (C) 2001-2005 Refractions Research Inc.
22  *
23  **********************************************************************/
24 
25 
26 #include "postgres.h"
27 #include "funcapi.h"
28 #include "fmgr.h"
29 #include "liblwgeom.h"
30 #include "liblwgeom_internal.h" /* For FP comparators. */
31 #include "lwgeom_pg.h"
32 #include "math.h"
33 #include "lwgeom_rtree.h"
35 
36 
37 #include "access/htup_details.h"
38 
39 /* Prototypes */
40 Datum LWGEOM_simplify2d(PG_FUNCTION_ARGS);
41 Datum LWGEOM_SetEffectiveArea(PG_FUNCTION_ARGS);
42 Datum LWGEOM_line_interpolate_point(PG_FUNCTION_ARGS);
43 Datum ST_LineCrossingDirection(PG_FUNCTION_ARGS);
44 Datum ST_MinimumBoundingRadius(PG_FUNCTION_ARGS);
45 Datum ST_MinimumBoundingCircle(PG_FUNCTION_ARGS);
46 Datum ST_GeometricMedian(PG_FUNCTION_ARGS);
47 Datum ST_IsPolygonCCW(PG_FUNCTION_ARGS);
48 Datum ST_IsPolygonCW(PG_FUNCTION_ARGS);
49 
50 
51 static double determineSide(const POINT2D *seg1, const POINT2D *seg2, const POINT2D *point);
52 static int isOnSegment(const POINT2D *seg1, const POINT2D *seg2, const POINT2D *point);
53 static int point_in_ring(POINTARRAY *pts, const POINT2D *point);
54 static int point_in_ring_rtree(RTREE_NODE *root, const POINT2D *point);
55 
56 /***********************************************************************
57  * Simple Douglas-Peucker line simplification.
58  * No checks are done to avoid introduction of self-intersections.
59  * No topology relations are considered.
60  *
61  * --strk@kbt.io;
62  ***********************************************************************/
63 
65 Datum LWGEOM_simplify2d(PG_FUNCTION_ARGS)
66 {
67  GSERIALIZED *geom = PG_GETARG_GSERIALIZED_P(0);
68  double dist = PG_GETARG_FLOAT8(1);
69  GSERIALIZED *result;
70  int type = gserialized_get_type(geom);
71  LWGEOM *in, *out;
72  bool preserve_collapsed = false;
73 
74  /* Handle optional argument to preserve collapsed features */
75  if ( PG_NARGS() > 2 && ! PG_ARGISNULL(2) )
76  preserve_collapsed = true;
77 
78  /* Can't simplify points! */
79  if ( type == POINTTYPE || type == MULTIPOINTTYPE )
80  PG_RETURN_POINTER(geom);
81 
82  in = lwgeom_from_gserialized(geom);
83 
84  out = lwgeom_simplify(in, dist, preserve_collapsed);
85  if ( ! out ) PG_RETURN_NULL();
86 
87  /* COMPUTE_BBOX TAINTING */
88  if ( in->bbox ) lwgeom_add_bbox(out);
89 
90  result = geometry_serialize(out);
91  lwgeom_free(out);
92  PG_FREE_IF_COPY(geom, 0);
93  PG_RETURN_POINTER(result);
94 }
95 
97 Datum LWGEOM_SetEffectiveArea(PG_FUNCTION_ARGS)
98 {
99  GSERIALIZED *geom = PG_GETARG_GSERIALIZED_P(0);
100  GSERIALIZED *result;
101  int type = gserialized_get_type(geom);
102  LWGEOM *in;
103  LWGEOM *out;
104  double area=0;
105  int set_area=0;
106 
107  if ( type == POINTTYPE || type == MULTIPOINTTYPE )
108  PG_RETURN_POINTER(geom);
109 
110  if ( (PG_NARGS()>1) && (!PG_ARGISNULL(1)) )
111  area = PG_GETARG_FLOAT8(1);
112 
113  if ( (PG_NARGS()>2) && (!PG_ARGISNULL(2)) )
114  set_area = PG_GETARG_INT32(2);
115 
116  in = lwgeom_from_gserialized(geom);
117 
118  out = lwgeom_set_effective_area(in,set_area, area);
119  if ( ! out ) PG_RETURN_NULL();
120 
121  /* COMPUTE_BBOX TAINTING */
122  if ( in->bbox ) lwgeom_add_bbox(out);
123 
124  result = geometry_serialize(out);
125  lwgeom_free(out);
126  PG_FREE_IF_COPY(geom, 0);
127  PG_RETURN_POINTER(result);
128 }
129 
131 Datum LWGEOM_ChaikinSmoothing(PG_FUNCTION_ARGS)
132 {
133  GSERIALIZED *geom = PG_GETARG_GSERIALIZED_P(0);
134  GSERIALIZED *result;
135  int type = gserialized_get_type(geom);
136  LWGEOM *in;
137  LWGEOM *out;
138  int preserve_endpoints=1;
139  int n_iterations=1;
140 
141  if ( type == POINTTYPE || type == MULTIPOINTTYPE )
142  PG_RETURN_POINTER(geom);
143 
144  if ( (PG_NARGS()>1) && (!PG_ARGISNULL(1)) )
145  n_iterations = PG_GETARG_INT32(1);
146 
147  if (n_iterations< 1 || n_iterations>5)
148  elog(ERROR,"Number of iterations must be between 1 and 5 : %s", __func__);
149 
150  if ( (PG_NARGS()>2) && (!PG_ARGISNULL(2)) )
151  {
152  if(PG_GETARG_BOOL(2))
153  preserve_endpoints = 1;
154  else
155  preserve_endpoints = 0;
156  }
157 
158  in = lwgeom_from_gserialized(geom);
159 
160  out = lwgeom_chaikin(in, n_iterations, preserve_endpoints);
161  if ( ! out ) PG_RETURN_NULL();
162 
163  /* COMPUTE_BBOX TAINTING */
164  if ( in->bbox ) lwgeom_add_bbox(out);
165 
166  result = geometry_serialize(out);
167  lwgeom_free(out);
168  PG_FREE_IF_COPY(geom, 0);
169  PG_RETURN_POINTER(result);
170 }
171 
172 
173 /***********************************************************************
174  * --strk@kbt.io;
175  ***********************************************************************/
176 
177 /***********************************************************************
178  * Interpolate a point along a line, useful for Geocoding applications
179  * SELECT line_interpolate_point( 'LINESTRING( 0 0, 2 2'::geometry, .5 )
180  * returns POINT( 1 1 ).
181  ***********************************************************************/
183 Datum LWGEOM_line_interpolate_point(PG_FUNCTION_ARGS)
184 {
185  GSERIALIZED *gser = PG_GETARG_GSERIALIZED_P(0);
186  GSERIALIZED *result;
187  double distance_fraction = PG_GETARG_FLOAT8(1);
188  int repeat = PG_NARGS() > 2 && PG_GETARG_BOOL(2);
189  int srid = gserialized_get_srid(gser);
190  LWLINE* lwline;
191  LWGEOM* lwresult;
192  POINTARRAY* opa;
193 
194  if ( distance_fraction < 0 || distance_fraction > 1 )
195  {
196  elog(ERROR,"line_interpolate_point: 2nd arg isn't within [0,1]");
197  PG_FREE_IF_COPY(gser, 0);
198  PG_RETURN_NULL();
199  }
200 
201  if ( gserialized_get_type(gser) != LINETYPE )
202  {
203  elog(ERROR,"line_interpolate_point: 1st arg isn't a line");
204  PG_FREE_IF_COPY(gser, 0);
205  PG_RETURN_NULL();
206  }
207 
209  opa = lwline_interpolate_points(lwline, distance_fraction, repeat);
210 
211  lwgeom_free(lwline_as_lwgeom(lwline));
212  PG_FREE_IF_COPY(gser, 0);
213 
214  if (opa->npoints <= 1)
215  {
216  lwresult = lwpoint_as_lwgeom(lwpoint_construct(srid, NULL, opa));
217  } else {
218  lwresult = lwmpoint_as_lwgeom(lwmpoint_construct(srid, opa));
219  }
220 
221  result = geometry_serialize(lwresult);
222  lwgeom_free(lwresult);
223 
224  PG_RETURN_POINTER(result);
225 }
226 /***********************************************************************
227  * --jsunday@rochgrp.com;
228  ***********************************************************************/
229 
230 /***********************************************************************
231  *
232  * Grid application function for postgis.
233  *
234  * WHAT IS
235  * -------
236  *
237  * This is a grid application function for postgis.
238  * You use it to stick all of a geometry points to
239  * a custom grid defined by its origin and cell size
240  * in geometry units.
241  *
242  * Points reduction is obtained by collapsing all
243  * consecutive points falling on the same grid node and
244  * removing all consecutive segments S1,S2 having
245  * S2.startpoint = S1.endpoint and S2.endpoint = S1.startpoint.
246  *
247  * ISSUES
248  * ------
249  *
250  * Only 2D is contemplated in grid application.
251  *
252  * Consecutive coincident segments removal does not work
253  * on first/last segments (They are not considered consecutive).
254  *
255  * Grid application occurs on a geometry subobject basis, thus no
256  * points reduction occurs for multipoint geometries.
257  *
258  * USAGE TIPS
259  * ----------
260  *
261  * Run like this:
262  *
263  * SELECT SnapToGrid(<geometry>, <originX>, <originY>, <sizeX>, <sizeY>);
264  *
265  * Grid cells are not visible on a screen as long as the screen
266  * point size is equal or bigger then the grid cell size.
267  * This assumption may be used to reduce the number of points in
268  * a map for a given display scale.
269  *
270  * Keeping multiple resolution versions of the same data may be used
271  * in conjunction with MINSCALE/MAXSCALE keywords of mapserv to speed
272  * up rendering.
273  *
274  * Check also the use of DP simplification to reduce grid visibility.
275  * I'm still researching about the relationship between grid size and
276  * DP epsilon values - please tell me if you know more about this.
277  *
278  *
279  * --strk@kbt.io;
280  *
281  ***********************************************************************/
282 
283 #define CHECK_RING_IS_CLOSE
284 #define SAMEPOINT(a,b) ((a)->x==(b)->x&&(a)->y==(b)->y)
285 
286 
287 
288 /* Forward declarations */
289 Datum LWGEOM_snaptogrid(PG_FUNCTION_ARGS);
290 Datum LWGEOM_snaptogrid_pointoff(PG_FUNCTION_ARGS);
291 
292 
293 
295 Datum LWGEOM_snaptogrid(PG_FUNCTION_ARGS)
296 {
297  LWGEOM *in_lwgeom;
298  GSERIALIZED *out_geom = NULL;
299  LWGEOM *out_lwgeom;
300  gridspec grid;
301 
302  GSERIALIZED *in_geom = PG_GETARG_GSERIALIZED_P(0);
303 
304  /* Set grid values to zero to start */
305  memset(&grid, 0, sizeof(gridspec));
306 
307  grid.ipx = PG_GETARG_FLOAT8(1);
308  grid.ipy = PG_GETARG_FLOAT8(2);
309  grid.xsize = PG_GETARG_FLOAT8(3);
310  grid.ysize = PG_GETARG_FLOAT8(4);
311 
312  /* Return input geometry if input geometry is empty */
313  if ( gserialized_is_empty(in_geom) )
314  {
315  PG_RETURN_POINTER(in_geom);
316  }
317 
318  /* Return input geometry if input grid is meaningless */
319  if ( grid.xsize==0 && grid.ysize==0 && grid.zsize==0 && grid.msize==0 )
320  {
321  PG_RETURN_POINTER(in_geom);
322  }
323 
324  in_lwgeom = lwgeom_from_gserialized(in_geom);
325 
326  POSTGIS_DEBUGF(3, "SnapToGrid got a %s", lwtype_name(in_lwgeom->type));
327 
328  out_lwgeom = lwgeom_grid(in_lwgeom, &grid);
329  if ( out_lwgeom == NULL ) PG_RETURN_NULL();
330 
331  /* COMPUTE_BBOX TAINTING */
332  if ( in_lwgeom->bbox )
333  lwgeom_add_bbox(out_lwgeom);
334 
335 
336  POSTGIS_DEBUGF(3, "SnapToGrid made a %s", lwtype_name(out_lwgeom->type));
337 
338  out_geom = geometry_serialize(out_lwgeom);
339 
340  PG_RETURN_POINTER(out_geom);
341 }
342 
343 
344 #if POSTGIS_DEBUG_LEVEL >= 4
345 /* Print grid using given reporter */
346 static void
347 grid_print(const gridspec *grid)
348 {
349  lwpgnotice("GRID(%g %g %g %g, %g %g %g %g)",
350  grid->ipx, grid->ipy, grid->ipz, grid->ipm,
351  grid->xsize, grid->ysize, grid->zsize, grid->msize);
352 }
353 #endif
354 
355 
357 Datum LWGEOM_snaptogrid_pointoff(PG_FUNCTION_ARGS)
358 {
359  GSERIALIZED *in_geom, *in_point;
360  LWGEOM *in_lwgeom;
361  LWPOINT *in_lwpoint;
362  GSERIALIZED *out_geom = NULL;
363  LWGEOM *out_lwgeom;
364  gridspec grid;
365  /* BOX3D box3d; */
366  POINT4D offsetpoint;
367 
368  in_geom = PG_GETARG_GSERIALIZED_P(0);
369 
370  /* Return input geometry if input geometry is empty */
371  if ( gserialized_is_empty(in_geom) )
372  {
373  PG_RETURN_POINTER(in_geom);
374  }
375 
376  in_point = PG_GETARG_GSERIALIZED_P(1);
377  in_lwpoint = lwgeom_as_lwpoint(lwgeom_from_gserialized(in_point));
378  if ( in_lwpoint == NULL )
379  {
380  lwpgerror("Offset geometry must be a point");
381  }
382 
383  grid.xsize = PG_GETARG_FLOAT8(2);
384  grid.ysize = PG_GETARG_FLOAT8(3);
385  grid.zsize = PG_GETARG_FLOAT8(4);
386  grid.msize = PG_GETARG_FLOAT8(5);
387 
388  /* Take offsets from point geometry */
389  getPoint4d_p(in_lwpoint->point, 0, &offsetpoint);
390  grid.ipx = offsetpoint.x;
391  grid.ipy = offsetpoint.y;
392  if (FLAGS_GET_Z(in_lwpoint->flags) ) grid.ipz = offsetpoint.z;
393  else grid.ipz=0;
394  if (FLAGS_GET_M(in_lwpoint->flags) ) grid.ipm = offsetpoint.m;
395  else grid.ipm=0;
396 
397 #if POSTGIS_DEBUG_LEVEL >= 4
398  grid_print(&grid);
399 #endif
400 
401  /* Return input geometry if input grid is meaningless */
402  if ( grid.xsize==0 && grid.ysize==0 && grid.zsize==0 && grid.msize==0 )
403  {
404  PG_RETURN_POINTER(in_geom);
405  }
406 
407  in_lwgeom = lwgeom_from_gserialized(in_geom);
408 
409  POSTGIS_DEBUGF(3, "SnapToGrid got a %s", lwtype_name(in_lwgeom->type));
410 
411  out_lwgeom = lwgeom_grid(in_lwgeom, &grid);
412  if ( out_lwgeom == NULL ) PG_RETURN_NULL();
413 
414  /* COMPUTE_BBOX TAINTING */
415  if ( in_lwgeom->bbox ) lwgeom_add_bbox(out_lwgeom);
416 
417  POSTGIS_DEBUGF(3, "SnapToGrid made a %s", lwtype_name(out_lwgeom->type));
418 
419  out_geom = geometry_serialize(out_lwgeom);
420 
421  PG_RETURN_POINTER(out_geom);
422 }
423 
424 
425 /*
426 ** crossingDirection(line1, line2)
427 **
428 ** Determines crossing direction of line2 relative to line1.
429 ** Only accepts LINESTRING as parameters!
430 */
432 Datum ST_LineCrossingDirection(PG_FUNCTION_ARGS)
433 {
434  int type1, type2, rv;
435  LWLINE *l1 = NULL;
436  LWLINE *l2 = NULL;
437  GSERIALIZED *geom1 = PG_GETARG_GSERIALIZED_P(0);
438  GSERIALIZED *geom2 = PG_GETARG_GSERIALIZED_P(1);
439 
441 
442  type1 = gserialized_get_type(geom1);
443  type2 = gserialized_get_type(geom2);
444 
445  if ( type1 != LINETYPE || type2 != LINETYPE )
446  {
447  elog(ERROR,"This function only accepts LINESTRING as arguments.");
448  PG_RETURN_NULL();
449  }
450 
453 
454  rv = lwline_crossing_direction(l1, l2);
455 
456  PG_FREE_IF_COPY(geom1, 0);
457  PG_FREE_IF_COPY(geom2, 1);
458 
459  PG_RETURN_INT32(rv);
460 
461 }
462 
463 
464 
465 /***********************************************************************
466  * --strk@kbt.io
467  ***********************************************************************/
468 
469 Datum LWGEOM_line_substring(PG_FUNCTION_ARGS);
470 
472 Datum LWGEOM_line_substring(PG_FUNCTION_ARGS)
473 {
474  GSERIALIZED *geom = PG_GETARG_GSERIALIZED_P(0);
475  double from = PG_GETARG_FLOAT8(1);
476  double to = PG_GETARG_FLOAT8(2);
477  LWGEOM *olwgeom;
478  POINTARRAY *ipa, *opa;
479  GSERIALIZED *ret;
480  int type = gserialized_get_type(geom);
481 
482  if ( from < 0 || from > 1 )
483  {
484  elog(ERROR,"line_interpolate_point: 2nd arg isn't within [0,1]");
485  PG_RETURN_NULL();
486  }
487 
488  if ( to < 0 || to > 1 )
489  {
490  elog(ERROR,"line_interpolate_point: 3rd arg isn't within [0,1]");
491  PG_RETURN_NULL();
492  }
493 
494  if ( from > to )
495  {
496  elog(ERROR, "2nd arg must be smaller then 3rd arg");
497  PG_RETURN_NULL();
498  }
499 
500  if ( type == LINETYPE )
501  {
503 
504  if ( lwgeom_is_empty((LWGEOM*)iline) )
505  {
506  /* TODO return empty line */
507  lwline_release(iline);
508  PG_FREE_IF_COPY(geom, 0);
509  PG_RETURN_NULL();
510  }
511 
512  ipa = iline->points;
513 
514  opa = ptarray_substring(ipa, from, to, 0);
515 
516  if ( opa->npoints == 1 ) /* Point returned */
517  olwgeom = (LWGEOM *)lwpoint_construct(iline->srid, NULL, opa);
518  else
519  olwgeom = (LWGEOM *)lwline_construct(iline->srid, NULL, opa);
520 
521  }
522  else if ( type == MULTILINETYPE )
523  {
524  LWMLINE *iline;
525  uint32_t i = 0, g = 0;
526  int homogeneous = LW_TRUE;
527  LWGEOM **geoms = NULL;
528  double length = 0.0, sublength = 0.0, minprop = 0.0, maxprop = 0.0;
529 
531 
532  if ( lwgeom_is_empty((LWGEOM*)iline) )
533  {
534  /* TODO return empty collection */
535  lwmline_release(iline);
536  PG_FREE_IF_COPY(geom, 0);
537  PG_RETURN_NULL();
538  }
539 
540  /* Calculate the total length of the mline */
541  for ( i = 0; i < iline->ngeoms; i++ )
542  {
543  LWLINE *subline = (LWLINE*)iline->geoms[i];
544  if ( subline->points && subline->points->npoints > 1 )
545  length += ptarray_length_2d(subline->points);
546  }
547 
548  geoms = lwalloc(sizeof(LWGEOM*) * iline->ngeoms);
549 
550  /* Slice each sub-geometry of the multiline */
551  for ( i = 0; i < iline->ngeoms; i++ )
552  {
553  LWLINE *subline = (LWLINE*)iline->geoms[i];
554  double subfrom = 0.0, subto = 0.0;
555 
556  if ( subline->points && subline->points->npoints > 1 )
557  sublength += ptarray_length_2d(subline->points);
558 
559  /* Calculate proportions for this subline */
560  minprop = maxprop;
561  maxprop = sublength / length;
562 
563  /* This subline doesn't reach the lowest proportion requested
564  or is beyond the highest proporton */
565  if ( from > maxprop || to < minprop )
566  continue;
567 
568  if ( from <= minprop )
569  subfrom = 0.0;
570  if ( to >= maxprop )
571  subto = 1.0;
572 
573  if ( from > minprop && from <= maxprop )
574  subfrom = (from - minprop) / (maxprop - minprop);
575 
576  if ( to < maxprop && to >= minprop )
577  subto = (to - minprop) / (maxprop - minprop);
578 
579 
580  opa = ptarray_substring(subline->points, subfrom, subto, 0);
581  if ( opa && opa->npoints > 0 )
582  {
583  if ( opa->npoints == 1 ) /* Point returned */
584  {
585  geoms[g] = (LWGEOM *)lwpoint_construct(SRID_UNKNOWN, NULL, opa);
586  homogeneous = LW_FALSE;
587  }
588  else
589  {
590  geoms[g] = (LWGEOM *)lwline_construct(SRID_UNKNOWN, NULL, opa);
591  }
592  g++;
593  }
594 
595 
596 
597  }
598  /* If we got any points, we need to return a GEOMETRYCOLLECTION */
599  if ( ! homogeneous )
600  type = COLLECTIONTYPE;
601 
602  olwgeom = (LWGEOM*)lwcollection_construct(type, iline->srid, NULL, g, geoms);
603  }
604  else
605  {
606  elog(ERROR,"line_substring: 1st arg isn't a line");
607  PG_RETURN_NULL();
608  }
609 
610  ret = geometry_serialize(olwgeom);
611  lwgeom_free(olwgeom);
612  PG_FREE_IF_COPY(geom, 0);
613  PG_RETURN_POINTER(ret);
614 
615 }
616 
617 /*******************************************************************************
618  * The following is based on the "Fast Winding Number Inclusion of a Point
619  * in a Polygon" algorithm by Dan Sunday.
620  * http://softsurfer.com/Archive/algorithm_0103/algorithm_0103.htm#Winding%20Number
621  ******************************************************************************/
622 
623 /*
624  * returns: >0 for a point to the left of the segment,
625  * <0 for a point to the right of the segment,
626  * 0 for a point on the segment
627  */
628 static double determineSide(const POINT2D *seg1, const POINT2D *seg2, const POINT2D *point)
629 {
630  return ((seg2->x-seg1->x)*(point->y-seg1->y)-(point->x-seg1->x)*(seg2->y-seg1->y));
631 }
632 
633 /*
634  * This function doesn't test that the point falls on the line defined by
635  * the two points. It assumes that that has already been determined
636  * by having determineSide return within the tolerance. It simply checks
637  * that if the point is on the line, it is within the endpoints.
638  *
639  * returns: 1 if the point is not outside the bounds of the segment
640  * 0 if it is
641  */
642 static int isOnSegment(const POINT2D *seg1, const POINT2D *seg2, const POINT2D *point)
643 {
644  double maxX;
645  double maxY;
646  double minX;
647  double minY;
648 
649  if (seg1->x > seg2->x)
650  {
651  maxX = seg1->x;
652  minX = seg2->x;
653  }
654  else
655  {
656  maxX = seg2->x;
657  minX = seg1->x;
658  }
659  if (seg1->y > seg2->y)
660  {
661  maxY = seg1->y;
662  minY = seg2->y;
663  }
664  else
665  {
666  maxY = seg2->y;
667  minY = seg1->y;
668  }
669 
670  POSTGIS_DEBUGF(3, "maxX minX/maxY minY: %.8f %.8f/%.8f %.8f", maxX, minX, maxY, minY);
671 
672  if (maxX < point->x || minX > point->x)
673  {
674  POSTGIS_DEBUGF(3, "X value %.8f falls outside the range %.8f-%.8f", point->x, minX, maxX);
675 
676  return 0;
677  }
678  else if (maxY < point->y || minY > point->y)
679  {
680  POSTGIS_DEBUGF(3, "Y value %.8f falls outside the range %.8f-%.8f", point->y, minY, maxY);
681 
682  return 0;
683  }
684  return 1;
685 }
686 
687 /*
688  * return -1 iff point is outside ring pts
689  * return 1 iff point is inside ring pts
690  * return 0 iff point is on ring pts
691  */
692 static int point_in_ring_rtree(RTREE_NODE *root, const POINT2D *point)
693 {
694  int wn = 0;
695  uint32_t i;
696  double side;
697  const POINT2D *seg1;
698  const POINT2D *seg2;
699  LWMLINE *lines;
700 
701  POSTGIS_DEBUG(2, "point_in_ring called.");
702 
703  lines = RTreeFindLineSegments(root, point->y);
704  if (!lines)
705  return -1;
706 
707  for (i=0; i<lines->ngeoms; i++)
708  {
709  seg1 = getPoint2d_cp(lines->geoms[i]->points, 0);
710  seg2 = getPoint2d_cp(lines->geoms[i]->points, 1);
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 increment.
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  */
775 static int point_in_ring(POINTARRAY *pts, const POINT2D *point)
776 {
777  int wn = 0;
778  uint32_t i;
779  double side;
780  const POINT2D* seg1;
781  const POINT2D* seg2;
782 
783  POSTGIS_DEBUG(2, "point_in_ring called.");
784 
785  seg2 = getPoint2d_cp(pts, 0);
786  for (i=0; i<pts->npoints-1; i++)
787  {
788  seg1 = seg2;
789  seg2 = getPoint2d_cp(pts, i+1);
790 
791  side = determineSide(seg1, seg2, point);
792 
793  POSTGIS_DEBUGF(3, "segment: (%.8f, %.8f),(%.8f, %.8f)", seg1->x, seg1->y, seg2->x, seg2->y);
794  POSTGIS_DEBUGF(3, "side result: %.8f", side);
795  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));
796 
797  /* zero length segments are ignored. */
798  if (((seg2->x - seg1->x)*(seg2->x - seg1->x) + (seg2->y - seg1->y)*(seg2->y - seg1->y)) < 1e-12*1e-12)
799  {
800  POSTGIS_DEBUG(3, "segment is zero length... ignoring.");
801 
802  continue;
803  }
804 
805  /* a point on the boundary of a ring is not contained. */
806  /* WAS: if (fabs(side) < 1e-12), see #852 */
807  if (side == 0.0)
808  {
809  if (isOnSegment(seg1, seg2, point) == 1)
810  {
811  POSTGIS_DEBUGF(3, "point on ring boundary between points %d, %d", i, i+1);
812 
813  return 0;
814  }
815  }
816 
817  /*
818  * If the point is to the left of the line, and it's rising,
819  * then the line is to the right of the point and
820  * circling counter-clockwise, so increment.
821  */
822  if (FP_CONTAINS_BOTTOM(seg1->y, point->y, seg2->y) && side>0)
823  {
824  POSTGIS_DEBUG(3, "incrementing winding number.");
825 
826  ++wn;
827  }
828  /*
829  * If the point is to the right of the line, and it's falling,
830  * then the line is to the right of the point and circling
831  * clockwise, so decrement.
832  */
833  else if (FP_CONTAINS_BOTTOM(seg2->y, point->y, seg1->y) && side<0)
834  {
835  POSTGIS_DEBUG(3, "decrementing winding number.");
836 
837  --wn;
838  }
839  }
840 
841  POSTGIS_DEBUGF(3, "winding number %d", wn);
842 
843  if (wn == 0)
844  return -1;
845  return 1;
846 }
847 
848 /*
849  * return 0 iff point outside polygon or on boundary
850  * return 1 iff point inside polygon
851  */
852 int point_in_polygon_rtree(RTREE_NODE **root, int ringCount, LWPOINT *point)
853 {
854  int i;
855  POINT2D pt;
856 
857  POSTGIS_DEBUGF(2, "point_in_polygon called for %p %d %p.", root, ringCount, point);
858 
859  getPoint2d_p(point->point, 0, &pt);
860  /* assume bbox short-circuit has already been attempted */
861 
862  if (point_in_ring_rtree(root[0], &pt) != 1)
863  {
864  POSTGIS_DEBUG(3, "point_in_polygon_rtree: outside exterior ring.");
865 
866  return 0;
867  }
868 
869  for (i=1; i<ringCount; i++)
870  {
871  if (point_in_ring_rtree(root[i], &pt) != -1)
872  {
873  POSTGIS_DEBUGF(3, "point_in_polygon_rtree: within hole %d.", i);
874 
875  return 0;
876  }
877  }
878  return 1;
879 }
880 
881 /*
882  * return -1 if point outside polygon
883  * return 0 if point on boundary
884  * return 1 if point inside polygon
885  *
886  * Expected **root order is each exterior ring followed by its holes, eg. EIIEIIEI
887  */
888 int point_in_multipolygon_rtree(RTREE_NODE **root, int polyCount, int *ringCounts, LWPOINT *point)
889 {
890  int i, p, r, in_ring;
891  POINT2D pt;
892  int result = -1;
893 
894  POSTGIS_DEBUGF(2, "point_in_multipolygon_rtree called for %p %d %p.", root, polyCount, point);
895 
896  getPoint2d_p(point->point, 0, &pt);
897  /* assume bbox short-circuit has already been attempted */
898 
899  i = 0; /* the current index into the root array */
900 
901  /* is the point inside any of the sub-polygons? */
902  for ( p = 0; p < polyCount; p++ )
903  {
904  in_ring = point_in_ring_rtree(root[i], &pt);
905  POSTGIS_DEBUGF(4, "point_in_multipolygon_rtree: exterior ring (%d), point_in_ring returned %d", p, in_ring);
906  if ( in_ring == -1 ) /* outside the exterior ring */
907  {
908  POSTGIS_DEBUG(3, "point_in_multipolygon_rtree: outside exterior ring.");
909  }
910  else if ( in_ring == 0 ) /* on the boundary */
911  {
912  POSTGIS_DEBUGF(3, "point_in_multipolygon_rtree: on edge of exterior ring %d", p);
913  return 0;
914  } else {
915  result = in_ring;
916 
917  for(r=1; r<ringCounts[p]; r++)
918  {
919  in_ring = point_in_ring_rtree(root[i+r], &pt);
920  POSTGIS_DEBUGF(4, "point_in_multipolygon_rtree: interior ring (%d), point_in_ring returned %d", r, in_ring);
921  if (in_ring == 1) /* inside a hole => outside the polygon */
922  {
923  POSTGIS_DEBUGF(3, "point_in_multipolygon_rtree: within hole %d of exterior ring %d", r, p);
924  result = -1;
925  break;
926  }
927  if (in_ring == 0) /* on the edge of a hole */
928  {
929  POSTGIS_DEBUGF(3, "point_in_multipolygon_rtree: on edge of hole %d of exterior ring %d", r, p);
930  return 0;
931  }
932  }
933  /* if we have a positive result, we can short-circuit and return it */
934  if ( result != -1)
935  {
936  return result;
937  }
938  }
939  /* increment the index by the total number of rings in the sub-poly */
940  /* we do this here in case we short-cutted out of the poly before looking at all the rings */
941  i += ringCounts[p];
942  }
943 
944  return result; /* -1 = outside, 0 = boundary, 1 = inside */
945 
946 }
947 
948 /*
949  * return -1 iff point outside polygon
950  * return 0 iff point on boundary
951  * return 1 iff point inside polygon
952  */
953 int point_in_polygon(LWPOLY *polygon, LWPOINT *point)
954 {
955  uint32_t i;
956  int 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  uint32_t i, j;
1000  int result, in_ring;
1001  POINT2D pt;
1002 
1003  POSTGIS_DEBUG(2, "point_in_polygon called.");
1004 
1005  getPoint2d_p(point->point, 0, &pt);
1006  /* assume bbox short-circuit has already been attempted */
1007 
1008  result = -1;
1009 
1010  for (j = 0; j < mpolygon->ngeoms; j++ )
1011  {
1012 
1013  LWPOLY *polygon = mpolygon->geoms[j];
1014 
1015  /* everything is outside of an empty polygon */
1016  if ( polygon->nrings == 0 ) continue;
1017 
1018  in_ring = point_in_ring(polygon->rings[0], &pt);
1019  if ( in_ring == -1) /* outside the exterior ring */
1020  {
1021  POSTGIS_DEBUG(3, "point_in_polygon: outside exterior ring.");
1022  continue;
1023  }
1024  if ( in_ring == 0 )
1025  {
1026  return 0;
1027  }
1028 
1029  result = in_ring;
1030 
1031  for (i=1; i<polygon->nrings; i++)
1032  {
1033  in_ring = point_in_ring(polygon->rings[i], &pt);
1034  if (in_ring == 1) /* inside a hole => outside the polygon */
1035  {
1036  POSTGIS_DEBUGF(3, "point_in_polygon: within hole %d.", i);
1037  result = -1;
1038  break;
1039  }
1040  if (in_ring == 0) /* on the edge of a hole */
1041  {
1042  POSTGIS_DEBUGF(3, "point_in_polygon: on edge of hole %d.", i);
1043  return 0;
1044  }
1045  }
1046  if ( result != -1)
1047  {
1048  return result;
1049  }
1050  }
1051  return result;
1052 }
1053 
1054 
1055 /*******************************************************************************
1056  * End of "Fast Winding Number Inclusion of a Point in a Polygon" derivative.
1057  ******************************************************************************/
1058 
1059 /**********************************************************************
1060  *
1061  * ST_MinimumBoundingRadius
1062  *
1063  **********************************************************************/
1064 
1066 Datum ST_MinimumBoundingRadius(PG_FUNCTION_ARGS)
1067 {
1068  GSERIALIZED* geom;
1069  LWGEOM* input;
1070  LWBOUNDINGCIRCLE* mbc = NULL;
1071  LWGEOM* lwcenter;
1072  GSERIALIZED* center;
1073  TupleDesc resultTupleDesc;
1074  HeapTuple resultTuple;
1075  Datum result;
1076  Datum result_values[2];
1077  bool result_is_null[2];
1078  double radius = 0;
1079 
1080  if (PG_ARGISNULL(0))
1081  PG_RETURN_NULL();
1082 
1083  geom = PG_GETARG_GSERIALIZED_P(0);
1084 
1085  /* Empty geometry? Return POINT EMPTY with zero radius */
1086  if (gserialized_is_empty(geom))
1087  {
1089  }
1090  else
1091  {
1092  input = lwgeom_from_gserialized(geom);
1093  mbc = lwgeom_calculate_mbc(input);
1094 
1095  if (!(mbc && mbc->center))
1096  {
1097  lwpgerror("Error calculating minimum bounding circle.");
1098  lwgeom_free(input);
1099  PG_RETURN_NULL();
1100  }
1101 
1102  lwcenter = (LWGEOM*) lwpoint_make2d(input->srid, mbc->center->x, mbc->center->y);
1103  radius = mbc->radius;
1104 
1106  lwgeom_free(input);
1107  }
1108 
1109  center = geometry_serialize(lwcenter);
1110  lwgeom_free(lwcenter);
1111 
1112  get_call_result_type(fcinfo, NULL, &resultTupleDesc);
1113  BlessTupleDesc(resultTupleDesc);
1114 
1115  result_values[0] = PointerGetDatum(center);
1116  result_is_null[0] = false;
1117  result_values[1] = Float8GetDatum(radius);
1118  result_is_null[1] = false;
1119 
1120  resultTuple = heap_form_tuple(resultTupleDesc, result_values, result_is_null);
1121 
1122  result = HeapTupleGetDatum(resultTuple);
1123 
1124  PG_RETURN_DATUM(result);
1125 }
1126 
1127 /**********************************************************************
1128  *
1129  * ST_MinimumBoundingCircle
1130  *
1131  **********************************************************************/
1132 
1134 Datum ST_MinimumBoundingCircle(PG_FUNCTION_ARGS)
1135 {
1136  GSERIALIZED* geom;
1137  LWGEOM* input;
1138  LWBOUNDINGCIRCLE* mbc = NULL;
1139  LWGEOM* lwcircle;
1140  GSERIALIZED* center;
1141  int segs_per_quarter;
1142 
1143  if (PG_ARGISNULL(0))
1144  PG_RETURN_NULL();
1145 
1146  geom = PG_GETARG_GSERIALIZED_P(0);
1147  segs_per_quarter = PG_GETARG_INT32(1);
1148 
1149  /* Empty geometry? Return POINT EMPTY */
1150  if (gserialized_is_empty(geom))
1151  {
1153  }
1154  else
1155  {
1156  input = lwgeom_from_gserialized(geom);
1157  mbc = lwgeom_calculate_mbc(input);
1158 
1159  if (!(mbc && mbc->center))
1160  {
1161  lwpgerror("Error calculating minimum bounding circle.");
1162  lwgeom_free(input);
1163  PG_RETURN_NULL();
1164  }
1165 
1166  /* Zero radius? Return a point. */
1167  if (mbc->radius == 0)
1168  lwcircle = lwpoint_as_lwgeom(lwpoint_make2d(input->srid, mbc->center->x, mbc->center->y));
1169  else
1170  lwcircle = lwpoly_as_lwgeom(lwpoly_construct_circle(input->srid, mbc->center->x, mbc->center->y, mbc->radius, segs_per_quarter, LW_TRUE));
1171 
1173  lwgeom_free(input);
1174  }
1175 
1176  center = geometry_serialize(lwcircle);
1177  lwgeom_free(lwcircle);
1178 
1179  PG_RETURN_POINTER(center);
1180 }
1181 
1182 /**********************************************************************
1183  *
1184  * ST_GeometricMedian
1185  *
1186  **********************************************************************/
1187 
1189 Datum ST_GeometricMedian(PG_FUNCTION_ARGS)
1190 {
1191  GSERIALIZED* geom;
1192  GSERIALIZED* result;
1193  LWGEOM* input;
1194  LWPOINT* lwresult;
1195  static const double min_default_tolerance = 1e-8;
1196  double tolerance = min_default_tolerance;
1197  bool compute_tolerance_from_box;
1198  bool fail_if_not_converged;
1199  int max_iter;
1200 
1201  /* Read and validate our input arguments */
1202  if (PG_ARGISNULL(0))
1203  PG_RETURN_NULL();
1204 
1205  compute_tolerance_from_box = PG_ARGISNULL(1);
1206 
1207  if (!compute_tolerance_from_box)
1208  {
1209  tolerance = PG_GETARG_FLOAT8(1);
1210  if (tolerance < 0)
1211  {
1212  lwpgerror("Tolerance must be positive.");
1213  PG_RETURN_NULL();
1214  }
1215  }
1216 
1217  max_iter = PG_ARGISNULL(2) ? -1 : PG_GETARG_INT32(2);
1218  fail_if_not_converged = PG_ARGISNULL(3) ? LW_FALSE : PG_GETARG_BOOL(3);
1219 
1220  if (max_iter < 0)
1221  {
1222  lwpgerror("Maximum iterations must be positive.");
1223  PG_RETURN_NULL();
1224  }
1225 
1226  /* OK, inputs are valid. */
1227  geom = PG_GETARG_GSERIALIZED_P(0);
1228  input = lwgeom_from_gserialized(geom);
1229 
1230  if (compute_tolerance_from_box)
1231  {
1232  /* Compute a default tolerance based on the smallest dimension
1233  * of the geometry's bounding box.
1234  */
1235  static const double tolerance_coefficient = 1e-6;
1236  const GBOX* box = lwgeom_get_bbox(input);
1237 
1238  if (box)
1239  {
1240  double min_dim = FP_MIN(box->xmax - box->xmin, box->ymax - box->ymin);
1241  if (lwgeom_has_z(input))
1242  min_dim = FP_MIN(min_dim, box->zmax - box->zmin);
1243 
1244  /* Apply a lower bound to the computed default tolerance to
1245  * avoid a tolerance of zero in the case of collinear
1246  * points.
1247  */
1248  tolerance = FP_MAX(min_default_tolerance, tolerance_coefficient * min_dim);
1249  }
1250  }
1251 
1252  lwresult = lwgeom_median(input, tolerance, max_iter, fail_if_not_converged);
1253  lwgeom_free(input);
1254 
1255  if(!lwresult)
1256  {
1257  lwpgerror("Error computing geometric median.");
1258  PG_RETURN_NULL();
1259  }
1260 
1261  result = geometry_serialize(lwpoint_as_lwgeom(lwresult));
1262 
1263  PG_RETURN_POINTER(result);
1264 }
1265 
1266 /**********************************************************************
1267  *
1268  * ST_IsPolygonCW
1269  *
1270  **********************************************************************/
1271 
1273 Datum ST_IsPolygonCW(PG_FUNCTION_ARGS)
1274 {
1275  GSERIALIZED* geom;
1276  LWGEOM* input;
1277  bool is_clockwise;
1278 
1279  if (PG_ARGISNULL(0))
1280  PG_RETURN_NULL();
1281 
1282  geom = PG_GETARG_GSERIALIZED_P(0);
1283  input = lwgeom_from_gserialized(geom);
1284 
1285  is_clockwise = lwgeom_is_clockwise(input);
1286 
1287  lwgeom_free(input);
1288  PG_FREE_IF_COPY(geom, 0);
1289 
1290  PG_RETURN_BOOL(is_clockwise);
1291 }
1292 
1293 /**********************************************************************
1294  *
1295  * ST_IsPolygonCCW
1296  *
1297  **********************************************************************/
1298 
1300 Datum ST_IsPolygonCCW(PG_FUNCTION_ARGS)
1301 {
1302  GSERIALIZED* geom;
1303  LWGEOM* input;
1304  bool is_ccw;
1305 
1306  if (PG_ARGISNULL(0))
1307  PG_RETURN_NULL();
1308 
1309  geom = PG_GETARG_GSERIALIZED_P_COPY(0);
1310  input = lwgeom_from_gserialized(geom);
1311  lwgeom_reverse_in_place(input);
1312  is_ccw = lwgeom_is_clockwise(input);
1313  lwgeom_free(input);
1314  PG_FREE_IF_COPY(geom, 0);
1315 
1316  PG_RETURN_BOOL(is_ccw);
1317 }
1318 
double x
Definition: liblwgeom.h:354
#define LINETYPE
Definition: liblwgeom.h:85
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:86
static int point_in_ring_rtree(RTREE_NODE *root, const POINT2D *point)
Datum LWGEOM_simplify2d(PG_FUNCTION_ARGS)
GBOX * bbox
Definition: liblwgeom.h:400
static int point_in_ring(POINTARRAY *pts, const POINT2D *point)
double m
Definition: liblwgeom.h:354
LWGEOM * lwgeom_grid(const LWGEOM *lwgeom, const gridspec *grid)
Definition: lwgeom.c:2243
LWCOLLECTION * lwcollection_construct(uint8_t type, int srid, GBOX *bbox, uint32_t ngeoms, LWGEOM **geoms)
Definition: lwcollection.c:43
char * r
Definition: cu_in_wkt.c:24
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
int lwline_crossing_direction(const LWLINE *l1, const LWLINE *l2)
Given two lines, characterize how (and if) they cross each other.
Definition: lwalgorithm.c:461
void lwgeom_reverse_in_place(LWGEOM *lwgeom)
Reverse vertex order of LWGEOM.
Definition: lwgeom.c:102
uint32_t ngeoms
Definition: liblwgeom.h:496
Datum ST_IsPolygonCW(PG_FUNCTION_ARGS)
The following struct and methods are used for a 1D RTree implementation, described at: http://lin-ear...
Definition: lwgeom_rtree.h:46
LWPOINT * lwpoint_make2d(int srid, double x, double y)
Definition: lwpoint.c:163
Datum area(PG_FUNCTION_ARGS)
double xmax
Definition: liblwgeom.h:295
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1144
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:1053
LWMPOINT * lwmpoint_construct(int srid, const POINTARRAY *pa)
Definition: lwmpoint.c:52
#define MULTIPOINTTYPE
Definition: liblwgeom.h:87
double ptarray_length_2d(const POINTARRAY *pts)
Find the 2d length of the given POINTARRAY (even if it&#39;s 3d)
Definition: ptarray.c:1682
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:450
Datum ST_GeometricMedian(PG_FUNCTION_ARGS)
void error_if_srid_mismatch(int srid1, int srid2)
Definition: lwutil.c:338
LWPOINT * lwpoint_construct_empty(int srid, char hasz, char hasm)
Definition: lwpoint.c:151
LWPOINT * lwgeom_as_lwpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:161
Datum LWGEOM_SetEffectiveArea(PG_FUNCTION_ARGS)
#define FP_MIN(A, B)
POINTARRAY * point
Definition: liblwgeom.h:413
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:930
Datum ST_IsPolygonCCW(PG_FUNCTION_ARGS)
int32_t srid
Definition: liblwgeom.h:401
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:320
uint32_t nrings
Definition: liblwgeom.h:457
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:178
Datum LWGEOM_line_substring(PG_FUNCTION_ARGS)
unsigned int uint32_t
Definition: uthash.h:78
double x
Definition: liblwgeom.h:330
Datum LWGEOM_snaptogrid_pointoff(PG_FUNCTION_ARGS)
LWBOUNDINGCIRCLE * lwgeom_calculate_mbc(const LWGEOM *g)
double zmax
Definition: liblwgeom.h:299
double ymin
Definition: liblwgeom.h:296
POINTARRAY * lwline_interpolate_points(const LWLINE *line, double length_fraction, char repeat)
Interpolate one or more points along a line.
Definition: lwline.c:542
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:218
LWGEOM * geom
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:330
double xmin
Definition: liblwgeom.h:294
#define LW_FALSE
Definition: liblwgeom.h:76
LWPOLY ** geoms
Definition: liblwgeom.h:498
static int is_clockwise(int num_points, double *x, double *y, double *z)
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:42
LWGEOM * lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed)
Definition: lwgeom.c:1842
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:187
LWMLINE * lwgeom_as_lwmline(const LWGEOM *lwgeom)
Definition: lwgeom.c:242
const GBOX * lwgeom_get_bbox(const LWGEOM *lwgeom)
Get a non-empty geometry bounding box, computing and caching it if not already there.
Definition: lwgeom.c:734
POINTARRAY ** rings
Definition: liblwgeom.h:459
double ymax
Definition: liblwgeom.h:297
double y
Definition: liblwgeom.h:330
#define FLAGS_GET_Z(flags)
Macros for manipulating the &#39;flags&#39; byte.
Definition: liblwgeom.h:139
double z
Definition: liblwgeom.h:354
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:170
int getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *point)
Definition: lwgeom_api.c:113
POINT2D * center
Definition: liblwgeom.h:1649
uint8_t flags
Definition: liblwgeom.h:410
uint32_t ngeoms
Definition: liblwgeom.h:483
LWLINE ** geoms
Definition: liblwgeom.h:485
int getPoint2d_p(const POINTARRAY *pa, uint32_t n, POINT2D *point)
Definition: lwgeom_api.c:338
int point_in_polygon(LWPOLY *polygon, LWPOINT *point)
static int isOnSegment(const POINT2D *seg1, const POINT2D *seg2, const POINT2D *point)
LWGEOM * lwgeom_chaikin(const LWGEOM *igeom, int n_iterations, int preserve_endpoint)
Definition: lwchaikins.c:182
GSERIALIZED * geometry_serialize(LWGEOM *lwgeom)
Datum LWGEOM_line_interpolate_point(PG_FUNCTION_ARGS)
double zmin
Definition: liblwgeom.h:298
Datum ST_LineCrossingDirection(PG_FUNCTION_ARGS)
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
void lwgeom_add_bbox(LWGEOM *lwgeom)
Compute a bbox if not already computed.
Definition: lwgeom.c:686
#define FLAGS_GET_M(flags)
Definition: liblwgeom.h:140
uint8_t type
Definition: liblwgeom.h:398
type
Definition: ovdump.py:41
static double determineSide(const POINT2D *seg1, const POINT2D *seg2, const POINT2D *point)
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:335
Datum LWGEOM_ChaikinSmoothing(PG_FUNCTION_ARGS)
LWPOLY * lwpoly_construct_circle(int srid, double x, double y, double radius, uint32_t segments_per_quarter, char exterior)
Definition: lwpoly.c:120
int point_in_polygon_rtree(RTREE_NODE **root, int ringCount, LWPOINT *point)
int lwgeom_is_clockwise(LWGEOM *lwgeom)
Ensure the outer ring is clockwise oriented and all inner rings are counter-clockwise.
Definition: lwgeom.c:65
LWPOINT * lwgeom_median(const LWGEOM *g, double tol, uint32_t maxiter, char fail_if_not_converged)
LWPOINT * lwpoint_construct(int srid, GBOX *bbox, POINTARRAY *point)
Definition: lwpoint.c:129
void * lwalloc(size_t size)
Definition: lwutil.c:229
Datum ST_MinimumBoundingCircle(PG_FUNCTION_ARGS)
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:1393
void lwmline_release(LWMLINE *lwline)
Definition: lwmline.c:32
#define FP_CONTAINS_BOTTOM(A, X, B)
PG_FUNCTION_INFO_V1(LWGEOM_simplify2d)
double y
Definition: liblwgeom.h:354
#define MULTILINETYPE
Definition: liblwgeom.h:88
void lwline_release(LWLINE *lwline)
Definition: lwline.c:134
Datum ST_MinimumBoundingRadius(PG_FUNCTION_ARGS)
void lwboundingcircle_destroy(LWBOUNDINGCIRCLE *c)
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:99
if(!(yy_init))
#define COLLECTIONTYPE
Definition: liblwgeom.h:90
This library is the generic geometry handling section of PostGIS.
#define FP_MAX(A, B)
const POINT2D * getPoint2d_cp(const POINTARRAY *pa, uint32_t n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:364
POINTARRAY * points
Definition: liblwgeom.h:424
LWGEOM * lwmpoint_as_lwgeom(const LWMPOINT *obj)
Definition: lwgeom.c:295
Snap to grid.
uint32_t npoints
Definition: liblwgeom.h:373