#3074, add in lwgeom backend for subdivision

git-svn-id: http://svn.osgeo.org/postgis/trunk@13323 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/liblwgeom/cunit/cu_geos.c b/liblwgeom/cunit/cu_geos.c
index 877b7bb51d8c99e720bad5b3d8b8d0934b37a1bd..8986ce9319e830c5dee08fb5b0acefc00842c91a 100644 (file)
--- a/liblwgeom/cunit/cu_geos.c
+++ b/liblwgeom/cunit/cu_geos.c
@@ -18,6 +18,8 @@
 #include "lwgeom_geos.h"
 #include "cu_tester.h"
 
+#include "liblwgeom_internal.h"
+
 static void test_geos_noop(void)
 {
        int i;
@@ -65,6 +67,38 @@ static void test_geos_noop(void)
 }
 
 
+static void test_geos_subdivide(void)
+{
+#if POSTGIS_GEOS_VERSION < 35
+       printf("%d\n", POSTGIS_GEOS_VERSION);
+       return;
+#else
+       char *ewkt = "MULTILINESTRING((0 0, 0 100))";
+       char *out_ewkt;
+       LWGEOM *geom1 = lwgeom_from_wkt(ewkt, LW_PARSER_CHECK_NONE);
+
+       LWGEOM *geom2 = lwgeom_segmentize2d(geom1, 1.0);
+       LWCOLLECTION *geom3 = lwgeom_subdivide(geom2, 80);
+       out_ewkt = lwgeom_to_ewkt((LWGEOM*)geom3);
+       // printf("\n--------\n%s\n--------\n", out_ewkt);
+       CU_ASSERT_EQUAL(2, geom3->ngeoms);
+       lwfree(out_ewkt);
+       lwcollection_free(geom3);
+       lwgeom_free(geom2);
+
+       geom2 = lwgeom_segmentize2d(geom1, 1.0);
+       geom3 = lwgeom_subdivide(geom2, 20);
+       out_ewkt = lwgeom_to_ewkt((LWGEOM*)geom3);
+       // printf("\n--------\n%s\n--------\n", out_ewkt);
+       CU_ASSERT_EQUAL(8, geom3->ngeoms);
+       lwfree(out_ewkt);
+       lwcollection_free(geom3);
+       lwgeom_free(geom2);
+
+       lwgeom_free(geom1);
+#endif
+}
+
 /*
 ** Used by test harness to register the tests in this file.
 */
@@ -73,4 +107,5 @@ void geos_suite_setup(void)
 {
        CU_pSuite suite = CU_add_suite("GEOS", NULL, NULL);
        PG_ADD_TEST(suite, test_geos_noop);
+       PG_ADD_TEST(suite, test_geos_subdivide);
 }

diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in
index cf7c52afe9234898a21060a92430afc7d2468159..44fe32b4764a93079e526d8fd3d4ba5b178cac75 100644 (file)
--- a/liblwgeom/liblwgeom.h.in
+++ b/liblwgeom/liblwgeom.h.in
@@ -2014,6 +2014,7 @@ LWGEOM *lwgeom_difference(const LWGEOM *geom1, const LWGEOM *geom2);
 LWGEOM *lwgeom_symdifference(const LWGEOM* geom1, const LWGEOM* geom2);
 LWGEOM *lwgeom_union(const LWGEOM *geom1, const LWGEOM *geom2);
 LWGEOM *lwgeom_clip_by_rect(const LWGEOM *geom1, double x0, double y0, double x1, double y1);
