source: trunk/postgis/lwgeom_functions_lrs.c @ 4745

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

Make non-M attempts to run LRS functions error out instead of return NULL (#113)

  • 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: 12.3 KB
Line 
1/**********************************************************************
2 * $Id: lwgeom_functions_lrs.c 4745 2009-11-05 00:37:32Z 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
19
20Datum LWGEOM_locate_between_m(PG_FUNCTION_ARGS);
21
22typedef struct
23{
24        POINTARRAY **ptarrays;
25        uint32 nptarrays;
26}
27POINTARRAYSET;
28
29static POINTARRAYSET ptarray_locate_between_m(
30    POINTARRAY *ipa, double m0, double m1);
31
32static LWGEOM *lwcollection_locate_between_m(
33    LWCOLLECTION *lwcoll, double m0, double m1);
34
35static LWGEOM *lwgeom_locate_between_m(
36    LWGEOM *lwin, double m0, double m1);
37
38static LWGEOM *lwline_locate_between_m(
39    LWLINE *lwline_in, double m0, double m1);
40
41static LWGEOM *lwpoint_locate_between_m(
42    LWPOINT *lwpoint, double m0, double m1);
43
44static int clip_seg_by_m_range(
45    POINT4D *p1, POINT4D *p2, double m0, double m1);
46
47
48/*
49 * Clip a segment by a range of measures.
50 * Z and M values are interpolated in case of clipping.
51 *
52 * Returns a bitfield, flags being:
53 *    0x0001 : segment intersects the range
54 *    0x0010 : first point is modified
55 *    0x0100 : second point is modified
56 *
57 * Values:
58 *     - 0 segment fully outside the range, no modifications
59 *     - 1 segment fully inside the range, no modifications
60 *     - 7 segment crosses the range, both points modified.
61 *     - 3 first point out, second in, first point modified
62 *     - 5 first point in, second out, second point modified
63 */
64static int
65clip_seg_by_m_range(POINT4D *p1, POINT4D *p2, double m0, double m1)
66{
67        double dM0, dM1, dX, dY, dZ;
68        POINT4D *tmp;
69        int swapped=0;
70        int ret=0;
71
72        LWDEBUGF(3, "m0: %g m1: %g", m0, m1);
73
74        /* Handle corner case of m values being the same */
75        if ( p1->m == p2->m )
76        {
77                /* out of range, no clipping */
78                if ( p1->m < m0 || p1->m > m1 )
79                        return 0;
80
81                /* inside range, no clipping */
82                return 1;
83        }
84
85        /*
86         * Order points so that p1 has the smaller M
87         */
88        if ( p1->m > p2->m )
89        {
90                tmp=p2;
91                p2=p1;
92                p1=tmp;
93                swapped=1;
94        }
95
96        /*
97         * The M range is not intersected, segment
98         * fully out of range, no clipping.
99         */
100        if ( p2->m < m0 || p1->m > m1 )
101                return 0;
102
103        /*
104         * The segment is fully inside the range,
105         * no clipping.
106         */
107        if ( p1->m >= m0 && p2->m <= m1 )
108                return 1;
109
110        /*
111         * Segment intersects range, lets compute
112         * the proportional location of the two
113         * measures wrt p1/p2 m range.
114         *
115         * if p1 and p2 have the same measure
116         * this should never be reached (either
117         * both inside or both outside)
118         *
119         */
120        dM0=(m0-p1->m)/(p2->m-p1->m); /* delta-M0 */
121        dM1=(m1-p2->m)/(p2->m-p1->m); /* delta-M1 */
122        dX=p2->x-p1->x;
123        dY=p2->y-p1->y;
124        dZ=p2->z-p1->z;
125
126        LWDEBUGF(3, "dM0:%g dM1:%g", dM0, dM1);
127        LWDEBUGF(3, "dX:%g dY:%g dZ:%g", dX, dY, dZ);
128        LWDEBUGF(3, "swapped: %d", swapped);
129
130        /*
131         * First point out of range, project
132         * it on the range
133         */
134        if ( p1->m < m0 )
135        {
136                /*
137                 * To prevent rounding errors, then if m0==m1 and p2 lies within the range, copy
138                 * p1 as a direct copy of p2
139                 */
140                if (m0 == m1 && p2->m <= m1)
141                {
142                        memcpy(p1, p2, sizeof(POINT4D));
143
144                        LWDEBUG(3, "Projected p1 on range (as copy of p2)");
145                }
146                else
147                {
148                        /* Otherwise interpolate coordinates */
149                        p1->x += (dX*dM0);
150                        p1->y += (dY*dM0);
151                        p1->z += (dZ*dM0);
152                        p1->m = m0;
153
154                        LWDEBUG(3, "Projected p1 on range");
155                }
156
157                if ( swapped ) ret |= 0x0100;
158                else ret |= 0x0010;
159        }
160
161        /*
162         * Second point out of range, project
163         * it on the range
164         */
165        if ( p2->m > m1 )
166        {
167                /*
168                 * To prevent rounding errors, then if m0==m1 and p1 lies within the range, copy
169                 * p2 as a direct copy of p1
170                 */
171                if (m0 == m1 && p1->m >= m0)
172                {
173                        memcpy(p2, p1, sizeof(POINT4D));
174
175                        LWDEBUG(3, "Projected p2 on range (as copy of p1)");
176                }
177                else
178                {
179                        /* Otherwise interpolate coordinates */
180                        p2->x += (dX*dM1);
181                        p2->y += (dY*dM1);
182                        p2->z += (dZ*dM1);
183                        p2->m = m1;
184
185                        LWDEBUG(3, "Projected p2 on range");
186                }
187
188                if ( swapped ) ret |= 0x0010;
189                else ret |= 0x0100;
190        }
191
192        /* Clipping occurred */
193        return ret;
194
195}
196
197static POINTARRAYSET
198ptarray_locate_between_m(POINTARRAY *ipa, double m0, double m1)
199{
200        POINTARRAYSET ret;
201        DYNPTARRAY *dpa=NULL;
202        int i;
203
204        ret.nptarrays=0;
205
206        /*
207         * We allocate space for as many pointarray as
208         * segments in the input POINTARRAY, as worst
209         * case is for each segment to cross the M range
210         * window.
211         * TODO: rework this to reduce used memory
212         */
213        ret.ptarrays=lwalloc(sizeof(POINTARRAY *)*ipa->npoints-1);
214
215        LWDEBUGF(2, "ptarray_locate...: called for pointarray %x, m0:%g, m1:%g",
216                 ipa, m0, m1);
217
218
219        for (i=1; i<ipa->npoints; i++)
220        {
221                POINT4D p1, p2;
222                int clipval;
223
224                getPoint4d_p(ipa, i-1, &p1);
225                getPoint4d_p(ipa, i, &p2);
226
227                LWDEBUGF(3, " segment %d-%d [ %g %g %g %g -  %g %g %g %g ]",
228                         i-1, i,
229                         p1.x, p1.y, p1.z, p1.m,
230                         p2.x, p2.y, p2.z, p2.m);
231
232                clipval = clip_seg_by_m_range(&p1, &p2, m0, m1);
233
234                /* segment completely outside, nothing to do */
235                if (! clipval ) continue;
236
237                LWDEBUGF(3, " clipped to: [ %g %g %g %g - %g %g %g %g ]   clipval: %x", p1.x, p1.y, p1.z, p1.m,
238                         p2.x, p2.y, p2.z, p2.m, clipval);
239
240                /* If no points have been accumulated so far, then if clipval != 0 the first point must be the
241                   start of the intersection */
242                if (dpa == NULL)
243                {
244                        LWDEBUGF(3, " 1 creating new POINTARRAY with first point %g,%g,%g,%g", p1.x, p1.y, p1.z, p1.m);
245
246                        dpa = dynptarray_create(ipa->npoints-i, ipa->dims);
247                        dynptarray_addPoint4d(dpa, &p1, 1);
248                }
249
250                /* Otherwise always add the next point, avoiding duplicates */
251                if (dpa)
252                        dynptarray_addPoint4d(dpa, &p2, 0);
253
254                /*
255                 * second point has been clipped
256                 */
257                if ( clipval & 0x0100 || i == ipa->npoints-1 )
258                {
259                        LWDEBUGF(3, " closing pointarray %x with %d points", dpa->pa, dpa->pa->npoints);
260
261                        ret.ptarrays[ret.nptarrays++] = dpa->pa;
262                        lwfree(dpa);
263                        dpa=NULL;
264                }
265        }
266
267        /*
268         * if dpa!=NULL it means we didn't close it yet.
269         * this should never happen.
270         */
271        if ( dpa != NULL ) lwerror("Something wrong with algorithm");
272
273        return ret;
274}
275
276/*
277 * Point is assumed to have an M value.
278 * Return NULL if point is not in the given range (inclusive)
279 * Return an LWPOINT *copy* otherwise.
280 */
281static LWGEOM *
282lwpoint_locate_between_m(LWPOINT *lwpoint, double m0, double m1)
283{
284        POINT3DM p3dm;
285
286        LWDEBUGF(2, "lwpoint_locate_between called for lwpoint %x", lwpoint);
287
288        lwpoint_getPoint3dm_p(lwpoint, &p3dm);
289        if ( p3dm.m >= m0 && p3dm.m <= m1)
290        {
291                LWDEBUG(3, " lwpoint... returning a clone of input");
292
293                return (LWGEOM *)lwpoint_clone(lwpoint);
294        }
295        else
296        {
297                LWDEBUG(3, " lwpoint... returning a clone of input");
298
299                return NULL;
300        }
301}
302
303/*
304 * Line is assumed to have an M value.
305 *
306 * Return NULL if no parts of the line are in the given range (inclusive)
307 *
308 * Return an LWCOLLECTION with LWLINES and LWPOINT being consecutive
309 * and isolated points on the line falling in the range.
310 *
311 * X,Y and Z (if present) ordinates are interpolated.
312 *
313 */
314static LWGEOM *
315lwline_locate_between_m(LWLINE *lwline_in, double m0, double m1)
316{
317        POINTARRAY *ipa=lwline_in->points;
318        int i;
319        LWGEOM **geoms;
320        int ngeoms;
321        int outtype;
322        int typeflag=0; /* see flags below */
323        const int pointflag=0x01;
324        const int lineflag=0x10;
325        POINTARRAYSET paset=ptarray_locate_between_m(ipa, m0, m1);
326
327        LWDEBUGF(2, "lwline_locate_between called for lwline %x", lwline_in);
328
329        LWDEBUGF(3, " ptarray_locate... returned %d pointarrays",
330                 paset.nptarrays);
331
332        if ( paset.nptarrays == 0 )
333        {
334                return NULL;
335        }
336
337        ngeoms=paset.nptarrays;
338        /* TODO: rework this to reduce used memory */
339        geoms=lwalloc(sizeof(LWGEOM *)*ngeoms);
340        for (i=0; i<ngeoms; i++)
341        {
342                LWPOINT *lwpoint;
343                LWLINE *lwline;
344
345                POINTARRAY *pa=paset.ptarrays[i];
346
347                /* This is a point */
348                if ( pa->npoints == 1 )
349                {
350                        lwpoint=lwalloc(sizeof(LWPOINT));
351                        lwpoint->type=lwgeom_makeType_full(
352                                          TYPE_HASZ(pa->dims),
353                                          TYPE_HASM(pa->dims),
354                                          lwline_in->SRID,
355                                          POINTTYPE,
356                                          0);
357                        lwpoint->SRID=lwline_in->SRID;
358                        lwpoint->bbox=NULL;
359                        lwpoint->point=pa;
360                        geoms[i]=(LWGEOM *)lwpoint;
361                        typeflag|=pointflag;
362                }
363
364                /* This is a line */
365                else if ( pa->npoints > 1 )
366                {
367                        lwline=lwalloc(sizeof(LWLINE));
368                        lwline->type=lwgeom_makeType_full(
369                                         TYPE_HASZ(pa->dims),
370                                         TYPE_HASM(pa->dims),
371                                         lwline_in->SRID,
372                                         LINETYPE,
373                                         0);
374                        lwline->SRID=lwline_in->SRID;
375                        lwline->bbox=NULL;
376                        lwline->points=pa;
377                        geoms[i]=(LWGEOM *)lwline;
378                        typeflag|=lineflag;
379                }
380
381                /* This is a bug */
382                else
383                {
384                        lwerror("ptarray_locate_between_m returned a POINARRAY set containing POINTARRAY with 0 points");
385                }
386
387        }
388
389        if ( ngeoms == 1 )
390        {
391                return geoms[0];
392        }
393        else
394        {
395                /* Choose best type */
396                if ( typeflag == 1 ) outtype=MULTIPOINTTYPE;
397                else if ( typeflag == 2 ) outtype=MULTILINETYPE;
398                else outtype = COLLECTIONTYPE;
399
400                return (LWGEOM *)lwcollection_construct(outtype,
401                                                        lwline_in->SRID, NULL, ngeoms, geoms);
402        }
403}
404
405/*
406 * Return a fully new allocated LWCOLLECTION
407 * always tagged as COLLECTIONTYPE.
408 */
409static LWGEOM *
410lwcollection_locate_between_m(LWCOLLECTION *lwcoll, double m0, double m1)
411{
412        int i;
413        int ngeoms=0;
414        LWGEOM **geoms;
415
416        LWDEBUGF(2, "lwcollection_locate_between_m called for lwcoll %x", lwcoll);
417
418        geoms=lwalloc(sizeof(LWGEOM *)*lwcoll->ngeoms);
419        for (i=0; i<lwcoll->ngeoms; i++)
420        {
421                LWGEOM *sub=lwgeom_locate_between_m(lwcoll->geoms[i],
422                                                    m0, m1);
423                if ( sub != NULL )
424                        geoms[ngeoms++] = sub;
425        }
426
427        if ( ngeoms == 0 ) return NULL;
428
429        return (LWGEOM *)lwcollection_construct(COLLECTIONTYPE,
430                                                lwcoll->SRID, NULL, ngeoms, geoms);
431}
432
433/*
434 * Return a fully allocated LWGEOM containing elements
435 * intersected/interpolated with the given M range.
436 * Return NULL if none of the elements fall in the range.
437 *
438 * m0 is assumed to be less-or-equal to m1.
439 * input LWGEOM is assumed to contain an M value.
440 *
441 */
442static LWGEOM *
443lwgeom_locate_between_m(LWGEOM *lwin, double m0, double m1)
444{
445        LWDEBUGF(2, "lwgeom_locate_between called for lwgeom %x", lwin);
446
447        switch (TYPE_GETTYPE(lwin->type))
448        {
449        case POINTTYPE:
450                return lwpoint_locate_between_m(
451                           (LWPOINT *)lwin, m0, m1);
452        case LINETYPE:
453                return lwline_locate_between_m(
454                           (LWLINE *)lwin, m0, m1);
455
456        case MULTIPOINTTYPE:
457        case MULTILINETYPE:
458        case COLLECTIONTYPE:
459                return lwcollection_locate_between_m(
460                           (LWCOLLECTION *)lwin, m0, m1);
461
462                /* Polygon types are not supported */
463        case POLYGONTYPE:
464        case MULTIPOLYGONTYPE:
465                lwerror("Areal geometries are not supported by locate_between_measures");
466                return NULL;
467        }
468
469        lwerror("Unkonwn geometry type (%s:%d)", __FILE__, __LINE__);
470        return NULL;
471}
472
473/*
474 * Return a derived geometry collection value with elements that match
475 * the specified range of measures inclusively.
476 *
477 * Implements SQL/MM ST_LocateBetween(measure, measure) method.
478 * See ISO/IEC CD 13249-3:200x(E)
479 *
480 */
481PG_FUNCTION_INFO_V1(LWGEOM_locate_between_m);
482Datum LWGEOM_locate_between_m(PG_FUNCTION_ARGS)
483{
484        PG_LWGEOM *gin = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
485        PG_LWGEOM *gout;
486        double start_measure = PG_GETARG_FLOAT8(1);
487        double end_measure = PG_GETARG_FLOAT8(2);
488        LWGEOM *lwin, *lwout;
489        int type;
490
491        if ( end_measure < start_measure )
492        {
493                lwerror("locate_between_m: 2nd arg must be bigger then 1st arg");
494                PG_RETURN_NULL();
495        }
496
497        /*
498         * Return error if input doesn't have a measure
499         */
500        if ( ! lwgeom_hasM(gin->type) )
501        {
502                lwerror("Geometry argument does not have an 'M' ordinate");
503                PG_RETURN_NULL();
504        }
505
506        /*
507         * Raise an error if input is a polygon, a multipolygon
508         * or a collection
509         */
510        type=lwgeom_getType(gin->type);
511               
512        if ( type == POLYGONTYPE || type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE )
513        {
514                lwerror("Areal or Collection types are not supported");
515                PG_RETURN_NULL();
516        }
517
518        lwin = pglwgeom_deserialize(gin);
519
520        lwout = lwgeom_locate_between_m(lwin,
521                                        start_measure, end_measure);
522
523        lwgeom_release(lwin);
524
525        if ( lwout == NULL )
526        {
527                lwout = (LWGEOM *)lwcollection_construct_empty(
528                            pglwgeom_getSRID(gin), lwgeom_hasZ(gin->type),
529                            lwgeom_hasM(gin->type));
530        }
531
532        gout = pglwgeom_serialize(lwout);
533        lwgeom_release(lwout);
534
535        PG_RETURN_POINTER(gout);
536}
537
Note: See TracBrowser for help on using the repository browser.