From: Nicklas Avén Date: Sun, 11 Jan 2015 20:13:12 +0000 (+0000) Subject: Add function ST_EffectiveArea, Visvalingam’s algorithm simplification #2227 X-Git-Tag: 2.2.0rc1~710 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=fe5423a1d892476877f9973c183a3a074ff40508;p=postgis Add function ST_EffectiveArea, Visvalingam’s algorithm simplification #2227 git-svn-id: http://svn.osgeo.org/postgis/trunk@13174 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/doc/reference_processing.xml b/doc/reference_processing.xml index 164892928..357340ac6 100644 --- a/doc/reference_processing.xml +++ b/doc/reference_processing.xml @@ -2670,6 +2670,68 @@ FROM (SELECT ST_Buffer('POINT(1 3)', 10,12) As the_geom) As foo; + + + + ST_SetEffectiveArea + Sets for each vertex point it's effective area, + and can by filtring on this area return a simplified geometry + + + + + + geometry ST_SetEffectiveArea + geometry geomA + float threashold = 0 + + + + + + Description + Sets for each vertex point it's effective area from Visvalingam’s algorithm. + The effective area is stored as the M-value of the geomtries. + If the second optional parameter is used, the resulting geometriy will be build obnly on vertex points with an effective area + greater than or equal to that threashold value. That will be a simplified geometry. + + This function can be used for server side simplification by using the threashold. Another option is to not give any threashold value. + Then you get the full geometry back, but with effective areas as M-values wich can be used by the client to simplify very fast. + + Will actually do something only with + (multi)lines and (multi)polygons but you can safely call it with + any kind of geometry. Since simplification occurs on a + object-by-object basis you can also feed a GeometryCollection to + this function. + + This function handles 3D and the third dimmension will affect the effective area. + + Note that returned geometry might loose its + simplicity (see ) + Note topology may not be preserved and may result in invalid geometries. Use (see ) to preserve topology. + The output geoemtry will loose all previous information in the M-values + Availability: 2.2.0 + + + + Examples + A linestring that get the efffective area calculated. All points is returned since we give 0 as themin area threashold + + +select ST_AStext(ST_SetEffectiveArea(geom)) all_pts, ST_AStext(ST_SetEffectiveArea(geom,30) ) thrshld_30 +FROM (SELECT 'LINESTRING(5 2, 3 8, 6 20, 7 25, 10 10)'::geometry geom) As foo; +-result + all_pts | thrshld_30 +-----------+-------------------+ +LINESTRING M (5 2 1.79769313486232e+308,3 8 29,6 20 1.5,7 25 49.5,10 10 1.79769313486232e+308) | LINESTRING M (5 2 1.79769313486232e+308,7 25 49.5,10 10 1.79769313486232e+308) + + + + + See Also + , , Topology + + diff --git a/liblwgeom/Makefile.in b/liblwgeom/Makefile.in index 08140d740..fc27456c2 100644 --- a/liblwgeom/Makefile.in +++ b/liblwgeom/Makefile.in @@ -88,6 +88,7 @@ SA_OBJS = \ lwgeom_geos_node.o \ lwgeom_geos_split.o \ lwgeom_transform.o \ + effectivearea.o \ varint.o NM_OBJS = \ diff --git a/liblwgeom/effectivearea.c b/liblwgeom/effectivearea.c new file mode 100644 index 000000000..b7ea1e0e9 --- /dev/null +++ b/liblwgeom/effectivearea.c @@ -0,0 +1,334 @@ +/********************************************************************** + * + * PostGIS - Spatial Types for PostgreSQL + * http://postgis.net + * Copyright 2014 Nicklas Avén + * + * This is free software; you can redistribute and/or modify it under + * the terms of the GNU General Public Licence. See the COPYING file. + * + **********************************************************************/ + + #include "liblwgeom_internal.h" + #include "lwgeom_log.h" + +typedef struct +{ + double area; + int prev; + int next; +} areanode; + +/** + +Structure to hold pointarray and it's arealist +*/ +typedef struct +{ + const POINTARRAY *inpts; + areanode *initial_arealist; + double *res_arealist; +} EFFECTIVE_AREAS; + +/** + +Calculate the area of a triangle in 2d +*/ +static double triarea2d(const double *P1, const double *P2, const double *P3) +{ + LWDEBUG(2, "Entered triarea2d"); + double ax,bx,ay,by, area; + + ax=P1[0]-P2[0]; + bx=P3[0]-P2[0]; + ay=P1[1]-P2[1]; + by=P3[1]-P2[1]; + + area= fabs(0.5*(ax*by-ay*bx)); + + LWDEBUGF(2, "Calculated area is %f",(float) area); + return area; +} + +/** + +Calculate the area of a triangle in 3d space +*/ +static double triarea3d(const double *P1, const double *P2, const double *P3) +{ + LWDEBUG(2, "Entered triarea3d"); + double ax,bx,ay,by,az,bz,cx,cy,cz, area; + + ax=P1[0]-P2[0]; + bx=P3[0]-P2[0]; + ay=P1[1]-P2[1]; + by=P3[1]-P2[1]; + az=P1[2]-P2[2]; + bz=P3[2]-P2[2]; + + cx = ay*bz - az*by; + cy = az*bx - ax*bz; + cz = ax*by - ay*bx; + + area = fabs(0.5*(sqrt(cx*cx+cy*cy+cz*cz))); + + LWDEBUGF(2, "Calculated area is %f",(float) area); + return area; +} + +/** + +To get the effective area, we have to check what area a point results in when all smaller areas are eliminated +*/ +static void tune_areas(EFFECTIVE_AREAS ea) +{ + LWDEBUG(2, "Entered tune_areas"); + const double *P1; + const double *P2; + const double *P3; + double area; + + int npoints=ea.inpts->npoints; + int i; + int current, prevcurrent, nextcurrent; + int is3d = FLAGS_GET_Z(ea.inpts->flags); + + do + { + /*i is our cursor*/ + i=0; + + /*current is the point that currently gives the smallest area*/ + current=0; + + do + { + i=ea.initial_arealist[i].next; + if(ea.initial_arealist[i].area0) + { + + LWDEBUGF(2, "calculate previous, %d",prevcurrent); + P1= (double*)getPoint_internal(ea.inpts, ea.initial_arealist[prevcurrent].prev); + if(is3d) + area=triarea3d(P1, P2, P3); + else + area=triarea2d(P1, P2, P3); + ea.initial_arealist[prevcurrent].area = FP_MAX(area,ea.res_arealist[current]); + } + if(nextcurrent0); + return; +} + + +static void ptarray_calc_areas(EFFECTIVE_AREAS ea) +{ + LWDEBUG(2, "Entered ptarray_set_effective_area2d"); + int i; + int npoints=ea.inpts->npoints; + int is3d = FLAGS_GET_Z(ea.inpts->flags); + + const double *P1; + const double *P2; + const double *P3; + + P1 = (double*)getPoint_internal(ea.inpts, 0); + P2 = (double*)getPoint_internal(ea.inpts, 1); + + /*The first and last point shall always have the maximum effective area.*/ + ea.initial_arealist[0].area=ea.initial_arealist[npoints-1].area=DBL_MAX; + ea.res_arealist[0]=ea.res_arealist[npoints-1]=DBL_MAX; + + ea.initial_arealist[0].next=1; + ea.initial_arealist[0].prev=0; + + LWDEBUGF(2, "Time to iterate the pointlist with %d points",npoints); + for (i=1;i<(npoints)-1;i++) + { + ea.initial_arealist[i].next=i+1; + ea.initial_arealist[i].prev=i-1; + LWDEBUGF(2, "i=%d",i); + P3 = (double*)getPoint_internal(ea.inpts, i+1); + + if(is3d) + ea.initial_arealist[i].area=triarea3d(P1, P2, P3); + else + ea.initial_arealist[i].area=triarea2d(P1, P2, P3); + + LWDEBUGF(2, "area = %f",ea.initial_arealist[i].area); + P1=P2; + P2=P3; + + } + ea.initial_arealist[npoints-1].next=npoints-1; + ea.initial_arealist[npoints-1].prev=npoints-2; + LWDEBUG(2, "Time to go on"); + + for (i=1;i<(npoints)-1;i++) + { + ea.res_arealist[i]=-1; + } + + tune_areas(ea); + return ; +} + + +static POINTARRAY * ptarray_set_effective_area(POINTARRAY *inpts,double trshld) +{ + LWDEBUG(2, "Entered ptarray_set_effective_area"); + int p; + POINT4D pt; + EFFECTIVE_AREAS ea; + + ea.initial_arealist = lwalloc(inpts->npoints*sizeof(areanode)); + LWDEBUGF(2, "sizeof(areanode)=%d",sizeof(areanode)); + ea.res_arealist = lwalloc(inpts->npoints*sizeof(double)); + ea.inpts=inpts; + POINTARRAY *opts = ptarray_construct_empty(FLAGS_GET_Z(inpts->flags), 1, inpts->npoints); + + ptarray_calc_areas(ea); + + /*Only return points with an effective area above the threashold*/ + for (p=0;pnpoints;p++) + { + if(ea.res_arealist[p]>=trshld) + { + pt=getPoint4d(ea.inpts, p); + pt.m=ea.res_arealist[p]; + ptarray_append_point(opts, &pt, LW_FALSE); + } + } + + lwfree(ea.initial_arealist); + lwfree(ea.res_arealist); + return opts; + +} + +static LWLINE* lwline_set_effective_area(const LWLINE *iline, double trshld) +{ + LWDEBUG(2, "Entered lwline_set_effective_area"); + + LWLINE *oline = lwline_construct_empty(iline->srid, FLAGS_GET_Z(iline->flags), 1); + + /* Skip empty case */ + if( lwline_is_empty(iline) ) + return lwline_clone(iline); + + oline = lwline_construct(iline->srid, NULL, ptarray_set_effective_area(iline->points,trshld)); + + oline->type = iline->type; + return oline; + +} + + +static LWPOLY* lwpoly_set_effective_area(const LWPOLY *ipoly, double trshld) +{ + LWDEBUG(2, "Entered lwpoly_set_effective_area"); + int i; + LWPOLY *opoly = lwpoly_construct_empty(ipoly->srid, FLAGS_GET_Z(ipoly->flags), 1); + + + if( lwpoly_is_empty(ipoly) ) + return opoly; /* should we return NULL instead ? */ + + for (i = 0; i < ipoly->nrings; i++) + { + /* Add ring to simplified polygon */ + if( lwpoly_add_ring(opoly, ptarray_set_effective_area(ipoly->rings[i],trshld)) == LW_FAILURE ) + return NULL; + } + + + opoly->type = ipoly->type; + + if( lwpoly_is_empty(opoly) ) + return NULL; + + return opoly; + +} + + +static LWCOLLECTION* lwcollection_set_effective_area(const LWCOLLECTION *igeom, double trshld) +{ + LWDEBUG(2, "Entered lwcollection_set_effective_area"); + int i; + LWCOLLECTION *out = lwcollection_construct_empty(igeom->type, igeom->srid, FLAGS_GET_Z(igeom->flags), 1); + + if( lwcollection_is_empty(igeom) ) + return out; /* should we return NULL instead ? */ + + for( i = 0; i < igeom->ngeoms; i++ ) + { + LWGEOM *ngeom = lwgeom_set_effective_area(igeom->geoms[i],trshld); + if ( ngeom ) out = lwcollection_add_lwgeom(out, ngeom); + } + + return out; +} + + +LWGEOM* lwgeom_set_effective_area(const LWGEOM *igeom, double trshld) +{ + LWDEBUG(2, "Entered lwgeom_set_effective_area"); + switch (igeom->type) + { + case POINTTYPE: + case MULTIPOINTTYPE: + return lwgeom_clone(igeom); + case LINETYPE: + return (LWGEOM*)lwline_set_effective_area((LWLINE*)igeom, trshld); + case POLYGONTYPE: + return (LWGEOM*)lwpoly_set_effective_area((LWPOLY*)igeom, trshld); + case MULTILINETYPE: + case MULTIPOLYGONTYPE: + case COLLECTIONTYPE: + return (LWGEOM*)lwcollection_set_effective_area((LWCOLLECTION *)igeom, trshld); + default: + lwerror("lwgeom_simplify: unsupported geometry type: %s",lwtype_name(igeom->type)); + } + return NULL; +} + \ No newline at end of file diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in index f756f0a98..d9bc9dd6e 100644 --- a/liblwgeom/liblwgeom.h.in +++ b/liblwgeom/liblwgeom.h.in @@ -928,6 +928,7 @@ extern LWGEOM* lwgeom_force_3dm(const LWGEOM *geom); extern LWGEOM* lwgeom_force_4d(const LWGEOM *geom); extern LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist); +extern LWGEOM* lwgeom_set_effective_area(const LWGEOM *igeom, double area); /* * Force to use SFS 1.1 geometry type diff --git a/postgis/lwgeom_functions_analytic.c b/postgis/lwgeom_functions_analytic.c index c1f999001..60b32c29e 100644 --- a/postgis/lwgeom_functions_analytic.c +++ b/postgis/lwgeom_functions_analytic.c @@ -31,6 +31,7 @@ /* Prototypes */ Datum LWGEOM_simplify2d(PG_FUNCTION_ARGS); +Datum LWGEOM_SetEffectiveArea(PG_FUNCTION_ARGS); Datum ST_LineCrossingDirection(PG_FUNCTION_ARGS); double determineSide(POINT2D *seg1, POINT2D *seg2, POINT2D *point); @@ -67,6 +68,35 @@ Datum LWGEOM_simplify2d(PG_FUNCTION_ARGS) PG_RETURN_POINTER(result); } +PG_FUNCTION_INFO_V1(LWGEOM_SetEffectiveArea); +Datum LWGEOM_SetEffectiveArea(PG_FUNCTION_ARGS) +{ + GSERIALIZED *geom = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + GSERIALIZED *result; + int type = gserialized_get_type(geom); + LWGEOM *in; + LWGEOM *out; + double area; + + if ( type == POINTTYPE || type == MULTIPOINTTYPE ) + PG_RETURN_POINTER(geom); + + area = PG_GETARG_FLOAT8(1); + in = lwgeom_from_gserialized(geom); + + out = lwgeom_set_effective_area(in, area); + if ( ! out ) PG_RETURN_NULL(); + + /* COMPUTE_BBOX TAINTING */ + if ( in->bbox ) lwgeom_add_bbox(out); + + result = geometry_serialize(out); + lwgeom_free(out); + PG_FREE_IF_COPY(geom, 0); + PG_RETURN_POINTER(result); +} + + /*********************************************************************** * --strk@keybit.net; ***********************************************************************/ diff --git a/postgis/postgis.sql.in b/postgis/postgis.sql.in index 2c35bc06c..6f3cba824 100644 --- a/postgis/postgis.sql.in +++ b/postgis/postgis.sql.in @@ -2841,6 +2841,12 @@ CREATE OR REPLACE FUNCTION ST_Simplify(geometry, float8) RETURNS geometry AS 'MODULE_PATHNAME', 'LWGEOM_simplify2d' LANGUAGE 'c' IMMUTABLE STRICT; + +-- Availability: 2.2.0 +CREATE OR REPLACE FUNCTION ST_SetEffectiveArea(geometry, float8 default 0) + RETURNS geometry + AS 'MODULE_PATHNAME', 'LWGEOM_SetEffectiveArea' + LANGUAGE 'c' IMMUTABLE STRICT; -- ST_SnapToGrid(input, xoff, yoff, xsize, ysize) -- Availability: 1.2.2