From: Sandro Santilli Date: Wed, 25 Aug 2004 12:32:12 +0000 (+0000) Subject: Added perimeter,perimeter2d,perimeter3d. X-Git-Tag: pgis_0_9_1~48 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=5244c5bcd49fcf0c30a007da45a2bb9d26f146b5;p=postgis Added perimeter,perimeter2d,perimeter3d. Modified length,length2d,length3d semantic. Added force_2d. git-svn-id: http://svn.osgeo.org/postgis/trunk@746 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/lwgeom/MISSING_OBJECTS b/lwgeom/MISSING_OBJECTS index 4fdbfc657..8928427ee 100644 --- a/lwgeom/MISSING_OBJECTS +++ b/lwgeom/MISSING_OBJECTS @@ -55,9 +55,6 @@ FUNC: KEEPING FUNCTION: [geomcollfromtext(geometry)] FUNC: KEEPING FUNCTION: [isempty(geometry)] FUNC: KEEPING FUNCTION: [issimple(geometry)] FUNC: KEEPING FUNCTION: [equals(geometry, geometry)] -FUNC: KEEPING FUNCTION: [length_spheroid(geometry, spheroid)] -FUNC: KEEPING FUNCTION: [length3d_spheroid(geometry, spheroid)] -FUNC: KEEPING FUNCTION: [distance_spheroid(geometry, geometry, spheroid)] FUNC: KEEPING FUNCTION: [xmin(box3d)] FUNC: KEEPING FUNCTION: [ymin(box3d)] FUNC: KEEPING FUNCTION: [zmin(box3d)] @@ -67,10 +64,9 @@ FUNC: KEEPING FUNCTION: [zmax(box3d)] FUNC: KEEPING FUNCTION: [box3dtobox(box3d)] FUNC: KEEPING FUNCTION: [combine_bbox(box3d, geometry)] -FUNC: KEEPING FUNCTION: [length3d(geometry)] -FUNC: KEEPING FUNCTION: [length(geometry)] -FUNC: KEEPING FUNCTION: [perimeter3d(geometry)] -FUNC: KEEPING FUNCTION: [perimeter(geometry)] +FUNC: KEEPING FUNCTION: [length_spheroid(geometry, spheroid)] +FUNC: KEEPING FUNCTION: [length3d_spheroid(geometry, spheroid)] +FUNC: KEEPING FUNCTION: [distance_spheroid(geometry, geometry, spheroid)] FUNC: KEEPING FUNCTION: [point_inside_circle(geometry, double precision, double precision, double precision)] @@ -79,8 +75,9 @@ FUNC: KEEPING FUNCTION: [centroid(geometry)] FUNC: KEEPING FUNCTION: [isring(geometry)] FUNC: KEEPING FUNCTION: [pointonsurface(geometry)] FUNC: KEEPING FUNCTION: [collector(geometry, geometry)] -FUNC: KEEPING FUNCTION: [force_2d(geometry)] + FUNC: KEEPING FUNCTION: [force_3d(geometry)] + FUNC: KEEPING FUNCTION: [force_collection(geometry)] FUNC: KEEPING FUNCTION: [intersection(geometry, geometry)] FUNC: KEEPING FUNCTION: [buffer(geometry, double precision)] diff --git a/lwgeom/lwgeom_functions_basic.c b/lwgeom/lwgeom_functions_basic.c index 98ecf3bc4..ee15122aa 100644 --- a/lwgeom/lwgeom_functions_basic.c +++ b/lwgeom/lwgeom_functions_basic.c @@ -20,7 +20,7 @@ -#define DEBUG +//#define DEBUG #include "wktparse.h" @@ -28,16 +28,135 @@ Datum combine_box2d(PG_FUNCTION_ARGS); Datum LWGEOM_mem_size(PG_FUNCTION_ARGS); Datum LWGEOM_summary(PG_FUNCTION_ARGS); Datum LWGEOM_npoints(PG_FUNCTION_ARGS); -Datum LWGEOM_area2d(PG_FUNCTION_ARGS); +Datum LWGEOM_area_polygon(PG_FUNCTION_ARGS); Datum postgis_uses_stats(PG_FUNCTION_ARGS); Datum postgis_scripts_released(PG_FUNCTION_ARGS); Datum postgis_lib_version(PG_FUNCTION_ARGS); - +Datum LWGEOM_length2d_linestring(PG_FUNCTION_ARGS); +Datum LWGEOM_length3d_linestring(PG_FUNCTION_ARGS); +Datum LWGEOM_length_linestring(PG_FUNCTION_ARGS); +Datum LWGEOM_perimeter2d_poly(PG_FUNCTION_ARGS); +Datum LWGEOM_perimeter_poly(PG_FUNCTION_ARGS); +Datum LWGEOM_force_2d(PG_FUNCTION_ARGS); // internal char * lwgeom_summary_recursive(char *serialized, int offset); int32 lwgeom_npoints_recursive(char *serialized); -double lwgeom_polygon_area2d(LWPOLY *poly); +double lwgeom_polygon_area(LWPOLY *poly); +double lwgeom_polygon_perimeter(LWPOLY *poly); +double lwgeom_polygon_perimeter2d(LWPOLY *poly); +double lwgeom_line_length2d(LWLINE *line); +double lwgeom_pointarray_length2d(POINTARRAY *pts); +double lwgeom_pointarray_length(POINTARRAY *pts); +void lwgeom_force2d_recursive(char *serialized, char *optr, int *retsize); + +/*------------------------------------------------------------------*/ + +//find the 2d length of the given POINTARRAY (even if it's 3d) +double lwgeom_pointarray_length2d(POINTARRAY *pts) +{ + double dist = 0.0; + int i; + + if ( pts->npoints < 2 ) return 0.0; + for (i=0; inpoints-1;i++) + { + POINT2D *frm = (POINT2D *)getPoint(pts, i); + POINT2D *to = (POINT2D *)getPoint(pts, i+1); + dist += sqrt( ( (frm->x - to->x)*(frm->x - to->x) ) + + ((frm->y - to->y)*(frm->y - to->y) ) ); + } + return dist; +} + +//find the 3d/2d length of the given POINTARRAY (depending on its dimensions) +double lwgeom_pointarray_length(POINTARRAY *pts) +{ + double dist = 0.0; + int i; + + if ( pts->npoints < 2 ) return 0.0; + + // compute 2d length if 3d is not available + if ( pts->ndims < 3 ) return lwgeom_pointarray_length2d(pts); + + for (i=0; inpoints-1;i++) + { + POINT3D *frm = (POINT3D *)getPoint(pts, i); + POINT3D *to = (POINT3D *)getPoint(pts, i+1); + dist += sqrt( ( (frm->x - to->x)*(frm->x - to->x) ) + + ((frm->y - to->y)*(frm->y - to->y) ) + + ((frm->z - to->z)*(frm->z - to->z) ) ); + } + + return dist; +} + +//find the area of the outer ring - sum (area of inner rings) +// Could use a more numerically stable calculator... +double lwgeom_polygon_area(LWPOLY *poly) +{ + double poly_area=0.0; + int i; + +//elog(NOTICE,"in lwgeom_polygon_area (%d rings)", poly->nrings); + + for (i=0; inrings; i++) + { + int j; + POINTARRAY *ring = poly->rings[i]; + double ringarea = 0.0; + +//elog(NOTICE," rings %d has %d points", i, ring->npoints); + for (j=0; jnpoints-1; j++) + { + POINT2D *p1 = (POINT2D *)getPoint(ring, j); + POINT2D *p2 = (POINT2D *)getPoint(ring, j+1); + ringarea += ( p1->x * p2->y ) - ( p1->y * p2->x ); + } + + ringarea /= 2.0; +//elog(NOTICE," ring 1 has area %lf",ringarea); + ringarea = fabs(ringarea ); + if (i != 0) //outer + ringarea = -1.0*ringarea ; // its a hole + + poly_area += ringarea; + } + + return poly_area; +} + +// Compute the sum of polygon rings length. +// Could use a more numerically stable calculator... +double lwgeom_polygon_perimeter(LWPOLY *poly) +{ + double result=0.0; + int i; + +//elog(NOTICE,"in lwgeom_polygon_perimeter (%d rings)", poly->nrings); + + for (i=0; inrings; i++) + result += lwgeom_pointarray_length(poly->rings[i]); + + return result; +} + +// Compute the sum of polygon rings length (forcing 2d computation). +// Could use a more numerically stable calculator... +double lwgeom_polygon_perimeter2d(LWPOLY *poly) +{ + double result=0.0; + int i; + +//elog(NOTICE,"in lwgeom_polygon_perimeter (%d rings)", poly->nrings); + + for (i=0; inrings; i++) + result += lwgeom_pointarray_length2d(poly->rings[i]); + + return result; +} + /*------------------------------------------------------------------*/ @@ -314,65 +433,280 @@ Datum LWGEOM_npoints(PG_FUNCTION_ARGS) PG_RETURN_INT32(npoints); } -//find the 2d area of the outer ring - sum (area 2d of inner rings) -// Could use a more numerically stable calculator... -double lwgeom_polygon_area2d(LWPOLY *poly) +// Calculate the area of all the subobj in a polygon +// area(point) = 0 +// area (line) = 0 +// area(polygon) = find its 2d area +PG_FUNCTION_INFO_V1(LWGEOM_area_polygon); +Datum LWGEOM_area_polygon(PG_FUNCTION_ARGS) { - double poly_area=0.0; + LWGEOM *geom = (LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + LWGEOM_INSPECTED *inspected = lwgeom_inspect(SERIALIZED_FORM(geom)); + LWPOLY *poly; + double area = 0.0; int i; -//elog(NOTICE,"in lwgeom_polygon_area2d (%d rings)", poly->nrings); +//elog(NOTICE, "in LWGEOM_area_polygon"); - for (i=0; inrings; i++) + for (i=0; ingeometries; i++) { - int j; - POINTARRAY *ring = poly->rings[i]; - double ringarea = 0.0; + poly = lwgeom_getpoly_inspected(inspected, i); + if ( poly == NULL ) continue; + area += lwgeom_polygon_area(poly); +//elog(NOTICE, " LWGEOM_area_polygon found a poly (%f)", area); + } + + pfree_inspected(inspected); -//elog(NOTICE," rings %d has %d points", i, ring->npoints); - for (j=0; jnpoints-1; j++) - { - POINT2D *p1 = (POINT2D *)getPoint(ring, j); - POINT2D *p2 = (POINT2D *)getPoint(ring, j+1); - ringarea += ( p1->x * p2->y ) - ( p1->y * p2->x ); - } + PG_RETURN_FLOAT8(area); +} - ringarea /= 2.0; -//elog(NOTICE," ring 1 has area %lf",ringarea); - ringarea = fabs(ringarea ); - if (i != 0) //outer - ringarea = -1.0*ringarea ; // its a hole +//find the "length of a geometry" +// length2d(point) = 0 +// length2d(line) = length of line +// length2d(polygon) = 0 -- could make sense to return sum(ring perimeter) +// uses euclidian 2d length (even if input is 3d) +PG_FUNCTION_INFO_V1(LWGEOM_length2d_linestring); +Datum LWGEOM_length2d_linestring(PG_FUNCTION_ARGS) +{ + LWGEOM *geom = (LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + LWGEOM_INSPECTED *inspected = lwgeom_inspect(SERIALIZED_FORM(geom)); + LWLINE *line; + double dist = 0.0; + int i; - poly_area += ringarea; +//elog(NOTICE, "in LWGEOM_length2d"); + + for (i=0; ingeometries; i++) + { + line = lwgeom_getline_inspected(inspected, i); + if ( line == NULL ) continue; + dist += lwgeom_pointarray_length2d(line->points); +//elog(NOTICE, " LWGEOM_length2d found a line (%f)", dist); } - return poly_area; + pfree_inspected(inspected); + + PG_RETURN_FLOAT8(dist); } -// Calculate the area of all the subobj in a polygon -// area(point) = 0 -// area (line) = 0 -// area(polygon) = find its 2d area -PG_FUNCTION_INFO_V1(LWGEOM_area2d); -Datum LWGEOM_area2d(PG_FUNCTION_ARGS) +//find the "length of a geometry" +// length(point) = 0 +// length(line) = length of line +// length(polygon) = 0 -- could make sense to return sum(ring perimeter) +// uses euclidian 3d/2d length depending on input dimensions. +PG_FUNCTION_INFO_V1(LWGEOM_length_linestring); +Datum LWGEOM_length_linestring(PG_FUNCTION_ARGS) { LWGEOM *geom = (LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); LWGEOM_INSPECTED *inspected = lwgeom_inspect(SERIALIZED_FORM(geom)); - LWPOLY *poly; - double area = 0.0; + LWLINE *line; + double dist = 0.0; int i; -//elog(NOTICE, "in LWGEOM_area2d"); + for (i=0; ingeometries; i++) + { + line = lwgeom_getline_inspected(inspected, i); + if ( line == NULL ) continue; + dist += lwgeom_pointarray_length(line->points); + } + + pfree_inspected(inspected); + + PG_RETURN_FLOAT8(dist); +} + +// find the "perimeter of a geometry" +// perimeter(point) = 0 +// perimeter(line) = 0 +// perimeter(polygon) = sum of ring perimeters +// uses euclidian 3d/2d computation depending on input dimension. +PG_FUNCTION_INFO_V1(LWGEOM_perimeter_poly); +Datum LWGEOM_perimeter_poly(PG_FUNCTION_ARGS) +{ + LWGEOM *geom = (LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + LWGEOM_INSPECTED *inspected = lwgeom_inspect(SERIALIZED_FORM(geom)); + double ret = 0.0; + int i; for (i=0; ingeometries; i++) { + LWPOLY *poly; poly = lwgeom_getpoly_inspected(inspected, i); if ( poly == NULL ) continue; - area += lwgeom_polygon_area2d(poly); -//elog(NOTICE, " LWGEOM_area2d found a poly (%f)", area); + ret += lwgeom_polygon_perimeter(poly); } - + pfree_inspected(inspected); - PG_RETURN_FLOAT8(area); + PG_RETURN_FLOAT8(ret); +} + +// find the "perimeter of a geometry" +// perimeter(point) = 0 +// perimeter(line) = 0 +// perimeter(polygon) = sum of ring perimeters +// uses euclidian 2d computation even if input is 3d +PG_FUNCTION_INFO_V1(LWGEOM_perimeter2d_poly); +Datum LWGEOM_perimeter2d_poly(PG_FUNCTION_ARGS) +{ + LWGEOM *geom = (LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + LWGEOM_INSPECTED *inspected = lwgeom_inspect(SERIALIZED_FORM(geom)); + double ret = 0.0; + int i; + + for (i=0; ingeometries; i++) + { + LWPOLY *poly; + poly = lwgeom_getpoly_inspected(inspected, i); + if ( poly == NULL ) continue; + ret += lwgeom_polygon_perimeter2d(poly); + } + + pfree_inspected(inspected); + + PG_RETURN_FLOAT8(ret); +} + + +/* + * Write to already allocated memory 'optr' a 2d version of + * the given serialized form. + * Return number bytes written in given int pointer. + */ +void +lwgeom_force2d_recursive(char *serialized, char *optr, int *retsize) +{ + LWGEOM_INSPECTED *inspected; + int i; + int totsize=0; + int size=0; + int type; + LWPOINT *point = NULL; + LWLINE *line = NULL; + LWPOLY *poly = NULL; + char *loc; + + +#ifdef DEBUG + elog(NOTICE, "lwgeom_force2d_recursive: call"); +#endif + + type = lwgeom_getType(serialized[0]); + + if ( type == POINTTYPE ) + { + point = lwpoint_deserialize(serialized); + point->ndims = 2; + lwpoint_serialize_buf(point, optr, retsize); +#ifdef DEBUG +elog(NOTICE, "lwgeom_force2d_recursive: it's a point, size:%d", *retsize); +#endif + return; + } + + if ( type == LINETYPE ) + { + line = lwline_deserialize(serialized); + line->ndims = 2; + lwline_serialize_buf(line, optr, retsize); +#ifdef DEBUG +elog(NOTICE, "lwgeom_force2d_recursive: it's a line, size:%d", *retsize); +#endif + return; + } + + if ( type == POLYGONTYPE ) + { + poly = lwpoly_deserialize(serialized); + poly->ndims = 2; + lwpoly_serialize_buf(poly, optr, retsize); +#ifdef DEBUG +elog(NOTICE, "lwgeom_force2d_recursive: it's a poly, size:%d", *retsize); +#endif + return; + } + + // OK, this is a collection, so we write down its metadata + // first and then call us again + +#ifdef DEBUG +elog(NOTICE, "lwgeom_force2d_recursive: it's a collection (type:%d)", type); +#endif + + // Add type + *optr = lwgeom_makeType_full(2, lwgeom_hasSRID(serialized[0]), + type, lwgeom_hasBBOX(serialized[0])); + optr++; + totsize++; + loc=serialized+1; + + // Add BBOX if any + if (lwgeom_hasBBOX(serialized[0])) + { + memcpy(optr, loc, sizeof(BOX2DFLOAT4)); + optr += sizeof(BOX2DFLOAT4); + totsize += sizeof(BOX2DFLOAT4); + loc += sizeof(BOX2DFLOAT4); + } + + // Add SRID if any + if (lwgeom_hasSRID(serialized[0])) + { + memcpy(optr, loc, 4); + optr += 4; + totsize += 4; + loc += 4; + } + + // Add numsubobjects + memcpy(optr, loc, 4); + optr += 4; + totsize += 4; + +#ifdef DEBUG +elog(NOTICE, " collection header size:%d", totsize); +#endif + + // Now recurse for each suboject + inspected = lwgeom_inspect(serialized); + for (i=0; ingeometries; i++) + { + char *subgeom = lwgeom_getsubgeometry_inspected(inspected, i); + lwgeom_force2d_recursive(subgeom, optr, &size); + totsize += size; + optr += size; +#ifdef DEBUG +elog(NOTICE, " elem %d size: %d (tot: %d)", i, size, totsize); +#endif + } + pfree_inspected(inspected); + + *retsize = totsize; +} + +// transform input geometry to 2d if not 2d already +PG_FUNCTION_INFO_V1(LWGEOM_force_2d); +Datum LWGEOM_force_2d(PG_FUNCTION_ARGS) +{ + LWGEOM *geom = (LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + LWGEOM *result; + int32 size = 0; + + // already 2d + if ( lwgeom_ndims(geom->type) == 2 ) PG_RETURN_POINTER(geom); + + // allocate a larger for safety and simplicity + result = (LWGEOM *) palloc(geom->size); + + lwgeom_force2d_recursive(SERIALIZED_FORM(geom), + SERIALIZED_FORM(result), &size); + + // we can safely avoid this... memory will be freed at + // end of query processing anyway. + //result = repalloc(result, size+4); + + result->size = size+4; + + PG_RETURN_POINTER(result); } diff --git a/lwgeom/lwpostgis.sql.in b/lwgeom/lwpostgis.sql.in index b963ae52c..bfda37af2 100644 --- a/lwgeom/lwpostgis.sql.in +++ b/lwgeom/lwpostgis.sql.in @@ -10,60 +10,6 @@ -- the terms of the GNU General Public Licence. See the COPYING file. -- -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - --- $Log$ --- Revision 1.16 2004/08/25 07:29:32 strk --- Moved some OGC functions from lwgeom_inout.c to lwgeom_ogc.c. --- Added area2d (area) to lwgeom_functions_basic.c --- --- Revision 1.15 2004/08/24 16:33:41 strk --- Added StartPoint(), EndPoint() --- --- Revision 1.14 2004/08/24 16:20:10 strk --- Added X(), Y() and Z() funx --- --- Revision 1.13 2004/08/24 15:50:04 strk --- PointN() ported. --- --- Revision 1.12 2004/08/24 15:05:34 strk --- Added NumInteriorRings() and InteriorRingN() --- --- Revision 1.11 2004/08/24 14:48:58 strk --- Added dimension() and exteriorring() --- --- Revision 1.10 2004/08/24 10:01:16 strk --- OGC functions (not implemented by GEOS) moved to lwgeom_ogc.c. --- Renamed PG-exposed functions to start with LWGEOM --- --- Revision 1.9 2004/08/24 09:34:33 strk --- Added npoints,numpoints,numgeometries,geometryn --- --- Revision 1.8 2004/08/23 15:57:56 strk --- versioning functions completed --- --- Revision 1.7 2004/08/23 15:37:16 strk --- Changed SCRIPTS_VERSION to 0.0.1 --- --- Revision 1.6 2004/08/23 08:32:14 strk --- Removed custom allocator from lwgeom_api. --- Added summary(geometry). --- Some indentation. --- --- Revision 1.5 2004/08/20 16:36:22 strk --- transform() support --- --- Revision 1.4 2004/08/20 14:54:35 strk --- gist operators support functions renamed to allow for finer control by postgis_restore.pl --- --- Revision 1.3 2004/08/20 14:08:41 strk --- Added Geom{etry,}FromWkb(,[]) funx. --- Added LWGEOM typedef and SERIALIZED_FORM(LWGEOM) macro. --- Made lwgeom_setSRID an API function. --- Added LWGEOM_setAllocator(). --- --- Revision 1.2 2004/08/19 13:54:15 strk --- cpp checks updated to use 80 instead of 75 for USE_VERSION --- --- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #if USE_VERSION > 71 #define CREATEFUNCTION CREATE OR REPLACE FUNCTION @@ -940,17 +886,62 @@ CREATEFUNCTION GeometryFromText(geometry, int4) -- Misures ------------------------------------------------------------------------ +-- this is a fake (for back-compatibility) +-- uses 3d if 3d is available, 2d otherwise +CREATEFUNCTION length3d(geometry) + RETURNS FLOAT8 + AS '@MODULE_FILENAME@', 'LWGEOM_length_linestring' + LANGUAGE 'C' WITH (isstrict); + +CREATEFUNCTION length2d(geometry) + RETURNS FLOAT8 + AS '@MODULE_FILENAME@', 'LWGEOM_length2d_linestring' + LANGUAGE 'C' WITH (isstrict); + +CREATEFUNCTION length(geometry) + RETURNS FLOAT8 + AS '@MODULE_FILENAME@', 'LWGEOM_length_linestring' + LANGUAGE 'C' WITH (isstrict); + +-- this is a fake (for back-compatibility) +-- uses 3d if 3d is available, 2d otherwise +CREATEFUNCTION perimeter3d(geometry) + RETURNS FLOAT8 + AS '@MODULE_FILENAME@', 'LWGEOM_perimeter_poly' + LANGUAGE 'C' WITH (isstrict); + +CREATEFUNCTION perimeter2d(geometry) + RETURNS FLOAT8 + AS '@MODULE_FILENAME@', 'LWGEOM_perimeter2d_poly' + LANGUAGE 'C' WITH (isstrict); + +CREATEFUNCTION perimeter(geometry) + RETURNS FLOAT8 + AS '@MODULE_FILENAME@', 'LWGEOM_perimeter_poly' + LANGUAGE 'C' WITH (isstrict); + +-- this is an alias for 'area(geometry)' +-- there is nothing such an 'area3d'... CREATEFUNCTION area2d(geometry) RETURNS FLOAT8 - AS '@MODULE_FILENAME@', 'LWGEOM_area2d' + AS '@MODULE_FILENAME@', 'LWGEOM_area_polygon' LANGUAGE 'C' WITH (isstrict); CREATEFUNCTION area(geometry) RETURNS FLOAT8 - AS '@MODULE_FILENAME@','LWGEOM_area2d' + AS '@MODULE_FILENAME@','LWGEOM_area_polygon' LANGUAGE 'C' WITH (isstrict); +------------------------------------------------------------------------ +-- MISC +------------------------------------------------------------------------ + +CREATEFUNCTION force_2d(geometry) + RETURNS geometry + AS '@MODULE_FILENAME@', 'LWGEOM_force_2d' + LANGUAGE 'C' WITH (isstrict); + ------------------------------------------------------------------------ --