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