]> granicus.if.org Git - postgis/commitdiff
GEOS support -first run
authorDavid Blasby <dblasby@gmail.com>
Wed, 9 Apr 2003 18:34:43 +0000 (18:34 +0000)
committerDavid Blasby <dblasby@gmail.com>
Wed, 9 Apr 2003 18:34:43 +0000 (18:34 +0000)
git-svn-id: http://svn.osgeo.org/postgis/trunk@256 b70326c6-7e19-0410-871a-916f4a2858ee

Attic/postgis_sql_common.sql.in
Makefile
postgis_geos.c [new file with mode: 0644]
postgis_geos_wrapper.cpp [new file with mode: 0644]

index 796dc5da68c8b1a8a66049da1d16c319f1070d51..b9023ed247eb6befdad173b5c5e86bc3655abc40 100644 (file)
@@ -815,3 +815,55 @@ CREATE OPERATOR > (
    RESTRICT = contsel, JOIN = contjoinsel
 );
 
+
+--- GEOS FUNCTIONS
+CREATE FUNCTION relate(geometry,geometry)
+   RETURNS text
+   AS '@MODULE_FILENAME@','relate_full'
+   LANGUAGE 'c' with (isstrict);
+
+CREATE FUNCTION relate(geometry,geometry,text)
+   RETURNS boolean
+   AS '@MODULE_FILENAME@','relate_pattern'
+   LANGUAGE 'c' with (isstrict);
+
+CREATE FUNCTION disjoint(geometry,geometry)
+   RETURNS boolean
+   AS '@MODULE_FILENAME@'
+   LANGUAGE 'c' with (isstrict);
+
+CREATE FUNCTION touches(geometry,geometry)
+   RETURNS boolean
+   AS '@MODULE_FILENAME@'
+   LANGUAGE 'c' with (isstrict);
+
+CREATE FUNCTION intersects(geometry,geometry)
+   RETURNS boolean
+   AS '@MODULE_FILENAME@'
+   LANGUAGE 'c' with (isstrict);
+
+CREATE FUNCTION crosses(geometry,geometry)
+   RETURNS boolean
+   AS '@MODULE_FILENAME@'
+   LANGUAGE 'c' with (isstrict);
+
+CREATE FUNCTION within(geometry,geometry)
+   RETURNS boolean
+   AS '@MODULE_FILENAME@'
+   LANGUAGE 'c' with (isstrict);
+
+CREATE FUNCTION contains(geometry,geometry)
+   RETURNS boolean
+   AS '@MODULE_FILENAME@'
+   LANGUAGE 'c' with (isstrict);
+
+CREATE FUNCTION overlaps(geometry,geometry)
+   RETURNS boolean
+   AS '@MODULE_FILENAME@'
+   LANGUAGE 'c' with (isstrict);
+
+  CREATE FUNCTION isvalid(geometry)
+   RETURNS boolean
+   AS '@MODULE_FILENAME@'
+   LANGUAGE 'c' with (isstrict);
+   
\ No newline at end of file
index 3fd37ed7da0d6350bf7c9652099db4216d8576da..9defe7bc5df56ef828b6e81f7a3d65b7af15df0d 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -3,9 +3,16 @@
 #---------------------------------------------------------------
 # Set USE_PROJ to 1 for Proj4 reprojection support
 #
-USE_PROJ=0
+USE_PROJ=1
 PROJ_DIR=/usr/local
 
+#---------------------------------------------------------------
+# Set USE_GEOS to 1 for GEOS spatial predicate and operator
+# support
+#
+USE_GEOS=1
+GEOS_DIR=/usr/local
+
 #---------------------------------------------------------------
 # Set USE_STATS to 1 for new GiST statistics collection support
 # Note that this support requires additional columns in 
@@ -66,8 +73,13 @@ ifeq ($(USE_PROJ),1)
 else
        override CPPFLAGS := -g -I$(srcdir) $(CPPFLAGS) -DFRONTEND -DSYSCONFDIR='"$(sysconfdir)"' -DUSE_VERSION=$(USE_VERSION)
 endif
+ifeq ($(USE_GEOS),1)
+       override CPPFLAGS := $(CPPFLAGS) -I$(GEOS_DIR)/include/geos -DUSE_GEOS
+endif
 override DLLLIBS := $(BE_DLLLIBS) $(DLLLIBS)
 
+CXXFLAGS=-I$(PROJ_DIR)/include -I$(GEOS_DIR)/include/geos
+
 #---------------------------------------------------------------
 # Add index selectivity to C flags
 #
@@ -86,29 +98,36 @@ else
        GIST_ESTIMATE=postgis_estimate.o
 endif
 
-OBJS=postgis_debug.o postgis_ops.o postgis_fn.o postgis_inout.o postgis_proj.o postgis_chip.o postgis_transform.o postgis_gist_$(GIST_SUPPORT).o $(GIST_ESTIMATE)
+OBJS=postgis_debug.o postgis_ops.o postgis_fn.o postgis_inout.o postgis_proj.o postgis_chip.o postgis_transform.o postgis_geos.o postgis_geos_wrapper.o postgis_gist_$(GIST_SUPPORT).o $(GIST_ESTIMATE)
 
 #---------------------------------------------------------------
 # Add libraries that libpq depends (or might depend) on into the
 # shared library link.  (The order in which you list them here doesn't
 # matter.)
