<pre>Hello!<br><br>Add GEOS simplify functions to PostGIS <br>(testing only with PostGIS 1.1.6 and geos-3.0.0rc2)<br><br>With them, you can use big tolerance without get null result value...<br><br>SQL : simplify_geos(geometry,float8) = GEOS DouglasPeuckerSimplifier::simplify(g1, tolerance)
<br>SQL : simplifytopologypreserve_geos(geometry,float8) = GEOS TopologyPreservingSimplifier::simplify(g1, tolerance)<br><br><br>Thanks a lots for PostGIS!!!<br><br>François<br><br>diff -r -c /postgis-1.1.6/lwgeom/lwgeom_geos.c /new/postgis-
1.1.6/lwgeom/lwgeom_geos.c<br>*** /postgis-1.1.6/lwgeom/lwgeom_geos.c Mon Dec 12 09:40:50 2005<br>--- /new/postgis-<br>1.1.6/lwgeom/lwgeom_geos.c Fri Nov 17 08:30:15 2006<br>***************<br>*** 66,71 ****<br>--- 66,75 ----
<br> Datum linemerge(PG_FUNCTION_ARGS);<br> Datum LWGEOM_buildarea(PG_FUNCTION_ARGS);<br> <br>+<br>+ Datum simplify_geos(PG_FUNCTION_ARGS);<br><br>+ Datum simplifytopologypreserve_geos(PG_FUNCTION_ARGS);<br>+<br> <br>
typedef struct Geometry Geometry;<br> <br>***************<br>*** 101,106 ****<br>--- 105,115 ----<br> extern Geometry *GEOSGetCentroid(Geometry *g, int *failure);<br><br> extern bool GEOSHasZ(Geometry *g1);<br>
<br>+<br>+ extern Geometry *GEOSSimplify(Geometry *g1, double tolerance);<br>+ extern Geometry *GEOSTopologyPreserveSimplify(Geometry *g1, double tolerance);<br>+<br>+ <br> extern void GEOSSetSRID(Geometry *g, int SRID);
<br><br> extern void GEOSdeleteChar(char *a);<br> extern void GEOSdeleteGeometry(Geometry *a);<br>***************<br>*** 2896,2898 ****<br>--- 2905,3089 ----<br> PG_RETURN_NULL();<br> }<br> <br>+<br>+ PG_FUNCTION_INFO_V1(simplify_geos);
<br><br>+ Datum simplify_geos(PG_FUNCTION_ARGS)<br>+ {<br>+ PG_LWGEOM *geom1;<br>+ double tolerance;<br>+ GEOSGeom g1,g2;<br>+ PG_LWGEOM *result;<br>+ int is3d;<br>+ int SRID;<br>+ <br>+ #ifdef PROFILE<br>+ profstart(PROF_QRUN);
<br><br>+ #endif<br>+ <br>+ geom1 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));<br>+ tolerance = PG_GETARG_FLOAT8(1);<br>+ <br>+ is3d = ( TYPE_HASZ(geom1->type) );<br>+ SRID = pglwgeom_getSRID(geom1);<br>
+ initGEOS(lwnotice, lwnotice);<br><br>+ <br>+ #ifdef PGIS_DEBUG<br>+ elog(NOTICE,"simplify starting");<br>+ #endif<br>+ <br>+ #ifdef PROFILE<br>+ profstart(PROF_P2G1);<br>+ #endif<br>+ g1 = POSTGIS2GEOS(geom1);
<br>+ #ifdef PROFILE<br>+ profstop(PROF_P2G1);<br><br>+ #endif<br>+ <br>+ #ifdef PGIS_DEBUG<br>+ elog(NOTICE," constructed geometry - calling geos");<br>+ elog(NOTICE," g1 = %s", GEOSGeomToWKT(g1));
<br>+ /*elog(NOTICE,"g1 is valid = %i",GEOSisvalid(g1)); */<br><br>+ #endif<br>+ <br>+ #ifdef PROFILE<br>+ profstart(PROF_GRUN);<br>+ #endif <br>+ g2 = GEOSSimplify(g1,tolerance);<br>+ #ifdef PROFILE<br>+ profstop(PROF_GRUN);
<br>+ #endif<br>+ <br>+ #ifdef PGIS_DEBUG<br>+ elog(NOTICE,"simplify finished");<br><br>+ #endif<br>+ <br>+ if (g2 == NULL)<br>+ {<br>+ elog(ERROR,"GEOS Simplify() threw an error!");<br>+ GEOSGeom_destroy(g1);
<br>+ PG_RETURN_NULL(); /* never get here */<br>+ }<br>+ #ifdef PGIS_DEBUG<br><br>+ elog(NOTICE,"result: %s", GEOSGeomToWKT(g2) ) ;<br>+ #endif<br>+ GEOSSetSRID(g2, SRID);<br>+ #ifdef PROFILE<br>+ profstart(PROF_G2P);
<br>+ #endif<br>+ result = GEOS2POSTGIS(g2, is3d);<br>+ #ifdef PROFILE<br><br>+ profstop(PROF_G2P);<br>+ #endif<br>+ <br>+ if (result == NULL)<br>+ {<br>+ GEOSGeom_destroy(g1);<br>+ GEOSGeom_destroy(g2);<br>+ elog(ERROR,"GEOS Simplify() threw an error (result postgis geometry formation)!");
<br><br>+ PG_RETURN_NULL(); /* never get here */<br>+ }<br>+ GEOSGeom_destroy(g1);<br>+ GEOSGeom_destroy(g2);<br>+ <br>+ /* compressType(result); */<br>+ <br>+ #ifdef PROFILE<br>+ profstop(PROF_QRUN);<br>+ profreport("geos",geom1, result);
<br><br>+ #endif<br>+ PG_FREE_IF_COPY(geom1, 0);<br>+ <br>+ PG_RETURN_POINTER(result);<br>+ }<br>+ <br>+ PG_FUNCTION_INFO_V1(simplifytopologypreserve_geos);<br>+ Datum simplifytopologypreserve_geos(PG_FUNCTION_ARGS)<br>
+ {<br><br>+ PG_LWGEOM *geom1;<br>+ double tolerance;<br>+ GEOSGeom g1,g2;<br>+ PG_LWGEOM *result;<br>+ int is3d;<br>+ int SRID;<br>+ <br>+ #ifdef PROFILE<br>+ profstart(PROF_QRUN);<br>+ #endif<br>+ <br>+ geom1 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
<br><br>+ tolerance = PG_GETARG_FLOAT8(1);<br>+ <br>+ is3d = ( TYPE_HASZ(geom1->type) );<br>+ SRID = pglwgeom_getSRID(geom1);<br>+ initGEOS(lwnotice, lwnotice);<br>+ <br>+ #ifdef PGIS_DEBUG<br>+ elog(NOTICE,"simplify starting");
<br><br>+ #endif<br>+ <br>+ #ifdef PROFILE<br>+ profstart(PROF_P2G1);<br>+ #endif<br>+ g1 = POSTGIS2GEOS(geom1);<br>+ #ifdef PROFILE<br>+ profstop(PROF_P2G1);<br>+ #endif<br>+ <br>+ #ifdef PGIS_DEBUG<br>+ elog(NOTICE," constructed geometry - calling geos");
<br><br>+ elog(NOTICE," g1 = %s", GEOSGeomToWKT(g1));<br>+ /*elog(NOTICE,"g1 is valid = %i",GEOSisvalid(g1)); */<br>+ #endif<br>+ <br>+ #ifdef PROFILE<br>+ profstart(PROF_GRUN);<br>+ #endif <br>+ g2 = GEOSTopologyPreserveSimplify(g1,tolerance);
<br><br>+ #ifdef PROFILE<br>+ profstop(PROF_GRUN);<br>+ #endif<br>+ <br>+ #ifdef PGIS_DEBUG<br>+ elog(NOTICE,"simplify finished");<br>+ #endif<br>+ <br>+ if (g2 == NULL)<br>+ {<br>+ elog(ERROR,"GEOS TopologyPreserveSimplify() threw an error!");
<br><br>+ GEOSGeom_destroy(g1);<br>+ PG_RETURN_NULL(); /* never get here */<br>+ }<br>+ #ifdef PGIS_DEBUG<br>+ elog(NOTICE,"result: %s", GEOSGeomToWKT(g2) ) ;<br>+ #endif<br>+ GEOSSetSRID(g2, SRID);<br>+ #ifdef PROFILE
<br><br>+ profstart(PROF_G2P);<br>+ #endif<br>+ result = GEOS2POSTGIS(g2, is3d);<br>+ #ifdef PROFILE<br>+ profstop(PROF_G2P);<br>+ #endif<br>+ <br>+ if (result == NULL)<br>+ {<br>+ GEOSGeom_destroy(g1);<br>+ GEOSGeom_destroy(g2);
<br><br>+ elog(ERROR,"GEOS TopologyPreserveSimplify() threw an error (result postgis geometry formation)!");<br>+ PG_RETURN_NULL(); /* never get here */<br>+ }<br>+ GEOSGeom_destroy(g1);<br>+ GEOSGeom_destroy(g2);
<br><br>+ <br>+ /* compressType(result); */<br>+ <br>+ #ifdef PROFILE<br>+ profstop(PROF_QRUN);<br>+ profreport("geos",geom1, result);<br>+ #endif<br>+ PG_FREE_IF_COPY(geom1, 0);<br>+ <br>+ PG_RETURN_POINTER(result);
<br><br>+ }<br>+<br>+ <br>diff -r -c /postgis-1.1.6/lwgeom/lwgeom_geos_c.c /new/postgis-1.1.6/lwgeom/lwgeom_geos_c.c<br>*** /postgis-1.1.6/lwgeom/lwgeom_geos_c.c Wed Oct 25 03:23:00 2006<br>--- /new/postgis-1.1.6/lwgeom/lwgeom_geos_c.c Fri Nov 17 08:15:59 2006
<br><br>***************<br>*** 67,72 ****<br>--- 67,76 ----<br> Datum LWGEOM_buildarea(PG_FUNCTION_ARGS);<br> Datum linemerge(PG_FUNCTION_ARGS);<br> <br>+<br>+ Datum simplify_geos(PG_FUNCTION_ARGS);<br>+ Datum simplifytopologypreserve_geos(PG_FUNCTION_ARGS);
<br><br>+<br> <br> LWGEOM *GEOS2LWGEOM(GEOSGeom geom, char want3d);<br> PG_LWGEOM *GEOS2POSTGIS(GEOSGeom geom, char want3d);<br>***************<br>*** 2974,2976 ****<br>--- 2978,3162 ----<br> PG_RETURN_POINTER(result);
<br><br> <br> }<br>+ <br>+<br>+ PG_FUNCTION_INFO_V1(simplify_geos);<br>+ Datum simplify_geos(PG_FUNCTION_ARGS)<br>+ {<br>+ PG_LWGEOM *geom1;<br>+ double tolerance;<br>+ GEOSGeom g1,g2;<br>+ PG_LWGEOM *result;<br>+ int is3d;
<br><br>+ int SRID;<br>+ <br>+ #ifdef PROFILE<br>+ profstart(PROF_QRUN);<br>+ #endif<br>+ <br>+ geom1 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));<br>+ tolerance = PG_GETARG_FLOAT8(1);<br>+ <br>+ is3d = ( TYPE_HASZ(geom1->type) );
<br><br>+ SRID = pglwgeom_getSRID(geom1);<br>+ initGEOS(lwnotice, lwnotice);<br>+ <br>+ #ifdef PGIS_DEBUG<br>+ elog(NOTICE,"simplify starting");<br>+ #endif<br>+ <br>+ #ifdef PROFILE<br>+ profstart(PROF_P2G1);
<br><br>+ #endif<br>+ g1 = POSTGIS2GEOS(geom1);<br>+ #ifdef PROFILE<br>+ profstop(PROF_P2G1);<br>+ #endif<br>+ <br>+ #ifdef PGIS_DEBUG<br>+ elog(NOTICE," constructed geometry - calling geos");<br>+ elog(NOTICE," g1 = %s", GEOSGeomToWKT(g1));
<br><br>+ /*elog(NOTICE,"g1 is valid = %i",GEOSisvalid(g1)); */<br>+ #endif<br>+ <br>+ #ifdef PROFILE<br>+ profstart(PROF_GRUN);<br>+ #endif <br>+ g2 = GEOSSimplify(g1,tolerance);<br>+ #ifdef PROFILE<br>+ profstop(PROF_GRUN);
<br><br>+ #endif<br>+ <br>+ #ifdef PGIS_DEBUG<br>+ elog(NOTICE,"simplify finished");<br>+ #endif<br>+ <br>+ if (g2 == NULL)<br>+ {<br>+ elog(ERROR,"GEOS Simplify() threw an error!");<br>+ GEOSGeom_destroy(g1);
<br><br>+ PG_RETURN_NULL(); /* never get here */<br>+ }<br>+ #ifdef PGIS_DEBUG<br>+ elog(NOTICE,"result: %s", GEOSGeomToWKT(g2) ) ;<br>+ #endif<br>+ GEOSSetSRID(g2, SRID);<br>+ #ifdef PROFILE<br>+ profstart(PROF_G2P);
<br><br>+ #endif<br>+ result = GEOS2POSTGIS(g2, is3d);<br>+ #ifdef PROFILE<br>+ profstop(PROF_G2P);<br>+ #endif<br>+ <br>+ if (result == NULL)<br>+ {<br>+ GEOSGeom_destroy(g1);<br>+ GEOSGeom_destroy(g2);<br>+ elog(ERROR,"GEOS Simplify() threw an error (result postgis geometry formation)!");
<br><br>+ PG_RETURN_NULL(); /* never get here */<br>+ }<br>+ GEOSGeom_destroy(g1);<br>+ GEOSGeom_destroy(g2);<br>+ <br>+ /* compressType(result); */<br>+ <br>+ #ifdef PROFILE<br>+ profstop(PROF_QRUN);<br>+ profreport("geos",geom1, result);
<br><br>+ #endif<br>+ PG_FREE_IF_COPY(geom1, 0);<br>+ <br>+ PG_RETURN_POINTER(result);<br>+ }<br>+ <br>+ PG_FUNCTION_INFO_V1(simplifytopologypreserve_geos);<br>+ Datum simplifytopologypreserve_geos(PG_FUNCTION_ARGS)<br>
+ {<br><br>+ PG_LWGEOM *geom1;<br>+ double tolerance;<br>+ GEOSGeom g1,g2;<br>+ PG_LWGEOM *result;<br>+ int is3d;<br>+ int SRID;<br>+ <br>+ #ifdef PROFILE<br>+ profstart(PROF_QRUN);<br>+ #endif<br>+ <br>+ geom1 = (PG_LWGEOM *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
<br><br>+ tolerance = PG_GETARG_FLOAT8(1);<br>+ <br>+ is3d = ( TYPE_HASZ(geom1->type) );<br>+ SRID = pglwgeom_getSRID(geom1);<br>+ initGEOS(lwnotice, lwnotice);<br>+ <br>+ #ifdef PGIS_DEBUG<br>+ elog(NOTICE,"simplify starting");
<br><br>+ #endif<br>+ <br>+ #ifdef PROFILE<br>+ profstart(PROF_P2G1);<br>+ #endif<br>+ g1 = POSTGIS2GEOS(geom1);<br>+ #ifdef PROFILE<br>+ profstop(PROF_P2G1);<br>+ #endif<br>+ <br>+ #ifdef PGIS_DEBUG<br>+ elog(NOTICE," constructed geometry - calling geos");
<br><br>+ elog(NOTICE," g1 = %s", GEOSGeomToWKT(g1));<br>+ /*elog(NOTICE,"g1 is valid = %i",GEOSisvalid(g1)); */<br>+ #endif<br>+ <br>+ #ifdef PROFILE<br>+ profstart(PROF_GRUN);<br>+ #endif <br>+ g2 = GEOSTopologyPreserveSimplify(g1,tolerance);
<br><br>+ #ifdef PROFILE<br>+ profstop(PROF_GRUN);<br>+ #endif<br>+ <br>+ #ifdef PGIS_DEBUG<br>+ elog(NOTICE,"simplify finished");<br>+ #endif<br>+ <br>+ if (g2 == NULL)<br>+ {<br>+ elog(ERROR,"GEOS TopologyPreserveSimplify() threw an error!");
<br><br>+ GEOSGeom_destroy(g1);<br>+ PG_RETURN_NULL(); /* never get here */<br>+ }<br>+ #ifdef PGIS_DEBUG<br>+ elog(NOTICE,"result: %s", GEOSGeomToWKT(g2) ) ;<br>+ #endif<br>+ GEOSSetSRID(g2, SRID);<br>+ #ifdef PROFILE
<br><br>+ profstart(PROF_G2P);<br>+ #endif<br>+ result = GEOS2POSTGIS(g2, is3d);<br>+ #ifdef PROFILE<br>+ profstop(PROF_G2P);<br>+ #endif<br>+ <br>+ if (result == NULL)<br>+ {<br>+ GEOSGeom_destroy(g1);<br>+ GEOSGeom_destroy(g2);
<br><br>+ elog(ERROR,"GEOS TopologyPreserveSimplify() threw an error (result postgis geometry formation)!");<br>+ PG_RETURN_NULL(); /* never get here */<br>+ }<br>+ GEOSGeom_destroy(g1);<br>+ GEOSGeom_destroy(g2);
<br><br>+ <br>+ /* compressType(result); */<br>+ <br>+ #ifdef PROFILE<br>+ profstop(PROF_QRUN);<br>+ profreport("geos",geom1, result);<br>+ #endif<br>+ PG_FREE_IF_COPY(geom1, 0);<br>+ <br>+ PG_RETURN_POINTER(result);
<br><br>+ }<br>+<br>diff -r -c /postgis-1.1.6/lwgeom/lwgeom_geos_wrapper.cpp /new/postgis-1.1.6/lwgeom/lwgeom_geos_wrapper.cpp<br>*** /postgis-1.1.6/lwgeom/lwgeom_geos_wrapper.cpp Fri Dec 16 00:56:15 2005<br>--- /new/postgis-
<br>1.1.6/lwgeom/lwgeom_geos_wrapper.cpp Fri Nov 17 07:54:29 2006<br>***************<br>*** 10,15 ****<br>--- 10,21 ----<br> #include "postgis_geos_version.h"<br> #include "geos/geom.h"<br> #include "geos/util.h"
<br><br>+ <br>+<br>+ //#include <geos/simplify/DouglasPeuckerSimplifier.h><br>+ //#include <geos/simplify/TopologyPreservingSimplifier.h><br>+<br>+ <br> #if GEOS_FIRST_INTERFACE <= 3 && GEOS_LAST_INTERFACE >= 3
<br><br> #include "geos/opValid.h"<br> #include "geos/opPolygonize.h"<br>***************<br>*** 140,145 ****<br>--- 146,155 ----<br> extern "C" Geometry *GEOSSymDifference(Geometry *g1,Geometry *g2);
<br><br> extern "C" Geometry *GEOSUnion(Geometry *g1,Geometry *g2);<br> <br>+<br>+ extern "C" Geometry *GEOSSimplify(Geometry *g1, double tolerance);<br>+ extern "C" Geometry *GEOSTopologyPreserveSimplify(Geometry *g1, double tolerance);
<br><br>+<br> <br> extern "C" POINT3D *GEOSGetCoordinate(Geometry *g1);<br> extern "C" POINT3D *GEOSGetCoordinates_Polygon(Polygon *g1);<br>***************<br>*** 1857,1859 ****<br>--- 1867,1917 ----
<br><br> return NULL;<br> }<br> }<br>+ <br>+<br>+ Geometry *<br>+ GEOSSimplify(const Geometry *g1, double tolerance)<br>+ {<br>+ using namespace geos::simplify;<br>+ <br>+ try<br>+ {<br>+ Geometry::AutoPtr g(DouglasPeuckerSimplifier::simplify(g1, tolerance));
<br><br>+ return g.release();<br>+ }<br>+ catch (const std::exception &e)<br>+ {<br>+ ERROR_MESSAGE("%s", e.what());<br>+ return NULL;<br>+ }<br>+ <br>+ catch (...)<br>+ {<br>+ ERROR_MESSAGE("Unknown exception thrown");
<br><br>+ return NULL;<br>+ }<br>+ }<br>+ <br>+ Geometry *<br>+ GEOSTopologyPreserveSimplify(const Geometry *g1, double tolerance)<br>+ {<br>+ using namespace geos::simplify;<br>+ <br>+ try<br>+ {<br>+ Geometry::AutoPtr g(TopologyPreservingSimplifier::simplify(g1, tolerance));
<br><br>+ return g.release();<br>+ }<br>+ catch (const std::exception &e)<br>+ {<br>+ ERROR_MESSAGE("%s", e.what());<br>+ return NULL;<br>+ }<br>+ <br>+ catch (...)<br>+ {<br>+ ERROR_MESSAGE("Unknown exception thrown");
<br><br>+ return NULL;<br>+ }<br>+ }<br>+<br>diff -r -c /postgis-1.1.6/lwgeom/lwpostgis.sql.in /new/postgis-1.1.6/lwgeom/lwpostgis.sql.in<br>*** /postgis-1.1.6/lwgeom/lwpostgis.sql.in Fri Jul 7 06:56:52 2006<br>--- /new/postgis-
<br>1.1.6/lwgeom/lwpostgis.sql.in Fri Nov 17 08:55:30 2006<br>***************<br>*** 2876,2881 ****<br>--- 2876,2893 ----<br> -- GEOS<br> ---------------------------------------------------------------<br> <br>+<br>+ CREATEFUNCTION simplify_geos(geometry,float8)
<br><br>+ RETURNS geometry<br>+ AS '@MODULE_FILENAME@','simplify_geos'<br>+ LANGUAGE 'C' _IMMUTABLE_STRICT; -- WITH (isstrict,iscachable)<br>+ <br>+ CREATEFUNCTION simplifytopologypreserve_geos(geometry,float8)<br>
+ RETURNS geometry<br><br>+ AS '@MODULE_FILENAME@','simplifytopologypreserve_geos'<br>+ LANGUAGE 'C' _IMMUTABLE_STRICT; -- WITH (isstrict,iscachable)<br>+<br>+ <br> CREATEFUNCTION intersection(geometry,geometry)
<br> RETURNS geometry<br><br> AS '@MODULE_FILENAME@','intersection'<br></pre>