source: trunk/postgis/lwgeom_functions_analytic.c @ 4168

Last change on this file since 4168 was 4168, checked in by pramsey, 7 years ago

astyle --style=ansi --indent=tab (#133)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
  • Property svn:mime-type set to text/plain
File size: 39.9 KB
Line 
1/**********************************************************************
2 * $Id: lwgeom_functions_analytic.c 4168 2009-06-11 16:44:03Z pramsey $
3 *
4 * PostGIS - Spatial Types for PostgreSQL
5 * http://postgis.refractions.net
6 * Copyright 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 "lwgeom_pg.h"
17#include "math.h"
18#include "lwgeom_rtree.h"
19#include "lwalgorithm.h"
20
21
22/***********************************************************************
23 * Simple Douglas-Peucker line simplification.
24 * No checks are done to avoid introduction of self-intersections.
25 * No topology relations are considered.
26 *
27 * --strk@keybit.net;
28 ***********************************************************************/
29
30
31/* Prototypes */
32void DP_findsplit2d(POINTARRAY *pts, int p1, int p2, int *split, double *dist);
33POINTARRAY *DP_simplify2d(POINTARRAY *inpts, double epsilon);
34LWLINE *simplify2d_lwline(const LWLINE *iline, double dist);
35LWPOLY *simplify2d_lwpoly(const LWPOLY *ipoly, double dist);
36LWCOLLECTION *simplify2d_collection(const LWCOLLECTION *igeom, double dist);
37LWGEOM *simplify2d_lwgeom(const LWGEOM *igeom, double dist);
38Datum LWGEOM_simplify2d(PG_FUNCTION_ARGS);
39Datum ST_LineCrossingDirection(PG_FUNCTION_ARGS);
40Datum ST_LocateBetweenElevations(PG_FUNCTION_ARGS);
41
42double determineSide(POINT2D *seg1, POINT2D *seg2, POINT2D *point);
43int isOnSegment(POINT2D *seg1, POINT2D *seg2, POINT2D *point);
44int point_in_ring(POINTARRAY *pts, POINT2D *point);
45int point_in_polygon(LWPOLY *polygon, LWPOINT *point);
46int point_in_multipolygon(LWMPOLY *mpolygon, LWPOINT *point);
47int point_in_ring_rtree(RTREE_NODE *root, POINT2D *point);
48int point_in_polygon_rtree(RTREE_NODE **root, int ringCount, LWPOINT *point);
49int point_in_multipolygon_rtree(RTREE_NODE **root, int polyCount, int ringCount, LWPOINT *point);
50
51
52/*
53 * Search farthest point from segment p1-p2
54 * returns distance in an int pointer
55 */
56void
57DP_findsplit2d(POINTARRAY *pts, int p1, int p2, int *split, double *dist)
58{
59        int k;
60        POINT2D pa, pb, pk;
61        double tmp;
62
63        LWDEBUG(4, "DP_findsplit called");
64
65        *dist = -1;
66        *split = p1;
67
68        if (p1 + 1 < p2)
69        {
70
71                getPoint2d_p(pts, p1, &pa);
72                getPoint2d_p(pts, p2, &pb);
73
74                LWDEBUGF(4, "DP_findsplit: P%d(%f,%f) to P%d(%f,%f)",
75                         p1, pa.x, pa.y, p2, pb.x, pb.y);
76
77                for (k=p1+1; k<p2; k++)
78                {
79                        getPoint2d_p(pts, k, &pk);
80
81                        LWDEBUGF(4, "DP_findsplit: P%d(%f,%f)", k, pk.x, pk.y);
82
83                        /* distance computation */
84                        tmp = distance2d_pt_seg(&pk, &pa, &pb);
85
86                        if (tmp > *dist)
87                        {
88                                *dist = tmp;    /* record the maximum */
89                                *split = k;
90
91                                LWDEBUGF(4, "DP_findsplit: P%d is farthest (%g)", k, *dist);
92                        }
93                }
94
95        } /* length---should be redone if can == 0 */
96
97        else
98        {
99                LWDEBUG(3, "DP_findsplit: segment too short, no split/no dist");
100        }
101
102}
103
104
105POINTARRAY *
106DP_simplify2d(POINTARRAY *inpts, double epsilon)
107{
108        int *stack;                     /* recursion stack */
109        int sp=-1;                      /* recursion stack pointer */
110        int p1, split;
111        double dist;
112        POINTARRAY *outpts;
113        int ptsize = pointArray_ptsize(inpts);
114
115        /* Allocate recursion stack */
116        stack = lwalloc(sizeof(int)*inpts->npoints);
117
118        p1 = 0;
119        stack[++sp] = inpts->npoints-1;
120
121        LWDEBUGF(2, "DP_simplify called input has %d pts and %d dims (ptsize: %d)", inpts->npoints, inpts->dims, ptsize);
122
123        /* allocate space for output POINTARRAY */
124        outpts = palloc(sizeof(POINTARRAY));
125        outpts->dims = inpts->dims;
126        outpts->npoints=1;
127        outpts->serialized_pointlist = palloc(ptsize*inpts->npoints);
128        memcpy(getPoint_internal(outpts, 0), getPoint_internal(inpts, 0),
129               ptsize);
130
131        LWDEBUG(3, "DP_simplify: added P0 to simplified point array (size 1)");
132
133        do
134        {
135
136                DP_findsplit2d(inpts, p1, stack[sp], &split, &dist);
137
138                LWDEBUGF(3, "DP_simplify: farthest point from P%d-P%d is P%d (dist. %g)", p1, stack[sp], split, dist);
139
140                if (dist > epsilon)
141                {
142                        stack[++sp] = split;
143                }
144                else
145                {
146                        outpts->npoints++;
147                        memcpy(getPoint_internal(outpts, outpts->npoints-1),
148                               getPoint_internal(inpts, stack[sp]),
149                               ptsize);
150
151                        LWDEBUGF(4, "DP_simplify: added P%d to simplified point array (size: %d)", stack[sp], outpts->npoints);
152
153                        p1 = stack[sp--];
154                }
155
156                LWDEBUGF(4, "stack pointer = %d", sp);
157        }
158        while (! (sp<0) );
159
160        /*
161         * If we have reduced the number of points realloc
162         * outpoints array to free up some memory.
163         * Might be turned on and off with a SAVE_MEMORY define ...
164         */
165        if ( outpts->npoints < inpts->npoints )
166        {
167                outpts->serialized_pointlist = repalloc(
168                                                   outpts->serialized_pointlist,
169                                                   ptsize*outpts->npoints);
170                if ( outpts->serialized_pointlist == NULL )
171                {
172                        elog(ERROR, "Out of virtual memory");
173                }
174        }
175
176        lwfree(stack);
177        return outpts;
178}
179
180LWLINE *
181simplify2d_lwline(const LWLINE *iline, double dist)
182{
183        POINTARRAY *ipts;
184        POINTARRAY *opts;
185        LWLINE *oline;
186
187        LWDEBUG(2, "simplify2d_lwline called");
188
189        ipts = iline->points;
190        opts = DP_simplify2d(ipts, dist);
191        oline = lwline_construct(iline->SRID, NULL, opts);
192
193        return oline;
194}
195
196LWPOLY *
197simplify2d_lwpoly(const LWPOLY *ipoly, double dist)
198{
199        POINTARRAY *ipts;
200        POINTARRAY **orings = NULL;
201        LWPOLY *opoly;
202        int norings=0, ri;
203
204        LWDEBUGF(2, "simplify_polygon3d: simplifying polygon with %d rings", ipoly->nrings);
205
206        orings = (POINTARRAY **)palloc(sizeof(POINTARRAY *)*ipoly->nrings);
207
208        for (ri=0; ri<ipoly->nrings; ri++)
209        {
210                POINTARRAY *opts;
211
212                ipts = ipoly->rings[ri];
213
214                opts = DP_simplify2d(ipts, dist);
215
216
217                if ( opts->npoints < 2 )
218                {
219                        /* There as to be an error in DP_simplify */
220                        elog(NOTICE, "DP_simplify returned a <2 pts array");
221                        pfree(opts);
222                        continue;
223                }
224
225                if ( opts->npoints < 4 )
226                {
227                        pfree(opts);
228
229                        LWDEBUGF(3, "simplify_polygon3d: ring%d skipped ( <4 pts )", ri);
230
231                        if ( ri ) continue;
232                        else break;
233                }
234
235
236                LWDEBUGF(3, "simplify_polygon3d: ring%d simplified from %d to %d points", ri, ipts->npoints, opts->npoints);
237
238
239                /*
240                 * Add ring to simplified ring array
241                 * (TODO: dinamic allocation of pts_per_ring)
242                 */
243                orings[norings] = opts;
244                norings++;
245
246        }
247
248        LWDEBUGF(3, "simplify_polygon3d: simplified polygon with %d rings", norings);
249
250        if ( ! norings ) return NULL;
251
252        opoly = lwpoly_construct(ipoly->SRID, NULL, norings, orings);
253
254        return opoly;
255}
256
257LWCOLLECTION *
258simplify2d_collection(const LWCOLLECTION *igeom, double dist)
259{
260        unsigned int i;
261        unsigned int ngeoms=0;
262        LWGEOM **geoms = lwalloc(sizeof(LWGEOM *)*igeom->ngeoms);
263        LWCOLLECTION *out;
264
265        for (i=0; i<igeom->ngeoms; i++)
266        {
267                LWGEOM *ngeom = simplify2d_lwgeom(igeom->geoms[i], dist);
268                if ( ngeom ) geoms[ngeoms++] = ngeom;
269        }
270
271        out = lwcollection_construct(TYPE_GETTYPE(igeom->type), igeom->SRID,
272                                     NULL, ngeoms, geoms);
273
274        return out;
275}
276
277LWGEOM *
278simplify2d_lwgeom(const LWGEOM *igeom, double dist)
279{
280        switch (TYPE_GETTYPE(igeom->type))
281        {
282        case POINTTYPE:
283        case MULTIPOINTTYPE:
284                return lwgeom_clone(igeom);
285        case LINETYPE:
286                return (LWGEOM *)simplify2d_lwline(
287                           (LWLINE *)igeom, dist);
288        case POLYGONTYPE:
289                return (LWGEOM *)simplify2d_lwpoly(
290                           (LWPOLY *)igeom, dist);
291        case MULTILINETYPE:
292        case MULTIPOLYGONTYPE:
293        case COLLECTIONTYPE:
294                return (LWGEOM *)simplify2d_collection(
295                           (LWCOLLECTION *)igeom, dist);
296        default:
297                lwerror("simplify2d_lwgeom: unknown geometry type: %d",
298                        TYPE_GETTYPE(igeom->type));
299        }
300        return NULL;
301}
302
303PG_FUNCTION_INFO_V1(LWGEOM_simplify2d);
304Datum LWGEOM_simplify2d(PG_FUNCTION_ARGS)
305{
306        PG_LWGEOM *geom = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
307        LWGEOM *in = lwgeom_deserialize(SERIALIZED_FORM(geom));
308        LWGEOM *out;
309        PG_LWGEOM *result;
310        double dist = PG_GETARG_FLOAT8(1);
311
312        out = simplify2d_lwgeom(in, dist);
313        if ( ! out ) PG_RETURN_NULL();
314
315        /* COMPUTE_BBOX TAINTING */
316        if ( in->bbox ) lwgeom_add_bbox(out);
317
318        result = pglwgeom_serialize(out);
319
320        PG_RETURN_POINTER(result);
321}
322
323/***********************************************************************
324 * --strk@keybit.net;
325 ***********************************************************************/
326
327/***********************************************************************
328 * Interpolate a point along a line, useful for Geocoding applications
329 * SELECT line_interpolate_point( 'LINESTRING( 0 0, 2 2'::geometry, .5 )
330 * returns POINT( 1 1 ).
331 * Works in 2d space only.
332 *
333 * Initially written by: jsunday@rochgrp.com
334 * Ported to LWGEOM by: strk@refractions.net
335 ***********************************************************************/
336
337Datum LWGEOM_line_interpolate_point(PG_FUNCTION_ARGS);
338
339PG_FUNCTION_INFO_V1(LWGEOM_line_interpolate_point);
340Datum LWGEOM_line_interpolate_point(PG_FUNCTION_ARGS)
341{
342        PG_LWGEOM *geom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
343        double distance = PG_GETARG_FLOAT8(1);
344        LWLINE *line;
345        LWPOINT *point;
346        POINTARRAY *ipa, *opa;
347        POINT4D pt;
348        uchar *srl;
349        int nsegs, i;
350        double length, slength, tlength;
351
352        if ( distance < 0 || distance > 1 )
353        {
354                elog(ERROR,"line_interpolate_point: 2nd arg isnt within [0,1]");
355                PG_RETURN_NULL();
356        }
357
358        if ( lwgeom_getType(geom->type) != LINETYPE )
359        {
360                elog(ERROR,"line_interpolate_point: 1st arg isnt a line");
361                PG_RETURN_NULL();
362        }
363
364        line = lwline_deserialize(SERIALIZED_FORM(geom));
365        ipa = line->points;
366
367        /* If distance is one of the two extremes, return the point on that
368         * end rather than doing any expensive computations
369         */
370        if ( distance == 0.0 || distance == 1.0 )
371        {
372                if ( distance == 0.0 )
373                        getPoint4d_p(ipa, 0, &pt);
374                else
375                        getPoint4d_p(ipa, ipa->npoints-1, &pt);
376
377                opa = pointArray_construct((uchar *)&pt,
378                                           TYPE_HASZ(line->type),
379                                           TYPE_HASM(line->type),
380                                           1);
381                point = lwpoint_construct(line->SRID, 0, opa);
382                srl = lwpoint_serialize(point);
383                /* We shouldn't need this, the memory context is getting freed on the next line.
384                lwpoint_free(point); */
385                PG_RETURN_POINTER(PG_LWGEOM_construct(srl, line->SRID, 0));
386        }
387
388        /* Interpolate a point on the line */
389        nsegs = ipa->npoints - 1;
390        length = lwgeom_pointarray_length2d(ipa);
391        tlength = 0;
392        for ( i = 0; i < nsegs; i++ )
393        {
394                POINT4D p1, p2;
395                POINT4D *p1ptr=&p1, *p2ptr=&p2; /* don't break
396                                                                                 * strict-aliasing rules
397                                                                                 */
398
399                getPoint4d_p(ipa, i, &p1);
400                getPoint4d_p(ipa, i+1, &p2);
401
402                /* Find the relative length of this segment */
403                slength = distance2d_pt_pt((POINT2D*)p1ptr, (POINT2D*)p2ptr)/length;
404
405                /* If our target distance is before the total length we've seen
406                 * so far. create a new point some distance down the current
407                 * segment.
408                 */
409                if ( distance < tlength + slength )
410                {
411                        double dseg = (distance - tlength) / slength;
412                        interpolate_point4d(&p1, &p2, &pt, dseg);
413                        opa = pointArray_construct((uchar *)&pt,
414                                                   TYPE_HASZ(line->type),
415                                                   TYPE_HASM(line->type),
416                                                   1);
417                        point = lwpoint_construct(line->SRID, 0, opa);
418                        srl = lwpoint_serialize(point);
419                        /* We shouldn't need this, the memory context is getting freed on the next line
420                        lwpoint_free(point); */
421                        PG_RETURN_POINTER(PG_LWGEOM_construct(srl, line->SRID, 0));
422                }
423                tlength += slength;
424        }
425
426        /* Return the last point on the line. This shouldn't happen, but
427         * could if there's some floating point rounding errors. */
428        getPoint4d_p(ipa, ipa->npoints-1, &pt);
429        opa = pointArray_construct((uchar *)&pt,
430                                   TYPE_HASZ(line->type),
431                                   TYPE_HASM(line->type),
432                                   1);
433        point = lwpoint_construct(line->SRID, 0, opa);
434        srl = lwpoint_serialize(point);
435        /* We shouldn't need this, the memory context is getting freed on the next line
436        lwpoint_free(point); */
437        PG_RETURN_POINTER(PG_LWGEOM_construct(srl, line->SRID, 0));
438}
439/***********************************************************************
440 * --jsunday@rochgrp.com;
441 ***********************************************************************/
442
443/***********************************************************************
444 *
445 *  Grid application function for postgis.
446 *
447 *  WHAT IS
448 *  -------
449 *
450 *  This is a grid application function for postgis.
451 *  You use it to stick all of a geometry points to
452 *  a custom grid defined by its origin and cell size
453 *  in geometry units.
454 *
455 *  Points reduction is obtained by collapsing all
456 *  consecutive points falling on the same grid node and
457 *  removing all consecutive segments S1,S2 having
458 *  S2.startpoint = S1.endpoint and S2.endpoint = S1.startpoint.
459 *
460 *  ISSUES
461 *  ------
462 *
463 *  Only 2D is contemplated in grid application.
464 *
465 *  Consecutive coincident segments removal does not work
466 *  on first/last segments (They are not considered consecutive).
467 *
468 *  Grid application occurs on a geometry subobject basis, thus no
469 *  points reduction occurs for multipoint geometries.
470 *
471 *  USAGE TIPS
472 *  ----------
473 *
474 *  Run like this:
475 *
476 *     SELECT SnapToGrid(<geometry>, <originX>, <originY>, <sizeX>, <sizeY>);
477 *
478 *     Grid cells are not visible on a screen as long as the screen
479 *     point size is equal or bigger then the grid cell size.
480 *     This assumption may be used to reduce the number of points in
481 *     a map for a given display scale.
482 *
483 *     Keeping multiple resolution versions of the same data may be used
484 *     in conjunction with MINSCALE/MAXSCALE keywords of mapserv to speed
485 *     up rendering.
486 *
487 *     Check also the use of DP simplification to reduce grid visibility.
488 *     I'm still researching about the relationship between grid size and
489 *     DP epsilon values - please tell me if you know more about this.
490 *
491 *
492 * --strk@keybit.net;
493 *
494 ***********************************************************************/
495
496#define CHECK_RING_IS_CLOSE
497#define SAMEPOINT(a,b) ((a)->x==(b)->x&&(a)->y==(b)->y)
498
499typedef struct gridspec_t
500{
501        double ipx;
502        double ipy;
503        double ipz;
504        double ipm;
505        double xsize;
506        double ysize;
507        double zsize;
508        double msize;
509}
510gridspec;
511
512
513/* Forward declarations */
514LWGEOM *lwgeom_grid(LWGEOM *lwgeom, gridspec *grid);
515LWCOLLECTION *lwcollection_grid(LWCOLLECTION *coll, gridspec *grid);
516LWPOINT * lwpoint_grid(LWPOINT *point, gridspec *grid);
517LWPOLY * lwpoly_grid(LWPOLY *poly, gridspec *grid);
518LWLINE *lwline_grid(LWLINE *line, gridspec *grid);
519POINTARRAY *ptarray_grid(POINTARRAY *pa, gridspec *grid);
520Datum LWGEOM_snaptogrid(PG_FUNCTION_ARGS);
521Datum LWGEOM_snaptogrid_pointoff(PG_FUNCTION_ARGS);
522static int grid_isNull(const gridspec *grid);
523#if POSTGIS_DEBUG_LEVEL > 0
524static void grid_print(const gridspec *grid);
525#endif
526
527/* A NULL grid is a grid in which size in all dimensions is 0 */
528static int
529grid_isNull(const gridspec *grid)
530{
531        if ( grid->xsize==0 &&
532                grid->ysize==0 &&
533                grid->zsize==0 &&
534                grid->msize==0 ) return 1;
535        else return 0;
536}
537
538#if POSTGIS_DEBUG_LEVEL > 0
539/* Print grid using given reporter */
540static void
541grid_print(const gridspec *grid)
542{
543        lwnotice("GRID(%g %g %g %g, %g %g %g %g)",
544                 grid->ipx, grid->ipy, grid->ipz, grid->ipm,
545                 grid->xsize, grid->ysize, grid->zsize, grid->msize);
546}
547#endif
548
549/*
550 * Stick an array of points to the given gridspec.
551 * Return "gridded" points in *outpts and their number in *outptsn.
552 *
553 * Two consecutive points falling on the same grid cell are collapsed
554 * into one single point.
555 *
556 */
557POINTARRAY *
558ptarray_grid(POINTARRAY *pa, gridspec *grid)
559{
560        POINT4D pbuf;
561        int ipn, opn; /* point numbers (input/output) */
562        DYNPTARRAY *dpa;
563        POINTARRAY *opa;
564
565        LWDEBUGF(2, "ptarray_grid called on %p", pa);
566
567        dpa=dynptarray_create(pa->npoints, pa->dims);
568
569        for (ipn=0, opn=0; ipn<pa->npoints; ++ipn)
570        {
571
572                getPoint4d_p(pa, ipn, &pbuf);
573
574                if ( grid->xsize )
575                        pbuf.x = rint((pbuf.x - grid->ipx)/grid->xsize) *
576                                 grid->xsize + grid->ipx;
577
578                if ( grid->ysize )
579                        pbuf.y = rint((pbuf.y - grid->ipy)/grid->ysize) *
580                                 grid->ysize + grid->ipy;
581
582                if ( TYPE_HASZ(pa->dims) && grid->zsize )
583                        pbuf.z = rint((pbuf.z - grid->ipz)/grid->zsize) *
584                                 grid->zsize + grid->ipz;
585
586                if ( TYPE_HASM(pa->dims) && grid->msize )
587                        pbuf.m = rint((pbuf.m - grid->ipm)/grid->msize) *
588                                 grid->msize + grid->ipm;
589
590                dynptarray_addPoint4d(dpa, &pbuf, 0);
591
592        }
593
594        opa = dpa->pa;
595        lwfree(dpa);
596
597        return opa;
598}
599
600LWLINE *
601lwline_grid(LWLINE *line, gridspec *grid)
602{
603        LWLINE *oline;
604        POINTARRAY *opa;
605
606        opa = ptarray_grid(line->points, grid);
607
608        /* Skip line3d with less then 2 points */
609        if ( opa->npoints < 2 ) return NULL;
610
611        /* TODO: grid bounding box... */
612        oline = lwline_construct(line->SRID, NULL, opa);
613
614        return oline;
615}
616
617LWPOLY *
618lwpoly_grid(LWPOLY *poly, gridspec *grid)
619{
620        LWPOLY *opoly;
621        int ri;
622        POINTARRAY **newrings = NULL;
623        int nrings = 0;
624        double minvisiblearea;
625
626        /*
627         * TODO: control this assertion
628         * it is assumed that, since the grid size will be a pixel,
629         * a visible ring should show at least a white pixel inside,
630         * thus, for a square, that would be grid_xsize*grid_ysize
631         */
632        minvisiblearea = grid->xsize * grid->ysize;
633
634        nrings = 0;
635
636        LWDEBUGF(3, "grid_polygon3d: applying grid to polygon with %d rings",
637                 poly->nrings);
638
639        for (ri=0; ri<poly->nrings; ri++)
640        {
641                POINTARRAY *ring = poly->rings[ri];
642                POINTARRAY *newring;
643
644#if POSTGIS_DEBUG_LEVEL >= 4
645                POINT2D p1, p2;
646                getPoint2d_p(ring, 0, &p1);
647                getPoint2d_p(ring, ring->npoints-1, &p2);
648                if ( ! SAMEPOINT(&p1, &p2) )
649                        LWDEBUG(4, "Before gridding: first point != last point");
650#endif
651
652                newring = ptarray_grid(ring, grid);
653
654                /* Skip ring if not composed by at least 4 pts (3 segments) */
655                if ( newring->npoints < 4 )
656                {
657                        pfree(newring);
658
659                        LWDEBUGF(3, "grid_polygon3d: ring%d skipped ( <4 pts )", ri);
660
661                        if ( ri ) continue;
662                        else break; /* this is the external ring, no need to work on holes */
663                }
664
665#if POSTGIS_DEBUG_LEVEL >= 4
666                getPoint2d_p(newring, 0, &p1);
667                getPoint2d_p(newring, newring->npoints-1, &p2);
668                if ( ! SAMEPOINT(&p1, &p2) )
669                        LWDEBUG(4, "After gridding: first point != last point");
670#endif
671
672                LWDEBUGF(3, "grid_polygon3d: ring%d simplified from %d to %d points", ri,
673                         ring->npoints, newring->npoints);
674
675                /*
676                 * Add ring to simplified ring array
677                 * (TODO: dinamic allocation of pts_per_ring)
678                 */
679                if ( ! nrings )
680                {
681                        newrings = palloc(sizeof(POINTARRAY *));
682                }
683                else
684                {
685                        newrings = repalloc(newrings, sizeof(POINTARRAY *)*(nrings+1));
686                }
687                if ( ! newrings )
688                {
689                        elog(ERROR, "Out of virtual memory");
690                        return NULL;
691                }
692                newrings[nrings++] = newring;
693        }
694
695        LWDEBUGF(3, "grid_polygon3d: simplified polygon with %d rings", nrings);
696
697        if ( ! nrings ) return NULL;
698
699        opoly = lwpoly_construct(poly->SRID, NULL, nrings, newrings);
700        return opoly;
701}
702
703LWPOINT *
704lwpoint_grid(LWPOINT *point, gridspec *grid)
705{
706        LWPOINT *opoint;
707        POINTARRAY *opa;
708
709        opa = ptarray_grid(point->point, grid);
710
711        /* TODO: grid bounding box ? */
712        opoint = lwpoint_construct(point->SRID, NULL, opa);
713
714        LWDEBUG(2, "lwpoint_grid called");
715
716        return opoint;
717}
718
719LWCOLLECTION *
720lwcollection_grid(LWCOLLECTION *coll, gridspec *grid)
721{
722        unsigned int i;
723        LWGEOM **geoms;
724        unsigned int ngeoms=0;
725
726        geoms = palloc(coll->ngeoms * sizeof(LWGEOM *));
727
728        for (i=0; i<coll->ngeoms; i++)
729        {
730                LWGEOM *g = lwgeom_grid(coll->geoms[i], grid);
731                if ( g ) geoms[ngeoms++] = g;
732        }
733
734        if ( ! ngeoms ) return lwcollection_construct_empty(coll->SRID, 0, 0);
735
736        return lwcollection_construct(TYPE_GETTYPE(coll->type), coll->SRID,
737                                      NULL, ngeoms, geoms);
738}
739
740LWGEOM *
741lwgeom_grid(LWGEOM *lwgeom, gridspec *grid)
742{
743        switch (TYPE_GETTYPE(lwgeom->type))
744        {
745        case POINTTYPE:
746                return (LWGEOM *)lwpoint_grid((LWPOINT *)lwgeom, grid);
747        case LINETYPE:
748                return (LWGEOM *)lwline_grid((LWLINE *)lwgeom, grid);
749        case POLYGONTYPE:
750                return (LWGEOM *)lwpoly_grid((LWPOLY *)lwgeom, grid);
751        case MULTIPOINTTYPE:
752        case MULTILINETYPE:
753        case MULTIPOLYGONTYPE:
754        case COLLECTIONTYPE:
755                return (LWGEOM *)lwcollection_grid((LWCOLLECTION *)lwgeom, grid);
756        default:
757                elog(ERROR, "lwgeom_grid: Unsupported geometry type: %s",
758                     lwgeom_typename(TYPE_GETTYPE(lwgeom->type)));
759                return NULL;
760        }
761}
762
763PG_FUNCTION_INFO_V1(LWGEOM_snaptogrid);
764Datum LWGEOM_snaptogrid(PG_FUNCTION_ARGS)
765{
766        Datum datum;
767        PG_LWGEOM *in_geom;
768        LWGEOM *in_lwgeom;
769        PG_LWGEOM *out_geom = NULL;
770        LWGEOM *out_lwgeom;
771        gridspec grid;
772        /* BOX3D box3d; */
773
774        if ( PG_ARGISNULL(0) ) PG_RETURN_NULL();
775        datum = PG_GETARG_DATUM(0);
776        in_geom = (PG_LWGEOM *)PG_DETOAST_DATUM(datum);
777
778        if ( PG_ARGISNULL(1) ) PG_RETURN_NULL();
779        grid.ipx = PG_GETARG_FLOAT8(1);
780
781        if ( PG_ARGISNULL(2) ) PG_RETURN_NULL();
782        grid.ipy = PG_GETARG_FLOAT8(2);
783
784        if ( PG_ARGISNULL(3) ) PG_RETURN_NULL();
785        grid.xsize = PG_GETARG_FLOAT8(3);
786
787        if ( PG_ARGISNULL(4) ) PG_RETURN_NULL();
788        grid.ysize = PG_GETARG_FLOAT8(4);
789
790        /* Do not support gridding Z and M values for now */
791        grid.ipz=grid.ipm=grid.zsize=grid.msize=0;
792
793        /* Return input geometry if grid is null */
794        if ( grid_isNull(&grid) )
795        {
796                PG_RETURN_POINTER(in_geom);
797        }
798
799        in_lwgeom = lwgeom_deserialize(SERIALIZED_FORM(in_geom));
800
801        POSTGIS_DEBUGF(3, "SnapToGrid got a %s", lwgeom_typename(TYPE_GETTYPE(in_lwgeom->type)));
802
803        out_lwgeom = lwgeom_grid(in_lwgeom, &grid);
804        if ( out_lwgeom == NULL ) PG_RETURN_NULL();
805
806        /* COMPUTE_BBOX TAINTING */
807        if ( in_lwgeom->bbox ) lwgeom_add_bbox(out_lwgeom);
808
809#if 0
810        /*
811         * COMPUTE_BBOX WHEN SIMPLE
812         *
813         * WARNING: this is not SIMPLE, as snapping
814         * an existing bbox to a grid does not
815         * give the same result as computing a
816         * new bounding box on the snapped coordinates.
817         *
818         * This bug has been fixed in postgis-1.1.2
819         */
820        if ( in_lwgeom->bbox )
821        {
822                box2df_to_box3d_p(in_lwgeom->bbox, &box3d);
823
824                box3d.xmin = rint((box3d.xmin - grid.ipx)/grid.xsize)
825                             * grid.xsize + grid.ipx;
826                box3d.xmax = rint((box3d.xmax - grid.ipx)/grid.xsize)
827                             * grid.xsize + grid.ipx;
828                box3d.ymin = rint((box3d.ymin - grid.ipy)/grid.ysize)
829                             * grid.ysize + grid.ipy;
830                box3d.ymax = rint((box3d.ymax - grid.ipy)/grid.ysize)
831                             * grid.ysize + grid.ipy;
832
833                out_lwgeom->bbox = box3d_to_box2df(&box3d);
834        }
835#endif /* 0 */
836
837        POSTGIS_DEBUGF(3, "SnapToGrid made a %s", lwgeom_typename(TYPE_GETTYPE(out_lwgeom->type)));
838
839        out_geom = pglwgeom_serialize(out_lwgeom);
840
841        PG_RETURN_POINTER(out_geom);
842}
843
844PG_FUNCTION_INFO_V1(LWGEOM_snaptogrid_pointoff);
845Datum LWGEOM_snaptogrid_pointoff(PG_FUNCTION_ARGS)
846{
847        Datum datum;
848        PG_LWGEOM *in_geom, *in_point;
849        LWGEOM *in_lwgeom;
850        LWPOINT *in_lwpoint;
851        PG_LWGEOM *out_geom = NULL;
852        LWGEOM *out_lwgeom;
853        gridspec grid;
854        /* BOX3D box3d; */
855        POINT4D offsetpoint;
856
857        if ( PG_ARGISNULL(0) ) PG_RETURN_NULL();
858        datum = PG_GETARG_DATUM(0);
859        in_geom = (PG_LWGEOM *)PG_DETOAST_DATUM(datum);
860
861        if ( PG_ARGISNULL(1) ) PG_RETURN_NULL();
862        datum = PG_GETARG_DATUM(1);
863        in_point = (PG_LWGEOM *)PG_DETOAST_DATUM(datum);
864        in_lwpoint = lwpoint_deserialize(SERIALIZED_FORM(in_point));
865        if ( in_lwpoint == NULL )
866        {
867                lwerror("Offset geometry must be a point");
868        }
869
870        if ( PG_ARGISNULL(2) ) PG_RETURN_NULL();
871        grid.xsize = PG_GETARG_FLOAT8(2);
872
873        if ( PG_ARGISNULL(3) ) PG_RETURN_NULL();
874        grid.ysize = PG_GETARG_FLOAT8(3);
875
876        if ( PG_ARGISNULL(4) ) PG_RETURN_NULL();
877        grid.zsize = PG_GETARG_FLOAT8(4);
878
879        if ( PG_ARGISNULL(5) ) PG_RETURN_NULL();
880        grid.msize = PG_GETARG_FLOAT8(5);
881
882        /* Take offsets from point geometry */
883        getPoint4d_p(in_lwpoint->point, 0, &offsetpoint);
884        grid.ipx = offsetpoint.x;
885        grid.ipy = offsetpoint.y;
886        if (TYPE_HASZ(in_lwpoint->type) ) grid.ipz = offsetpoint.z;
887        else grid.ipz=0;
888        if (TYPE_HASM(in_lwpoint->type) ) grid.ipm = offsetpoint.m;
889        else grid.ipm=0;
890
891#if POSTGIS_DEBUG_LEVEL >= 4
892        grid_print(&grid);
893#endif
894
895        /* Return input geometry if grid is null */
896        if ( grid_isNull(&grid) )
897        {
898                PG_RETURN_POINTER(in_geom);
899        }
900
901        in_lwgeom = lwgeom_deserialize(SERIALIZED_FORM(in_geom));
902
903        POSTGIS_DEBUGF(3, "SnapToGrid got a %s", lwgeom_typename(TYPE_GETTYPE(in_lwgeom->type)));
904
905        out_lwgeom = lwgeom_grid(in_lwgeom, &grid);
906        if ( out_lwgeom == NULL ) PG_RETURN_NULL();
907
908        /* COMPUTE_BBOX TAINTING */
909        if ( in_lwgeom->bbox ) lwgeom_add_bbox(out_lwgeom);
910
911#if 0
912        /*
913         * COMPUTE_BBOX WHEN SIMPLE
914         *
915         * WARNING: this is not SIMPLE, as snapping
916         * an existing bbox to a grid does not
917         * give the same result as computing a
918         * new bounding box on the snapped coordinates.
919         *
920         * This bug has been fixed in postgis-1.1.2
921         */
922        if ( in_lwgeom->bbox )
923        {
924                box2df_to_box3d_p(in_lwgeom->bbox, &box3d);
925
926                box3d.xmin = rint((box3d.xmin - grid.ipx)/grid.xsize)
927                             * grid.xsize + grid.ipx;
928                box3d.xmax = rint((box3d.xmax - grid.ipx)/grid.xsize)
929                             * grid.xsize + grid.ipx;
930                box3d.ymin = rint((box3d.ymin - grid.ipy)/grid.ysize)
931                             * grid.ysize + grid.ipy;
932                box3d.ymax = rint((box3d.ymax - grid.ipy)/grid.ysize)
933                             * grid.ysize + grid.ipy;
934
935                out_lwgeom->bbox = box3d_to_box2df(&box3d);
936        }
937#endif /* 0 */
938
939        POSTGIS_DEBUGF(3, "SnapToGrid made a %s", lwgeom_typename(TYPE_GETTYPE(out_lwgeom->type)));
940
941        out_geom = pglwgeom_serialize(out_lwgeom);
942
943        PG_RETURN_POINTER(out_geom);
944}
945
946
947/*
948** crossingDirection(line1, line2)
949**
950** Determines crossing direction of line2 relative to line1.
951** Only accepts LINESTRING ass parameters!
952*/
953PG_FUNCTION_INFO_V1(ST_LineCrossingDirection);
954Datum ST_LineCrossingDirection(PG_FUNCTION_ARGS)
955{
956        int type1, type2, rv;
957        BOX2DFLOAT4 box1, box2;
958        LWLINE *l1 = NULL;
959        LWLINE *l2 = NULL;
960        PG_LWGEOM *geom1 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
961        PG_LWGEOM *geom2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
962
963        errorIfSRIDMismatch(pglwgeom_getSRID(geom1), pglwgeom_getSRID(geom2));
964
965        /*
966        ** If the bounding boxes don't interact, then there can't be any
967        ** crossing, return right away.
968        */
969        if ( getbox2d_p(SERIALIZED_FORM(geom1), &box1) &&
970                getbox2d_p(SERIALIZED_FORM(geom2), &box2) )
971        {
972                if ( ( box2.xmax < box1.xmin ) || ( box2.xmin > box1.xmax ) ||
973                        ( box2.ymax < box1.ymin ) || ( box2.ymin > box1.ymax ) )
974                {
975                        PG_RETURN_INT32(LINE_NO_CROSS);
976                }
977        }
978
979        type1 = lwgeom_getType((uchar)SERIALIZED_FORM(geom1)[0]);
980        type2 = lwgeom_getType((uchar)SERIALIZED_FORM(geom2)[0]);
981
982        if ( type1 != LINETYPE || type2 != LINETYPE )
983        {
984                elog(ERROR,"This function only accepts LINESTRING as arguments.");
985                PG_RETURN_NULL();
986        }
987
988        l1 = lwline_deserialize(SERIALIZED_FORM(geom1));
989        l2 = lwline_deserialize(SERIALIZED_FORM(geom2));
990
991        rv = lwline_crossing_direction(l1, l2);
992
993        PG_FREE_IF_COPY(geom1, 0);
994        PG_FREE_IF_COPY(geom2, 0);
995
996        PG_RETURN_INT32(rv);
997
998}
999
1000PG_FUNCTION_INFO_V1(ST_LocateBetweenElevations);
1001Datum ST_LocateBetweenElevations(PG_FUNCTION_ARGS)
1002{
1003        PG_LWGEOM *geom_in = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1004        double from = PG_GETARG_FLOAT8(1);
1005        double to = PG_GETARG_FLOAT8(2);
1006        LWCOLLECTION *geom_out = NULL;
1007        LWGEOM *line_in = NULL;
1008        uchar type = (uchar)SERIALIZED_FORM(geom_in)[0];
1009        char geomtype = TYPE_GETTYPE(type);
1010        char hasz = TYPE_HASZ(type);
1011        static int ordinate = 2; /* Z */
1012
1013        if ( ! ( geomtype == LINETYPE || geomtype == MULTILINETYPE ) )
1014        {
1015                elog(ERROR,"This function only accepts LINESTRING or MULTILINESTRING as arguments.");
1016                PG_RETURN_NULL();
1017        }
1018
1019        if ( ! hasz )
1020        {
1021                elog(ERROR,"This function only accepts LINESTRING or MULTILINESTRING with Z values as arguments.");
1022                PG_RETURN_NULL();
1023        }
1024
1025        line_in = lwgeom_deserialize(SERIALIZED_FORM(geom_in));
1026        if ( geomtype == LINETYPE )
1027        {
1028                geom_out = lwline_clip_to_ordinate_range((LWLINE*)line_in, ordinate, from, to);
1029        }
1030        else if ( geomtype == MULTILINETYPE )
1031        {
1032                geom_out = lwmline_clip_to_ordinate_range((LWMLINE*)line_in, ordinate, from, to);
1033        }
1034        lwgeom_free(line_in);
1035
1036        if ( ! geom_out )
1037        {
1038                elog(ERROR,"The lwline_clip_to_ordinate_range returned null.");
1039                PG_RETURN_NULL();
1040        }
1041
1042        PG_FREE_IF_COPY(geom_in, 0);
1043        PG_RETURN_POINTER(pglwgeom_serialize((LWGEOM*)geom_out));
1044}
1045
1046/***********************************************************************
1047 * --strk@keybit.net
1048 ***********************************************************************/
1049
1050Datum LWGEOM_line_substring(PG_FUNCTION_ARGS);
1051
1052PG_FUNCTION_INFO_V1(LWGEOM_line_substring);
1053Datum LWGEOM_line_substring(PG_FUNCTION_ARGS)
1054{
1055        PG_LWGEOM *geom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1056        double from = PG_GETARG_FLOAT8(1);
1057        double to = PG_GETARG_FLOAT8(2);
1058        LWLINE *iline;
1059        LWGEOM *olwgeom;
1060        POINTARRAY *ipa, *opa;
1061        PG_LWGEOM *ret;
1062
1063        if ( from < 0 || from > 1 )
1064        {
1065                elog(ERROR,"line_interpolate_point: 2nd arg isnt within [0,1]");
1066                PG_RETURN_NULL();
1067        }
1068
1069        if ( to < 0 || to > 1 )
1070        {
1071                elog(ERROR,"line_interpolate_point: 3rd arg isnt within [0,1]");
1072                PG_RETURN_NULL();
1073        }
1074
1075        if ( from > to )
1076        {
1077                elog(ERROR, "2nd arg must be smaller then 3rd arg");
1078                PG_RETURN_NULL();
1079        }
1080
1081        if ( lwgeom_getType(geom->type) != LINETYPE )
1082        {
1083                elog(ERROR,"line_interpolate_point: 1st arg isnt a line");
1084                PG_RETURN_NULL();
1085        }
1086
1087        iline = lwline_deserialize(SERIALIZED_FORM(geom));
1088        ipa = iline->points;
1089
1090        opa = ptarray_substring(ipa, from, to);
1091
1092        if ( opa->npoints == 1 ) /* Point returned */
1093                olwgeom = (LWGEOM *)lwpoint_construct(iline->SRID, NULL, opa);
1094        else
1095                olwgeom = (LWGEOM *)lwline_construct(iline->SRID, NULL, opa);
1096
1097        ret = pglwgeom_serialize(olwgeom);
1098        PG_FREE_IF_COPY(geom, 0);
1099        lwgeom_release(olwgeom);
1100        PG_RETURN_POINTER(ret);
1101}
1102
1103Datum LWGEOM_line_locate_point(PG_FUNCTION_ARGS);
1104
1105PG_FUNCTION_INFO_V1(LWGEOM_line_locate_point);
1106Datum LWGEOM_line_locate_point(PG_FUNCTION_ARGS)
1107{
1108        PG_LWGEOM *geom1 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1109        PG_LWGEOM *geom2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
1110        LWLINE *lwline;
1111        LWPOINT *lwpoint;
1112        POINTARRAY *pa;
1113        POINT2D p;
1114        double ret;
1115
1116        if ( lwgeom_getType(geom1->type) != LINETYPE )
1117        {
1118                elog(ERROR,"line_locate_point: 1st arg isnt a line");
1119                PG_RETURN_NULL();
1120        }
1121        if ( lwgeom_getType(geom2->type) != POINTTYPE )
1122        {
1123                elog(ERROR,"line_locate_point: 2st arg isnt a point");
1124                PG_RETURN_NULL();
1125        }
1126        if ( pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2) )
1127        {
1128                elog(ERROR, "Operation on two geometries with different SRIDs");
1129                PG_RETURN_NULL();
1130        }
1131
1132        lwline = lwline_deserialize(SERIALIZED_FORM(geom1));
1133        lwpoint = lwpoint_deserialize(SERIALIZED_FORM(geom2));
1134
1135        pa = lwline->points;
1136        lwpoint_getPoint2d_p(lwpoint, &p);
1137
1138        ret = ptarray_locate_point(pa, &p);
1139
1140        PG_RETURN_FLOAT8(ret);
1141}
1142
1143/*******************************************************************************
1144 * The following is based on the "Fast Winding Number Inclusion of a Point
1145 * in a Polygon" algorithm by Dan Sunday.
1146 * http://softsurfer.com/Archive/algorithm_0103/algorithm_0103.htm#Winding%20Number
1147 ******************************************************************************/
1148
1149/*
1150 * returns: >0 for a point to the left of the segment,
1151 *          <0 for a point to the right of the segment,
1152 *          0 for a point on the segment
1153 */
1154double determineSide(POINT2D *seg1, POINT2D *seg2, POINT2D *point)
1155{
1156        return ((seg2->x-seg1->x)*(point->y-seg1->y)-(point->x-seg1->x)*(seg2->y-seg1->y));
1157}
1158
1159/*
1160 * This function doesn't test that the point falls on the line defined by
1161 * the two points.  It assumes that that has already been determined
1162 * by having determineSide return within the tolerance.  It simply checks
1163 * that if the point is on the line, it is within the endpoints.
1164 *
1165 * returns: 1 if the point is not outside the bounds of the segment
1166 *          0 if it is
1167 */
1168int isOnSegment(POINT2D *seg1, POINT2D *seg2, POINT2D *point)
1169{
1170        double maxX;
1171        double maxY;
1172        double minX;
1173        double minY;
1174
1175        if (seg1->x > seg2->x)
1176        {
1177                maxX = seg1->x;
1178                minX = seg2->x;
1179        }
1180        else
1181        {
1182                maxX = seg2->x;
1183                minX = seg1->x;
1184        }
1185        if (seg1->y > seg2->y)
1186        {
1187                maxY = seg1->y;
1188                minY = seg2->y;
1189        }
1190        else
1191        {
1192                maxY = seg2->y;
1193                minY = seg1->y;
1194        }
1195
1196        LWDEBUGF(3, "maxX minX/maxY minY: %.8f %.8f/%.8f %.8f", maxX, minX, maxY, minY);
1197
1198        if (maxX < point->x || minX > point->x)
1199        {
1200                LWDEBUGF(3, "X value %.8f falls outside the range %.8f-%.8f", point->x, minX, maxX);
1201
1202                return 0;
1203        }
1204        else if (maxY < point->y || minY > point->y)
1205        {
1206                LWDEBUGF(3, "Y value %.8f falls outside the range %.8f-%.8f", point->y, minY, maxY);
1207
1208                return 0;
1209        }
1210        return 1;
1211}
1212
1213/*
1214 * return -1 iff point is outside ring pts
1215 * return 1 iff point is inside ring pts
1216 * return 0 iff point is on ring pts
1217 */
1218int point_in_ring_rtree(RTREE_NODE *root, POINT2D *point)
1219{
1220        int wn = 0;
1221        int i;
1222        double side;
1223        POINT2D seg1;
1224        POINT2D seg2;
1225        LWMLINE *lines;
1226
1227        LWDEBUG(2, "point_in_ring called.");
1228
1229        lines = findLineSegments(root, point->y);
1230        if (!lines)
1231                return -1;
1232
1233        for (i=0; i<lines->ngeoms; i++)
1234        {
1235                getPoint2d_p(lines->geoms[i]->points, 0, &seg1);
1236                getPoint2d_p(lines->geoms[i]->points, 1, &seg2);
1237
1238
1239                side = determineSide(&seg1, &seg2, point);
1240
1241                LWDEBUGF(3, "segment: (%.8f, %.8f),(%.8f, %.8f)", seg1.x, seg1.y, seg2.x, seg2.y);
1242                LWDEBUGF(3, "side result: %.8f", side);
1243                LWDEBUGF(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));
1244
1245                /* zero length segments are ignored. */
1246                if (((seg2.x-seg1.x)*(seg2.x-seg1.x)+(seg2.y-seg1.y)*(seg2.y-seg1.y)) < 1e-12*1e-12)
1247                {
1248                        LWDEBUG(3, "segment is zero length... ignoring.");
1249
1250                        continue;
1251                }
1252
1253                /* a point on the boundary of a ring is not contained. */
1254                if (fabs(side) < 1e-12)
1255                {
1256                        if (isOnSegment(&seg1, &seg2, point) == 1)
1257                        {
1258                                LWDEBUGF(3, "point on ring boundary between points %d, %d", i, i+1);
1259
1260                                return 0;
1261                        }
1262                }
1263                /*
1264                 * If the point is to the left of the line, and it's rising,
1265                 * then the line is to the right of the point and
1266                 * circling counter-clockwise, so incremement.
1267                 */
1268                else if (FP_CONTAINS_BOTTOM(seg1.y,point->y,seg2.y) && side>0)
1269                {
1270                        LWDEBUG(3, "incrementing winding number.");
1271
1272                        ++wn;
1273                }
1274                /*
1275                 * If the point is to the right of the line, and it's falling,
1276                 * then the line is to the right of the point and circling
1277                 * clockwise, so decrement.
1278                 */
1279                else if (FP_CONTAINS_BOTTOM(seg2.y,point->y,seg1.y) && side<0)
1280                {
1281                        LWDEBUG(3, "decrementing winding number.");
1282
1283                        --wn;
1284                }
1285        }
1286
1287        LWDEBUGF(3, "winding number %d", wn);
1288
1289        if (wn == 0)
1290                return -1;
1291        return 1;
1292}
1293
1294
1295/*
1296 * return -1 iff point is outside ring pts
1297 * return 1 iff point is inside ring pts
1298 * return 0 iff point is on ring pts
1299 */
1300int point_in_ring(POINTARRAY *pts, POINT2D *point)
1301{
1302        int wn = 0;
1303        int i;
1304        double side;
1305        POINT2D seg1;
1306        POINT2D seg2;
1307
1308        LWDEBUG(2, "point_in_ring called.");
1309
1310
1311        for (i=0; i<pts->npoints-1; i++)
1312        {
1313                getPoint2d_p(pts, i, &seg1);
1314                getPoint2d_p(pts, i+1, &seg2);
1315
1316
1317                side = determineSide(&seg1, &seg2, point);
1318
1319                LWDEBUGF(3, "segment: (%.8f, %.8f),(%.8f, %.8f)", seg1.x, seg1.y, seg2.x, seg2.y);
1320                LWDEBUGF(3, "side result: %.8f", side);
1321                LWDEBUGF(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));
1322
1323                /* zero length segments are ignored. */
1324                if (((seg2.x-seg1.x)*(seg2.x-seg1.x)+(seg2.y-seg1.y)*(seg2.y-seg1.y)) < 1e-12*1e-12)
1325                {
1326                        LWDEBUG(3, "segment is zero length... ignoring.");
1327
1328                        continue;
1329                }
1330
1331                /* a point on the boundary of a ring is not contained. */
1332                if (fabs(side) < 1e-12)
1333                {
1334                        if (isOnSegment(&seg1, &seg2, point) == 1)
1335                        {
1336                                LWDEBUGF(3, "point on ring boundary between points %d, %d", i, i+1);
1337
1338                                return 0;
1339                        }
1340                }
1341                /*
1342                 * If the point is to the left of the line, and it's rising,
1343                 * then the line is to the right of the point and
1344                 * circling counter-clockwise, so incremement.
1345                 */
1346                else if (FP_CONTAINS_BOTTOM(seg1.y,point->y,seg2.y) && side>0)
1347                {
1348                        LWDEBUG(3, "incrementing winding number.");
1349
1350                        ++wn;
1351                }
1352                /*
1353                 * If the point is to the right of the line, and it's falling,
1354                 * then the line is to the right of the point and circling
1355                 * clockwise, so decrement.
1356                 */
1357                else if (FP_CONTAINS_BOTTOM(seg2.y,point->y,seg1.y) && side<0)
1358                {
1359                        LWDEBUG(3, "decrementing winding number.");
1360
1361                        --wn;
1362                }
1363        }
1364
1365        LWDEBUGF(3, "winding number %d", wn);
1366
1367        if (wn == 0)
1368                return -1;
1369        return 1;
1370}
1371
1372/*
1373 * return 0 iff point outside polygon or on boundary
1374 * return 1 iff point inside polygon
1375 */
1376int point_in_polygon_rtree(RTREE_NODE **root, int ringCount, LWPOINT *point)
1377{
1378        int i;
1379        POINT2D pt;
1380
1381        LWDEBUGF(2, "point_in_polygon called for %p %d %p.", root, ringCount, point);
1382
1383        getPoint2d_p(point->point, 0, &pt);
1384        /* assume bbox short-circuit has already been attempted */
1385
1386        if (point_in_ring_rtree(root[0], &pt) != 1)
1387        {
1388                LWDEBUG(3, "point_in_polygon_rtree: outside exterior ring.");
1389
1390                return 0;
1391        }
1392
1393        for (i=1; i<ringCount; i++)
1394        {
1395                if (point_in_ring_rtree(root[i], &pt) != -1)
1396                {
1397                        LWDEBUGF(3, "point_in_polygon_rtree: within hole %d.", i);
1398
1399                        return 0;
1400                }
1401        }
1402        return 1;
1403}
1404
1405/*
1406 * return -1 if point outside polygon
1407 * return 0 if point on boundary
1408 * return 1 if point inside polygon
1409 *
1410 * Expected **root order is all the exterior rings first, then all the holes
1411 *
1412 * TODO: this could be made slightly more efficient by ordering the rings in
1413 * EIIIEIIIEIEI order (exterior/interior) and including list of exterior ring
1414 * positions on the cache object.
1415 */
1416int point_in_multipolygon_rtree(RTREE_NODE **root, int polyCount, int ringCount, LWPOINT *point)
1417{
1418        int i;
1419        POINT2D pt;
1420        int result = -1;
1421
1422        LWDEBUGF(2, "point_in_multipolygon_rtree called for %p %d %d %p.", root, polyCount, ringCount, point);
1423
1424        getPoint2d_p(point->point, 0, &pt);
1425        /* assume bbox short-circuit has already been attempted */
1426
1427        /* is the point inside (not outside) any of the exterior rings? */
1428        for ( i = 0; i < polyCount; i++ )
1429        {
1430                int in_ring = point_in_ring_rtree(root[i], &pt);
1431                LWDEBUGF(4, "point_in_multipolygon_rtree: exterior ring (%d), point_in_ring returned %d", i, in_ring);
1432                if ( in_ring != -1 ) /* not outside this ring */
1433                {
1434                        LWDEBUG(3, "point_in_multipolygon_rtree: inside exterior ring.");
1435                        result = in_ring;
1436                        break;
1437                }
1438        }
1439
1440        if ( result == -1 ) /* strictly outside all rings */
1441                return result;
1442
1443        /* ok, it's in a ring, but if it's in a hole it's still outside */
1444        for ( i = polyCount; i < ringCount; i++ )
1445        {
1446                int in_ring = point_in_ring_rtree(root[i], &pt);
1447                LWDEBUGF(4, "point_in_multipolygon_rtree: hole (%d), point_in_ring returned %d", i, in_ring);
1448                if ( in_ring == 1 ) /* completely inside hole */
1449                {
1450                        LWDEBUGF(3, "point_in_multipolygon_rtree: within hole %d.", i);
1451                        return -1;
1452                }
1453                if ( in_ring == 0 ) /* on the boundary of a hole */
1454                {
1455                        result = 0;
1456                }
1457        }
1458        return result; /* -1 = outside, 0 = boundary, 1 = inside */
1459
1460}
1461
1462/*
1463 * return -1 iff point outside polygon
1464 * return 0 iff point on boundary
1465 * return 1 iff point inside polygon
1466 */
1467int point_in_polygon(LWPOLY *polygon, LWPOINT *point)
1468{
1469        int i, result, in_ring;
1470        POINTARRAY *ring;
1471        POINT2D pt;
1472
1473        LWDEBUG(2, "point_in_polygon called.");
1474
1475        getPoint2d_p(point->point, 0, &pt);
1476        /* assume bbox short-circuit has already been attempted */
1477
1478        ring = polygon->rings[0];
1479        in_ring = point_in_ring(polygon->rings[0], &pt);
1480        if ( in_ring == -1) /* outside the exterior ring */
1481        {
1482                LWDEBUG(3, "point_in_polygon: outside exterior ring.");
1483                return -1;
1484        }
1485        result = in_ring;
1486
1487        for (i=1; i<polygon->nrings; i++)
1488        {
1489                ring = polygon->rings[i];
1490                in_ring = point_in_ring(polygon->rings[i], &pt);
1491                if (in_ring == 1) /* inside a hole => outside the polygon */
1492                {
1493                        LWDEBUGF(3, "point_in_polygon: within hole %d.", i);
1494                        return -1;
1495                }
1496                if (in_ring == 0) /* on the edge of a hole */
1497                {
1498                        LWDEBUGF(3, "point_in_polygon: on edge of hole %d.", i);
1499                        return 0;
1500                }
1501        }
1502        return result; /* -1 = outside, 0 = boundary, 1 = inside */
1503}
1504
1505/*
1506 * return -1 iff point outside multipolygon
1507 * return 0 iff point on multipolygon boundary
1508 * return 1 iff point inside multipolygon
1509 */
1510int point_in_multipolygon(LWMPOLY *mpolygon, LWPOINT *point)
1511{
1512        int i, j, result, in_ring;
1513        POINTARRAY *ring;
1514        POINT2D pt;
1515
1516        LWDEBUG(2, "point_in_polygon called.");
1517
1518        getPoint2d_p(point->point, 0, &pt);
1519        /* assume bbox short-circuit has already been attempted */
1520
1521        result = -1;
1522
1523        for (j = 0; j < mpolygon->ngeoms; j++ )
1524        {
1525
1526                LWPOLY *polygon = mpolygon->geoms[j];
1527                ring = polygon->rings[0];
1528                in_ring = point_in_ring(polygon->rings[0], &pt);
1529                if ( in_ring == -1) /* outside the exterior ring */
1530                {
1531                        LWDEBUG(3, "point_in_polygon: outside exterior ring.");
1532                        continue;
1533                }
1534                if ( in_ring == 0 )
1535                {
1536                        return 0;
1537                }
1538
1539                result = in_ring;
1540
1541                for (i=1; i<polygon->nrings; i++)
1542                {
1543                        ring = polygon->rings[i];
1544                        in_ring = point_in_ring(polygon->rings[i], &pt);
1545                        if (in_ring == 1) /* inside a hole => outside the polygon */
1546                        {
1547                                LWDEBUGF(3, "point_in_polygon: within hole %d.", i);
1548                                result = -1;
1549                                break;
1550                        }
1551                        if (in_ring == 0) /* on the edge of a hole */
1552                        {
1553                                LWDEBUGF(3, "point_in_polygon: on edge of hole %d.", i);
1554                                return 0;
1555                        }
1556                }
1557                if ( result != -1)
1558                {
1559                        return result;
1560                }
1561        }
1562        return result;
1563}
1564
1565
1566/*******************************************************************************
1567 * End of "Fast Winding Number Inclusion of a Point in a Polygon" derivative.
1568 ******************************************************************************/
1569
Note: See TracBrowser for help on using the repository browser.