From: Paul Ramsey Date: Mon, 9 Jan 2012 18:27:30 +0000 (+0000) Subject: Interim progress on LRS work. X-Git-Tag: 2.0.0alpha1~150 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=acc2cfbad56c6a33b604a4c1b6d293583ea294a6;p=postgis Interim progress on LRS work. git-svn-id: http://svn.osgeo.org/postgis/trunk@8726 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in index a084bc417..530cd44a1 100644 --- a/liblwgeom/liblwgeom.h.in +++ b/liblwgeom/liblwgeom.h.in @@ -750,7 +750,7 @@ extern int getPoint2d_p(const POINTARRAY *pa, int n, POINT2D *point); * will be extracted from it * */ -extern void ptarray_set_point4d(POINTARRAY *pa, int n, POINT4D *p4d); +extern void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d); /* * get a pointer to nth point of a POINTARRAY @@ -799,7 +799,7 @@ extern POINTARRAY* ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoi * If allow_duplicate is LW_TRUE, then a duplicate point will * not be added. */ -extern int ptarray_append_point(POINTARRAY *pa, POINT4D *pt, int allow_duplicates); +extern int ptarray_append_point(POINTARRAY *pa, const POINT4D *pt, int allow_duplicates); /** * Append a #POINTARRAY, pa2 to the end of an existing #POINTARRAY, pa1. @@ -812,7 +812,7 @@ extern int ptarray_append_ptarray(POINTARRAY *pa1, POINTARRAY *pa2, int splice_e * Insert a point into an existing #POINTARRAY. Zero * is the index of the start of the array. */ -extern int ptarray_insert_point(POINTARRAY *pa, POINT4D *p, int where); +extern int ptarray_insert_point(POINTARRAY *pa, const POINT4D *p, int where); /** * Remove a point from an existing #POINTARRAY. Zero @@ -1157,6 +1157,7 @@ POINTARRAY *ptarray_clone_deep(const POINTARRAY *ptarray); * from underneath the geometries. */ extern LWPOINT* lwpoint_construct(int srid, GBOX *bbox, POINTARRAY *point); +extern LWMPOINT *lwmpoint_construct(int srid, const POINTARRAY *pa); extern LWLINE* lwline_construct(int srid, GBOX *bbox, POINTARRAY *points); extern LWCIRCSTRING* lwcircstring_construct(int srid, GBOX *bbox, POINTARRAY *points); extern LWPOLY* lwpoly_construct(int srid, GBOX *bbox, uint32_t nrings, POINTARRAY **points); @@ -1185,7 +1186,8 @@ extern LWPOINT *lwpoint_make2d(int srid, double x, double y); extern LWPOINT *lwpoint_make3dz(int srid, double x, double y, double z); extern LWPOINT *lwpoint_make3dm(int srid, double x, double y, double m); extern LWPOINT *lwpoint_make4d(int srid, double x, double y, double z, double m); -extern LWLINE *lwline_from_lwpointarray(int srid, uint32_t npoints, LWPOINT **points); +extern LWPOINT *lwpoint_make(int srid, int hasz, int hasm, const POINT4D *p); +extern LWLINE *lwline_from_ptarray(int srid, uint32_t npoints, LWPOINT **points); extern LWLINE *lwline_from_lwmpoint(int srid, LWMPOINT *mpoint); extern LWLINE *lwline_addpoint(LWLINE *line, LWPOINT *point, uint32_t where); extern LWLINE *lwline_removepoint(LWLINE *line, uint32_t which); @@ -1229,10 +1231,19 @@ extern int lwgeom_ndims(const LWGEOM *geom); */ extern double ptarray_locate_point(POINTARRAY *, POINT2D *, double *); - +/** +* Add a measure dimension to a line, interpolating linearly from the start +* to the end value. +*/ extern LWLINE *lwline_measured_from_lwline(const LWLINE *lwline, double m_start, double m_end); extern LWMLINE* lwmline_measured_from_lwmline(const LWMLINE *lwmline, double m_start, double m_end); +/** +* Determine the location(s) along a measured line where m occurs and +* return as a multipoint. Offset to left (negative) or right (positive). +*/ +extern LWGEOM* lwgeom_locate_along(const LWGEOM *lwin, double m, double offset); + /* * Ensure every segment is at most 'dist' long. * Returned LWGEOM might is unchanged if a POINT. diff --git a/liblwgeom/liblwgeom_internal.h b/liblwgeom/liblwgeom_internal.h index 9cc23d435..20a6a5f3e 100644 --- a/liblwgeom/liblwgeom_internal.h +++ b/liblwgeom/liblwgeom_internal.h @@ -274,6 +274,8 @@ void ptarray_affine(POINTARRAY *pa, const AFFINE *affine); * PointArray */ char ptarray_isccw(const POINTARRAY *pa); +int ptarray_has_z(const POINTARRAY *pa); +int ptarray_has_m(const POINTARRAY *pa); /* * Clone support diff --git a/liblwgeom/lwgeom_api.c b/liblwgeom/lwgeom_api.c index c53899ab4..78edb496b 100644 --- a/liblwgeom/lwgeom_api.c +++ b/liblwgeom/lwgeom_api.c @@ -455,7 +455,7 @@ getPoint2d_p(const POINTARRAY *pa, int n, POINT2D *point) * */ void -ptarray_set_point4d(POINTARRAY *pa, int n, POINT4D *p4d) +ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d) { uint8_t *ptr=getPoint_internal(pa, n); switch ( FLAGS_GET_ZM(pa->flags) ) diff --git a/liblwgeom/lwgeom_geos.c b/liblwgeom/lwgeom_geos.c index ccb024918..9371ae62f 100644 --- a/liblwgeom/lwgeom_geos.c +++ b/liblwgeom/lwgeom_geos.c @@ -1052,10 +1052,11 @@ lwgeom_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyl #else GEOSGeometry *g1, *g3; LWGEOM *lwgeom_result; + LWGEOM *lwgeom_in = lwline_as_lwgeom(lwline); initGEOS(lwnotice, lwgeom_geos_error); - g1 = (GEOSGeometry *)LWGEOM2GEOS(lwline_as_lwgeom(lwline)); + g1 = (GEOSGeometry *)LWGEOM2GEOS(lwgeom_in); if ( ! g1 ) { lwerror("lwgeom_offsetcurve: Geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg); @@ -1082,9 +1083,8 @@ lwgeom_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyl LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3)); - GEOSSetSRID(g3, lwgeom_get_srid(lwline_as_lwgeom(lwline))); - - lwgeom_result = GEOS2LWGEOM(g3, lwgeom_has_z(lwline_as_lwgeom(lwline))); + GEOSSetSRID(g3, lwgeom_get_srid(lwgeom_in)); + lwgeom_result = GEOS2LWGEOM(g3, lwgeom_has_z(lwgeom_in)); GEOSGeom_destroy(g3); if (lwgeom_result == NULL) diff --git a/liblwgeom/lwline.c b/liblwgeom/lwline.c index bfaa45f64..24057e445 100644 --- a/liblwgeom/lwline.c +++ b/liblwgeom/lwline.c @@ -146,7 +146,7 @@ lwline_same(const LWLINE *l1, const LWLINE *l2) * LWLINE dimensions are large enough to host all input dimensions. */ LWLINE * -lwline_from_lwpointarray(int srid, uint32_t npoints, LWPOINT **points) +lwline_from_ptarray(int srid, uint32_t npoints, LWPOINT **points) { int i; int hasz = LW_FALSE; @@ -162,7 +162,7 @@ lwline_from_lwpointarray(int srid, uint32_t npoints, LWPOINT **points) { if ( points[i]->type != POINTTYPE ) { - lwerror("lwline_from_lwpointarray: invalid input type: %s", lwtype_name(points[i]->type)); + lwerror("lwline_from_ptarray: invalid input type: %s", lwtype_name(points[i]->type)); return NULL; } if ( FLAGS_GET_Z(points[i]->flags) ) hasz = LW_TRUE; diff --git a/liblwgeom/lwlinearreferencing.c b/liblwgeom/lwlinearreferencing.c index d7c9762c8..bbe4096ca 100644 --- a/liblwgeom/lwlinearreferencing.c +++ b/liblwgeom/lwlinearreferencing.c @@ -14,6 +14,167 @@ #include "lwgeom_log.h" +static int +segment_locate_along(const POINT4D *p1, const POINT4D *p2, double m, double offset, POINT4D *pn) +{ + double m1 = p1->m; + double m2 = p2->m; + double mprop; + + /* M is out of range, no new point generated. */ + if ( (m < FP_MIN(m1,m2)) || (m > FP_MAX(m1,m2)) ) + { + return LW_FALSE; + } + + /* We'll just can out on this degenerate case for now. + Correct behavior is probably an mprop of 0.5? + Still would have to deal with case of true p1==p2. */ + if( m1 == m2 ) + { + lwerror("Zero measure-length line encountered!"); + } + + /* M is in range, new point to be generated. */ + mprop = (m - m1) / (m2 - m1); + pn->x = p1->x + (p2->x - p1->x) * mprop; + pn->y = p1->y + (p2->y - p1->y) * mprop; + pn->z = p1->z + (p2->z - p1->z) * mprop; + pn->m = m; + + /* Offset to the left or right, if necessary. */ + if ( offset != 0.0 ) + { + double theta = atan2(p2->y - p1->y, p2->x - p1->x); + pn->x -= sin(theta) * offset; + pn->y += cos(theta) * offset; + } + + return LW_TRUE; +} + + +static POINTARRAY* +ptarray_locate_along(const POINTARRAY *pa, double m, double offset) +{ + int i; + POINT4D p1, p2, pn; + POINTARRAY *dpa = NULL; + + /* Can't do anything with degenerate point arrays */ + if ( ! pa || pa->npoints < 2 ) return NULL; + + /* Walk through each segment in the point array */ + for ( i = 1; i < pa->npoints; i++ ) + { + getPoint4d_p(pa, i-1, &p1); + getPoint4d_p(pa, i, &p2); + + /* No derived point? Move to next segment. */ + if ( segment_locate_along(&p1, &p2, m, offset, &pn) == LW_FALSE ) + continue; + + /* No pointarray, make a fresh one */ + if ( dpa == NULL ) + dpa = ptarray_construct_empty(ptarray_has_z(pa), ptarray_has_m(pa), 8); + + /* Add our new point to the array */ + ptarray_append_point(dpa, &pn, 0); + } + + return dpa; +} + +static LWMPOINT* +lwline_locate_along(const LWLINE *lwline, double m, double offset) +{ + POINTARRAY *opa = NULL; + LWMPOINT *mp = NULL; + LWGEOM *lwg = lwline_as_lwgeom(lwline); + int hasz, hasm, srid; + + /* Return degenerates upwards */ + if ( ! lwline ) return NULL; + + /* Create empty return shell */ + srid = lwgeom_get_srid(lwg); + hasz = lwgeom_has_z(lwg); + hasm = lwgeom_has_m(lwg); + + if ( hasm ) + { + /* Find points along */ + opa = ptarray_locate_along(lwline->points, m, offset); + } + else + { + LWLINE *lwline_measured = lwline_measured_from_lwline(lwline, 0.0, 1.0); + opa = ptarray_locate_along(lwline_measured->points, m, offset); + lwline_free(lwline_measured); + } + + /* Return NULL as EMPTY */ + if ( ! opa ) + return lwmpoint_construct_empty(srid, hasz, hasm); + + /* Convert pointarray into a multipoint */ + mp = lwmpoint_construct(srid, opa); + ptarray_free(opa); + return mp; +} + +static LWMPOINT* +lwmline_locate_along(const LWMLINE *lwmline, double m, double offset) +{ + LWMPOINT *lwmpoint = NULL; + LWGEOM *lwg = lwmline_as_lwgeom(lwmline); + int i, j; + + /* Return degenerates upwards */ + if ( (!lwmline) || (lwmline->ngeoms < 1) ) return NULL; + + /* Construct return */ + lwmpoint = lwmpoint_construct_empty(lwgeom_get_srid(lwg), lwgeom_has_z(lwg), lwgeom_has_m(lwg)); + + /* Locate along each sub-line */ + for ( i = 0; i < lwmline->ngeoms; i++ ) + { + LWMPOINT *along = lwline_locate_along(lwmline->geoms[i], m, offset); + if ( along != NULL ) + { + for ( j = 0; j < along->ngeoms; j++ ) + { + lwmpoint_add_lwpoint(lwmpoint, along->geoms[i]); + } + /* Free the containing geometry, but leave the sub-geometries around */ + if ( along->bbox ) lwfree(along->bbox); + lwfree(along); + } + } + return lwmpoint; +} + +LWGEOM* +lwgeom_locate_along(const LWGEOM *lwin, double m, double offset) +{ + if ( ! lwin ) return NULL; + + switch (lwin->type) + { + case LINETYPE: + return (LWGEOM*)lwline_locate_along((LWLINE*)lwin, m, offset); + case MULTILINETYPE: + return (LWGEOM*)lwmline_locate_along((LWMLINE*)lwin, m, offset); + /* Only line types supported right now */ + /* TO DO: CurveString, CompoundCurve, MultiCurve */ + /* TO DO: Point, MultiPoint */ + default: + lwerror("Only linear geometries are supported, %s provided.",lwtype_name(lwin->type)); + return NULL; + } + return NULL; +} + /** * Given a POINT4D and an ordinate number, return * the value of the ordinate. diff --git a/liblwgeom/lwmpoint.c b/liblwgeom/lwmpoint.c index ff4e54bab..72c44284c 100644 --- a/liblwgeom/lwmpoint.c +++ b/liblwgeom/lwmpoint.c @@ -29,15 +29,33 @@ lwmpoint_construct_empty(int srid, char hasz, char hasm) return ret; } - - - LWMPOINT* lwmpoint_add_lwpoint(LWMPOINT *mobj, const LWPOINT *obj) { LWDEBUG(4, "Called"); return (LWMPOINT*)lwcollection_add_lwgeom((LWCOLLECTION*)mobj, (LWGEOM*)obj); } +LWMPOINT * +lwmpoint_construct(int srid, const POINTARRAY *pa) +{ + int i; + int hasz = ptarray_has_z(pa); + int hasm = ptarray_has_m(pa); + LWMPOINT *ret = (LWMPOINT*)lwcollection_construct_empty(MULTIPOINTTYPE, srid, hasz, hasm); + + for ( i = 0; i < pa->npoints; i++ ) + { + LWPOINT *lwp; + POINT4D p; + getPoint4d_p(pa, i, &p); + lwp = lwpoint_make(srid, hasz, hasm, &p); + lwmpoint_add_lwpoint(ret, lwp); + } + + return ret; +} + + void lwmpoint_free(LWMPOINT *mpt) { int i; diff --git a/liblwgeom/lwpoint.c b/liblwgeom/lwpoint.c index 8aaf4a259..c56a4cf6d 100644 --- a/liblwgeom/lwpoint.c +++ b/liblwgeom/lwpoint.c @@ -169,6 +169,14 @@ lwpoint_make4d(int srid, double x, double y, double z, double m) return lwpoint_construct(srid, NULL, pa); } +LWPOINT * +lwpoint_make(int srid, int hasz, int hasm, const POINT4D *p) +{ + POINTARRAY *pa = ptarray_construct_empty(hasz, hasm, 1); + ptarray_append_point(pa, p, LW_TRUE); + return lwpoint_construct(srid, NULL, pa); +} + void lwpoint_free(LWPOINT *pt) { if ( pt->bbox ) diff --git a/liblwgeom/measures.c b/liblwgeom/measures.c index 439020706..190628c77 100644 --- a/liblwgeom/measures.c +++ b/liblwgeom/measures.c @@ -64,7 +64,7 @@ lw_dist2d_distanceline(LWGEOM *lw1, LWGEOM *lw2,int srid,int mode) lwpoints[0] = lwpoint_make2d(srid, x1, y1); lwpoints[1] = lwpoint_make2d(srid, x2, y2); - result = (LWGEOM *)lwline_from_lwpointarray(srid, 2, lwpoints); + result = (LWGEOM *)lwline_from_ptarray(srid, 2, lwpoints); } return result; } diff --git a/liblwgeom/measures3d.c b/liblwgeom/measures3d.c index ef67bcbeb..d1de278bf 100644 --- a/liblwgeom/measures3d.c +++ b/liblwgeom/measures3d.c @@ -66,7 +66,7 @@ lw_dist3d_distanceline(LWGEOM *lw1, LWGEOM *lw2,int srid,int mode) lwpoints[0] = lwpoint_make3dz(srid, x1, y1, z1); lwpoints[1] = lwpoint_make3dz(srid, x2, y2, z2); - result = (LWGEOM *)lwline_from_lwpointarray(srid, 2, lwpoints); + result = (LWGEOM *)lwline_from_ptarray(srid, 2, lwpoints); } return result; diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c index c5f8834f1..51da16992 100644 --- a/liblwgeom/ptarray.c +++ b/liblwgeom/ptarray.c @@ -17,6 +17,20 @@ /*#define POSTGIS_DEBUG_LEVEL 4*/ #include "lwgeom_log.h" +int +ptarray_has_z(const POINTARRAY *pa) +{ + if ( ! pa ) return LW_FALSE; + return FLAGS_GET_Z(pa->flags); +} + +int +ptarray_has_m(const POINTARRAY *pa) +{ + if ( ! pa ) return LW_FALSE; + return FLAGS_GET_M(pa->flags); +} + /* * Size of point represeneted in the POINTARRAY * 16 for 2d, 24 for 3d, 32 for 4d @@ -65,7 +79,7 @@ ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints) * pointarray supports. */ int -ptarray_insert_point(POINTARRAY *pa, POINT4D *p, int where) +ptarray_insert_point(POINTARRAY *pa, const POINT4D *p, int where) { size_t point_size = ptarray_point_size(pa); LWDEBUGF(5,"pa = %p; p = %p; where = %d", pa, p, where); @@ -122,7 +136,7 @@ ptarray_insert_point(POINTARRAY *pa, POINT4D *p, int where) } int -ptarray_append_point(POINTARRAY *pa, POINT4D *pt, int repeated_points) +ptarray_append_point(POINTARRAY *pa, const POINT4D *pt, int repeated_points) { /* Check for pathology */ diff --git a/postgis/lwgeom_functions_basic.c b/postgis/lwgeom_functions_basic.c index dea38ea3c..f526391b1 100644 --- a/postgis/lwgeom_functions_basic.c +++ b/postgis/lwgeom_functions_basic.c @@ -1430,7 +1430,7 @@ Datum LWGEOM_makeline_garray(PG_FUNCTION_ARGS) POSTGIS_DEBUGF(3, "LWGEOM_makeline_garray: point elements: %d", npoints); - outlwg = (LWGEOM *)lwline_from_lwpointarray(srid, npoints, lwpoints); + outlwg = (LWGEOM *)lwline_from_ptarray(srid, npoints, lwpoints); result = geometry_serialize(outlwg); @@ -1467,7 +1467,7 @@ Datum LWGEOM_makeline(PG_FUNCTION_ARGS) lwpoints[0] = lwgeom_as_lwpoint(lwgeom_from_gserialized(pglwg1)); lwpoints[1] = lwgeom_as_lwpoint(lwgeom_from_gserialized(pglwg2)); - outline = lwline_from_lwpointarray(lwpoints[0]->srid, 2, lwpoints); + outline = lwline_from_ptarray(lwpoints[0]->srid, 2, lwpoints); result = geometry_serialize((LWGEOM *)outline); diff --git a/postgis/lwgeom_geos.c b/postgis/lwgeom_geos.c index 621cd3e98..f81e8edc3 100644 --- a/postgis/lwgeom_geos.c +++ b/postgis/lwgeom_geos.c @@ -1246,7 +1246,7 @@ Datum ST_OffsetCurve(PG_FUNCTION_ARGS) pfree(paramstr); /* alloc'ed in text2cstring */ } - lwgeom_result = lwline_as_lwgeom(lwgeom_offsetcurve(lwgeom_as_lwline(lwgeom_input), size, quadsegs, joinStyle, mitreLimit)); + lwgeom_result = lwgeom_offsetcurve(lwgeom_as_lwline(lwgeom_input), size, quadsegs, joinStyle, mitreLimit); if (lwgeom_result == NULL) lwerror("ST_OffsetCurve: lwgeom_offsetcurve returned NULL");