[postgis-devel] [patch] GEOS simplify functions
François Weykmans
fanfoi at gmail.com
Wed Dec 13 09:26:06 PST 2006
Hello!
This patch add GEOS simplify functions to PostGIS
(tested only with PostGIS 1.1.6 and geos-3.0.0rc2 but I can build path
to PostGIS 1.2.0)
With it, you can use large tolerance values without getting null
result values...
SQL : simplify_geos(geometry,float8) = GEOS
DouglasPeuckerSimplifier::simplify(g1, tolerance)
SQL : simplifytopologypreserve_geos(geometry,float8) = GEOS
TopologyPreservingSimplifier::simplify(g1, tolerance)
Example : I use simplify_geos into online SVG Map application (demo
here : http://siti.cmm.qc.ca/)
PostGIS 1.1.6/Geos-3.0.0rc2/PostgreSQL 8.2
Could this patch be included by default into PostGIS?
Thanks a lots for PostGIS!!!
François
diff -r -c /postgis-1.1.6/lwgeom/lwgeom_geos.c
/new/postgis-1.1.6/lwgeom/lwgeom_geos.c
*** /postgis-1.1.6/lwgeom/lwgeom_geos.c Mon Dec 12 09:40:50 2005
--- /new/postgis-1.1.6/lwgeom/lwgeom_geos.c Fri Nov 17 08:30:15 2006
***************
*** 66,71 ****
--- 66,75 ----
Datum linemerge(PG_FUNCTION_ARGS);
Datum LWGEOM_buildarea(PG_FUNCTION_ARGS);
+
+ Datum simplify_geos(PG_FUNCTION_ARGS);
+ Datum simplifytopologypreserve_geos(PG_FUNCTION_ARGS);
+
typedef struct Geometry Geometry;
***************
*** 101,106 ****
--- 105,115 ----
extern Geometry *GEOSGetCentroid(Geometry *g, int *failure);
extern bool GEOSHasZ(Geometry *g1);
+
+ extern Geometry *GEOSSimplify(Geometry *g1, double tolerance);
+ extern Geometry *GEOSTopologyPreserveSimplify(Geometry *g1, double tolerance);
+
+
extern void GEOSSetSRID(Geometry *g, int SRID);
extern void GEOSdeleteChar(char *a);
extern void GEOSdeleteGeometry(Geometry *a);
***************
*** 2896,2898 ****
--- 2905,3089 ----
PG_RETURN_NULL();
}
+
+ PG_FUNCTION_INFO_V1(simplify_geos);
+ Datum simplify_geos(PG_FUNCTION_ARGS)
+ {
+ PG_LWGEOM *geom1;
+ double tolerance;
+ GEOSGeom g1,g2;
+ PG_LWGEOM *result;
+ int is3d;
+ int SRID;
+
+ #ifdef PROFILE
+ profstart(PROF_QRUN);
+ #endif
+
+ geom1 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ tolerance = PG_GETARG_FLOAT8(1);
+
+ is3d = ( TYPE_HASZ(geom1->type) );
+ SRID = pglwgeom_getSRID(geom1);
+ initGEOS(lwnotice, lwnotice);
+
+ #ifdef PGIS_DEBUG
+ elog(NOTICE,"simplify starting");
+ #endif
+
+ #ifdef PROFILE
+ profstart(PROF_P2G1);
+ #endif
+ g1 = POSTGIS2GEOS(geom1);
+ #ifdef PROFILE
+ profstop(PROF_P2G1);
+ #endif
+
+ #ifdef PGIS_DEBUG
+ elog(NOTICE," constructed geometry - calling geos");
+ elog(NOTICE," g1 = %s", GEOSGeomToWKT(g1));
+ /*elog(NOTICE,"g1 is valid = %i",GEOSisvalid(g1)); */
+ #endif
+
+ #ifdef PROFILE
+ profstart(PROF_GRUN);
+ #endif
+ g2 = GEOSSimplify(g1,tolerance);
+ #ifdef PROFILE
+ profstop(PROF_GRUN);
+ #endif
+
+ #ifdef PGIS_DEBUG
+ elog(NOTICE,"simplify finished");
+ #endif
+
+ if (g2 == NULL)
+ {
+ elog(ERROR,"GEOS Simplify() threw an error!");
+ GEOSGeom_destroy(g1);
+ PG_RETURN_NULL(); /* never get here */
+ }
+ #ifdef PGIS_DEBUG
+ elog(NOTICE,"result: %s", GEOSGeomToWKT(g2) ) ;
+ #endif
+ GEOSSetSRID(g2, SRID);
+ #ifdef PROFILE
+ profstart(PROF_G2P);
+ #endif
+ result = GEOS2POSTGIS(g2, is3d);
+ #ifdef PROFILE
+ profstop(PROF_G2P);
+ #endif
+
+ if (result == NULL)
+ {
+ GEOSGeom_destroy(g1);
+ GEOSGeom_destroy(g2);
+ elog(ERROR,"GEOS Simplify() threw an error (result postgis
geometry formation)!");
+ PG_RETURN_NULL(); /* never get here */
+ }
+ GEOSGeom_destroy(g1);
+ GEOSGeom_destroy(g2);
+
+ /* compressType(result); */
+
+ #ifdef PROFILE
+ profstop(PROF_QRUN);
+ profreport("geos",geom1, result);
+ #endif
+ PG_FREE_IF_COPY(geom1, 0);
+
+ PG_RETURN_POINTER(result);
+ }
+
+ PG_FUNCTION_INFO_V1(simplifytopologypreserve_geos);
+ Datum simplifytopologypreserve_geos(PG_FUNCTION_ARGS)
+ {
+ PG_LWGEOM *geom1;
+ double tolerance;
+ GEOSGeom g1,g2;
+ PG_LWGEOM *result;
+ int is3d;
+ int SRID;
+
+ #ifdef PROFILE
+ profstart(PROF_QRUN);
+ #endif
+
+ geom1 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ tolerance = PG_GETARG_FLOAT8(1);
+
+ is3d = ( TYPE_HASZ(geom1->type) );
+ SRID = pglwgeom_getSRID(geom1);
+ initGEOS(lwnotice, lwnotice);
+
+ #ifdef PGIS_DEBUG
+ elog(NOTICE,"simplify starting");
+ #endif
+
+ #ifdef PROFILE
+ profstart(PROF_P2G1);
+ #endif
+ g1 = POSTGIS2GEOS(geom1);
+ #ifdef PROFILE
+ profstop(PROF_P2G1);
+ #endif
+
+ #ifdef PGIS_DEBUG
+ elog(NOTICE," constructed geometry - calling geos");
+ elog(NOTICE," g1 = %s", GEOSGeomToWKT(g1));
+ /*elog(NOTICE,"g1 is valid = %i",GEOSisvalid(g1)); */
+ #endif
+
+ #ifdef PROFILE
+ profstart(PROF_GRUN);
+ #endif
+ g2 = GEOSTopologyPreserveSimplify(g1,tolerance);
+ #ifdef PROFILE
+ profstop(PROF_GRUN);
+ #endif
+
+ #ifdef PGIS_DEBUG
+ elog(NOTICE,"simplify finished");
+ #endif
+
+ if (g2 == NULL)
+ {
+ elog(ERROR,"GEOS TopologyPreserveSimplify() threw an error!");
+ GEOSGeom_destroy(g1);
+ PG_RETURN_NULL(); /* never get here */
+ }
+ #ifdef PGIS_DEBUG
+ elog(NOTICE,"result: %s", GEOSGeomToWKT(g2) ) ;
+ #endif
+ GEOSSetSRID(g2, SRID);
+ #ifdef PROFILE
+ profstart(PROF_G2P);
+ #endif
+ result = GEOS2POSTGIS(g2, is3d);
+ #ifdef PROFILE
+ profstop(PROF_G2P);
+ #endif
+
+ if (result == NULL)
+ {
+ GEOSGeom_destroy(g1);
+ GEOSGeom_destroy(g2);
+ elog(ERROR,"GEOS TopologyPreserveSimplify() threw an error (result
postgis geometry formation)!");
+ PG_RETURN_NULL(); /* never get here */
+ }
+ GEOSGeom_destroy(g1);
+ GEOSGeom_destroy(g2);
+
+ /* compressType(result); */
+
+ #ifdef PROFILE
+ profstop(PROF_QRUN);
+ profreport("geos",geom1, result);
+ #endif
+ PG_FREE_IF_COPY(geom1, 0);
+
+ PG_RETURN_POINTER(result);
+ }
+
+
diff -r -c /postgis-1.1.6/lwgeom/lwgeom_geos_c.c
/new/postgis-1.1.6/lwgeom/lwgeom_geos_c.c
*** /postgis-1.1.6/lwgeom/lwgeom_geos_c.c Wed Oct 25 03:23:00 2006
--- /new/postgis-1.1.6/lwgeom/lwgeom_geos_c.c Fri Nov 17 08:15:59 2006
***************
*** 67,72 ****
--- 67,76 ----
Datum LWGEOM_buildarea(PG_FUNCTION_ARGS);
Datum linemerge(PG_FUNCTION_ARGS);
+
+ Datum simplify_geos(PG_FUNCTION_ARGS);
+ Datum simplifytopologypreserve_geos(PG_FUNCTION_ARGS);
+
LWGEOM *GEOS2LWGEOM(GEOSGeom geom, char want3d);
PG_LWGEOM *GEOS2POSTGIS(GEOSGeom geom, char want3d);
***************
*** 2974,2976 ****
--- 2978,3162 ----
PG_RETURN_POINTER(result);
}
+
+
+ PG_FUNCTION_INFO_V1(simplify_geos);
+ Datum simplify_geos(PG_FUNCTION_ARGS)
+ {
+ PG_LWGEOM *geom1;
+ double tolerance;
+ GEOSGeom g1,g2;
+ PG_LWGEOM *result;
+ int is3d;
+ int SRID;
+
+ #ifdef PROFILE
+ profstart(PROF_QRUN);
+ #endif
+
+ geom1 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ tolerance = PG_GETARG_FLOAT8(1);
+
+ is3d = ( TYPE_HASZ(geom1->type) );
+ SRID = pglwgeom_getSRID(geom1);
+ initGEOS(lwnotice, lwnotice);
+
+ #ifdef PGIS_DEBUG
+ elog(NOTICE,"simplify starting");
+ #endif
+
+ #ifdef PROFILE
+ profstart(PROF_P2G1);
+ #endif
+ g1 = POSTGIS2GEOS(geom1);
+ #ifdef PROFILE
+ profstop(PROF_P2G1);
+ #endif
+
+ #ifdef PGIS_DEBUG
+ elog(NOTICE," constructed geometry - calling geos");
+ elog(NOTICE," g1 = %s", GEOSGeomToWKT(g1));
+ /*elog(NOTICE,"g1 is valid = %i",GEOSisvalid(g1)); */
+ #endif
+
+ #ifdef PROFILE
+ profstart(PROF_GRUN);
+ #endif
+ g2 = GEOSSimplify(g1,tolerance);
+ #ifdef PROFILE
+ profstop(PROF_GRUN);
+ #endif
+
+ #ifdef PGIS_DEBUG
+ elog(NOTICE,"simplify finished");
+ #endif
+
+ if (g2 == NULL)
+ {
+ elog(ERROR,"GEOS Simplify() threw an error!");
+ GEOSGeom_destroy(g1);
+ PG_RETURN_NULL(); /* never get here */
+ }
+ #ifdef PGIS_DEBUG
+ elog(NOTICE,"result: %s", GEOSGeomToWKT(g2) ) ;
+ #endif
+ GEOSSetSRID(g2, SRID);
+ #ifdef PROFILE
+ profstart(PROF_G2P);
+ #endif
+ result = GEOS2POSTGIS(g2, is3d);
+ #ifdef PROFILE
+ profstop(PROF_G2P);
+ #endif
+
+ if (result == NULL)
+ {
+ GEOSGeom_destroy(g1);
+ GEOSGeom_destroy(g2);
+ elog(ERROR,"GEOS Simplify() threw an error (result postgis
geometry formation)!");
+ PG_RETURN_NULL(); /* never get here */
+ }
+ GEOSGeom_destroy(g1);
+ GEOSGeom_destroy(g2);
+
+ /* compressType(result); */
+
+ #ifdef PROFILE
+ profstop(PROF_QRUN);
+ profreport("geos",geom1, result);
+ #endif
+ PG_FREE_IF_COPY(geom1, 0);
+
+ PG_RETURN_POINTER(result);
+ }
+
+ PG_FUNCTION_INFO_V1(simplifytopologypreserve_geos);
+ Datum simplifytopologypreserve_geos(PG_FUNCTION_ARGS)
+ {
+ PG_LWGEOM *geom1;
+ double tolerance;
+ GEOSGeom g1,g2;
+ PG_LWGEOM *result;
+ int is3d;
+ int SRID;
+
+ #ifdef PROFILE
+ profstart(PROF_QRUN);
+ #endif
+
+ geom1 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ tolerance = PG_GETARG_FLOAT8(1);
+
+ is3d = ( TYPE_HASZ(geom1->type) );
+ SRID = pglwgeom_getSRID(geom1);
+ initGEOS(lwnotice, lwnotice);
+
+ #ifdef PGIS_DEBUG
+ elog(NOTICE,"simplify starting");
+ #endif
+
+ #ifdef PROFILE
+ profstart(PROF_P2G1);
+ #endif
+ g1 = POSTGIS2GEOS(geom1);
+ #ifdef PROFILE
+ profstop(PROF_P2G1);
+ #endif
+
+ #ifdef PGIS_DEBUG
+ elog(NOTICE," constructed geometry - calling geos");
+ elog(NOTICE," g1 = %s", GEOSGeomToWKT(g1));
+ /*elog(NOTICE,"g1 is valid = %i",GEOSisvalid(g1)); */
+ #endif
+
+ #ifdef PROFILE
+ profstart(PROF_GRUN);
+ #endif
+ g2 = GEOSTopologyPreserveSimplify(g1,tolerance);
+ #ifdef PROFILE
+ profstop(PROF_GRUN);
+ #endif
+
+ #ifdef PGIS_DEBUG
+ elog(NOTICE,"simplify finished");
+ #endif
+
+ if (g2 == NULL)
+ {
+ elog(ERROR,"GEOS TopologyPreserveSimplify() threw an error!");
+ GEOSGeom_destroy(g1);
+ PG_RETURN_NULL(); /* never get here */
+ }
+ #ifdef PGIS_DEBUG
+ elog(NOTICE,"result: %s", GEOSGeomToWKT(g2) ) ;
+ #endif
+ GEOSSetSRID(g2, SRID);
+ #ifdef PROFILE
+ profstart(PROF_G2P);
+ #endif
+ result = GEOS2POSTGIS(g2, is3d);
+ #ifdef PROFILE
+ profstop(PROF_G2P);
+ #endif
+
+ if (result == NULL)
+ {
+ GEOSGeom_destroy(g1);
+ GEOSGeom_destroy(g2);
+ elog(ERROR,"GEOS TopologyPreserveSimplify() threw an error (result
postgis geometry formation)!");
+ PG_RETURN_NULL(); /* never get here */
+ }
+ GEOSGeom_destroy(g1);
+ GEOSGeom_destroy(g2);
+
+ /* compressType(result); */
+
+ #ifdef PROFILE
+ profstop(PROF_QRUN);
+ profreport("geos",geom1, result);
+ #endif
+ PG_FREE_IF_COPY(geom1, 0);
+
+ PG_RETURN_POINTER(result);
+ }
+
diff -r -c /postgis-1.1.6/lwgeom/lwgeom_geos_wrapper.cpp
/new/postgis-1.1.6/lwgeom/lwgeom_geos_wrapper.cpp
*** /postgis-1.1.6/lwgeom/lwgeom_geos_wrapper.cpp Fri Dec 16 00:56:15 2005
--- /new/postgis-1.1.6/lwgeom/lwgeom_geos_wrapper.cpp Fri Nov 17 07:54:29 2006
***************
*** 10,15 ****
--- 10,21 ----
#include "postgis_geos_version.h"
#include "geos/geom.h"
#include "geos/util.h"
+
+
+ //#include <geos/simplify/DouglasPeuckerSimplifier.h>
+ //#include <geos/simplify/TopologyPreservingSimplifier.h>
+
+
#if GEOS_FIRST_INTERFACE <= 3 && GEOS_LAST_INTERFACE >= 3
#include "geos/opValid.h"
#include "geos/opPolygonize.h"
***************
*** 140,145 ****
--- 146,155 ----
extern "C" Geometry *GEOSSymDifference(Geometry *g1,Geometry *g2);
extern "C" Geometry *GEOSUnion(Geometry *g1,Geometry *g2);
+
+ extern "C" Geometry *GEOSSimplify(Geometry *g1, double tolerance);
+ extern "C" Geometry *GEOSTopologyPreserveSimplify(Geometry *g1,
double tolerance);
+
extern "C" POINT3D *GEOSGetCoordinate(Geometry *g1);
extern "C" POINT3D *GEOSGetCoordinates_Polygon(Polygon *g1);
***************
*** 1857,1859 ****
--- 1867,1917 ----
return NULL;
}
}
+
+
+ Geometry *
+ GEOSSimplify(const Geometry *g1, double tolerance)
+ {
+ using namespace geos::simplify;
+
+ try
+ {
+ Geometry::AutoPtr g(DouglasPeuckerSimplifier::simplify(g1, tolerance));
+ return g.release();
+ }
+ catch (const std::exception &e)
+ {
+ ERROR_MESSAGE("%s", e.what());
+ return NULL;
+ }
+
+ catch (...)
+ {
+ ERROR_MESSAGE("Unknown exception thrown");
+ return NULL;
+ }
+ }
+
+ Geometry *
+ GEOSTopologyPreserveSimplify(const Geometry *g1, double tolerance)
+ {
+ using namespace geos::simplify;
+
+ try
+ {
+ Geometry::AutoPtr g(TopologyPreservingSimplifier::simplify(g1, tolerance));
+ return g.release();
+ }
+ catch (const std::exception &e)
+ {
+ ERROR_MESSAGE("%s", e.what());
+ return NULL;
+ }
+
+ catch (...)
+ {
+ ERROR_MESSAGE("Unknown exception thrown");
+ return NULL;
+ }
+ }
+
diff -r -c /postgis-1.1.6/lwgeom/lwpostgis.sql.in
/new/postgis-1.1.6/lwgeom/lwpostgis.sql.in
*** /postgis-1.1.6/lwgeom/lwpostgis.sql.in Fri Jul 7 06:56:52 2006
--- /new/postgis-1.1.6/lwgeom/lwpostgis.sql.in Fri Nov 17 08:55:30 2006
***************
*** 2876,2881 ****
--- 2876,2893 ----
-- GEOS
---------------------------------------------------------------
+
+ CREATEFUNCTION simplify_geos(geometry,float8)
+ RETURNS geometry
+ AS '@MODULE_FILENAME@','simplify_geos'
+ LANGUAGE 'C' _IMMUTABLE_STRICT; -- WITH (isstrict,iscachable)
+
+ CREATEFUNCTION simplifytopologypreserve_geos(geometry,float8)
+ RETURNS geometry
+ AS '@MODULE_FILENAME@','simplifytopologypreserve_geos'
+ LANGUAGE 'C' _IMMUTABLE_STRICT; -- WITH (isstrict,iscachable)
+
+
CREATEFUNCTION intersection(geometry,geometry)
RETURNS geometry
AS '@MODULE_FILENAME@','intersection'
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-devel/attachments/20061213/998ce04a/attachment.html>
More information about the postgis-devel
mailing list