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