From: Nicklas Avén Date: Tue, 20 Mar 2018 13:11:17 +0000 (+0000) Subject: Add ChaikinSmoothing #4050 X-Git-Tag: 2.5.0alpha~52 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=744863ab8cd861471abb9e435175ac0bda878db1;p=postgis Add ChaikinSmoothing #4050 git-svn-id: http://svn.osgeo.org/postgis/trunk@16479 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/NEWS b/NEWS index 396f982e9..14b2f08e6 100644 --- a/NEWS +++ b/NEWS @@ -2,6 +2,7 @@ PostGIS 2.5.0 2018/xx/xx * New Features * + - #4050, ST_ChaikinSmoothing (Nicklas Avén) - #3989, ST_Buffer single sided option (Stephen Knox) - #3876, ST_Angle function (Rémi Cura) - #3564, ST_LineInterpolatePoints (Dan Baston) diff --git a/doc/reference_processing.xml b/doc/reference_processing.xml index 07416df8c..86cf8327b 100644 --- a/doc/reference_processing.xml +++ b/doc/reference_processing.xml @@ -3081,7 +3081,7 @@ FROM (SELECT ST_Buffer('POINT(1 3)', 10,12) As the_geom) As foo; - + ST_SimplifyVW Returns a "simplified" version of the given geometry using the Visvalingam-Whyatt algorithm @@ -3129,7 +3129,63 @@ LINESTRING(5 2,7 25,10 10) , , , Topology - + + + + ST_ChaikinSmoothing + Returns a "smoothed" version of the given geometry using the Chaikin algorithm + + + + + + geometry ST_ChaikinSmoothing + geometry geomA + integer nIterations + boolean preserveEndPoints + + + + + + Description + Returns a "smoothed" version of the given geometry using the Chaikin algorithm. + See Chaikins-Algorithm for an explanation of the process. + For each iteration the number of vertex points will double. + The function puts new vertex points at 1/4 of the line before and after each point and removes the original point. + To reduce the number of points use one of the simplification functions on the result. + The new points gets interpolated values for all included dimensions, also z and m. + + Note that returned geometry will get more points than the original. + To reduce the number of points again use one of the simplification functions on the result. + (see and ) + Second argument, number of iterations is limited to max 5 iterations + Note third argument is only valid for polygons, and will be ignored for linestrings + This function handles 3D and the third dimension will affect the result. + Availability: 2.5.0 + + + + Examples + A triangle is smoothed + + +select ST_AsText(ST_ChaikinSmoothing(geom)) smoothed +FROM (SELECT 'LINESTRING(0 0, 8 8, 0 16)'::geometry geom) As foo; +-result + smoothed +-----------+-------------------+ +LINESTRING(0 0,6 6,6 10,0 16) + + + + + See Also + , + + + + ST_SetEffectiveArea diff --git a/liblwgeom/Makefile.in b/liblwgeom/Makefile.in index d14dc6b36..d2649f15f 100644 --- a/liblwgeom/Makefile.in +++ b/liblwgeom/Makefile.in @@ -114,6 +114,7 @@ SA_OBJS = \ lwgeom_wrapx.o \ lwunionfind.o \ effectivearea.o \ + lwchaikins.o \ lwkmeans.o \ varint.o diff --git a/liblwgeom/cunit/Makefile.in b/liblwgeom/cunit/Makefile.in index 8ff18ad06..10553ec8a 100644 --- a/liblwgeom/cunit/Makefile.in +++ b/liblwgeom/cunit/Makefile.in @@ -36,6 +36,7 @@ OBJS= \ cu_tree.o \ cu_measures.o \ cu_effectivearea.o \ + cu_chaikin.o \ cu_node.o \ cu_clip_by_rect.o \ cu_libgeom.o \ diff --git a/liblwgeom/cunit/cu_chaikin.c b/liblwgeom/cunit/cu_chaikin.c new file mode 100644 index 000000000..b5d52243f --- /dev/null +++ b/liblwgeom/cunit/cu_chaikin.c @@ -0,0 +1,66 @@ +/********************************************************************** + * + * PostGIS - Spatial Types for PostgreSQL + * http://postgis.net + * + * Copyright 2012 Sandro Santilli + * + * This is free software; you can redistribute and/or modify it under + * the terms of the GNU General Public Licence. See the COPYING file. + * + **********************************************************************/ + +#include +#include +#include +#include "CUnit/Basic.h" + +#include "liblwgeom_internal.h" +#include "cu_tester.h" + +static void do_test_chaikin(char *geom_txt,char *expected, int n_iterations, int preserve_end_points) +{ + LWGEOM *geom_in, *geom_out; + char *out_txt; + geom_in = lwgeom_from_wkt(geom_txt, LW_PARSER_CHECK_NONE); + geom_out = lwgeom_chaikin(geom_in, n_iterations, preserve_end_points); + out_txt = lwgeom_to_wkt(geom_out, WKT_EXTENDED, 3, NULL); + if(strcmp(expected, out_txt)) + printf("%s is not equal to %s\n", expected, out_txt); + CU_ASSERT_STRING_EQUAL(expected, out_txt) + lwfree(out_txt); + lwgeom_free(geom_in); + lwgeom_free(geom_out); + return; +} + +static void do_test_chaikin_lines(void) +{ + /*Simpliest test*/ + do_test_chaikin("LINESTRING(0 0,8 8,16 0)","LINESTRING(0 0,6 6,10 6,16 0)", 1, 1); + /*2 iterations*/ + do_test_chaikin("LINESTRING(0 0,8 8,16 0)","LINESTRING(0 0,4.5 4.5,7 6,9 6,11.5 4.5,16 0)", 2, 1); + /*check so it really ignores preserve_end_points set to off for linestrings*/ + do_test_chaikin("LINESTRING(0 0,8 8,16 0)","LINESTRING(0 0,4.5 4.5,7 6,9 6,11.5 4.5,16 0)", 2, 0); + return; +} + +static void do_test_chaikin_polygons(void) +{ + /*Simpliest test*/ + do_test_chaikin("POLYGON((0 0,8 8,16 0,0 0))","POLYGON((0 0,6 6,10 6,14 2,12 0,0 0))", 1, 1); + /*2 iterations*/ + do_test_chaikin("POLYGON((0 0,8 8,16 0,0 0))","POLYGON((0 0,4.5 4.5,7 6,9 6,11 5,13 3,13.5 1.5,12.5 0.5,9 0,0 0))", 2, 1); + /*2 iterations without preserving end points*/ + do_test_chaikin("POLYGON((0 0,8 8,16 0,0 0))","POLYGON((3 3,5 5,7 6,9 6,11 5,13 3,13.5 1.5,12.5 0.5,10 0,6 0,3.5 0.5,2.5 1.5,3 3))", 2, 0); + return; +} + + +void chaikin_suite_setup(void); +void chaikin_suite_setup(void) +{ + CU_pSuite suite = CU_add_suite("chaikin",NULL,NULL); + PG_ADD_TEST(suite, do_test_chaikin_lines); + PG_ADD_TEST(suite, do_test_chaikin_polygons); +} diff --git a/liblwgeom/cunit/cu_tester.c b/liblwgeom/cunit/cu_tester.c index 28d4df9c2..f44bb7876 100644 --- a/liblwgeom/cunit/cu_tester.c +++ b/liblwgeom/cunit/cu_tester.c @@ -47,6 +47,7 @@ extern void libgeom_suite_setup(void); extern void lwstroke_suite_setup(void); extern void measures_suite_setup(void); extern void effectivearea_suite_setup(void); +extern void chaikin_suite_setup(void); extern void minimum_bounding_circle_suite_setup(void); extern void misc_suite_setup(void); extern void node_suite_setup(void); @@ -97,6 +98,7 @@ PG_SuiteSetup setupfuncs[] = lwstroke_suite_setup, measures_suite_setup, effectivearea_suite_setup, + chaikin_suite_setup, minimum_bounding_circle_suite_setup, misc_suite_setup, node_suite_setup, diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in index 290908e2f..4b54cc41a 100644 --- a/liblwgeom/liblwgeom.h.in +++ b/liblwgeom/liblwgeom.h.in @@ -990,6 +990,7 @@ extern LWGEOM* lwgeom_force_3dm(const LWGEOM *geom); extern LWGEOM* lwgeom_force_4d(const LWGEOM *geom); extern LWGEOM* lwgeom_set_effective_area(const LWGEOM *igeom, int set_area, double area); +extern LWGEOM* lwgeom_chaikin(const LWGEOM *igeom, int n_iterations, int preserve_endpoint); /* * Force to use SFS 1.1 geometry type diff --git a/liblwgeom/lwchaikins.c b/liblwgeom/lwchaikins.c new file mode 100644 index 000000000..b69c0dd1b --- /dev/null +++ b/liblwgeom/lwchaikins.c @@ -0,0 +1,202 @@ + +/********************************************************************** + * + * PostGIS - Spatial Types for PostgreSQL + * http://postgis.net + * + * PostGIS is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 2 of the License, or + * (at your option) any later version. + * + * PostGIS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with PostGIS. If not, see . + * + ********************************************************************** + * + * Copyright 2018 Nicklas Avén + * + **********************************************************************/ + + + #include "liblwgeom_internal.h" + + + +static POINTARRAY * ptarray_chaikin(POINTARRAY *inpts, int preserve_endpoint, int isclosed) +{ + uint32_t p, i, n_out_points=0, p1_set=0, p2_set=0; + POINT4D p1, p2; + POINTARRAY *opts; + double *dlist; + double deltaval, quarter_delta, val1, val2; + int ndims = 2 + ptarray_has_z(inpts) + ptarray_has_m(inpts); + int new_npoints = inpts->npoints * 2; + opts = ptarray_construct_empty(FLAGS_GET_Z(inpts->flags), FLAGS_GET_M(inpts->flags), new_npoints); + + dlist = (double*)(opts->serialized_pointlist); + + p1 = getPoint4d(inpts, 0); + /*if first point*/ + if(preserve_endpoint) + { + ptarray_append_point(opts, &p1, LW_TRUE); + n_out_points++; + } + + for (p=1;pnpoints;p++) + { + memcpy(&p2, &p1, sizeof(POINT4D)); + p1 = getPoint4d(inpts, p); + if(p>0) + { + p1_set = p2_set = 0; + for (i=0;i 1) + { + dlist[n_out_points * ndims + i] = val2 + quarter_delta; + p1_set = 1; + } + if(!preserve_endpoint || p < inpts->npoints - 1) + { + dlist[(n_out_points + p1_set) * ndims + i] = val1 - quarter_delta; + p2_set = 1; + } + } + n_out_points+=(p1_set + p2_set); + } + } + + /*if last point*/ + if(preserve_endpoint) + { + opts->npoints = n_out_points; + ptarray_append_point(opts, &p1, LW_TRUE); + n_out_points++; + } + + if(isclosed && !preserve_endpoint) + { + opts->npoints = n_out_points; + POINT4D first_point = getPoint4d(opts, 0); + ptarray_append_point(opts, &first_point, LW_TRUE); + n_out_points++; + } + opts->npoints = n_out_points; + + return opts; + +} + +static LWLINE* lwline_chaikin(const LWLINE *iline, int n_iterations) +{ + POINTARRAY *pa, *pa_new; + int j; + + if( lwline_is_empty(iline)) + return lwline_clone(iline); + + + LWLINE *oline = lwline_construct_empty(iline->srid, FLAGS_GET_Z(iline->flags), FLAGS_GET_M(iline->flags)); + pa = iline->points; + for (j=0;j0) + ptarray_free(pa); + pa=pa_new; + } + oline = lwline_construct(iline->srid, NULL, pa); + + oline->type = iline->type; + return oline; + +} + + +static LWPOLY* lwpoly_chaikin(const LWPOLY *ipoly, int n_iterations, int preserve_endpoint) +{ + uint32_t i, j; + POINTARRAY *pa, *pa_new; + LWPOLY *opoly = lwpoly_construct_empty(ipoly->srid, FLAGS_GET_Z(ipoly->flags), FLAGS_GET_M(ipoly->flags)); + + if( lwpoly_is_empty(ipoly) ) + return opoly; + for (i = 0; i < ipoly->nrings; i++) + { + pa = ipoly->rings[i]; + for(j=0;j0) + ptarray_free(pa); + pa=pa_new; + } + if(pa->npoints>=4) + { + if( lwpoly_add_ring(opoly,pa ) == LW_FAILURE ) + return NULL; + } + } + + opoly->type = ipoly->type; + + if( lwpoly_is_empty(opoly) ) + return NULL; + + return opoly; + +} + + +static LWCOLLECTION* lwcollection_chaikin(const LWCOLLECTION *igeom, int n_iterations, int preserve_endpoint) +{ + LWDEBUG(2, "Entered lwcollection_set_effective_area"); + uint32_t 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_chaikin(igeom->geoms[i],n_iterations,preserve_endpoint); + if ( ngeom ) out = lwcollection_add_lwgeom(out, ngeom); + } + + return out; +} + + +LWGEOM* lwgeom_chaikin(const LWGEOM *igeom, int n_iterations, int preserve_endpoint) +{ + LWDEBUG(2, "Entered lwgeom_set_effective_area"); + switch (igeom->type) + { + case POINTTYPE: + case MULTIPOINTTYPE: + return lwgeom_clone(igeom); + case LINETYPE: + return (LWGEOM*)lwline_chaikin((LWLINE*)igeom, n_iterations); + case POLYGONTYPE: + return (LWGEOM*)lwpoly_chaikin((LWPOLY*)igeom,n_iterations,preserve_endpoint); + case MULTILINETYPE: + case MULTIPOLYGONTYPE: + case COLLECTIONTYPE: + return (LWGEOM*)lwcollection_chaikin((LWCOLLECTION *)igeom,n_iterations,preserve_endpoint); + default: + lwerror("lwgeom_chaikin: unsupported geometry type: %s",lwtype_name(igeom->type)); + } + return NULL; +} diff --git a/postgis/lwgeom_functions_analytic.c b/postgis/lwgeom_functions_analytic.c index 5e18473f4..18dde2a11 100644 --- a/postgis/lwgeom_functions_analytic.c +++ b/postgis/lwgeom_functions_analytic.c @@ -127,6 +127,50 @@ Datum LWGEOM_SetEffectiveArea(PG_FUNCTION_ARGS) PG_RETURN_POINTER(result); } +PG_FUNCTION_INFO_V1(LWGEOM_ChaikinSmoothing); +Datum LWGEOM_ChaikinSmoothing(PG_FUNCTION_ARGS) +{ + GSERIALIZED *geom = PG_GETARG_GSERIALIZED_P(0); + GSERIALIZED *result; + int type = gserialized_get_type(geom); + LWGEOM *in; + LWGEOM *out; + int preserve_endpoints=1; + int n_iterations=1; + + if ( type == POINTTYPE || type == MULTIPOINTTYPE ) + PG_RETURN_POINTER(geom); + + if ( (PG_NARGS()>1) && (!PG_ARGISNULL(1)) ) + n_iterations = PG_GETARG_INT32(1); + + if (n_iterations>5) + elog(ERROR,"Not more than 5 iterations please"); + if (n_iterations< 1) + elog(ERROR,"Number of iterations must be between 1 and 5"); + + if ( (PG_NARGS()>2) && (!PG_ARGISNULL(2)) ) + { + if(PG_GETARG_BOOL(2)) + preserve_endpoints = 1; + else + preserve_endpoints = 0; + } + + in = lwgeom_from_gserialized(geom); + + out = lwgeom_chaikin(in, n_iterations, preserve_endpoints); + if ( ! out ) PG_RETURN_NULL(); + + /* COMPUTE_BBOX TAINTING */ + if ( in->bbox ) lwgeom_add_bbox(out); + + result = geometry_serialize(out); + lwgeom_free(out); + PG_FREE_IF_COPY(geom, 0); + PG_RETURN_POINTER(result); +} + /*********************************************************************** * --strk@kbt.io; diff --git a/postgis/postgis.sql.in b/postgis/postgis.sql.in index ec456ff26..020715dcd 100644 --- a/postgis/postgis.sql.in +++ b/postgis/postgis.sql.in @@ -3144,7 +3144,14 @@ CREATE OR REPLACE FUNCTION ST_SetEffectiveArea(geometry, float8 default -1, int AS 'MODULE_PATHNAME', 'LWGEOM_SetEffectiveArea' LANGUAGE 'c' IMMUTABLE STRICT _PARALLEL COST 1; -- reset cost, see #3675 - + +-- Availability: 2.5.0 +CREATE OR REPLACE FUNCTION ST_ChaikinSmoothing(geometry, integer default 1, boolean default true) + RETURNS geometry + AS 'MODULE_PATHNAME', 'LWGEOM_ChaikinSmoothing' + LANGUAGE 'c' IMMUTABLE STRICT _PARALLEL + COST 1; -- reset cost, see #3675 + -- ST_SnapToGrid(input, xoff, yoff, xsize, ysize) -- Availability: 1.2.2 CREATE OR REPLACE FUNCTION ST_SnapToGrid(geometry, float8, float8, float8, float8) diff --git a/regress/Makefile.in b/regress/Makefile.in index 0e76194e1..e5a71399e 100644 --- a/regress/Makefile.in +++ b/regress/Makefile.in @@ -80,6 +80,7 @@ TESTS = \ bestsrid \ binary \ boundary \ + chaikin \ cluster \ concave_hull\ ctors \ diff --git a/regress/chaikin.sql b/regress/chaikin.sql new file mode 100644 index 000000000..f799c71fb --- /dev/null +++ b/regress/chaikin.sql @@ -0,0 +1,6 @@ +SELECT '1', ST_astext(ST_ChaikinSmoothing('LINESTRING(0 0, 8 8, 0 16)')); +SELECT '2', ST_astext(ST_ChaikinSmoothing('LINESTRING(0 0, 8 8, 0 16)',10)); +SELECT '3', ST_astext(ST_ChaikinSmoothing('LINESTRING(0 0, 8 8, 0 16)',0)); +SELECT '4', ST_astext(ST_ChaikinSmoothing('LINESTRING(0 0, 8 8, 0 16)',2)); +SELECT '5', ST_astext(ST_ChaikinSmoothing('POINT(0 0)')); +SELECT '6', ST_astext(ST_ChaikinSmoothing('GEOMETRYCOLLECTION(POINT(1 1), LINESTRING(1 1, 1 3, 1 5), POLYGON((5 5, 5 10, 10 10, 10 5, 5 5), (6 6, 6 7, 7 7, 7 6, 6 6 )))', 2, 't')); diff --git a/regress/chaikin_expected b/regress/chaikin_expected new file mode 100644 index 000000000..737d6b7bf --- /dev/null +++ b/regress/chaikin_expected @@ -0,0 +1,6 @@ +1|LINESTRING(0 0,6 6,6 10,0 16) +ERROR: Not more than 5 iterations please +ERROR: Number of iterations must be between 1 and 5 +4|LINESTRING(0 0,4.5 4.5,6 7,6 9,4.5 11.5,0 16) +5|POINT(0 0) +6|GEOMETRYCOLLECTION(POINT(1 1),LINESTRING(1 1,1 2.125,1 2.75,1 3.25,1 3.875,1 5),POLYGON((5 5,5 7.8125,5.3125 9.0625,5.9375 9.6875,6.875 10,8.125 10,9.0625 9.6875,9.6875 9.0625,10 8.125,10 6.875,9.6875 5.9375,9.0625 5.3125,7.8125 5,5 5),(6 6,6 6.5625,6.0625 6.8125,6.1875 6.9375,6.375 7,6.625 7,6.8125 6.9375,6.9375 6.8125,7 6.625,7 6.375,6.9375 6.1875,6.8125 6.0625,6.5625 6,6 6)))