From: Paul Ramsey Date: Thu, 5 Oct 2017 19:34:16 +0000 (+0000) Subject: Move simplify code to an in_place basis, X-Git-Tag: 2.5.0alpha~452 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=aa2c187f827dc20f8fab9111e02868a010621398;p=postgis Move simplify code to an in_place basis, references #3877 git-svn-id: http://svn.osgeo.org/postgis/trunk@15908 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/liblwgeom/cunit/cu_algorithm.c b/liblwgeom/cunit/cu_algorithm.c index c73fb8104..ea6612350 100644 --- a/liblwgeom/cunit/cu_algorithm.c +++ b/liblwgeom/cunit/cu_algorithm.c @@ -927,59 +927,59 @@ static void test_geohash_point_as_int(void) static void test_lwgeom_simplify(void) { - LWGEOM *l; - LWGEOM *g; - char *ewkt; - - /* Simplify but only so far... */ - g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE); - l = lwgeom_simplify(g, 10, LW_TRUE); - ewkt = lwgeom_to_ewkt(l); - CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,0 0)"); - lwgeom_free(g); - lwgeom_free(l); - lwfree(ewkt); - - /* Simplify but only so far... */ - g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE); - l = lwgeom_simplify(g, 10, LW_TRUE); - ewkt = lwgeom_to_ewkt(l); - CU_ASSERT_STRING_EQUAL(ewkt, "POLYGON((0 0,1 0,1 1,0 0))"); - lwgeom_free(g); - lwgeom_free(l); - lwfree(ewkt); - - /* Simplify and collapse */ - g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE); - l = lwgeom_simplify(g, 10, LW_FALSE); - CU_ASSERT_EQUAL(l, NULL); - lwgeom_free(g); - lwgeom_free(l); - - /* Simplify and collapse */ - g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE); - l = lwgeom_simplify(g, 10, LW_FALSE); - CU_ASSERT_EQUAL(l, NULL); - lwgeom_free(g); - lwgeom_free(l); - - /* Not simplifiable */ - g = lwgeom_from_wkt("LINESTRING(0 0, 50 1.00001, 100 0)", LW_PARSER_CHECK_NONE); - l = lwgeom_simplify(g, 1.0, LW_FALSE); - ewkt = lwgeom_to_ewkt(l); - CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,50 1.00001,100 0)"); - lwgeom_free(g); - lwgeom_free(l); - lwfree(ewkt); - - /* Simplifiable */ - g = lwgeom_from_wkt("LINESTRING(0 0,50 0.99999,100 0)", LW_PARSER_CHECK_NONE); - l = lwgeom_simplify(g, 1.0, LW_FALSE); - ewkt = lwgeom_to_ewkt(l); - CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,100 0)"); - lwgeom_free(g); - lwgeom_free(l); - lwfree(ewkt); + LWGEOM *l; + LWGEOM *g; + char *ewkt; + + /* Simplify but only so far... */ + g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE); + l = lwgeom_simplify(g, 10, LW_TRUE); + ewkt = lwgeom_to_ewkt(l); + CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,0 0)"); + lwgeom_free(g); + lwgeom_free(l); + lwfree(ewkt); + + /* Simplify but only so far... */ + g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE); + l = lwgeom_simplify(g, 10, LW_TRUE); + ewkt = lwgeom_to_ewkt(l); + CU_ASSERT_STRING_EQUAL(ewkt, "POLYGON((0 0,1 0,1 1,0 0))"); + lwgeom_free(g); + lwgeom_free(l); + lwfree(ewkt); + + /* Simplify and collapse */ + g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE); + l = lwgeom_simplify(g, 10, LW_FALSE); + CU_ASSERT_EQUAL(l, NULL); + lwgeom_free(g); + lwgeom_free(l); + + /* Simplify and collapse */ + g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE); + l = lwgeom_simplify(g, 10, LW_FALSE); + CU_ASSERT_EQUAL(l, NULL); + lwgeom_free(g); + lwgeom_free(l); + + /* Not simplifiable */ + g = lwgeom_from_wkt("LINESTRING(0 0, 50 1.00001, 100 0)", LW_PARSER_CHECK_NONE); + l = lwgeom_simplify(g, 1.0, LW_FALSE); + ewkt = lwgeom_to_ewkt(l); + CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,50 1.00001,100 0)"); + lwgeom_free(g); + lwgeom_free(l); + lwfree(ewkt); + + /* Simplifiable */ + g = lwgeom_from_wkt("LINESTRING(0 0,50 0.99999,100 0)", LW_PARSER_CHECK_NONE); + l = lwgeom_simplify(g, 1.0, LW_FALSE); + ewkt = lwgeom_to_ewkt(l); + CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,100 0)"); + lwgeom_free(g); + lwgeom_free(l); + lwfree(ewkt); } diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in index 76d88c192..0e902a7e4 100644 --- a/liblwgeom/liblwgeom.h.in +++ b/liblwgeom/liblwgeom.h.in @@ -1224,6 +1224,8 @@ extern int lwgeom_is_clockwise(LWGEOM *lwgeom); +extern LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed); + /**************************************************************** * READ/WRITE FUNCTIONS * @@ -1237,16 +1239,14 @@ extern int lwgeom_is_clockwise(LWGEOM *lwgeom); extern void lwgeom_reverse_in_place(LWGEOM *lwgeom); extern void lwgeom_force_clockwise(LWGEOM *lwgeom); extern void lwgeom_longitude_shift(LWGEOM *lwgeom); -extern LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed); +extern void lwgeom_simplify_in_place(LWGEOM *igeom, double dist, int preserve_collapsed); extern void lwgeom_affine(LWGEOM *geom, const AFFINE *affine); extern void lwgeom_scale(LWGEOM *geom, const POINT4D *factors); +extern LWGEOM* lwgeom_remove_repeated_points(const LWGEOM *in, double tolerance); -/****************************************************************/ - - /** * @brief wrap geometry on given cut x value * @@ -1656,7 +1656,6 @@ extern LWBOUNDINGCIRCLE* lwgeom_calculate_mbc(const LWGEOM* g); /** * Remove repeated points! */ -extern LWGEOM* lwgeom_remove_repeated_points(const LWGEOM *in, double tolerance); /** * Swap ordinate values in every vertex of the geometry. diff --git a/liblwgeom/liblwgeom_internal.h b/liblwgeom/liblwgeom_internal.h index f7a9877e0..86eefc552 100644 --- a/liblwgeom/liblwgeom_internal.h +++ b/liblwgeom/liblwgeom_internal.h @@ -201,10 +201,7 @@ int lwcollection_count_vertices(LWCOLLECTION *col); /** * @param minpts minimun number of points to retain, if possible. */ -POINTARRAY* ptarray_simplify(POINTARRAY *inpts, double epsilon, unsigned int minpts); -LWLINE* lwline_simplify(const LWLINE *iline, double dist, int preserve_collapsed); -LWPOLY* lwpoly_simplify(const LWPOLY *ipoly, double dist, int preserve_collapsed); -LWCOLLECTION* lwcollection_simplify(const LWCOLLECTION *igeom, double dist, int preserve_collapsed); +void ptarray_simplify_in_place(POINTARRAY *pa, double epsilon, unsigned int minpts); /* * The possible ways a pair of segments can interact. Returned by lw_segment_intersects @@ -374,6 +371,12 @@ int point4d_transform(POINT4D *pt, projPJ srcpj, projPJ dstpj); */ void ptarray_longitude_shift(POINTARRAY *pa); +/* +* Support for in place modification of point arrays, fast +* function to move coordinate values around +*/ +void ptarray_copy_point(POINTARRAY *pa, int from, int to); + /* * Reverse */ diff --git a/liblwgeom/lwcollection.c b/liblwgeom/lwcollection.c index d6d36e0a8..b0e8d080e 100644 --- a/liblwgeom/lwcollection.c +++ b/liblwgeom/lwcollection.c @@ -523,22 +523,6 @@ int lwcollection_count_vertices(LWCOLLECTION *col) return v; } -LWCOLLECTION* lwcollection_simplify(const LWCOLLECTION *igeom, double dist, int preserve_collapsed) -{ - int i; - LWCOLLECTION *out = lwcollection_construct_empty(igeom->type, igeom->srid, FLAGS_GET_Z(igeom->flags), FLAGS_GET_M(igeom->flags)); - - if( lwcollection_is_empty(igeom) ) - return out; /* should we return NULL instead ? */ - - for( i = 0; i < igeom->ngeoms; i++ ) - { - LWGEOM *ngeom = lwgeom_simplify(igeom->geoms[i], dist, preserve_collapsed); - if ( ngeom ) out = lwcollection_add_lwgeom(out, ngeom); - } - - return out; -} int lwcollection_allows_subtype(int collectiontype, int subtype) { diff --git a/liblwgeom/lwgeom.c b/liblwgeom/lwgeom.c index ef2aa75a7..2ce8ddf97 100644 --- a/liblwgeom/lwgeom.c +++ b/liblwgeom/lwgeom.c @@ -1630,27 +1630,116 @@ void lwgeom_set_srid(LWGEOM *geom, int32_t srid) } } -LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed) +/**************************************************************/ + + +void +lwgeom_simplify_in_place(LWGEOM *geom, double epsilon, int preserve_collapsed) { - switch (igeom->type) + switch (geom->type) { - case POINTTYPE: - case MULTIPOINTTYPE: - return lwgeom_clone(igeom); - case LINETYPE: - return (LWGEOM*)lwline_simplify((LWLINE*)igeom, dist, preserve_collapsed); - case POLYGONTYPE: - return (LWGEOM*)lwpoly_simplify((LWPOLY*)igeom, dist, preserve_collapsed); - case MULTILINETYPE: - case MULTIPOLYGONTYPE: - case COLLECTIONTYPE: - return (LWGEOM*)lwcollection_simplify((LWCOLLECTION *)igeom, dist, preserve_collapsed); - default: - lwerror("%s: unsupported geometry type: %s", __func__, lwtype_name(igeom->type)); + /* No-op! Cannot simplify points */ + case POINTTYPE: + return; + case LINETYPE: + { + LWLINE *g = (LWLINE*)(geom); + POINTARRAY *pa = g->points; + ptarray_simplify_in_place(pa, epsilon, 2); + /* Invalid output */ + if (pa->npoints == 1 && pa->maxpoints > 1) + { + /* Use first point as last point */ + if (preserve_collapsed) + { + pa->npoints = 2; + ptarray_copy_point(pa, 0, 1); + } + /* Finish the collapse process */ + else + { + pa->npoints = 0; + } + } + /* Duped output, force collapse */ + if (pa->npoints == 2 && !preserve_collapsed) + { + if (p2d_same(getPoint2d_cp(pa, 0), getPoint2d_cp(pa, 1))) + pa->npoints = 0; + } + break; + } + case POLYGONTYPE: + { + int i, j = 0; + LWPOLY *g = (LWPOLY*)(geom); + for (i = 0; i < g->nrings; i++) + { + POINTARRAY *pa = g->rings[i]; + /* Only stop collapse on first ring */ + int minpoints = (preserve_collapsed && i == 0) ? 4 : 0; + /* Skip zero'ed out rings */ + if(!pa) + continue; + ptarray_simplify_in_place(pa, epsilon, minpoints); + /* Drop collapsed rings */ + if(pa->npoints < 4) + { + ptarray_free(pa); + continue; + } + g->rings[j++] = pa; + } + /* Update ring count */ + g->nrings = j; + break; + } + /* Can process all multi* types as generic collection */ + case MULTIPOINTTYPE: + case MULTILINETYPE: + case MULTIPOLYGONTYPE: + case COLLECTIONTYPE: + { + int i, j = 0; + LWCOLLECTION *col = (LWCOLLECTION*)geom; + for (i = 0; i < col->ngeoms; i++) + { + LWGEOM *g = col->geoms[i]; + if (!g) continue; + lwgeom_simplify_in_place(g, epsilon, preserve_collapsed); + /* Drop zero'ed out geometries */ + if(lwgeom_is_empty(g)) + { + lwgeom_free(g); + continue; + } + col->geoms[j++] = g; + } + /* Update geometry count */ + col->ngeoms = j; + break; + } + default: + { + lwerror("%s: unsupported geometry type: %s", __func__, lwtype_name(geom->type)); + break; + } } - return NULL; + return; } + +/**************************************************************/ + +LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed) +{ + LWGEOM *lwgeom_out = lwgeom_clone_deep(igeom); + lwgeom_simplify_in_place(lwgeom_out, dist, preserve_collapsed); + if (lwgeom_is_empty(lwgeom_out)) return NULL; + return lwgeom_out; +} + + double lwgeom_area(const LWGEOM *geom) { int type = geom->type; diff --git a/liblwgeom/lwgeom_api.c b/liblwgeom/lwgeom_api.c index b1b888fad..1c7112c34 100644 --- a/liblwgeom/lwgeom_api.c +++ b/liblwgeom/lwgeom_api.c @@ -448,6 +448,42 @@ ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d) } } +void +ptarray_copy_point(POINTARRAY *pa, int from, int to) +{ + int ndims = FLAGS_NDIMS(pa->flags); + switch (ndims) + { + case 2: + { + POINT2D *p_from = (POINT2D*)(getPoint_internal(pa, from)); + POINT2D *p_to = (POINT2D*)(getPoint_internal(pa, to)); + *p_to = *p_from; + return; + } + case 3: + { + POINT3D *p_from = (POINT3D*)(getPoint_internal(pa, from)); + POINT3D *p_to = (POINT3D*)(getPoint_internal(pa, to)); + *p_to = *p_from; + return; + } + case 4: + { + POINT4D *p_from = (POINT4D*)(getPoint_internal(pa, from)); + POINT4D *p_to = (POINT4D*)(getPoint_internal(pa, to)); + *p_to = *p_from; + return; + } + default: + { + lwerror("%s: unsupported number of dimensions - %d", __func__, ndims); + return; + } + } + return; +} + /************************************************ * debugging routines diff --git a/liblwgeom/lwline.c b/liblwgeom/lwline.c index 9d34b6a65..eec96cab5 100644 --- a/liblwgeom/lwline.c +++ b/liblwgeom/lwline.c @@ -532,44 +532,6 @@ int lwline_count_vertices(LWLINE *line) return line->points->npoints; } -LWLINE* lwline_simplify(const LWLINE *iline, double dist, int preserve_collapsed) -{ - static const int minvertices = 2; /* TODO: allow setting this */ - LWLINE *oline; - POINTARRAY *pa; - - LWDEBUG(2, "function called"); - - /* Skip empty case */ - if( lwline_is_empty(iline) ) - return NULL; - - pa = ptarray_simplify(iline->points, dist, minvertices); - if ( ! pa ) return NULL; - - /* Make sure single-point collapses have two points */ - if ( pa->npoints == 1 ) - { - /* Make sure single-point collapses have two points */ - if ( preserve_collapsed ) - { - POINT4D pt; - getPoint4d_p(pa, 0, &pt); - ptarray_append_point(pa, &pt, LW_TRUE); - } - /* Return null for collapse */ - else - { - ptarray_free(pa); - return NULL; - } - } - - oline = lwline_construct(iline->srid, NULL, pa); - oline->type = iline->type; - return oline; -} - double lwline_length(const LWLINE *line) { if ( lwline_is_empty(line) ) diff --git a/liblwgeom/lwpoly.c b/liblwgeom/lwpoly.c index e32bcbeb2..643d67557 100644 --- a/liblwgeom/lwpoly.c +++ b/liblwgeom/lwpoly.c @@ -450,62 +450,6 @@ int lwpoly_count_vertices(LWPOLY *poly) return v; } -LWPOLY* lwpoly_simplify(const LWPOLY *ipoly, double dist, int preserve_collapsed) -{ - int i; - LWPOLY *opoly = lwpoly_construct_empty(ipoly->srid, FLAGS_GET_Z(ipoly->flags), FLAGS_GET_M(ipoly->flags)); - - LWDEBUGF(2, "%s: simplifying polygon with %d rings", __func__, ipoly->nrings); - - if ( lwpoly_is_empty(ipoly) ) - { - lwpoly_free(opoly); - return NULL; - } - - for ( i = 0; i < ipoly->nrings; i++ ) - { - POINTARRAY *opts; - int minvertices = 0; - - /* We'll still let holes collapse, but if we're preserving */ - /* and this is a shell, we ensure it is kept */ - if ( preserve_collapsed && i == 0 ) - minvertices = 4; - - opts = ptarray_simplify(ipoly->rings[i], dist, minvertices); - - LWDEBUGF(3, "ring%d simplified from %d to %d points", i, ipoly->rings[i]->npoints, opts->npoints); - - /* Less points than are needed to form a closed ring, we can't use this */ - if ( opts->npoints < 4 ) - { - LWDEBUGF(3, "ring%d skipped (% pts)", i, opts->npoints); - ptarray_free(opts); - if ( i ) continue; - else break; /* Don't scan holes if shell is collapsed */ - } - - /* Add ring to simplified polygon */ - if( lwpoly_add_ring(opoly, opts) == LW_FAILURE ) - { - lwpoly_free(opoly); - return NULL; - } - } - - LWDEBUGF(3, "simplified polygon with %d rings", ipoly->nrings); - opoly->type = ipoly->type; - - if( lwpoly_is_empty(opoly) ) - { - lwpoly_free(opoly); - return NULL; - } - - return opoly; -} - /** * Find the area of the outer ring - sum (area of inner rings). */ diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c index 5f2ba1571..f0d46c13a 100644 --- a/liblwgeom/ptarray.c +++ b/liblwgeom/ptarray.c @@ -1517,8 +1517,10 @@ ptarray_remove_repeated_points(const POINTARRAY *in, double tolerance) return ptarray_remove_repeated_points_minpoints(in, tolerance, 2); } +/************************************************************************/ + static void -ptarray_dp_findsplit(POINTARRAY *pts, int p1, int p2, int *split, double *dist) +ptarray_dp_findsplit_in_place(const POINTARRAY *pts, int p1, int p2, int *split, double *dist) { int k; const POINT2D *pk, *pa, *pb; @@ -1556,75 +1558,105 @@ ptarray_dp_findsplit(POINTARRAY *pts, int p1, int p2, int *split, double *dist) } } *dist = d; - - } /* length---should be redone if can == 0 */ + } else { LWDEBUG(3, "segment too short, no split/no dist"); *dist = -1; } +} +static int +int_cmp(const void *a, const void *b) +{ + /* casting pointer types */ + const int *ia = (const int *)a; + const int *ib = (const int *)b; + /* returns negative if b > a and positive if a > b */ + return *ia - *ib; } -POINTARRAY * -ptarray_simplify(POINTARRAY *inpts, double epsilon, unsigned int minpts) +void +ptarray_simplify_in_place(POINTARRAY *pa, double epsilon, unsigned int minpts) { - int *stack; /* recursion stack */ - int sp=-1; /* recursion stack pointer */ + static size_t stack_size = 256; + int *stack, *outlist; /* recursion stack */ + int stack_static[stack_size]; + int outlist_static[stack_size]; + int sp = -1; /* recursion stack pointer */ int p1, split; + int outn = 0; + int pai = 0; + int i; double dist; - POINTARRAY *outpts; - POINT4D pt; - double eps_sqr = epsilon * epsilon; - /* Allocate recursion stack */ - stack = lwalloc(sizeof(int)*inpts->npoints); - - p1 = 0; - stack[++sp] = inpts->npoints-1; - - LWDEBUGF(2, "Input has %d pts and %d dims", inpts->npoints, - FLAGS_NDIMS(inpts->flags)); + /* Do not try to simplify really short things */ + if (pa->npoints < 3) return; - /* Allocate output POINTARRAY, and add first point. */ - outpts = ptarray_construct_empty(FLAGS_GET_Z(inpts->flags), FLAGS_GET_M(inpts->flags), inpts->npoints); - getPoint4d_p(inpts, 0, &pt); - ptarray_append_point(outpts, &pt, LW_FALSE); + /* Only heap allocate book-keeping arrays if necessary */ + if (pa->npoints > stack_size) + { + stack = lwalloc(sizeof(int) * pa->npoints); + outlist = lwalloc(sizeof(int) * pa->npoints); + } + else + { + stack = stack_static; + outlist = outlist_static; + } - LWDEBUG(3, "Added P0 to simplified point array (size 1)"); + p1 = 0; + stack[++sp] = pa->npoints-1; + /* Add first point to output list */ + outlist[outn++] = 0; do { + ptarray_dp_findsplit_in_place(pa, p1, stack[sp], &split, &dist); - ptarray_dp_findsplit(inpts, p1, stack[sp], &split, &dist); - - LWDEBUGF(3, "Farthest point from P%d-P%d is P%d (dist. %g)", p1, stack[sp], split, dist); - - if (dist > eps_sqr || ( outpts->npoints+sp+1 < minpts && dist >= 0 ) ) + if ((dist > eps_sqr) || ((outn + sp+1 < minpts) && (dist >= 0))) { - LWDEBUGF(4, "Added P%d to stack (outpts:%d)", split, sp); stack[++sp] = split; } else { - getPoint4d_p(inpts, stack[sp], &pt); - LWDEBUGF(4, "npoints , minpoints %d %d", outpts->npoints, minpts); - ptarray_append_point(outpts, &pt, LW_FALSE); - - LWDEBUGF(4, "Added P%d to simplified point array (size: %d)", stack[sp], outpts->npoints); - + outlist[outn++] = stack[sp]; p1 = stack[sp--]; } + } + while (!(sp<0)); - LWDEBUGF(4, "stack pointer = %d", sp); + /* Put list of retained points into order */ + qsort(outlist, outn, sizeof(int), int_cmp); + /* Copy retained points to front of array */ + for (i = 0; i < outn; i++) + { + int j = outlist[i]; + /* Indexes the same, means no copy required */ + if (j == pai) + { + pai++; + continue; + } + /* Indexes different, copy value down */ + ptarray_copy_point(pa, j, pai++); } - while (! (sp<0) ); - lwfree(stack); - return outpts; + /* Adjust point count on array */ + pa->npoints = outn; + + /* Only free if arrays are on heap */ + if (stack != stack_static) + lwfree(stack); + if (outlist != outlist_static) + lwfree(outlist); + + return; } +/************************************************************************/ + /** * Find the 2d length of the given #POINTARRAY, using circular * arc interpolation between each coordinate triple.