# [postgis-users] [patch] GEOS simplify functions

François Weykmans fanfoi at gmail.com
Sun Nov 19 09:25:16 PST 2006

```Hello!

Add GEOS simplify functions to PostGIS
(testing only with PostGIS 1.1.6 and geos-3.0.0rc2)

With them, you can use big tolerance without get null result value...
SQL : simplify_geos(geometry,float8) = GEOS
DouglasPeuckerSimplifier::simplify(g1, tolerance)
SQL : simplifytopologypreserve_geos(geometry,float8) = GEOS
TopologyPreservingSimplifier::simplify(g1, tolerance)

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-users/attachments/20061119/7bd855c8/attachment.html>
```