PostGIS  2.5.0dev-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 
130 
131 /***********************************************************************
132  * --strk@kbt.io;
133  ***********************************************************************/
134 
135 /***********************************************************************
136  * Interpolate a point along a line, useful for Geocoding applications
137  * SELECT line_interpolate_point( 'LINESTRING( 0 0, 2 2'::geometry, .5 )
138  * returns POINT( 1 1 ).
139  ***********************************************************************/
141 Datum LWGEOM_line_interpolate_point(PG_FUNCTION_ARGS)
142 {
143  GSERIALIZED *gser = PG_GETARG_GSERIALIZED_P(0);
144  GSERIALIZED *result;
145  double distance_fraction = PG_GETARG_FLOAT8(1);
146  int repeat = PG_NARGS() > 2 && PG_GETARG_BOOL(2);
147  int srid = gserialized_get_srid(gser);
148  LWLINE* lwline;
149  LWGEOM* lwresult;
150  POINTARRAY* opa;
151 
152  if ( distance_fraction < 0 || distance_fraction > 1 )
153  {
154  elog(ERROR,"line_interpolate_point: 2nd arg isn't within [0,1]");
155  PG_FREE_IF_COPY(gser, 0);
156  PG_RETURN_NULL();
157  }
158 
159  if ( gserialized_get_type(gser) != LINETYPE )
160  {
161  elog(ERROR,"line_interpolate_point: 1st arg isn't a line");
162  PG_FREE_IF_COPY(gser, 0);
163  PG_RETURN_NULL();
164  }
165 
167  opa = lwline_interpolate_points(lwline, distance_fraction, repeat);
168 
169  lwgeom_free(lwline_as_lwgeom(lwline));
170  PG_FREE_IF_COPY(gser, 0);
171 
172  if (opa->npoints <= 1)
173  {
174  lwresult = lwpoint_as_lwgeom(lwpoint_construct(srid, NULL, opa));
175  } else {
176  lwresult = lwmpoint_as_lwgeom(lwmpoint_construct(srid, opa));
177  }
178 
179  result = geometry_serialize(lwresult);
180  lwgeom_free(lwresult);
181 
182  PG_RETURN_POINTER(result);
183 }
184 /***********************************************************************
185  * --jsunday@rochgrp.com;
186  ***********************************************************************/
187 
188 /***********************************************************************
189  *
190  * Grid application function for postgis.
191  *
192  * WHAT IS
193  * -------
194  *
195  * This is a grid application function for postgis.
196  * You use it to stick all of a geometry points to
197  * a custom grid defined by its origin and cell size
198  * in geometry units.
199  *
200  * Points reduction is obtained by collapsing all
201  * consecutive points falling on the same grid node and
202  * removing all consecutive segments S1,S2 having
203  * S2.startpoint = S1.endpoint and S2.endpoint = S1.startpoint.
204  *
205  * ISSUES
206  * ------
207  *
208  * Only 2D is contemplated in grid application.
209  *
210  * Consecutive coincident segments removal does not work
211  * on first/last segments (They are not considered consecutive).
212  *
213  * Grid application occurs on a geometry subobject basis, thus no
214  * points reduction occurs for multipoint geometries.
215  *
216  * USAGE TIPS
217  * ----------
218  *
219  * Run like this:
220  *
221  * SELECT SnapToGrid(<geometry>, <originX>, <originY>, <sizeX>, <sizeY>);
222  *
223  * Grid cells are not visible on a screen as long as the screen
224  * point size is equal or bigger then the grid cell size.
225  * This assumption may be used to reduce the number of points in
226  * a map for a given display scale.
227  *
228  * Keeping multiple resolution versions of the same data may be used
229  * in conjunction with MINSCALE/MAXSCALE keywords of mapserv to speed
230  * up rendering.
231  *
232  * Check also the use of DP simplification to reduce grid visibility.
233  * I'm still researching about the relationship between grid size and
234  * DP epsilon values - please tell me if you know more about this.
235  *
236  *
237  * --strk@kbt.io;
238  *
239  ***********************************************************************/
240 
241 #define CHECK_RING_IS_CLOSE
242 #define SAMEPOINT(a,b) ((a)->x==(b)->x&&(a)->y==(b)->y)
243 
244 
245 
246 /* Forward declarations */
247 Datum LWGEOM_snaptogrid(PG_FUNCTION_ARGS);
248 Datum LWGEOM_snaptogrid_pointoff(PG_FUNCTION_ARGS);
249 
250 
251 
253 Datum LWGEOM_snaptogrid(PG_FUNCTION_ARGS)
254 {
255  LWGEOM *in_lwgeom;
256  GSERIALIZED *out_geom = NULL;
257  LWGEOM *out_lwgeom;
258  gridspec grid;
259 
260  GSERIALIZED *in_geom = PG_GETARG_GSERIALIZED_P(0);
261 
262  /* Set grid values to zero to start */
263  memset(&grid, 0, sizeof(gridspec));
264 
265  grid.ipx = PG_GETARG_FLOAT8(1);
266  grid.ipy = PG_GETARG_FLOAT8(2);
267  grid.xsize = PG_GETARG_FLOAT8(3);
268  grid.ysize = PG_GETARG_FLOAT8(4);
269 
270  /* Return input geometry if input geometry is empty */
271  if ( gserialized_is_empty(in_geom) )
272  {
273  PG_RETURN_POINTER(in_geom);
274  }
275 
276  /* Return input geometry if input grid is meaningless */
277  if ( grid.xsize==0 && grid.ysize==0 && grid.zsize==0 && grid.msize==0 )
278  {
279  PG_RETURN_POINTER(in_geom);
280  }
281 
282  in_lwgeom = lwgeom_from_gserialized(in_geom);
283 
284  POSTGIS_DEBUGF(3, "SnapToGrid got a %s", lwtype_name(in_lwgeom->type));
285 
286  out_lwgeom = lwgeom_grid(in_lwgeom, &grid);
287  if ( out_lwgeom == NULL ) PG_RETURN_NULL();
288 
289  /* COMPUTE_BBOX TAINTING */
290  if ( in_lwgeom->bbox )
291  lwgeom_add_bbox(out_lwgeom);
292 
293 
294  POSTGIS_DEBUGF(3, "SnapToGrid made a %s", lwtype_name(out_lwgeom->type));
295 
296  out_geom = geometry_serialize(out_lwgeom);
297 
298  PG_RETURN_POINTER(out_geom);
299 }
300 
301 
302 #if POSTGIS_DEBUG_LEVEL >= 4
303 /* Print grid using given reporter */
304 static void
305 grid_print(const gridspec *grid)
306 {
307  lwpgnotice("GRID(%g %g %g %g, %g %g %g %g)",
308  grid->ipx, grid->ipy, grid->ipz, grid->ipm,
309  grid->xsize, grid->ysize, grid->zsize, grid->msize);
310 }
311 #endif
312 
313 
315 Datum LWGEOM_snaptogrid_pointoff(PG_FUNCTION_ARGS)
316 {
317  GSERIALIZED *in_geom, *in_point;
318  LWGEOM *in_lwgeom;
319  LWPOINT *in_lwpoint;
320  GSERIALIZED *out_geom = NULL;
321  LWGEOM *out_lwgeom;
322  gridspec grid;
323  /* BOX3D box3d; */
324  POINT4D offsetpoint;
325 
326  in_geom = PG_GETARG_GSERIALIZED_P(0);
327 
328  /* Return input geometry if input geometry is empty */
329  if ( gserialized_is_empty(in_geom) )
330  {
331  PG_RETURN_POINTER(in_geom);
332  }
333 
334  in_point = PG_GETARG_GSERIALIZED_P(1);
335  in_lwpoint = lwgeom_as_lwpoint(lwgeom_from_gserialized(in_point));
336  if ( in_lwpoint == NULL )
337  {
338  lwpgerror("Offset geometry must be a point");
339  }
340 
341  grid.xsize = PG_GETARG_FLOAT8(2);
342  grid.ysize = PG_GETARG_FLOAT8(3);
343  grid.zsize = PG_GETARG_FLOAT8(4);
344  grid.msize = PG_GETARG_FLOAT8(5);
345 
346  /* Take offsets from point geometry */
347  getPoint4d_p(in_lwpoint->point, 0, &offsetpoint);
348  grid.ipx = offsetpoint.x;
349  grid.ipy = offsetpoint.y;
350  if (FLAGS_GET_Z(in_lwpoint->flags) ) grid.ipz = offsetpoint.z;
351  else grid.ipz=0;
352  if (FLAGS_GET_M(in_lwpoint->flags) ) grid.ipm = offsetpoint.m;
353  else grid.ipm=0;
354 
355 #if POSTGIS_DEBUG_LEVEL >= 4
356  grid_print(&grid);
357 #endif
358 
359  /* Return input geometry if input grid is meaningless */
360  if ( grid.xsize==0 && grid.ysize==0 && grid.zsize==0 && grid.msize==0 )
361  {
362  PG_RETURN_POINTER(in_geom);
363  }
364 
365  in_lwgeom = lwgeom_from_gserialized(in_geom);
366 
367  POSTGIS_DEBUGF(3, "SnapToGrid got a %s", lwtype_name(in_lwgeom->type));
368 
369  out_lwgeom = lwgeom_grid(in_lwgeom, &grid);
370  if ( out_lwgeom == NULL ) PG_RETURN_NULL();
371 
372  /* COMPUTE_BBOX TAINTING */
373  if ( in_lwgeom->bbox ) lwgeom_add_bbox(out_lwgeom);
374 
375  POSTGIS_DEBUGF(3, "SnapToGrid made a %s", lwtype_name(out_lwgeom->type));
376 
377  out_geom = geometry_serialize(out_lwgeom);
378 
379  PG_RETURN_POINTER(out_geom);
380 }
381 
382 
383 /*
384 ** crossingDirection(line1, line2)
385 **
386 ** Determines crossing direction of line2 relative to line1.
387 ** Only accepts LINESTRING as parameters!
388 */
390 Datum ST_LineCrossingDirection(PG_FUNCTION_ARGS)
391 {
392  int type1, type2, rv;
393  LWLINE *l1 = NULL;
394  LWLINE *l2 = NULL;
395  GSERIALIZED *geom1 = PG_GETARG_GSERIALIZED_P(0);
396  GSERIALIZED *geom2 = PG_GETARG_GSERIALIZED_P(1);
397 
399 
400  type1 = gserialized_get_type(geom1);
401  type2 = gserialized_get_type(geom2);
402 
403  if ( type1 != LINETYPE || type2 != LINETYPE )
404  {
405  elog(ERROR,"This function only accepts LINESTRING as arguments.");
406  PG_RETURN_NULL();
407  }
408 
411 
412  rv = lwline_crossing_direction(l1, l2);
413 
414  PG_FREE_IF_COPY(geom1, 0);
415  PG_FREE_IF_COPY(geom2, 1);
416 
417  PG_RETURN_INT32(rv);
418 
419 }
420 
421 
422 
423 /***********************************************************************
424  * --strk@kbt.io
425  ***********************************************************************/
426 
427 Datum LWGEOM_line_substring(PG_FUNCTION_ARGS);
428 
430 Datum LWGEOM_line_substring(PG_FUNCTION_ARGS)
431 {
432  GSERIALIZED *geom = PG_GETARG_GSERIALIZED_P(0);
433  double from = PG_GETARG_FLOAT8(1);
434  double to = PG_GETARG_FLOAT8(2);
435  LWGEOM *olwgeom;
436  POINTARRAY *ipa, *opa;
437  GSERIALIZED *ret;
438  int type = gserialized_get_type(geom);
439 
440  if ( from < 0 || from > 1 )
441  {
442  elog(ERROR,"line_interpolate_point: 2nd arg isn't within [0,1]");
443  PG_RETURN_NULL();
444  }
445 
446  if ( to < 0 || to > 1 )
447  {
448  elog(ERROR,"line_interpolate_point: 3rd arg isn't within [0,1]");
449  PG_RETURN_NULL();
450  }
451 
452  if ( from > to )
453  {
454  elog(ERROR, "2nd arg must be smaller then 3rd arg");
455  PG_RETURN_NULL();
456  }
457 
458  if ( type == LINETYPE )
459  {
461 
462  if ( lwgeom_is_empty((LWGEOM*)iline) )
463  {
464  /* TODO return empty line */
465  lwline_release(iline);
466  PG_FREE_IF_COPY(geom, 0);
467  PG_RETURN_NULL();
468  }
469 
470  ipa = iline->points;
471 
472  opa = ptarray_substring(ipa, from, to, 0);
473 
474  if ( opa->npoints == 1 ) /* Point returned */
475  olwgeom = (LWGEOM *)lwpoint_construct(iline->srid, NULL, opa);
476  else
477  olwgeom = (LWGEOM *)lwline_construct(iline->srid, NULL, opa);
478 
479  }
480  else if ( type == MULTILINETYPE )
481  {
482  LWMLINE *iline;
483  uint32_t i = 0, g = 0;
484  int homogeneous = LW_TRUE;
485  LWGEOM **geoms = NULL;
486  double length = 0.0, sublength = 0.0, minprop = 0.0, maxprop = 0.0;
487 
489 
490  if ( lwgeom_is_empty((LWGEOM*)iline) )
491  {
492  /* TODO return empty collection */
493  lwmline_release(iline);
494  PG_FREE_IF_COPY(geom, 0);
495  PG_RETURN_NULL();
496  }
497 
498  /* Calculate the total length of the mline */
499  for ( i = 0; i < iline->ngeoms; i++ )
500  {
501  LWLINE *subline = (LWLINE*)iline->geoms[i];
502  if ( subline->points && subline->points->npoints > 1 )
503  length += ptarray_length_2d(subline->points);
504  }
505 
506  geoms = lwalloc(sizeof(LWGEOM*) * iline->ngeoms);
507 
508  /* Slice each sub-geometry of the multiline */
509  for ( i = 0; i < iline->ngeoms; i++ )
510  {
511  LWLINE *subline = (LWLINE*)iline->geoms[i];
512  double subfrom = 0.0, subto = 0.0;
513 
514  if ( subline->points && subline->points->npoints > 1 )
515  sublength += ptarray_length_2d(subline->points);
516 
517  /* Calculate proportions for this subline */
518  minprop = maxprop;
519  maxprop = sublength / length;
520 
521  /* This subline doesn't reach the lowest proportion requested
522  or is beyond the highest proporton */
523  if ( from > maxprop || to < minprop )
524  continue;
525 
526  if ( from <= minprop )
527  subfrom = 0.0;
528  if ( to >= maxprop )
529  subto = 1.0;
530 
531  if ( from > minprop && from <= maxprop )
532  subfrom = (from - minprop) / (maxprop - minprop);
533 
534  if ( to < maxprop && to >= minprop )
535  subto = (to - minprop) / (maxprop - minprop);
536 
537 
538  opa = ptarray_substring(subline->points, subfrom, subto, 0);
539  if ( opa && opa->npoints > 0 )
540  {
541  if ( opa->npoints == 1 ) /* Point returned */
542  {
543  geoms[g] = (LWGEOM *)lwpoint_construct(SRID_UNKNOWN, NULL, opa);
544  homogeneous = LW_FALSE;
545  }
546  else
547  {
548  geoms[g] = (LWGEOM *)lwline_construct(SRID_UNKNOWN, NULL, opa);
549  }
550  g++;
551  }
552 
553 
554 
555  }
556  /* If we got any points, we need to return a GEOMETRYCOLLECTION */
557  if ( ! homogeneous )
558  type = COLLECTIONTYPE;
559 
560  olwgeom = (LWGEOM*)lwcollection_construct(type, iline->srid, NULL, g, geoms);
561  }
562  else
563  {
564  elog(ERROR,"line_substring: 1st arg isn't a line");
565  PG_RETURN_NULL();
566  }
567 
568  ret = geometry_serialize(olwgeom);
569  lwgeom_free(olwgeom);
570  PG_FREE_IF_COPY(geom, 0);
571  PG_RETURN_POINTER(ret);
572 
573 }
574 
575 /*******************************************************************************
576  * The following is based on the "Fast Winding Number Inclusion of a Point
577  * in a Polygon" algorithm by Dan Sunday.
578  * http://softsurfer.com/Archive/algorithm_0103/algorithm_0103.htm#Winding%20Number
579  ******************************************************************************/
580 
581 /*
582  * returns: >0 for a point to the left of the segment,
583  * <0 for a point to the right of the segment,
584  * 0 for a point on the segment
585  */
586 static double determineSide(const POINT2D *seg1, const POINT2D *seg2, const POINT2D *point)
587 {
588  return ((seg2->x-seg1->x)*(point->y-seg1->y)-(point->x-seg1->x)*(seg2->y-seg1->y));
589 }
590 
591 /*
592  * This function doesn't test that the point falls on the line defined by
593  * the two points. It assumes that that has already been determined
594  * by having determineSide return within the tolerance. It simply checks
595  * that if the point is on the line, it is within the endpoints.
596  *
597  * returns: 1 if the point is not outside the bounds of the segment
598  * 0 if it is
599  */
600 static int isOnSegment(const POINT2D *seg1, const POINT2D *seg2, const POINT2D *point)
601 {
602  double maxX;
603  double maxY;
604  double minX;
605  double minY;
606 
607  if (seg1->x > seg2->x)
608  {
609  maxX = seg1->x;
610  minX = seg2->x;
611  }
612  else
613  {
614  maxX = seg2->x;
615  minX = seg1->x;
616  }
617  if (seg1->y > seg2->y)
618  {
619  maxY = seg1->y;
620  minY = seg2->y;
621  }
622  else
623  {
624  maxY = seg2->y;
625  minY = seg1->y;
626  }
627 
628  POSTGIS_DEBUGF(3, "maxX minX/maxY minY: %.8f %.8f/%.8f %.8f", maxX, minX, maxY, minY);
629 
630  if (maxX < point->x || minX > point->x)
631  {
632  POSTGIS_DEBUGF(3, "X value %.8f falls outside the range %.8f-%.8f", point->x, minX, maxX);
633 
634  return 0;
635  }
636  else if (maxY < point->y || minY > point->y)
637  {
638  POSTGIS_DEBUGF(3, "Y value %.8f falls outside the range %.8f-%.8f", point->y, minY, maxY);
639 
640  return 0;
641  }
642  return 1;
643 }
644 
645 /*
646  * return -1 iff point is outside ring pts
647  * return 1 iff point is inside ring pts
648  * return 0 iff point is on ring pts
649  */
650 static int point_in_ring_rtree(RTREE_NODE *root, const POINT2D *point)
651 {
652  int wn = 0;
653  uint32_t i;
654  double side;
655  const POINT2D *seg1;
656  const POINT2D *seg2;
657  LWMLINE *lines;
658 
659  POSTGIS_DEBUG(2, "point_in_ring called.");
660 
661  lines = RTreeFindLineSegments(root, point->y);
662  if (!lines)
663  return -1;
664 
665  for (i=0; i<lines->ngeoms; i++)
666  {
667  seg1 = getPoint2d_cp(lines->geoms[i]->points, 0);
668  seg2 = getPoint2d_cp(lines->geoms[i]->points, 1);
669 
670  side = determineSide(seg1, seg2, point);
671 
672  POSTGIS_DEBUGF(3, "segment: (%.8f, %.8f),(%.8f, %.8f)", seg1->x, seg1->y, seg2->x, seg2->y);
673  POSTGIS_DEBUGF(3, "side result: %.8f", side);
674  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));
675 
676  /* zero length segments are ignored. */
677  if (((seg2->x - seg1->x)*(seg2->x - seg1->x) + (seg2->y - seg1->y)*(seg2->y - seg1->y)) < 1e-12*1e-12)
678  {
679  POSTGIS_DEBUG(3, "segment is zero length... ignoring.");
680 
681  continue;
682  }
683 
684  /* a point on the boundary of a ring is not contained. */
685  /* WAS: if (fabs(side) < 1e-12), see #852 */
686  if (side == 0.0)
687  {
688  if (isOnSegment(seg1, seg2, point) == 1)
689  {
690  POSTGIS_DEBUGF(3, "point on ring boundary between points %d, %d", i, i+1);
691 
692  return 0;
693  }
694  }
695 
696  /*
697  * If the point is to the left of the line, and it's rising,
698  * then the line is to the right of the point and
699  * circling counter-clockwise, so incremement.
700  */
701  if (FP_CONTAINS_BOTTOM(seg1->y, point->y, seg2->y) && side>0)
702  {
703  POSTGIS_DEBUG(3, "incrementing winding number.");
704 
705  ++wn;
706  }
707  /*
708  * If the point is to the right of the line, and it's falling,
709  * then the line is to the right of the point and circling
710  * clockwise, so decrement.
711  */
712  else if (FP_CONTAINS_BOTTOM(seg2->y, point->y, seg1->y) && side<0)
713  {
714  POSTGIS_DEBUG(3, "decrementing winding number.");
715 
716  --wn;
717  }
718  }
719 
720  POSTGIS_DEBUGF(3, "winding number %d", wn);
721 
722  if (wn == 0)
723  return -1;
724  return 1;
725 }
726 
727 
728 /*
729  * return -1 iff point is outside ring pts
730  * return 1 iff point is inside ring pts
731  * return 0 iff point is on ring pts
732  */
733 static int point_in_ring(POINTARRAY *pts, const POINT2D *point)
734 {
735  int wn = 0;
736  uint32_t i;
737  double side;
738  const POINT2D* seg1;
739  const POINT2D* seg2;
740 
741  POSTGIS_DEBUG(2, "point_in_ring called.");
742 
743  seg2 = getPoint2d_cp(pts, 0);
744  for (i=0; i<pts->npoints-1; i++)
745  {
746  seg1 = seg2;
747  seg2 = getPoint2d_cp(pts, i+1);
748 
749  side = determineSide(seg1, seg2, point);
750 
751  POSTGIS_DEBUGF(3, "segment: (%.8f, %.8f),(%.8f, %.8f)", seg1->x, seg1->y, seg2->x, seg2->y);
752  POSTGIS_DEBUGF(3, "side result: %.8f", side);
753  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));
754 
755  /* zero length segments are ignored. */
756  if (((seg2->x - seg1->x)*(seg2->x - seg1->x) + (seg2->y - seg1->y)*(seg2->y - seg1->y)) < 1e-12*1e-12)
757  {
758  POSTGIS_DEBUG(3, "segment is zero length... ignoring.");
759 
760  continue;
761  }
762 
763  /* a point on the boundary of a ring is not contained. */
764  /* WAS: if (fabs(side) < 1e-12), see #852 */
765  if (side == 0.0)
766  {
767  if (isOnSegment(seg1, seg2, point) == 1)
768  {
769  POSTGIS_DEBUGF(3, "point on ring boundary between points %d, %d", i, i+1);
770 
771  return 0;
772  }
773  }
774 
775  /*
776  * If the point is to the left of the line, and it's rising,
777  * then the line is to the right of the point and
778  * circling counter-clockwise, so incremement.
779  */
780  if (FP_CONTAINS_BOTTOM(seg1->y, point->y, seg2->y) && side>0)
781  {
782  POSTGIS_DEBUG(3, "incrementing winding number.");
783 
784  ++wn;
785  }
786  /*
787  * If the point is to the right of the line, and it's falling,
788  * then the line is to the right of the point and circling
789  * clockwise, so decrement.
790  */
791  else if (FP_CONTAINS_BOTTOM(seg2->y, point->y, seg1->y) && side<0)
792  {
793  POSTGIS_DEBUG(3, "decrementing winding number.");
794 
795  --wn;
796  }
797  }
798 
799  POSTGIS_DEBUGF(3, "winding number %d", wn);
800 
801  if (wn == 0)
802  return -1;
803  return 1;
804 }
805 
806 /*
807  * return 0 iff point outside polygon or on boundary
808  * return 1 iff point inside polygon
809  */
810 int point_in_polygon_rtree(RTREE_NODE **root, int ringCount, LWPOINT *point)
811 {
812  int i;
813  POINT2D pt;
814 
815  POSTGIS_DEBUGF(2, "point_in_polygon called for %p %d %p.", root, ringCount, point);
816 
817  getPoint2d_p(point->point, 0, &pt);
818  /* assume bbox short-circuit has already been attempted */
819 
820  if (point_in_ring_rtree(root[0], &pt) != 1)
821  {
822  POSTGIS_DEBUG(3, "point_in_polygon_rtree: outside exterior ring.");
823 
824  return 0;
825  }
826 
827  for (i=1; i<ringCount; i++)
828  {
829  if (point_in_ring_rtree(root[i], &pt) != -1)
830  {
831  POSTGIS_DEBUGF(3, "point_in_polygon_rtree: within hole %d.", i);
832 
833  return 0;
834  }
835  }
836  return 1;
837 }
838 
839 /*
840  * return -1 if point outside polygon
841  * return 0 if point on boundary
842  * return 1 if point inside polygon
843  *
844  * Expected **root order is each exterior ring followed by its holes, eg. EIIEIIEI
845  */
846 int point_in_multipolygon_rtree(RTREE_NODE **root, int polyCount, int *ringCounts, LWPOINT *point)
847 {
848  int i, p, r, in_ring;
849  POINT2D pt;
850  int result = -1;
851 
852  POSTGIS_DEBUGF(2, "point_in_multipolygon_rtree called for %p %d %p.", root, polyCount, point);
853 
854  getPoint2d_p(point->point, 0, &pt);
855  /* assume bbox short-circuit has already been attempted */
856 
857  i = 0; /* the current index into the root array */
858 
859  /* is the point inside any of the sub-polygons? */
860  for ( p = 0; p < polyCount; p++ )
861  {
862  in_ring = point_in_ring_rtree(root[i], &pt);
863  POSTGIS_DEBUGF(4, "point_in_multipolygon_rtree: exterior ring (%d), point_in_ring returned %d", p, in_ring);
864  if ( in_ring == -1 ) /* outside the exterior ring */
865  {
866  POSTGIS_DEBUG(3, "point_in_multipolygon_rtree: outside exterior ring.");
867  }
868  else if ( in_ring == 0 ) /* on the boundary */
869  {
870  POSTGIS_DEBUGF(3, "point_in_multipolygon_rtree: on edge of exterior ring %d", p);
871  return 0;
872  } else {
873  result = in_ring;
874 
875  for(r=1; r<ringCounts[p]; r++)
876  {
877  in_ring = point_in_ring_rtree(root[i+r], &pt);
878  POSTGIS_DEBUGF(4, "point_in_multipolygon_rtree: interior ring (%d), point_in_ring returned %d", r, in_ring);
879  if (in_ring == 1) /* inside a hole => outside the polygon */
880  {
881  POSTGIS_DEBUGF(3, "point_in_multipolygon_rtree: within hole %d of exterior ring %d", r, p);
882  result = -1;
883  break;
884  }
885  if (in_ring == 0) /* on the edge of a hole */
886  {
887  POSTGIS_DEBUGF(3, "point_in_multipolygon_rtree: on edge of hole %d of exterior ring %d", r, p);
888  return 0;
889  }
890  }
891  /* if we have a positive result, we can short-circuit and return it */
892  if ( result != -1)
893  {
894  return result;
895  }
896  }
897  /* increment the index by the total number of rings in the sub-poly */
898  /* we do this here in case we short-cutted out of the poly before looking at all the rings */
899  i += ringCounts[p];
900  }
901 
902  return result; /* -1 = outside, 0 = boundary, 1 = inside */
903 
904 }
905 
906 /*
907  * return -1 iff point outside polygon
908  * return 0 iff point on boundary
909  * return 1 iff point inside polygon
910  */
911 int point_in_polygon(LWPOLY *polygon, LWPOINT *point)
912 {
913  uint32_t i;
914  int result, in_ring;
915  POINT2D pt;
916 
917  POSTGIS_DEBUG(2, "point_in_polygon called.");
918 
919  getPoint2d_p(point->point, 0, &pt);
920  /* assume bbox short-circuit has already been attempted */
921 
922  /* everything is outside of an empty polygon */
923  if ( polygon->nrings == 0 ) return -1;
924 
925  in_ring = point_in_ring(polygon->rings[0], &pt);
926  if ( in_ring == -1) /* outside the exterior ring */
927  {
928  POSTGIS_DEBUG(3, "point_in_polygon: outside exterior ring.");
929  return -1;
930  }
931  result = in_ring;
932 
933  for (i=1; i<polygon->nrings; i++)
934  {
935  in_ring = point_in_ring(polygon->rings[i], &pt);
936  if (in_ring == 1) /* inside a hole => outside the polygon */
937  {
938  POSTGIS_DEBUGF(3, "point_in_polygon: within hole %d.", i);
939  return -1;
940  }
941  if (in_ring == 0) /* on the edge of a hole */
942  {
943  POSTGIS_DEBUGF(3, "point_in_polygon: on edge of hole %d.", i);
944  return 0;
945  }
946  }
947  return result; /* -1 = outside, 0 = boundary, 1 = inside */
948 }
949 
950 /*
951  * return -1 iff point outside multipolygon
952  * return 0 iff point on multipolygon boundary
953  * return 1 iff point inside multipolygon
954  */
955 int point_in_multipolygon(LWMPOLY *mpolygon, LWPOINT *point)
956 {
957  uint32_t i, j;
958  int result, in_ring;
959  POINT2D pt;
960 
961  POSTGIS_DEBUG(2, "point_in_polygon called.");
962 
963  getPoint2d_p(point->point, 0, &pt);
964  /* assume bbox short-circuit has already been attempted */
965 
966  result = -1;
967 
968  for (j = 0; j < mpolygon->ngeoms; j++ )
969  {
970 
971  LWPOLY *polygon = mpolygon->geoms[j];
972 
973  /* everything is outside of an empty polygon */
974  if ( polygon->nrings == 0 ) continue;
975 
976  in_ring = point_in_ring(polygon->rings[0], &pt);
977  if ( in_ring == -1) /* outside the exterior ring */
978  {
979  POSTGIS_DEBUG(3, "point_in_polygon: outside exterior ring.");
980  continue;
981  }
982  if ( in_ring == 0 )
983  {
984  return 0;
985  }
986 
987  result = in_ring;
988 
989  for (i=1; i<polygon->nrings; i++)
990  {
991  in_ring = point_in_ring(polygon->rings[i], &pt);
992  if (in_ring == 1) /* inside a hole => outside the polygon */
993  {
994  POSTGIS_DEBUGF(3, "point_in_polygon: within hole %d.", i);
995  result = -1;
996  break;
997  }
998  if (in_ring == 0) /* on the edge of a hole */
999  {
1000  POSTGIS_DEBUGF(3, "point_in_polygon: on edge of hole %d.", i);
1001  return 0;
1002  }
1003  }
1004  if ( result != -1)
1005  {
1006  return result;
1007  }
1008  }
1009  return result;
1010 }
1011 
1012 
1013 /*******************************************************************************
1014  * End of "Fast Winding Number Inclusion of a Point in a Polygon" derivative.
1015  ******************************************************************************/
1016 
1017 /**********************************************************************
1018  *
1019  * ST_MinimumBoundingRadius
1020  *
1021  **********************************************************************/
1022 
1024 Datum ST_MinimumBoundingRadius(PG_FUNCTION_ARGS)
1025 {
1026  GSERIALIZED* geom;
1027  LWGEOM* input;
1028  LWBOUNDINGCIRCLE* mbc = NULL;
1029  LWGEOM* lwcenter;
1030  GSERIALIZED* center;
1031  TupleDesc resultTupleDesc;
1032  HeapTuple resultTuple;
1033  Datum result;
1034  Datum result_values[2];
1035  bool result_is_null[2];
1036  double radius = 0;
1037 
1038  if (PG_ARGISNULL(0))
1039  PG_RETURN_NULL();
1040 
1041  geom = PG_GETARG_GSERIALIZED_P(0);
1042 
1043  /* Empty geometry? Return POINT EMPTY with zero radius */
1044  if (gserialized_is_empty(geom))
1045  {
1047  }
1048  else
1049  {
1050  input = lwgeom_from_gserialized(geom);
1051  mbc = lwgeom_calculate_mbc(input);
1052 
1053  if (!(mbc && mbc->center))
1054  {
1055  lwpgerror("Error calculating minimum bounding circle.");
1056  lwgeom_free(input);
1057  PG_RETURN_NULL();
1058  }
1059 
1060  lwcenter = (LWGEOM*) lwpoint_make2d(input->srid, mbc->center->x, mbc->center->y);
1061  radius = mbc->radius;
1062 
1064  lwgeom_free(input);
1065  }
1066 
1067  center = geometry_serialize(lwcenter);
1068  lwgeom_free(lwcenter);
1069 
1070  get_call_result_type(fcinfo, NULL, &resultTupleDesc);
1071  BlessTupleDesc(resultTupleDesc);
1072 
1073  result_values[0] = PointerGetDatum(center);
1074  result_is_null[0] = false;
1075  result_values[1] = Float8GetDatum(radius);
1076  result_is_null[1] = false;
1077 
1078  resultTuple = heap_form_tuple(resultTupleDesc, result_values, result_is_null);
1079 
1080  result = HeapTupleGetDatum(resultTuple);
1081 
1082  PG_RETURN_DATUM(result);
1083 }
1084 
1085 /**********************************************************************
1086  *
1087  * ST_MinimumBoundingCircle
1088  *
1089  **********************************************************************/
1090 
1092 Datum ST_MinimumBoundingCircle(PG_FUNCTION_ARGS)
1093 {
1094  GSERIALIZED* geom;
1095  LWGEOM* input;
1096  LWBOUNDINGCIRCLE* mbc = NULL;
1097  LWGEOM* lwcircle;
1098  GSERIALIZED* center;
1099  int segs_per_quarter;
1100 
1101  if (PG_ARGISNULL(0))
1102  PG_RETURN_NULL();
1103 
1104  geom = PG_GETARG_GSERIALIZED_P(0);
1105  segs_per_quarter = PG_GETARG_INT32(1);
1106 
1107  /* Empty geometry? Return POINT EMPTY */
1108  if (gserialized_is_empty(geom))
1109  {
1111  }
1112  else
1113  {
1114  input = lwgeom_from_gserialized(geom);
1115  mbc = lwgeom_calculate_mbc(input);
1116 
1117  if (!(mbc && mbc->center))
1118  {
1119  lwpgerror("Error calculating minimum bounding circle.");
1120  lwgeom_free(input);
1121  PG_RETURN_NULL();
1122  }
1123 
1124  /* Zero radius? Return a point. */
1125  if (mbc->radius == 0)
1126  lwcircle = lwpoint_as_lwgeom(lwpoint_make2d(input->srid, mbc->center->x, mbc->center->y));
1127  else
1128  lwcircle = lwpoly_as_lwgeom(lwpoly_construct_circle(input->srid, mbc->center->x, mbc->center->y, mbc->radius, segs_per_quarter, LW_TRUE));
1129 
1131  lwgeom_free(input);
1132  }
1133 
1134  center = geometry_serialize(lwcircle);
1135  lwgeom_free(lwcircle);
1136 
1137  PG_RETURN_POINTER(center);
1138 }
1139 
1140 /**********************************************************************
1141  *
1142  * ST_GeometricMedian
1143  *
1144  **********************************************************************/
1145 
1147 Datum ST_GeometricMedian(PG_FUNCTION_ARGS)
1148 {
1149  GSERIALIZED* geom;
1150  GSERIALIZED* result;
1151  LWGEOM* input;
1152  LWPOINT* lwresult;
1153  static const double min_default_tolerance = 1e-8;
1154  double tolerance = min_default_tolerance;
1155  bool compute_tolerance_from_box;
1156  bool fail_if_not_converged;
1157  int max_iter;
1158 
1159  /* Read and validate our input arguments */
1160  if (PG_ARGISNULL(0))
1161  PG_RETURN_NULL();
1162 
1163  compute_tolerance_from_box = PG_ARGISNULL(1);
1164 
1165  if (!compute_tolerance_from_box)
1166  {
1167  tolerance = PG_GETARG_FLOAT8(1);
1168  if (tolerance < 0)
1169  {
1170  lwpgerror("Tolerance must be positive.");
1171  PG_RETURN_NULL();
1172  }
1173  }
1174 
1175  max_iter = PG_ARGISNULL(2) ? -1 : PG_GETARG_INT32(2);
1176  fail_if_not_converged = PG_ARGISNULL(3) ? LW_FALSE : PG_GETARG_BOOL(3);
1177 
1178  if (max_iter < 0)
1179  {
1180  lwpgerror("Maximum iterations must be positive.");
1181  PG_RETURN_NULL();
1182  }
1183 
1184  /* OK, inputs are valid. */
1185  geom = PG_GETARG_GSERIALIZED_P(0);
1186  input = lwgeom_from_gserialized(geom);
1187 
1188  if (compute_tolerance_from_box)
1189  {
1190  /* Compute a default tolerance based on the smallest dimension
1191  * of the geometry's bounding box.
1192  */
1193  static const double tolerance_coefficient = 1e-6;
1194  const GBOX* box = lwgeom_get_bbox(input);
1195 
1196  if (box)
1197  {
1198  double min_dim = FP_MIN(box->xmax - box->xmin, box->ymax - box->ymin);
1199  if (lwgeom_has_z(input))
1200  min_dim = FP_MIN(min_dim, box->zmax - box->zmin);
1201 
1202  /* Apply a lower bound to the computed default tolerance to
1203  * avoid a tolerance of zero in the case of collinear
1204  * points.
1205  */
1206  tolerance = FP_MAX(min_default_tolerance, tolerance_coefficient * min_dim);
1207  }
1208  }
1209 
1210  lwresult = lwgeom_median(input, tolerance, max_iter, fail_if_not_converged);
1211  lwgeom_free(input);
1212 
1213  if(!lwresult)
1214  {
1215  lwpgerror("Error computing geometric median.");
1216  PG_RETURN_NULL();
1217  }
1218 
1219  result = geometry_serialize(lwpoint_as_lwgeom(lwresult));
1220 
1221  PG_RETURN_POINTER(result);
1222 }
1223 
1224 /**********************************************************************
1225  *
1226  * ST_IsPolygonCW
1227  *
1228  **********************************************************************/
1229 
1231 Datum ST_IsPolygonCW(PG_FUNCTION_ARGS)
1232 {
1233  GSERIALIZED* geom;
1234  LWGEOM* input;
1235  bool is_clockwise;
1236 
1237  if (PG_ARGISNULL(0))
1238  PG_RETURN_NULL();
1239 
1240  geom = PG_GETARG_GSERIALIZED_P(0);
1241  input = lwgeom_from_gserialized(geom);
1242 
1243  is_clockwise = lwgeom_is_clockwise(input);
1244 
1245  lwgeom_free(input);
1246  PG_FREE_IF_COPY(geom, 0);
1247 
1248  PG_RETURN_BOOL(is_clockwise);
1249 }
1250 
1251 /**********************************************************************
1252  *
1253  * ST_IsPolygonCCW
1254  *
1255  **********************************************************************/
1256 
1258 Datum ST_IsPolygonCCW(PG_FUNCTION_ARGS)
1259 {
1260  GSERIALIZED* geom;
1261  LWGEOM* input;
1262  bool is_ccw;
1263 
1264  if (PG_ARGISNULL(0))
1265  PG_RETURN_NULL();
1266 
1267  geom = PG_GETARG_GSERIALIZED_P_COPY(0);
1268  input = lwgeom_from_gserialized(geom);
1269  lwgeom_reverse_in_place(input);
1270  is_ccw = lwgeom_is_clockwise(input);
1271  lwgeom_free(input);
1272  PG_FREE_IF_COPY(geom, 0);
1273 
1274  PG_RETURN_BOOL(is_ccw);
1275 }
double x
Definition: liblwgeom.h:351
#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:397
static int point_in_ring(POINTARRAY *pts, const POINT2D *point)
double m
Definition: liblwgeom.h:351
LWGEOM * lwgeom_grid(const LWGEOM *lwgeom, const gridspec *grid)
Definition: lwgeom.c:2236
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:460
void lwgeom_reverse_in_place(LWGEOM *lwgeom)
Reverse vertex order of LWGEOM.
Definition: lwgeom.c:101
uint32_t ngeoms
Definition: liblwgeom.h:493
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:45
LWPOINT * lwpoint_make2d(int srid, double x, double y)
Definition: lwpoint.c:163
Datum area(PG_FUNCTION_ARGS)
double xmax
Definition: liblwgeom.h:292
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1137
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'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:160
int32_t srid
Definition: liblwgeom.h:420
Datum LWGEOM_SetEffectiveArea(PG_FUNCTION_ARGS)
#define FP_MIN(A, B)
POINTARRAY * point
Definition: liblwgeom.h:410
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:923
Datum ST_IsPolygonCCW(PG_FUNCTION_ARGS)
int32_t srid
Definition: liblwgeom.h:398
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:319
uint32_t nrings
Definition: liblwgeom.h:454
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:327
Datum LWGEOM_snaptogrid_pointoff(PG_FUNCTION_ARGS)
LWBOUNDINGCIRCLE * lwgeom_calculate_mbc(const LWGEOM *g)
double zmax
Definition: liblwgeom.h:296
double ymin
Definition: liblwgeom.h:293
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:329
double xmin
Definition: liblwgeom.h:291
#define LW_FALSE
Definition: liblwgeom.h:76
LWPOLY ** geoms
Definition: liblwgeom.h:495
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:1835
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:187
LWMLINE * lwgeom_as_lwmline(const LWGEOM *lwgeom)
Definition: lwgeom.c:241
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:727
POINTARRAY ** rings
Definition: liblwgeom.h:456
double ymax
Definition: liblwgeom.h:294
double y
Definition: liblwgeom.h:327
#define FLAGS_GET_Z(flags)
Macros for manipulating the 'flags' byte.
Definition: liblwgeom.h:139
double z
Definition: liblwgeom.h:351
tuple x
Definition: pixval.py:53
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:169
int getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *point)
Definition: lwgeom_api.c:113
POINT2D * center
Definition: liblwgeom.h:1643
uint8_t flags
Definition: liblwgeom.h:407
uint32_t ngeoms
Definition: liblwgeom.h:480
LWLINE ** geoms
Definition: liblwgeom.h:482
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)
GSERIALIZED * geometry_serialize(LWGEOM *lwgeom)
Datum LWGEOM_line_interpolate_point(PG_FUNCTION_ARGS)
double zmin
Definition: liblwgeom.h:295
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:679
#define FLAGS_GET_M(flags)
Definition: liblwgeom.h:140
uint8_t type
Definition: liblwgeom.h:395
static double determineSide(const POINT2D *seg1, const POINT2D *seg2, const POINT2D *point)
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:334
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:64
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:1386
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:351
#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)
tuple y
Definition: pixval.py:54
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))
Definition: lwin_wkt_lex.c:888
#define COLLECTIONTYPE
Definition: liblwgeom.h:90
This library is the generic geometry handling section of PostGIS.
#define FP_MAX(A, B)
int32_t srid
Definition: liblwgeom.h:479
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:421
LWGEOM * lwmpoint_as_lwgeom(const LWMPOINT *obj)
Definition: lwgeom.c:294
Snap to grid.
uint32_t npoints
Definition: liblwgeom.h:370