-SHLIB_LINK=$(filter -L%, $(LDFLAGS)) 
-ifeq ($(USE_PROJ),1)
-       SHLIB_LINK=$(filter -L%, $(LDFLAGS)) $(BE_DLLLIBS) -L$(PROJ_DIR)/lib -lproj
-else
-       SHLIB_LINK=$(filter -L%, $(LDFLAGS)) $(BE_DLLLIBS)
-endif
+
+
+SHLIB_LINK=$(filter -L%, $(LDFLAGS)) $(BE_DLLLIBS) -lstdc++ -L$(PROJ_DIR)/lib -lproj -L$(GEOS_DIR)/lib -lgeos
+#ifeq ($(USE_PROJ),1)
+#      SHLIB_LNK=$(SHLIB_LNK) -L$(PROJ_DIR)/lib -lproj
+#endif
+#ifeq ($(USE_GEOS),1)
+#      SHLIB_LNK=$(SHLIB_LNK) -L$(GEOS_DIR)/lib -lgeos
+#endif
+
 
 #---------------------------------------------------------------
 # Makefile targets
 
+include $(top_srcdir)/src/Makefile.shlib
+
+postgis_geos_wrapper.o: postgis_geos_wrapper.cpp
+
 all: all-lib $(NAME).sql $(NAME).sql $(NAME)_undef.sql loaderdumper
 
 loaderdumper:
        $(MAKE) -C loader
 
 # Shared library stuff
-include $(top_srcdir)/src/Makefile.shlib
 
 $(NAME).sql: $(NAME)_sql_common.sql.in $(NAME)_sql_$(USE_VERSION)_end.sql.in $(NAME)_sql_$(USE_VERSION)_start.sql.in 
        cat $(NAME)_sql_$(USE_VERSION)_start.sql.in $(NAME)_sql_common.sql.in $(NAME)_sql_$(USE_VERSION)_end.sql.in | sed -e 's:@MODULE_FILENAME@:$(LPATH)/$(shlib):g;s:@POSTGIS_VERSION@:$(SO_MAJOR_VERSION).$(SO_MINOR_VERSION):g'  > $@ 
