]> granicus.if.org Git - postgis/commitdiff
Add ChaikinSmoothing #4050
authorNicklas Avén <nicklas.aven@jordogskog.no>
Tue, 20 Mar 2018 13:11:17 +0000 (13:11 +0000)
committerNicklas Avén <nicklas.aven@jordogskog.no>
Tue, 20 Mar 2018 13:11:17 +0000 (13:11 +0000)
git-svn-id: http://svn.osgeo.org/postgis/trunk@16479 b70326c6-7e19-0410-871a-916f4a2858ee

13 files changed:
NEWS
doc/reference_processing.xml
liblwgeom/Makefile.in
liblwgeom/cunit/Makefile.in
liblwgeom/cunit/cu_chaikin.c [new file with mode: 0644]
liblwgeom/cunit/cu_tester.c
liblwgeom/liblwgeom.h.in
liblwgeom/lwchaikins.c [new file with mode: 0644]
postgis/lwgeom_functions_analytic.c
postgis/postgis.sql.in
regress/Makefile.in
regress/chaikin.sql [new file with mode: 0644]
regress/chaikin_expected [new file with mode: 0644]

diff --git a/NEWS b/NEWS
index 396f982e9d8244f772887e25d58a5899af6aa8b3..14b2f08e619ea9efc4f8b36f4c914d0fbc81a4db 100644 (file)
--- 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)
index 07416df8c438f741ebd111ecf4ea0acff9610aff..86cf8327b41510450e57340e2b98eea547f7f619 100644 (file)
@@ -3081,7 +3081,7 @@ FROM (SELECT ST_Buffer('POINT(1 3)', 10,12) As the_geom) As foo;
                  </refsection>
        </refentry>
 
-<refentry id="ST_SimplifyVW">
+    <refentry id="ST_SimplifyVW">
          <refnamediv>
                <refname>ST_SimplifyVW</refname>
                <refpurpose>Returns a "simplified" version of the given geometry using the Visvalingam-Whyatt algorithm</refpurpose>
@@ -3129,7 +3129,63 @@ LINESTRING(5 2,7 25,10 10)
                        <para><xref linkend="ST_SetEffectiveArea" />, <xref linkend="ST_Simplify" />, <xref linkend="ST_SimplifyPreserveTopology" />, Topology <xref linkend="TP_ST_Simplify"/></para>
                  </refsection>
        </refentry>
-               <refentry id="ST_SetEffectiveArea">
+
+    <refentry id="ST_ChaikinSmoothing">
+         <refnamediv>
+               <refname>ST_ChaikinSmoothing</refname>
+               <refpurpose>Returns a "smoothed" version of the given geometry using the Chaikin algorithm</refpurpose>
+         </refnamediv>
+
+         <refsynopsisdiv>
+               <funcsynopsis>
+                 <funcprototype>
+                       <funcdef>geometry <function>ST_ChaikinSmoothing</function></funcdef>
+                       <paramdef><type>geometry</type> <parameter>geomA</parameter></paramdef>
+                       <paramdef><type>integer</type> <parameter>nIterations</parameter></paramdef>
+                       <paramdef><type>boolean</type> <parameter>preserveEndPoints</parameter></paramdef>
+                 </funcprototype>
+               </funcsynopsis>
+         </refsynopsisdiv>
+
+         <refsection>
+               <title>Description</title>
+               <para> Returns a "smoothed" version of the given geometry using the Chaikin algorithm.
+        See <ulink url="http://www.idav.ucdavis.edu/education/CAGDNotes/Chaikins-Algorithm/Chaikins-Algorithm.html">Chaikins-Algorithm</ulink> 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.</para>
+
+               <note><para>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 <xref linkend="ST_Simplify" /> and <xref linkend="ST_SimplifyVW" />)</para></note>
+               <note><para>Second argument, number of iterations is limited to max 5 iterations</para></note>
+               <note><para>Note third argument is only valid for polygons, and will be ignored for linestrings</para></note>
+               <note><para>This function handles 3D and the third dimension will affect the result.</para></note>
+               <para>Availability: 2.5.0</para>
+         </refsection>
+
+                 <refsection>
+                       <title>Examples</title>
+                       <para>A triangle is smoothed</para>
+                               <programlisting>
+
+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)
+
+                               </programlisting>
+                 </refsection>
+                 <refsection>
+                       <title>See Also</title>
+                       <para><xref linkend="ST_Simplify" />, <xref linkend="ST_SimplifyVW" /></para>
+                 </refsection>
+       </refentry>
+    
+    <refentry id="ST_SetEffectiveArea">
          <refnamediv>
                <refname>ST_SetEffectiveArea</refname>
                <refpurpose>
index d14dc6b3622a21db07fa85d94dd5014679e176ab..d2649f15f802a38a35ba1a04cee47de0c58cf23b 100644 (file)
@@ -114,6 +114,7 @@ SA_OBJS = \
        lwgeom_wrapx.o \
        lwunionfind.o \
        effectivearea.o \
+       lwchaikins.o \
        lwkmeans.o \
        varint.o
 
index 8ff18ad06a12b54536ac84376d3a332fee700076..10553ec8a369eb9292b3a7aa415c0aedd136bcdc 100644 (file)
@@ -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 (file)
index 0000000..b5d5224
--- /dev/null
@@ -0,0 +1,66 @@
+/**********************************************************************
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.net
+ *
+ * Copyright 2012 Sandro Santilli <strk@kbt.io>
+ *
+ * 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 <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#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);
+}
index 28d4df9c2c4fce66c9b2ab018727db77d54cb10c..f44bb78769f14bcb6db11f152d5a79fcd95f5bb9 100644 (file)
@@ -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,
index 290908e2f767fd5318df715fb4a37eb1bf9b383d..4b54cc41a79edfd57e8a868d8619f6e7b1756b56 100644 (file)
@@ -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 (file)
index 0000000..b69c0dd
--- /dev/null
@@ -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 <http://www.gnu.org/licenses/>.
+ *
+ **********************************************************************
+ *
+ * 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;p<inpts->npoints;p++)
+       {
+               memcpy(&p2, &p1, sizeof(POINT4D));
+               p1 = getPoint4d(inpts, p);
+               if(p>0)
+               {
+                       p1_set = p2_set = 0;
+                       for (i=0;i<ndims;i++)
+                       {
+                               val1 = ((double*) &(p1))[i];
+                               val2 = ((double*) &(p2))[i];
+                               deltaval =  val1 - val2;
+                               quarter_delta = deltaval * 0.25;
+                               if(!preserve_endpoint || p > 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;j<n_iterations;j++)
+       {
+               pa_new = ptarray_chaikin(pa,1,LW_FALSE);
+               if(j>0)
+                       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;j<n_iterations;j++)
+               {
+                       pa_new = ptarray_chaikin(pa,preserve_endpoint,LW_TRUE);
+                       if(j>0)
+                               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;
+}
index 5e18473f4387821b2c52bf9f40019ffa02d2ec8c..18dde2a1155adf0f33ce309584beb15389fcd022 100644 (file)
@@ -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;
index ec456ff26b0de5ab7367d0251b3e7c320bfff646..020715dcd5724ab480cccee421bf3d109a2a94a9 100644 (file)
@@ -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)
index 0e76194e16165c3df5d4b527fd330baffc854b71..e5a71399e76d57e06dd59e75f9d27f3d863d4539 100644 (file)
@@ -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 (file)
index 0000000..f799c71
--- /dev/null
@@ -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 (file)
index 0000000..737d6b7
--- /dev/null
@@ -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)))