1
/**********************************************************************
2
* $Id: lwgeom_functions_lrs.c 3358 2008-12-03 16:27:44Z mcayland $
4
* PostGIS - Spatial Types for PostgreSQL
5
* http://postgis.refractions.net
6
* Copyright 2001-2005 Refractions Research Inc.
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.
11
**********************************************************************/
15
#include "liblwgeom.h"
16
#include "lwgeom_pg.h"
20
#define DEBUG_INTERPOLATION 0
22
Datum LWGEOM_locate_between_m(PG_FUNCTION_ARGS);
25
POINTARRAY **ptarrays;
29
static POINTARRAYSET ptarray_locate_between_m(
30
POINTARRAY *ipa, double m0, double m1);
32
static LWGEOM *lwcollection_locate_between_m(
33
LWCOLLECTION *lwcoll, double m0, double m1);
35
static LWGEOM *lwgeom_locate_between_m(
36
LWGEOM *lwin, double m0, double m1);
38
static LWGEOM *lwline_locate_between_m(
39
LWLINE *lwline_in, double m0, double m1);
41
static LWGEOM *lwpoint_locate_between_m(
42
LWPOINT *lwpoint, double m0, double m1);
44
static int clip_seg_by_m_range(
45
POINT4D *p1, POINT4D *p2, double m0, double m1);
49
* Clip a segment by a range of measures.
50
* Z and M values are interpolated in case of clipping.
52
* Returns a bitfield, flags being:
53
* 0x0001 : segment intersects the range
54
* 0x0010 : first point is modified
55
* 0x0100 : second point is modified
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
65
clip_seg_by_m_range(POINT4D *p1, POINT4D *p2, double m0, double m1)
67
double dM0, dM1, dX, dY, dZ;
72
/* Handle corner case of m values being the same */
75
/* out of range, no clipping */
76
if ( p1->m < m0 || p1->m > m1 )
79
/* inside range, no clipping */
84
* Order points so that p1 has the smaller M
95
* The M range is not intersected, segment
96
* fully out of range, no clipping.
98
if ( p2->m < m0 || p1->m > m1 )
102
* The segment is fully inside the range,
105
if ( p1->m >= m0 && p2->m <= m1 )
109
* Segment intersects range, lets compute
110
* the proportional location of the two
111
* measures wrt p1/p2 m range.
113
* if p1 and p2 have the same measure
114
* this should never be reached (either
115
* both inside or both outside)
118
dM0=(m0-p1->m)/(p2->m-p1->m); /* delta-M0 */
119
dM1=(m1-p2->m)/(p2->m-p1->m); /* delta-M1 */
123
#if DEBUG_INTERPOLATION
124
lwnotice("dM0:%g dM1:%g", dM0, dM1);
125
lwnotice("dX:%g dY:%g dZ:%g", dX, dY, dZ);
130
* First point out of range, project
136
* To prevent rounding errors, then if m0==m1 and p2 lies within the range, copy
137
* p1 as a direct copy of p2
139
if (m0 == m1 && p2->m <= m1)
141
memcpy(p1, p2, sizeof(POINT4D));
143
#if DEBUG_INTERPOLATION
144
lwnotice("Projected p1 on range (as copy of p2)");
149
/* Otherwise interpolate coordinates */
155
#if DEBUG_INTERPOLATION
156
lwnotice("Projected p1 on range");
160
if ( swapped ) ret |= 0x0100;
165
* Second point out of range, project
171
* To prevent rounding errors, then if m0==m1 and p1 lies within the range, copy
172
* p2 as a direct copy of p1
174
if (m0 == m1 && p1->m >= m0)
176
memcpy(p2, p1, sizeof(POINT4D));
178
#if DEBUG_INTERPOLATION
179
lwnotice("Projected p2 on range (as copy of p1)");
184
/* Otherwise interpolate coordinates */
190
#if DEBUG_INTERPOLATION
191
lwnotice("Projected p2 on range");
195
if ( swapped ) ret |= 0x0010;
199
/* Clipping occurred */
205
ptarray_locate_between_m(POINTARRAY *ipa, double m0, double m1)
208
DYNPTARRAY *dpa=NULL;
214
* We allocate space for as many pointarray as
215
* segments in the input POINTARRAY, as worst
216
* case is for each segment to cross the M range
218
* TODO: rework this to reduce used memory
220
ret.ptarrays=lwalloc(sizeof(POINTARRAY *)*ipa->npoints-1);
223
lwnotice("ptarray_locate...: called for pointarray %x, m0:%g, m1:%g",
228
for(i=1; i<ipa->npoints; i++)
233
getPoint4d_p(ipa, i-1, &p1);
234
getPoint4d_p(ipa, i, &p2);
237
lwnotice(" segment %d-%d [ %g %g %g %g - %g %g %g %g ]",
239
p1.x, p1.y, p1.z, p1.m,
240
p2.x, p2.y, p2.z, p2.m);
244
clipval = clip_seg_by_m_range(&p1, &p2, m0, m1);
246
/* segment completely outside, nothing to do */
247
if ( ! clipval ) continue;
250
lwnotice(" clipped to: [ %g %g %g %g - %g %g %g %g ] clipval: %x", p1.x, p1.y, p1.z, p1.m,
251
p2.x, p2.y, p2.z, p2.m, clipval);
254
/* If no points have been accumulated so far, then if clipval != 0 the first point must be the
255
start of the intersection */
259
lwnotice(" 1 creating new POINTARRAY with first point %g,%g,%g,%g", p1.x, p1.y, p1.z, p1.m);
261
dpa = dynptarray_create(ipa->npoints-i, ipa->dims);
262
dynptarray_addPoint4d(dpa, &p1, 1);
265
/* Otherwise always add the next point, avoiding duplicates */
267
dynptarray_addPoint4d(dpa, &p2, 0);
270
* second point has been clipped
272
if ( clipval & 0x0100 || i == ipa->npoints-1 )
275
lwnotice(" closing pointarray %x with %d points", dpa->pa, dpa->pa->npoints);
277
ret.ptarrays[ret.nptarrays++] = dpa->pa;
278
lwfree(dpa); dpa=NULL;
283
* if dpa!=NULL it means we didn't close it yet.
284
* this should never happen.
286
if ( dpa != NULL ) lwerror("Something wrong with algorithm");
292
* Point is assumed to have an M value.
293
* Return NULL if point is not in the given range (inclusive)
294
* Return an LWPOINT *copy* otherwise.
297
lwpoint_locate_between_m(LWPOINT *lwpoint, double m0, double m1)
302
lwnotice("lwpoint_locate_between called for lwpoint %x", lwpoint);
305
lwpoint_getPoint3dm_p(lwpoint, &p3dm);
306
if ( p3dm.m >= m0 && p3dm.m <= m1)
309
lwnotice(" lwpoint... returning a clone of input");
311
return (LWGEOM *)lwpoint_clone(lwpoint);
316
lwnotice(" lwpoint... returning a clone of input");
323
* Line is assumed to have an M value.
325
* Return NULL if no parts of the line are in the given range (inclusive)
327
* Return an LWCOLLECTION with LWLINES and LWPOINT being consecutive
328
* and isolated points on the line falling in the range.
330
* X,Y and Z (if present) ordinates are interpolated.
334
lwline_locate_between_m(LWLINE *lwline_in, double m0, double m1)
336
POINTARRAY *ipa=lwline_in->points;
341
int typeflag=0; /* see flags below */
342
const int pointflag=0x01;
343
const int lineflag=0x10;
346
lwnotice("lwline_locate_between called for lwline %x", lwline_in);
349
POINTARRAYSET paset=ptarray_locate_between_m(ipa, m0, m1);
352
lwnotice(" ptarray_locate... returned %d pointarrays",
356
if ( paset.nptarrays == 0 )
361
ngeoms=paset.nptarrays;
362
/* TODO: rework this to reduce used memory */
363
geoms=lwalloc(sizeof(LWGEOM *)*ngeoms);
364
for (i=0; i<ngeoms; i++)
369
POINTARRAY *pa=paset.ptarrays[i];
371
/* This is a point */
372
if ( pa->npoints == 1 )
375
lwpoint=lwalloc(sizeof(LWPOINT));
376
lwpoint->type=lwgeom_makeType_full(
382
lwpoint->SRID=lwline_in->SRID;
385
geoms[i]=(LWGEOM *)lwpoint;
390
else if ( pa->npoints > 1 )
392
lwline=lwalloc(sizeof(LWLINE));
393
lwline->type=lwgeom_makeType_full(
399
lwline->SRID=lwline_in->SRID;
402
geoms[i]=(LWGEOM *)lwline;
409
lwerror("ptarray_locate_between_m returned a POINARRAY set containing POINTARRAY with 0 points");
422
/* Choose best type */
423
if ( typeflag == 1 ) outtype=MULTIPOINTTYPE;
424
else if ( typeflag == 2 ) outtype=MULTILINETYPE;
425
else outtype = COLLECTIONTYPE;
427
return (LWGEOM *)lwcollection_construct(outtype,
428
lwline_in->SRID, NULL, ngeoms, geoms);
433
* Return a fully new allocated LWCOLLECTION
434
* always tagged as COLLECTIONTYPE.
437
lwcollection_locate_between_m(LWCOLLECTION *lwcoll, double m0, double m1)
444
lwnotice("lwcollection_locate_between_m called for lwcoll %x", lwcoll);
447
geoms=lwalloc(sizeof(LWGEOM *)*lwcoll->ngeoms);
448
for (i=0; i<lwcoll->ngeoms; i++)
450
LWGEOM *sub=lwgeom_locate_between_m(lwcoll->geoms[i],
453
geoms[ngeoms++] = sub;
456
if ( ngeoms == 0 ) return NULL;
458
return (LWGEOM *)lwcollection_construct(COLLECTIONTYPE,
459
lwcoll->SRID, NULL, ngeoms, geoms);
463
* Return a fully allocated LWGEOM containing elements
464
* intersected/interpolated with the given M range.
465
* Return NULL if none of the elements fall in the range.
467
* m0 is assumed to be less-or-equal to m1.
468
* input LWGEOM is assumed to contain an M value.
472
lwgeom_locate_between_m(LWGEOM *lwin, double m0, double m1)
475
lwnotice("lwgeom_locate_between called for lwgeom %x", lwin);
477
switch (TYPE_GETTYPE(lwin->type))
480
return lwpoint_locate_between_m(
481
(LWPOINT *)lwin, m0, m1);
483
return lwline_locate_between_m(
484
(LWLINE *)lwin, m0, m1);
489
return lwcollection_locate_between_m(
490
(LWCOLLECTION *)lwin, m0, m1);
492
/* Polygon types are not supported */
494
case MULTIPOLYGONTYPE:
495
lwerror("Areal geometries are not supported by locate_between_measures");
499
lwerror("Unkonwn geometry type (%s:%d)", __FILE__, __LINE__);
504
* Return a derived geometry collection value with elements that match
505
* the specified range of measures inclusively.
507
* Implements SQL/MM ST_LocateBetween(measure, measure) method.
508
* See ISO/IEC CD 13249-3:200x(E)
511
PG_FUNCTION_INFO_V1(LWGEOM_locate_between_m);
512
Datum LWGEOM_locate_between_m(PG_FUNCTION_ARGS)
514
PG_LWGEOM *gin = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
516
double start_measure = PG_GETARG_FLOAT8(1);
517
double end_measure = PG_GETARG_FLOAT8(2);
518
LWGEOM *lwin, *lwout;
521
if( end_measure < start_measure )
523
lwerror("locate_between_m: 2nd arg must be bigger then 1st arg");
528
* Return NULL if input doesn't have a measure
530
if ( ! lwgeom_hasM(gin->type) )
536
* Raise an error if input is a polygon, a multipolygon
539
type=lwgeom_getType(gin->type);
540
if ( type == POLYGONTYPE || type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE )
542
lwerror("Areal or Collection types are not supported");
546
lwin = pglwgeom_deserialize(gin);
548
lwout = lwgeom_locate_between_m(lwin,
549
start_measure, end_measure);
551
lwgeom_release(lwin);
555
lwout = (LWGEOM *)lwcollection_construct_empty(
556
pglwgeom_getSRID(gin), lwgeom_hasZ(gin->type),
557
lwgeom_hasM(gin->type));
560
gout = pglwgeom_serialize(lwout);
561
lwgeom_release(lwout);
563
PG_RETURN_POINTER(gout);
added
removed
lwgeom/lwgeom_functions_lrs.c
