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