diff --git a/postgis_geos.c b/postgis_geos.c
new file mode 100644 (file)
index 0000000..5db4997
--- /dev/null
@@ -0,0 +1,615 @@
+/*
+drop function geos_thing(geometry,geometry);
+CREATE FUNCTION relate(geometry,geometry)
+   RETURNS text
+   AS '/data1/Refractions/Projects/PostGIS/work_dave/postgis/libpostgis.so.0.7','relate_full'
+   LANGUAGE 'c' with (isstrict);
+
+CREATE FUNCTION relate(geometry,geometry,text)
+   RETURNS boolean
+   AS '/data1/Refractions/Projects/PostGIS/work_dave/postgis/libpostgis.so.0.7','relate_pattern'
+   LANGUAGE 'c' with (isstrict);
+
+CREATE FUNCTION disjoint(geometry,geometry)
+   RETURNS boolean
+   AS '/data1/Refractions/Projects/PostGIS/work_dave/postgis/libpostgis.so.0.7'
+   LANGUAGE 'c' with (isstrict);
+
+CREATE FUNCTION touches(geometry,geometry)
+   RETURNS boolean
+   AS '/data1/Refractions/Projects/PostGIS/work_dave/postgis/libpostgis.so.0.7'
+   LANGUAGE 'c' with (isstrict);
+
+CREATE FUNCTION intersects(geometry,geometry)
+   RETURNS boolean
+   AS '/data1/Refractions/Projects/PostGIS/work_dave/postgis/libpostgis.so.0.7'
+   LANGUAGE 'c' with (isstrict);
+
+CREATE FUNCTION crosses(geometry,geometry)
+   RETURNS boolean
+   AS '/data1/Refractions/Projects/PostGIS/work_dave/postgis/libpostgis.so.0.7'
+   LANGUAGE 'c' with (isstrict);
+
+CREATE FUNCTION within(geometry,geometry)
+   RETURNS boolean
+   AS '/data1/Refractions/Projects/PostGIS/work_dave/postgis/libpostgis.so.0.7'
+   LANGUAGE 'c' with (isstrict);
+
+CREATE FUNCTION contains(geometry,geometry)
+   RETURNS boolean
+   AS '/data1/Refractions/Projects/PostGIS/work_dave/postgis/libpostgis.so.0.7'
+   LANGUAGE 'c' with (isstrict);
+
+CREATE FUNCTION overlaps(geometry,geometry)
+   RETURNS boolean
+   AS '/data1/Refractions/Projects/PostGIS/work_dave/postgis/libpostgis.so.0.7'
+   LANGUAGE 'c' with (isstrict);
+
+  CREATE FUNCTION isvalid(geometry)
+   RETURNS boolean
+   AS '/data1/Refractions/Projects/PostGIS/work_dave/postgis/libpostgis.so.0.7'
+   LANGUAGE 'c' with (isstrict);
+*/
+
+//--------------------------------------------------------------------------
+//  
+#ifdef USE_GEOS
+
+#include "postgres.h"
+
+
+#include <math.h>
+#include <float.h>
+#include <string.h>
+#include <stdio.h>
+#include <errno.h>
+
+#include "fmgr.h"
+
+#include "postgis.h"
+#include "utils/elog.h"
+
+#include "access/gist.h"
+#include "access/itup.h"
+#include "access/rtree.h"
+
+#include "utils/builtins.h"
+
+
+
+typedef  struct Geometry Geometry;
+
+
+extern const char * createGEOSPoint(POINT3D *pt);
+extern void initGEOS(int maxalign);
+extern char *GEOSrelate(Geometry *g1, Geometry*g2);
+extern bool GEOSrelatePattern(Geometry *g1, Geometry*g2,char *pat);
+extern  bool GEOSrelateDisjoint(Geometry *g1, Geometry*g2);
+extern  bool GEOSrelateTouches(Geometry *g1, Geometry*g2);
+extern  bool GEOSrelateIntersects(Geometry *g1, Geometry*g2);
+extern  bool GEOSrelateCrosses(Geometry *g1, Geometry*g2);
+extern  bool GEOSrelateWithin(Geometry *g1, Geometry*g2);
+extern  bool GEOSrelateContains(Geometry *g1, Geometry*g2);
+extern  bool GEOSrelateOverlaps(Geometry *g1, Geometry*g2);
+
+extern char *GEOSasText(Geometry *g1);
+
+
+extern void GEOSdeleteChar(char *a);
+extern void GEOSdeleteGeometry(Geometry *a);
+
+extern  Geometry *PostGIS2GEOS_point(POINT3D *point,int SRID, bool is3d);
+extern  Geometry *PostGIS2GEOS_linestring(LINE3D *line,int SRID, bool is3d);
+extern  Geometry *PostGIS2GEOS_polygon(POLYGON3D *polygon,int SRID, bool is3d);
+extern  Geometry *PostGIS2GEOS_multipolygon(POLYGON3D **polygons,int npolys, int SRID, bool is3d);
+extern Geometry *PostGIS2GEOS_multilinestring(LINE3D **lines,int nlines, int SRID, bool is3d);
+extern Geometry *PostGIS2GEOS_multipoint(POINT3D **points,int npoints, int SRID, bool is3d);
+extern Geometry *PostGIS2GEOS_box3d(BOX3D *box, int SRID);
+extern Geometry *PostGIS2GEOS_collection(Geometry **geoms, int ngeoms,int SRID, bool is3d);
+
+extern bool GEOSisvalid(Geometry *g1);
+
+
+Datum relate_full(PG_FUNCTION_ARGS);
+Datum relate_pattern(PG_FUNCTION_ARGS);
+Datum disjoint(PG_FUNCTION_ARGS);
+Datum touches(PG_FUNCTION_ARGS);
+Datum intersects(PG_FUNCTION_ARGS);
+Datum crosses(PG_FUNCTION_ARGS);
+Datum within(PG_FUNCTION_ARGS);
+Datum contains(PG_FUNCTION_ARGS);
+Datum overlaps(PG_FUNCTION_ARGS);
+
+
+Datum isvalid(PG_FUNCTION_ARGS);
+
+Geometry *POSTGIS2GEOS(GEOMETRY *g);
+void errorIfGeometryCollection(GEOMETRY *g1, GEOMETRY *g2);
+
+// return a GEOS Geometry from a POSTGIS GEOMETRY
+
+void errorIfGeometryCollection(GEOMETRY *g1, GEOMETRY *g2)
+{
+       if (  (g1->type == COLLECTIONTYPE) || (g2->type == COLLECTIONTYPE) )
+               elog(ERROR,"Relate Operation called with a GEOMETRYCOLLECTION type.  This is unsupported");
+}
+
+PG_FUNCTION_INFO_V1(isvalid);
+Datum isvalid(PG_FUNCTION_ARGS) 
+{
+               GEOMETRY                *geom1 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+               bool result;
+               Geometry *g1;
+
+               initGEOS(MAXIMUM_ALIGNOF);
+
+               g1 =    POSTGIS2GEOS(geom1 );
+               result = GEOSisvalid(g1);
+
+       GEOSdeleteGeometry(g1);
+
+       PG_RETURN_BOOL(result);
+}
+
+PG_FUNCTION_INFO_V1(overlaps);
+Datum overlaps(PG_FUNCTION_ARGS)
+{
+       GEOMETRY                *geom1 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+       GEOMETRY                *geom2 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+
+
+       Geometry *g1,*g2;
+       bool result;
+
+       errorIfGeometryCollection(geom1,geom2);
+       initGEOS(MAXIMUM_ALIGNOF);
+
+       g1 =    POSTGIS2GEOS(geom1 );
+       g2 =    POSTGIS2GEOS(geom2 );
+
+
+       result = GEOSrelateOverlaps(g1,g2);
+
+       GEOSdeleteGeometry(g1);
+       GEOSdeleteGeometry(g2);
+
+       PG_RETURN_BOOL(result);
+}
+
+
+
+PG_FUNCTION_INFO_V1(contains);
+Datum contains(PG_FUNCTION_ARGS)
+{
+       GEOMETRY                *geom1 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+       GEOMETRY                *geom2 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+
+
+       Geometry *g1,*g2;
+       bool result;
+
+       errorIfGeometryCollection(geom1,geom2);
+       initGEOS(MAXIMUM_ALIGNOF);
+
+       g1 =    POSTGIS2GEOS(geom1 );
+       g2 =    POSTGIS2GEOS(geom2 );
+
+
+       result = GEOSrelateContains(g1,g2);
+
+       GEOSdeleteGeometry(g1);
+       GEOSdeleteGeometry(g2);
+
+       PG_RETURN_BOOL(result);
+}
+
+
+PG_FUNCTION_INFO_V1(within);
+Datum within(PG_FUNCTION_ARGS)
+{
+       GEOMETRY                *geom1 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+       GEOMETRY                *geom2 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+
+
+       Geometry *g1,*g2;
+       bool result;
+
+       errorIfGeometryCollection(geom1,geom2);
+       initGEOS(MAXIMUM_ALIGNOF);
+
+       g1 =    POSTGIS2GEOS(geom1 );
+       g2 =    POSTGIS2GEOS(geom2 );
+
+
+       result = GEOSrelateWithin(g1,g2);
+
+       GEOSdeleteGeometry(g1);
+       GEOSdeleteGeometry(g2);
+
+       PG_RETURN_BOOL(result);
+}
+
+
+
+PG_FUNCTION_INFO_V1(crosses);
+Datum crosses(PG_FUNCTION_ARGS)
+{
+       GEOMETRY                *geom1 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+       GEOMETRY                *geom2 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+
+
+       Geometry *g1,*g2;
+       bool result;
+
+       errorIfGeometryCollection(geom1,geom2);
+       initGEOS(MAXIMUM_ALIGNOF);
+
+       g1 =    POSTGIS2GEOS(geom1 );
+       g2 =    POSTGIS2GEOS(geom2 );
+
+
+       result = GEOSrelateCrosses(g1,g2);
+
+       GEOSdeleteGeometry(g1);
+       GEOSdeleteGeometry(g2);
+
+       PG_RETURN_BOOL(result);
+}
+
+
+
+PG_FUNCTION_INFO_V1(intersects);
+Datum intersects(PG_FUNCTION_ARGS)
+{
+       GEOMETRY                *geom1 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+       GEOMETRY                *geom2 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+
+
+       Geometry *g1,*g2;
+       bool result;
+
+       errorIfGeometryCollection(geom1,geom2);
+       initGEOS(MAXIMUM_ALIGNOF);
+
+       g1 =    POSTGIS2GEOS(geom1 );
+       g2 =    POSTGIS2GEOS(geom2 );
+
+
+       result = GEOSrelateIntersects(g1,g2);
+
+       GEOSdeleteGeometry(g1);
+       GEOSdeleteGeometry(g2);
+
+       PG_RETURN_BOOL(result);
+}
+
+
+PG_FUNCTION_INFO_V1(touches);
+Datum touches(PG_FUNCTION_ARGS)
+{
+       GEOMETRY                *geom1 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+       GEOMETRY                *geom2 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+
+
+       Geometry *g1,*g2;
+       bool result;
+
+       errorIfGeometryCollection(geom1,geom2);
+       initGEOS(MAXIMUM_ALIGNOF);
+
+       g1 =    POSTGIS2GEOS(geom1 );
+       g2 =    POSTGIS2GEOS(geom2 );
+
+
+       result = GEOSrelateTouches(g1,g2);
+
+       GEOSdeleteGeometry(g1);
+       GEOSdeleteGeometry(g2);
+
+       PG_RETURN_BOOL(result);
+}
+
+
+PG_FUNCTION_INFO_V1(disjoint);
+Datum disjoint(PG_FUNCTION_ARGS)
+{
+       GEOMETRY                *geom1 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+       GEOMETRY                *geom2 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+
+
+       Geometry *g1,*g2;
+       bool result;
+
+       errorIfGeometryCollection(geom1,geom2);
+       initGEOS(MAXIMUM_ALIGNOF);
+
+       g1 =    POSTGIS2GEOS(geom1 );
+       g2 =    POSTGIS2GEOS(geom2 );
+
+
+       result = GEOSrelateDisjoint(g1,g2);
+
+       GEOSdeleteGeometry(g1);
+       GEOSdeleteGeometry(g2);
+
+       PG_RETURN_BOOL(result);
+}
+
+
+PG_FUNCTION_INFO_V1(relate_pattern);
+Datum relate_pattern(PG_FUNCTION_ARGS)
+{
+       GEOMETRY                *geom1 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+       GEOMETRY                *geom2 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+       char *patt;
+       bool result;
+
+       Geometry *g1,*g2;
+
+
+       errorIfGeometryCollection(geom1,geom2);
+       initGEOS(MAXIMUM_ALIGNOF);
+
+       g1 =    POSTGIS2GEOS(geom1 );
+       g2 =    POSTGIS2GEOS(geom2 );
+
+
+       patt =  DatumGetCString(DirectFunctionCall1(textout,
+                        PointerGetDatum(PG_GETARG_DATUM(2))));
+
+       result = GEOSrelatePattern(g1,g2,patt);
+
+
+       GEOSdeleteGeometry(g1);
+       GEOSdeleteGeometry(g2);
+
+       pfree(patt);
+       PG_RETURN_BOOL(result);
+
+
+}
+
+
+
+PG_FUNCTION_INFO_V1(relate_full);
+Datum relate_full(PG_FUNCTION_ARGS)
+{
+       GEOMETRY                *geom1 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+       GEOMETRY                *geom2 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+
+       Geometry *g1,*g2;
+       char    *relate_str;
+       int len;
+       char *result;
+
+       errorIfGeometryCollection(geom1,geom2);
+
+//elog(NOTICE,"in relate_full: \n");
+       initGEOS(MAXIMUM_ALIGNOF);
+
+//elog(NOTICE,"making GOES geometries: \n");
+       g1 =    POSTGIS2GEOS(geom1 );
+       g2 =    POSTGIS2GEOS(geom2 );
+//     elog(NOTICE,"geometries made, here's what they are: \n");
+//g1 = createGEOSFromText("POINT(1 1)");
+//g2 = createGEOSFromText("POINT(1 1)");
+//elog(NOTICE,GEOSasText(g1));
+//elog(NOTICE,GEOSasText(g2));
+
+//elog(NOTICE,"calling relate: \n");
+       relate_str = GEOSrelate(g1, g2);
+
+//elog(NOTICE,"relate finished \n");
+
+
+       len = strlen(relate_str) + 4;
+
+       result= palloc(len);
+       *((int *) result) = len;
+
+       memcpy(result +4, relate_str, len-4);
+
+       free(relate_str);
+//     GEOSdeleteGeometry(g1);
+//     GEOSdeleteGeometry(g2);
+
+       PG_RETURN_POINTER(result);
+}
+
+
+//BBOXONLYTYPE -> returns as a 2d polygon
+Geometry *POSTGIS2GEOS(GEOMETRY *g)
+{
+       POINT3D *pt;
+       LINE3D *line;
+       POLYGON3D *poly;
+       POLYGON3D **polys;
+       LINE3D **lines;
+       POINT3D **points;
+       Geometry **geoms;
+       Geometry *geos;
+       char     *obj;
+       int      obj_type;
+       int t;
+
+       int32 *offsets1 = (int32 *) ( ((char *) &(g->objType[0] ))+ sizeof(int32) * g->nobjs ) ;
+
+       switch(g->type)
+       {
+               case POINTTYPE: 
+                                                       pt = (POINT3D*) ((char *) g +offsets1[0]) ;
+                                                       return PostGIS2GEOS_point(pt,g->SRID,g->is3d);
+                                                       break;
+               case LINETYPE:          
+                                                       line = (LINE3D*) ((char *) g +offsets1[0]) ;
+                                                       return PostGIS2GEOS_linestring(line,g->SRID,g->is3d);
+                                                       break;
+               case POLYGONTYPE:   
+                                                       poly = (POLYGON3D*) ((char *) g +offsets1[0]) ;
+                                                       return PostGIS2GEOS_polygon(poly,g->SRID,g->is3d);
+                                                       break;
+               case MULTIPOLYGONTYPE:
+                                                                       //make an array of POLYGON3Ds
+                                                       polys = (POLYGON3D**) palloc(sizeof (POLYGON3D*) * g->nobjs);
+                                                       for (t=0;t<g->nobjs;t++)
+                                                       {
+                                                               polys[t] =      (POLYGON3D*) ((char *) g +offsets1[t]) ;
+                                                       }
+                                                       geos= PostGIS2GEOS_multipolygon(polys, g->nobjs, g->SRID,g->is3d);
+                                                       pfree(polys);
+                                                       return geos;
+                                                       break;
+               case MULTILINETYPE:
+                                                               //make an array of POLYGON3Ds
+                                                       lines = (LINE3D**) palloc(sizeof (LINE3D*) * g->nobjs);
+                                                       for (t=0;t<g->nobjs;t++)
+                                                       {
+                                                               lines[t] =      (LINE3D*) ((char *) g +offsets1[t]) ;
+                                                       }
+                                                       geos= PostGIS2GEOS_multilinestring(lines, g->nobjs, g->SRID,g->is3d);
+                                                       pfree(lines);
+                                                       return geos;
+                                                       break;
+               case MULTIPOINTTYPE:
+                                                               //make an array of POINT3Ds
+                                                       points = (POINT3D**) palloc(sizeof (POINT3D*) * g->nobjs);
+                                                       for (t=0;t<g->nobjs;t++)
+                                                       {
+                                                               points[t] =     (POINT3D*) ((char *) g +offsets1[t]) ;
+                                                       }
+                                                       geos= PostGIS2GEOS_multipoint(points, g->nobjs,g->SRID,g->is3d);
+                                                       pfree(points);
+                                                       return geos;
+                                                       break;
+               case BBOXONLYTYPE:
+                                                        return PostGIS2GEOS_box3d(&g->bvol, g->SRID);
+                                                        break;
+               case COLLECTIONTYPE:
+                                                               //make an array of GEOS Geometrys
+                                                       geoms = (Geometry**) palloc(sizeof (Geometry*) * g->nobjs);
+                                                       for (t=0;t<g->nobjs;t++)
+                                                       {
+                                                               obj = ((char *) g +offsets1[t]);
+                                                               obj_type =  g->objType[t];
+                                                               switch (obj_type)
+                                                               {
+                                                                       case POINTTYPE:
+                                                                                                       pt = (POINT3D*) obj ;
+                                                                                                       geoms[t] = PostGIS2GEOS_point(pt,g->SRID,g->is3d);
+                                                                                                       break;
+                                                                       case LINETYPE:
+                                                                                                       line = (LINE3D*) obj ;
+                                                                                                       geoms[t] = PostGIS2GEOS_linestring(line,g->SRID,g->is3d);
+                                                                                                       break;
+                                                                       case POLYGONTYPE:
+                                                                                                       poly = (POLYGON3D*) obj ;
+                                                                                                       geoms[t] = PostGIS2GEOS_polygon(poly,g->SRID,g->is3d);
+                                                                                                       break;
+                                                               }
+                                                       }
+                                                       geos= PostGIS2GEOS_collection(geoms,g->nobjs,g->SRID,g->is3d);
+                                                       pfree(geoms);
+                                                       return geos; 
+                                                       break;
+
+       }
+       return NULL;
+}
+
+
+
+
+//----------------------------------------------------------------------------
+// NULL implementation here
+// ---------------------------------------------------------------------------
+#else
+
+
+#include "postgres.h"
+
+
+#include <math.h>
+#include <float.h>
+#include <string.h>
+#include <stdio.h>
+#include <errno.h>
+
+#include "fmgr.h"
+
+
+Datum relate_full(PG_FUNCTION_ARGS);
+Datum relate_pattern(PG_FUNCTION_ARGS);
+Datum disjoint(PG_FUNCTION_ARGS);
+Datum touches(PG_FUNCTION_ARGS);
+Datum intersects(PG_FUNCTION_ARGS);
+Datum crosses(PG_FUNCTION_ARGS);
+Datum within(PG_FUNCTION_ARGS);
+Datum contains(PG_FUNCTION_ARGS);
+Datum overlaps(PG_FUNCTION_ARGS);
+Datum isvalid(PG_FUNCTION_ARGS);
+
+PG_FUNCTION_INFO_V1(relate_full);
+Datum relate_full(PG_FUNCTION_ARGS)
+{
+       elog(ERROR,"relate_full:: operation not implemented - compile PostGIS with GEOS support");
+       PG_RETURN_NULL(); // never get here
+}
+PG_FUNCTION_INFO_V1(relate_pattern);
+Datum relate_pattern(PG_FUNCTION_ARGS)
+{
+       elog(ERROR,"relate_pattern:: operation not implemented - compile PostGIS with GEOS support");
+       PG_RETURN_NULL(); // never get here
+}
+PG_FUNCTION_INFO_V1(disjoint);
+Datum disjoint(PG_FUNCTION_ARGS)
+{
+       elog(ERROR,"disjoint:: operation not implemented - compile PostGIS with GEOS support");
+       PG_RETURN_NULL(); // never get here
+}
+PG_FUNCTION_INFO_V1(intersects);
+Datum intersects(PG_FUNCTION_ARGS)
+{
+       elog(ERROR,"intersects:: operation not implemented - compile PostGIS with GEOS support");
+       PG_RETURN_NULL(); // never get here
+}
+PG_FUNCTION_INFO_V1(touches);
+Datum touches(PG_FUNCTION_ARGS)
+{
+       elog(ERROR,"touches:: operation not implemented - compile PostGIS with GEOS support");
+       PG_RETURN_NULL(); // never get here
+}
+PG_FUNCTION_INFO_V1(crosses);
+Datum crosses(PG_FUNCTION_ARGS)
+{
+       elog(ERROR,"crosses:: operation not implemented - compile PostGIS with GEOS support");
+       PG_RETURN_NULL(); // never get here
+}
+PG_FUNCTION_INFO_V1(within);
+Datum within(PG_FUNCTION_ARGS)
+{
+       elog(ERROR,"within:: operation not implemented - compile PostGIS with GEOS support");
+       PG_RETURN_NULL(); // never get here
+}
+PG_FUNCTION_INFO_V1(contains);
+Datum contains(PG_FUNCTION_ARGS)
+{
+       elog(ERROR,"contains:: operation not implemented - compile PostGIS with GEOS support");
+       PG_RETURN_NULL(); // never get here
+}
+PG_FUNCTION_INFO_V1(overlaps);
+Datum overlaps(PG_FUNCTION_ARGS)
+{
+       elog(ERROR,"overlaps:: operation not implemented - compile PostGIS with GEOS support");
+       PG_RETURN_NULL(); // never get here
+}
+PG_FUNCTION_INFO_V1(isvalid);
+Datum isvalid(PG_FUNCTION_ARGS)
+{
+       elog(ERROR,"isvalid:: operation not implemented - compile PostGIS with GEOS support");
+       PG_RETURN_NULL(); // never get here
+}
+
+#endif
+
diff --git a/postgis_geos_wrapper.cpp b/postgis_geos_wrapper.cpp
new file mode 100644 (file)
index 0000000..747a882
--- /dev/null
@@ -0,0 +1,401 @@
+//  g++ postgis_GEOSwrapper.cpp -c -I/usr/local/include  -I/usr/local/include/geos -I/usr/local/src/postgresql-7.2.3//src/include
+
+
+#include <stdio.h>
+
+#include <string>
+#include <iostream>
+#include <fstream>
+
+//#include "geom.h"
+//#include "util.h"
+//#include "graph.h"
+#include "io.h"
+//#include "opRelate.h"
+
+
+
+//WARNING THIS *MUST* BE SET CORRECTLY.
+int MAXIMUM_ALIGNOF = -999999;    // to be set during initialization - this will be either 4 (intel) or 8 (sparc)
+
+//for getting things to align properly  double are on 8byte align on solaris machines, and 4bytes on intel
+
+#define TYPEALIGN(ALIGNVAL,LEN) (((long)(LEN) + (ALIGNVAL-1)) & ~(ALIGNVAL-1))
+#define MAXALIGN(LEN)           TYPEALIGN(MAXIMUM_ALIGNOF, (LEN))
+
+typedef  int int32;
+
+typedef struct
+{
+       double          x,y,z;  //for lat/long   x=long, y=lat
+} POINT3D;
+
+typedef struct
+{
+       POINT3D         LLB,URT; /* corner POINT3Ds on long diagonal */
+} BOX3D;
+
+typedef struct
+{
+       int32   npoints; // how many points in the line
+       int32   junk;      // double-word alignment
+       POINT3D         points[1]; // array of actual points
+} LINE3D;
+
+
+typedef struct
+{
+       int32   nrings;  // how many rings in this polygon
+       int32           npoints[1]; //how many points in each ring
+       /* could be 4 byes of filler here to make sure points[] is
+         double-word aligned*/
+       POINT3D         points[1]; // array of actual points
+} POLYGON3D;
+
+typedef struct
+{
+       int32           size;           // postgres variable-length type requirement
+       int32           SRID;           // spatial reference system id
+       double  offsetX;        // for precision grid (future improvement)
+       double  offsetY;        // for precision grid (future improvement)
+       double  scale;  // for precision grid (future improvement)
+       int32           type;           // this type of geometry
+       bool            is3d;           // true if the points are 3d (only for output)
+       BOX3D           bvol;           // bounding volume of all the geo objects
+       int32           nobjs;  // how many sub-objects in this object
+       int32           objType[1];     // type of object
+       int32           objOffset[1];// offset (in bytes) into this structure where
+                                        // the object is located
+       char            objData[1];  // store for actual objects
+
+} GEOMETRY;
+
+
+//###########################################################
+
+extern "C" char *GEOSrelate(Geometry *g1, Geometry*g2);
+extern "C" void initGEOS(int maxalign);
+
+
+extern "C" void GEOSdeleteChar(char *a);
+extern "C" void GEOSdeleteGeometry(Geometry *a);
+extern "C" bool GEOSrelatePattern(Geometry *g1, Geometry*g2,char *pat);
+extern "C" bool GEOSrelateDisjoint(Geometry *g1, Geometry*g2);
+extern "C" bool GEOSrelateTouches(Geometry *g1, Geometry*g2);
+extern "C" bool GEOSrelateIntersects(Geometry *g1, Geometry*g2);
+extern "C" bool GEOSrelateCrosses(Geometry *g1, Geometry*g2);
+extern "C" bool GEOSrelateWithin(Geometry *g1, Geometry*g2);
+extern "C" bool GEOSrelateContains(Geometry *g1, Geometry*g2);
+extern "C" bool GEOSrelateOverlaps(Geometry *g1, Geometry*g2);
+
+
+extern "C" Geometry *PostGIS2GEOS_point(POINT3D *point,int SRID, bool is3d);
+extern "C" Geometry *PostGIS2GEOS_linestring(LINE3D *line,int SRID, bool is3d);
+extern "C" Geometry *PostGIS2GEOS_polygon(POLYGON3D *polygon,int SRID, bool is3d);
+extern "C" Geometry *PostGIS2GEOS_multipolygon(POLYGON3D **polygons,int npolys, int SRID, bool is3d);
+extern "C" Geometry *PostGIS2GEOS_multilinestring(LINE3D **lines,int nlines, int SRID, bool is3d);
+extern "C" Geometry *PostGIS2GEOS_multipoint(POINT3D **points,int npoints, int SRID, bool is3d);
+
+extern "C" Geometry *PostGIS2GEOS_box3d(BOX3D *box, int SRID);
+extern "C" Geometry *PostGIS2GEOS_collection(Geometry **geoms, int ngeoms,int SRID, bool is3d);
+
+extern "C" bool GEOSisvalid(Geometry *g1);
+
+
+//###########################################################
+
+GeometryFactory *geomFactory = NULL;
+
+void initGEOS (int maxalign)
+{
+       if (geomFactory == NULL)
+       {
+               geomFactory = new GeometryFactory( new PrecisionModel(), -1); // NOTE: SRID will have to be changed after geometry creation
+               MAXIMUM_ALIGNOF = maxalign;
+       }
+}
+
+               //note: you lose the 3d from this!
+Geometry *PostGIS2GEOS_box3d(BOX3D *box, int SRID)
+{
+       Geometry *g;
+       Envelope *envelope = new Envelope(box->LLB.x,box->URT.x ,box->LLB.y,box->URT.y);
+       g = geomFactory->toGeometry(envelope,geomFactory->getPrecisionModel(), SRID);
+       delete envelope;
+       return g;
+}
+
+Geometry *PostGIS2GEOS_collection(Geometry **geoms, int ngeoms,int SRID, bool is3d)
+{
+       Geometry *g;
+       int t;
+       vector<Geometry *> *subGeos=new vector<Geometry *>;
+
+       for (t =0; t< ngeoms; t++)
+       {
+               subGeos->push_back(geoms[t]);   
+       }
+       g = geomFactory->buildGeometry(subGeos);
+       g->setSRID(SRID);
+       return g;
+
+}
+
+Geometry *PostGIS2GEOS_point(POINT3D *point,int SRID, bool is3d)
+{
+       Coordinate *c;
+
+       if (is3d)
+               c = new Coordinate(point->x, point->y); 
+       else
+               c = new Coordinate(point->x, point->y, point->z);
+       Geometry *g = geomFactory->createPoint(*c);
+       g->setSRID(SRID);
+       return g;
+}
+
+
+Geometry *PostGIS2GEOS_linestring(LINE3D *line,int SRID, bool is3d)
+{
+       int t;
+       Coordinate c;
+
+       //build coordinatelist & pre-allocate space
+       BasicCoordinateList  *coords = new BasicCoordinateList(line->npoints);
+       if (is3d)
+       {
+               for (t=0;t<line->npoints;t++)
+               {
+                       c.x = line->points[t].x;
+                       c.y = line->points[t].y;
+                       c.z = line->points[t].z;
+                       coords->setAt( c ,t);
+               }
+
+       }
+       else  //make 2d points
+       {
+               for (t=0;t<line->npoints;t++)
+               {
+                       c.x = line->points[t].x;
+                       c.y = line->points[t].y;
+                       c.z = DoubleNotANumber;
+                       coords->setAt( c ,t);
+               }
+
+       }
+       Geometry *g = geomFactory->createLineString(coords);
+       g->setSRID(SRID);
+       return g;
+}
+
+       //polygons is an array of pointers to polygons
+Geometry *PostGIS2GEOS_multipolygon(POLYGON3D **polygons,int npolys, int SRID, bool is3d)
+{
+       int t;
+       vector<Geometry *> *subPolys=new vector<Geometry *>;
+       Geometry *g;
+
+       for (t =0; t< npolys; t++)
+       {
+               subPolys->push_back(PostGIS2GEOS_polygon(polygons[t], SRID,is3d ));     
+       }
+       g = geomFactory->createMultiPolygon(subPolys);
+       g->setSRID(SRID);
+       return g;
+}
+
+       //lines is an array of pointers to line3d
+Geometry *PostGIS2GEOS_multilinestring(LINE3D **lines,int nlines, int SRID, bool is3d)
+{
+       int t;
+       vector<Geometry *> *subLines =new vector<Geometry *>;
+       Geometry *g;
+
+       for (t =0; t< nlines; t++)
+       {
+               subLines->push_back(PostGIS2GEOS_linestring(lines[t], SRID,is3d ));     
+       }
+       g = geomFactory->createMultiLineString(subLines);
+       g->setSRID(SRID);
+       return g;
+}
+
+Geometry *PostGIS2GEOS_multipoint(POINT3D **points,int npoints, int SRID, bool is3d)
+{
+       int t;
+       vector<Geometry *> *subPoints =new vector<Geometry *>;
+       Geometry *g;
+
+       for (t =0; t< npoints; t++)
+       {
+               subPoints->push_back(PostGIS2GEOS_point(points[t], SRID,is3d ));        
+       }
+       g = geomFactory->createMultiPoint(subPoints);
+       g->setSRID(SRID);
+       return g;
+
+}
+
+
+Geometry *PostGIS2GEOS_polygon(POLYGON3D *polygon,int SRID, bool is3d)
+{
+               POINT3D *pts;
+               Coordinate c;
+               int  ring,t;
+               Geometry *g;
+               LinearRing *outerRing;
+               LinearRing *innerRing;
+               BasicCoordinateList *cl;
+               int pointOffset =0; // first point that we're looking at.  a POLYGON3D has all its points smooshed together
+               vector<Geometry *> *innerRings=new vector<Geometry *>;
+
+
+               pts = (POINT3D *) ( (char *)&(polygon->npoints[polygon->nrings] )  );
+               pts = (POINT3D *) MAXALIGN(pts);
+       
+                       // make outerRing
+                       cl = new BasicCoordinateList(polygon->npoints[0]);
+                       if (is3d)
+                       {
+                               for(t=0;t<polygon->npoints[0];t++)
+                               {
+                                               c.x = pts[t].x;
+                                               c.y = pts[t].y;
+                                               c.z = pts[t].z;
+                                               cl->setAt( c ,t);
+                               }
+                       }
+                       else
+                       {
+                               for(t=0;t<polygon->npoints[0];t++)
+                               {
+                                               c.x = pts[t].x;
+                                               c.y = pts[t].y;
+                                               c.z = DoubleNotANumber;
+                                               cl->setAt( c ,t);
+                               }
+                       }
+                       outerRing = geomFactory->createLinearRing(cl);
+                       outerRing->setSRID(SRID);
+                       pointOffset = polygon->npoints[0];
+
+               for(ring =1; ring< polygon->nrings; ring++)
+               {
+                       cl = new BasicCoordinateList(polygon->npoints[ring]);
+                       if (is3d)
+                       {
+                               for(t=0;t<polygon->npoints[ring];t++)
+                               {
+                                               c.x = pts[t+pointOffset].x;
+                                               c.y = pts[t+pointOffset].y;
+                                               c.z = pts[t+pointOffset].z;
+                                               cl->setAt( c ,t);
+                               }
+                       }
+                       else
+                       {
+                               for(t=0;t<polygon->npoints[ring];t++)
+                               {
+                                               c.x = pts[t+pointOffset].x;
+                                               c.y = pts[t+pointOffset].y;
+                                               c.z = DoubleNotANumber;
+                                               cl->setAt( c ,t);
+                               }
+                       }
+                       innerRing = geomFactory->createLinearRing(cl);
+                       innerRing->setSRID(SRID);
+                       innerRings->push_back(innerRing);
+                       pointOffset += polygon->npoints[ring];
+               }
+
+       g = geomFactory->createPolygon(outerRing, innerRings);
+       g->setSRID(SRID);
+       return g;
+}
+
+
+
+bool GEOSrelateDisjoint(Geometry *g1, Geometry*g2)
+{
+       return g1->disjoint(g2);
+}
+
+bool GEOSrelateTouches(Geometry *g1, Geometry*g2)
+{
+       return g1->touches(g2);
+}
+
+bool GEOSrelateIntersects(Geometry *g1, Geometry*g2)
+{
+       return g1->intersects(g2);
+}
+
+bool GEOSrelateCrosses(Geometry *g1, Geometry*g2)
+{
+       return g1->crosses(g2);
+}
+
+bool GEOSrelateWithin(Geometry *g1, Geometry*g2)
+{
+       return g1->within(g2);
+}
+
+bool GEOSrelateContains(Geometry *g1, Geometry*g2)
+{
+       return g1->contains(g2);
+}
+
+bool GEOSrelateOverlaps(Geometry *g1, Geometry*g2)
+{
+       return g1->overlaps(g2);
+}
+
+//BUG:: this leaks memory, but delete kills the PrecisionModel for ALL the geometries
+void GEOSdeleteGeometry(Geometry *a)
+{
+//     delete a;
+}
+
+void GEOSdeleteChar(char *a)
+{
+       free(a);
+}
+
+
+bool GEOSrelatePattern(Geometry *g1, Geometry*g2,char *pat)
+{
+       string s = pat;
+       return g1->relate(g2,pat);
+}
+
+
+char *GEOSrelate(Geometry *g1, Geometry*g2)
+{
+       IntersectionMatrix *im = g1->relate(g2);
+       string s;
+       char *result;
+
+       s= im->toString();
+       result = (char*) malloc( s.length() + 1);
+       strcpy(result, s.c_str() );
+
+
+       return result;
+}
+
+bool GEOSisvalid(Geometry *g1)
+{
+       return g1->isValid();
+}
+
+
+
+
+char *GEOSasText(Geometry *g1)
+{
+       string s = g1->toString();
+
+       return (char *) s.c_str() ;
+}
+