+LWCOLLECTION *lwgeom_subdivide(const LWGEOM *geom, int maxvertices);
 
 /**
  * Snap vertices and segments of a geometry to another using a given tolerance.

diff --git a/liblwgeom/lwgeom.c b/liblwgeom/lwgeom.c
index 186e8283c6cc88a79c8b788301c19ffb17caa136..e53c96acf676f4cee9d729d99e1cd0588862da53 100644 (file)
--- a/liblwgeom/lwgeom.c
+++ b/liblwgeom/lwgeom.c
@@ -1872,3 +1872,103 @@ lwgeom_npoints_in_rect(const LWGEOM *geom, const GBOX *gbox)
        return 0;
 }
 
+/* Prototype for recursion */
+static int lwgeom_subdivide_recursive(const LWGEOM *geom, int maxvertices, LWCOLLECTION *col, GBOX *clip);
+
+static int
+lwgeom_subdivide_recursive(const LWGEOM *geom, int maxvertices, LWCOLLECTION *col, GBOX *clip)
+{
+
+       /* Always just recurse into collections */
+       if ( lwgeom_is_collection(geom) )
+       {
+               LWCOLLECTION *gcol = (LWCOLLECTION*)geom;
+               int i, n = 0;
+               for ( i = 0; i < gcol->ngeoms; i++ )
+               {
+                       n += lwgeom_subdivide_recursive(gcol->geoms[i], maxvertices, col, clip);
+               }
+               return n;
+       }
+       
+       /* Points and single-point multipoints can be added right into the output result */
+       if ( (geom->type == POINTTYPE) || 
+            (geom->type == MULTIPOINTTYPE && lwgeom_count_vertices(geom) == 1) )
+       {
+               lwcollection_add_lwgeom(col, lwgeom_clone_deep(geom));
+               return 1;
+       }
+       
+       /* Haven't set up a clipping box yet? Make one then recurse again */
+       if ( ! clip )
+       {
+               double height, width;
+               const GBOX *box = lwgeom_get_bbox(geom);
+               GBOX square;
+               if ( ! box ) return 0;
+               width = box->xmax - box->xmin;
+               height = box->ymax - box->ymin;
+               if ( width > height )
+               {
+                       square.xmin = box->xmin;
+                       square.xmax = box->xmax;
+                       square.ymin = (box->ymin + box->ymax - width)/2;
+                       square.ymax = (box->ymin + box->ymax + width)/2;
+               }
+               else
+               {
+                       square.ymin = box->ymin;
+                       square.ymax = box->ymax;
+                       square.xmin = (box->xmin + box->xmax - height)/2;
+                       square.xmax = (box->xmin + box->xmax + height)/2;
+               }
+               return lwgeom_subdivide_recursive(geom, maxvertices, col, &square);
+       }
+       /* Everything is in place! Let's start subdividing! */
+       else
+       {
+               int npoints = lwgeom_npoints_in_rect(geom, clip);
+               /* Clipping rectangle didn't hit anything! No-op! */
+               if ( npoints == 0 )
+               {
+                       return 0;
+               }
+               /* Clipping rectangle has reached target size: apply it! */
+               else if ( npoints < maxvertices )
+               {
+                       lwcollection_add_lwgeom(col, lwgeom_clip_by_rect(geom, clip->xmin, clip->ymin, clip->xmax, clip->ymax));
+                       return npoints;
+               }
+               /* Clipping rectangle too small: subdivide more! */
+               else
+               {
+                       int i, n = 0;
+                       GBOX boxes[4];
+                       double width = FP_TOLERANCE + (clip->xmax - clip->xmin)/2;
+                       for( i = 0; i < 4; i++ )
+                       {
+                               int ix = i / 2;
+                               int iy = i % 2;
+                               boxes[i].xmin = clip->xmin + ix * width;
+                               boxes[i].xmax = clip->xmin + width + ix * width;
+                               boxes[i].ymin = clip->ymin + iy * width;
+                               boxes[i].ymax = clip->ymin + width + iy * width;
+                               n += lwgeom_subdivide_recursive(geom, maxvertices, col, &(boxes[i]));
+                       }
+                       return n;
+               }
+       }
+       return 0;
+}
+
+LWCOLLECTION *
+lwgeom_subdivide(const LWGEOM *geom, int maxvertices)
+{
+       int n = 0;
+       LWCOLLECTION *col;
+       col = lwcollection_construct_empty(COLLECTIONTYPE, geom->srid, lwgeom_has_z(geom), lwgeom_has_m(geom));
+       n = lwgeom_subdivide_recursive(geom, maxvertices, col, NULL);
+       return col;
+}
+
+