[SCM] PostGIS branch master updated. 3.6.0rc2-474-g0e794aa47
git at osgeo.org
git at osgeo.org
Fri May 1 10:42:26 PDT 2026
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "PostGIS".
The branch, master has been updated
via 0e794aa474544d0da986d4ba7a807103b97c129a (commit)
from 26012786958aedb55589c0696086f7c48fc9fc46 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 0e794aa474544d0da986d4ba7a807103b97c129a
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Fri May 1 10:41:07 2026 -0700
ST_CatmullSmoothing(geom, nSegments) "smoothes" a LineString or
LinearRing by adding vertices, and leaves all the input vertices
unchanged. The nSegments parameter controls how many vertices are
added between each input pair.
Closes #4398
diff --git a/NEWS b/NEWS
index 1a59e317b..b66809e9b 100644
--- a/NEWS
+++ b/NEWS
@@ -36,6 +36,7 @@ This version requires GEOS 3.10 or higher
- #5889, [topology] Include representative locations in topology build errors (Darafei Praliaskouski)
- #2614, Use GEOSPreparedDistance and caching to accelerate ST_DWithin (Paul Ramsey)
- GH-839, ST_Multi support for TIN and surfaces (Loïc Bartoletti)
+ - #4398, Add ST_CatmullSmoothing smoothes but retains original vertices (Paul Ramsey)
* Enhancements *
diff --git a/doc/html/images/Makefile.in b/doc/html/images/Makefile.in
index 203142611..288065530 100644
--- a/doc/html/images/Makefile.in
+++ b/doc/html/images/Makefile.in
@@ -85,6 +85,10 @@ GENERATED_IMAGES= \
st_centroid02.png \
st_centroid03.png \
st_centroid04.png \
+ st_catmullromsmoothing01.png \
+ st_catmullromsmoothing02.png \
+ st_catmullromsmoothing03.png \
+ st_catmullromsmoothing04.png \
st_chaikinsmoothing01.png \
st_chaikinsmoothing02.png \
st_chaikinsmoothing03.png \
diff --git a/doc/html/images/wkt/st_catmullromsmoothing01.wkt b/doc/html/images/wkt/st_catmullromsmoothing01.wkt
new file mode 100644
index 000000000..a2ecc12a5
--- /dev/null
+++ b/doc/html/images/wkt/st_catmullromsmoothing01.wkt
@@ -0,0 +1,2 @@
+ArgA-thin;POLYGON ((20 20, 60 90, 10 150, 100 190, 190 160, 130 120, 190 50, 140 70, 120 10, 90 60, 20 20))
+Result;POLYGON((20 20,22.4 27.12,32.4 40.56,45.2 57.44,56 74.88,60 90,53.52 102.96,39.76 115.68,24.24 127.92,12.48 139.44,10 150,19.04 160.4,35.92 170.8,57.28 180,79.76 186.8,100 190,120.4 188.64,143.2 183.52,164.8 176.08,181.6 167.76,190 160,185.68 153.12,171.04 146.16,152.56 138.64,136.72 130.08,130 120,136.08 106.48,150.64 89.84,168.16 72.96,183.12 58.72,190 50,186.56 49.52,176.48 55.36,163.12 63.44,149.84 69.68,140 70,134.24 61.36,130.32 46.48,127.28 29.92,124.16 16.24,120 10,115.28 14.4,110.64 26.4,105.36 41.2,98.72 54,90 60,76.8 56,59.6 45.2,42 32.4,27.6 22.4,20 20))
diff --git a/doc/html/images/wkt/st_catmullromsmoothing02.wkt b/doc/html/images/wkt/st_catmullromsmoothing02.wkt
new file mode 100644
index 000000000..eeeed2a86
--- /dev/null
+++ b/doc/html/images/wkt/st_catmullromsmoothing02.wkt
@@ -0,0 +1,2 @@
+ArgA-thin;POLYGON ((20 20, 60 90, 10 150, 100 190, 190 160, 130 120, 190 50, 140 70, 120 10, 90 60, 20 20))
+Result;POLYGON((20 20,19.95 22.59,22.4 27.12,26.75 33.23,32.4 40.56,38.75 48.75,45.2 57.44,51.15 66.27,56 74.88,59.15 82.91,60 90,58.015 96.495,53.52 102.96,47.205 109.365,39.76 115.68,31.875 121.875,24.24 127.92,17.545 133.785,12.48 139.44,9.735 144.855,10 150,13.33 155.125,19.04 160.4,26.71 165.675,35.92 170.8,46.25 175.625,57.28 180,68.59 183.775,79.76 186.8,90.37 188.925,100 190,109.675 189.88,120.4 188.64,131.725 186.46,143.2 183.52,154.375 180,164.8 176.08,174.025 171.94,181.6 167.76,187.075 163.72,190 160,189.535 156.54,185.68 153.12,179.245 149.68,171.04 146.16,161.875 142.5,152.56 138.64,143.905 134.52,136.72 130.08,131.815 125.26,130 120,131.635 113.81,136.08 106.48,142.645 98.37,150.64 89.84,159.375 81.25,168.16 72.96,176.305 65.33,183.12 58.72,187.915 53.49,190 50,189.32 48.715,186.56 49.52,182.14 51.905,176.48 55.36,170 59.375,163.12 63.44,156.26 67.045,149.84 69.68,144.28 70.835,140 70,136.83 66.745,134.24 61.36,132.11 54.415,130.32 46.48,128.75 38.125,127.28 29.92,125
.79 22.435,124.16 16.24,122.27 11.905,120 10,117.585 10.95,115.28 14.4,112.995 19.75,110.64 26.4,108.125 33.75,105.36 41.2,102.255 48.15,98.72 54,94.665 58.15,90 60,84.125 59.15,76.8 56,68.475 51.15,59.6 45.2,50.625 38.75,42 32.4,34.175 26.75,27.6 22.4,22.725 19.95,20 20))
diff --git a/doc/html/images/wkt/st_catmullromsmoothing03.wkt b/doc/html/images/wkt/st_catmullromsmoothing03.wkt
new file mode 100644
index 000000000..53817d248
--- /dev/null
+++ b/doc/html/images/wkt/st_catmullromsmoothing03.wkt
@@ -0,0 +1,2 @@
+ArgA-thin;LINESTRING (10 140, 80 130, 100 190, 190 150, 140 20, 120 120, 50 30, 30 100)
+Result;LINESTRING(10 140,24.8 136.88,40.4 132.64,55.6 128.96,69.2 127.52,80 130,86.08 139.12,88.24 153.76,89.36 169.84,92.32 183.28,100 190,115.76 189.84,137.68 185.52,160.72 177.28,179.84 165.36,190 150,188.48 126.08,178.64 93.44,164.56 59.76,150.32 32.72,140 20,134.88 28.32,132.24 52.56,130.16 82.64,126.72 108.48,120 120,108.4 111.6,93.2 90,76.8 63.6,61.6 40.8,50 30,42.8 33.76,38.4 46.48,35.6 64.32,33.2 83.44,30 100)
diff --git a/doc/html/images/wkt/st_catmullromsmoothing04.wkt b/doc/html/images/wkt/st_catmullromsmoothing04.wkt
new file mode 100644
index 000000000..1bfb3792f
--- /dev/null
+++ b/doc/html/images/wkt/st_catmullromsmoothing04.wkt
@@ -0,0 +1,2 @@
+ArgA-thin;LINESTRING (10 140, 80 130, 100 190, 190 150, 140 20, 120 120, 50 30, 30 100)
+Result;LINESTRING(10 140,17.225 138.685,24.8 136.88,32.575 134.795,40.4 132.64,48.125 130.625,55.6 128.96,62.675 127.855,69.2 127.52,75.025 128.165,80 130,83.71 133.615,86.08 139.12,87.47 146.005,88.24 153.76,88.75 161.875,89.36 169.84,90.43 177.145,92.32 183.28,95.39 187.735,100 190,106.795 190.455,115.76 189.84,126.265 188.185,137.68 185.52,149.375 181.875,160.72 177.28,171.085 171.765,179.84 165.36,186.355 158.095,190 150,190.535 139.61,188.48 126.08,184.345 110.37,178.64 93.44,171.875 76.25,164.56 59.76,157.205 44.93,150.32 32.72,144.415 24.09,140 20,137.01 21.54,134.88 28.32,133.37 39.08,132.24 52.56,131.25 67.5,130.16 82.64,128.73 96.72,126.72 108.48,123.89 116.66,120 120,114.8 117.975,108.4 111.6,101.1 101.925,93.2 90,85 76.875,76.8 63.6,68.9 51.225,61.6 40.8,55.2 33.375,50 30,45.975 30.52,42.8 33.76,40.325 39.24,38.4 46.48,36.875 55,35.6 64.32,34.425 73.96,33.2 83.44,31.775 92.28,30 100)
diff --git a/doc/reference_processing.xml b/doc/reference_processing.xml
index f4d7100ba..2f369b2ed 100644
--- a/doc/reference_processing.xml
+++ b/doc/reference_processing.xml
@@ -794,6 +794,140 @@ SELECT ST_ChaikinSmoothing(
</refentry>
+ <refentry xml:id="ST_CatmullRomSmoothing">
+ <refnamediv>
+ <refname>ST_CatmullRomSmoothing</refname>
+ <refpurpose>Returns a smoothed version of a geometry, using the Catmull-Rom spline algorithm</refpurpose>
+ </refnamediv>
+
+ <refsynopsisdiv>
+ <funcsynopsis>
+ <funcprototype>
+ <funcdef>geometry <function>ST_CatmullRomSmoothing</function></funcdef>
+ <paramdef><type>geometry</type> <parameter>geom</parameter></paramdef>
+ <paramdef choice="opt"><type>integer</type> <parameter>nSegments = 5</parameter></paramdef>
+ </funcprototype>
+ </funcsynopsis>
+ </refsynopsisdiv>
+
+ <refsection>
+ <title>Description</title>
+ <para>Smoothes a linear or polygonal geometry using the
+ <link xlink:href="https://en.wikipedia.org/wiki/Centripetal_Catmull%E2%80%93Rom_spline">Catmull-Rom spline</link>
+ algorithm. Unlike <xref linkend="ST_ChaikinSmoothing"/>, this is an <emphasis>interpolating</emphasis>
+ spline: the output curve passes through every original vertex. Between each pair of
+ consecutive original vertices, <varname>nSegments</varname> - 1 new intermediate
+ points are inserted.
+ </para>
+
+ <para>The <varname>nSegments</varname> parameter controls the density of the output.
+ With <varname>nSegments</varname> = 5 (the default), each original span is
+ divided into 5 sub-segments (inserting 4 new points per span). The minimum value is 2.
+ </para>
+
+ <para>At least 4 input vertices are required to apply smoothing;
+ geometries with fewer vertices are returned unchanged.
+ Points and multipoints are always returned unchanged.
+ </para>
+
+ <note><para>The output vertex count grows as
+ <literal>1 + (N-1) * nSegments</literal> for open lines,
+ and <literal>1 + N * nSegments</literal> for closed rings,
+ where N is the number of input vertices.
+ For large geometries, use a simplification function on the result
+ to reduce the number of points
+ (see <xref linkend="ST_Simplify"/>, <xref linkend="ST_SimplifyPreserveTopology"/>
+ and <xref linkend="ST_SimplifyVW"/>).
+ </para></note>
+
+ <para>The result has interpolated values for the Z and M dimensions when present.</para>
+
+ <para>&Z_support;</para>
+ <para role="availability" conformance="3.6.0">Availability: 3.6.0</para>
+ </refsection>
+
+ <refsection>
+ <title>Examples</title>
+ <para>Smoothing a 4-point collinear line with default nSegments=5:</para>
+ <programlisting>
+SELECT ST_AsText(ST_CatmullRomSmoothing('LINESTRING(0 0, 5 0, 10 0, 15 0)'));
+
+ LINESTRING(0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 0,10 0,11 0,12 0,13 0,14 0,15 0)
+</programlisting>
+
+ <para>Smoothing a Polygon using nSegments = 5 and 10:</para>
+ <informaltable>
+ <tgroup cols="2">
+ <tbody>
+ <row>
+ <entry><para><informalfigure>
+ <mediaobject>
+ <imageobject>
+ <imagedata fileref="images/st_catmullromsmoothing01.png"/>
+ </imageobject>
+ <caption><para>nSegments = 5</para></caption>
+ </mediaobject>
+ </informalfigure></para></entry>
+
+ <entry><para><informalfigure>
+ <mediaobject>
+ <imageobject>
+ <imagedata fileref="images/st_catmullromsmoothing02.png"/>
+ </imageobject>
+ <caption><para>nSegments = 10</para></caption>
+ </mediaobject>
+ </informalfigure></para></entry>
+ </row>
+ </tbody>
+ </tgroup>
+ </informaltable>
+ <programlisting>
+SELECT ST_CatmullRomSmoothing(
+ 'POLYGON ((20 20, 60 90, 10 150, 100 190, 190 160, 130 120, 190 50, 140 70, 120 10, 90 60, 20 20))',
+ generate_series(5, 10, 5) );
+</programlisting>
+
+ <para>Smoothing a LineString using nSegments = 5 and 10:</para>
+ <informaltable>
+ <tgroup cols="2">
+ <tbody>
+ <row>
+ <entry><para><informalfigure>
+ <mediaobject>
+ <imageobject>
+ <imagedata fileref="images/st_catmullromsmoothing03.png"/>
+ </imageobject>
+ <caption><para>nSegments = 5</para></caption>
+ </mediaobject>
+ </informalfigure></para></entry>
+
+ <entry><para><informalfigure>
+ <mediaobject>
+ <imageobject>
+ <imagedata fileref="images/st_catmullromsmoothing04.png"/>
+ </imageobject>
+ <caption><para>nSegments = 10</para></caption>
+ </mediaobject>
+ </informalfigure></para></entry>
+ </row>
+ </tbody>
+ </tgroup>
+ </informaltable>
+ <programlisting>
+SELECT ST_CatmullRomSmoothing(
+ 'LINESTRING (10 140, 80 130, 100 190, 190 150, 140 20, 120 120, 50 30, 30 100)',
+ generate_series(5, 10, 5) );
+</programlisting>
+ </refsection>
+
+ <refsection>
+ <title>See Also</title>
+ <para><xref linkend="ST_ChaikinSmoothing"/>, <xref linkend="ST_Simplify"/>,
+ <xref linkend="ST_SimplifyPreserveTopology"/>, <xref linkend="ST_SimplifyVW"/></para>
+ </refsection>
+ </refentry>
+
+
<refentry xml:id="ST_ConcaveHull">
<refnamediv>
<refname>ST_ConcaveHull</refname>
diff --git a/liblwgeom/Makefile.in b/liblwgeom/Makefile.in
index caed0abdb..f74c2d9a9 100644
--- a/liblwgeom/Makefile.in
+++ b/liblwgeom/Makefile.in
@@ -131,7 +131,7 @@ SA_OBJS = \
lwgeom_wrapx.o \
lwunionfind.o \
effectivearea.o \
- lwchaikins.o \
+ lwsmoothing.o \
lwmval.o \
lwkmeans.o \
varint.o \
diff --git a/liblwgeom/cunit/Makefile.in b/liblwgeom/cunit/Makefile.in
index f397414c9..791be060e 100644
--- a/liblwgeom/cunit/Makefile.in
+++ b/liblwgeom/cunit/Makefile.in
@@ -48,6 +48,7 @@ OBJS= \
cu_measures.o \
cu_effectivearea.o \
cu_chaikin.o \
+ cu_catmullrom.o \
cu_filterm.o \
cu_node.o \
cu_clip_by_rect.o \
diff --git a/liblwgeom/cunit/cu_catmullrom.c b/liblwgeom/cunit/cu_catmullrom.c
new file mode 100644
index 000000000..1d9ff4e0c
--- /dev/null
+++ b/liblwgeom/cunit/cu_catmullrom.c
@@ -0,0 +1,201 @@
+/**********************************************************************
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.net
+ *
+ * Copyright 2026 PostGIS contributors
+ *
+ * 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_catmull_rom(char *geom_txt, char *expected, int n_segments)
+{
+ LWGEOM *geom_in, *geom_out;
+ char *out_txt;
+ geom_in = lwgeom_from_wkt(geom_txt, LW_PARSER_CHECK_NONE);
+ geom_out = lwgeom_catmull_rom(geom_in, n_segments);
+ out_txt = lwgeom_to_wkt(geom_out, WKT_EXTENDED, 6, NULL);
+ if (strcmp(expected, out_txt))
+ printf("\nExpected: %s\n Got: %s\n", expected, out_txt);
+ ASSERT_STRING_EQUAL(expected, out_txt);
+ lwfree(out_txt);
+ lwgeom_free(geom_in);
+ lwgeom_free(geom_out);
+}
+
+/*
+ * A collinear 4-point line: Catmull-Rom on a straight line must produce
+ * a straight line with original vertices preserved at their positions.
+ * With n_segments=2 we get 1 + 3*2 = 7 output points.
+ */
+static void
+do_test_catmull_rom_lines(void)
+{
+ /* Collinear: output must be straight, original vertices preserved */
+ do_test_catmull_rom(
+ "LINESTRING(0 0, 1 0, 2 0, 3 0)",
+ "LINESTRING(0 0,0.5 0,1 0,1.5 0,2 0,2.5 0,3 0)",
+ 2);
+
+ /* Fewer than 4 points: return input unchanged */
+ do_test_catmull_rom(
+ "LINESTRING(0 0, 1 1, 2 0)",
+ "LINESTRING(0 0,1 1,2 0)",
+ 5);
+
+ /* Exactly 2 points: return input unchanged */
+ do_test_catmull_rom(
+ "LINESTRING(0 0, 10 10)",
+ "LINESTRING(0 0,10 10)",
+ 5);
+
+ /* Empty line: return empty */
+ do_test_catmull_rom(
+ "LINESTRING EMPTY",
+ "LINESTRING EMPTY",
+ 5);
+}
+
+/*
+ * Verify that original vertices are preserved by checking a symmetric
+ * curve where we can identify exactly where original points fall.
+ * LINESTRING(0 0, 0 10, 10 10, 10 0) with n_segments=1 is trivial
+ * (n_segments=1 is invalid per SQL, but at the C level we test =2).
+ */
+static void
+do_test_catmull_rom_vertex_preservation(void)
+{
+ LWGEOM *geom_in, *geom_out;
+ POINT4D pt;
+
+ /* With n_segments=3 on a 4-point line, output has 1+3*3=10 points.
+ * Original vertices must appear at indices 0, 3, 6, 9. */
+ geom_in = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 2 0, 3 0)", LW_PARSER_CHECK_NONE);
+ geom_out = lwgeom_catmull_rom(geom_in, 3);
+
+ LWLINE *oline = (LWLINE *)geom_out;
+ CU_ASSERT_EQUAL(oline->points->npoints, 10);
+
+ /* Check original vertices at expected positions */
+ pt = getPoint4d(oline->points, 0);
+ CU_ASSERT_DOUBLE_EQUAL(pt.x, 0.0, 1e-10);
+ CU_ASSERT_DOUBLE_EQUAL(pt.y, 0.0, 1e-10);
+
+ pt = getPoint4d(oline->points, 3);
+ CU_ASSERT_DOUBLE_EQUAL(pt.x, 1.0, 1e-10);
+ CU_ASSERT_DOUBLE_EQUAL(pt.y, 0.0, 1e-10);
+
+ pt = getPoint4d(oline->points, 6);
+ CU_ASSERT_DOUBLE_EQUAL(pt.x, 2.0, 1e-10);
+ CU_ASSERT_DOUBLE_EQUAL(pt.y, 0.0, 1e-10);
+
+ pt = getPoint4d(oline->points, 9);
+ CU_ASSERT_DOUBLE_EQUAL(pt.x, 3.0, 1e-10);
+ CU_ASSERT_DOUBLE_EQUAL(pt.y, 0.0, 1e-10);
+
+ lwgeom_free(geom_in);
+ lwgeom_free(geom_out);
+}
+
+/*
+ * Polygon ring smoothing: a 5-point square ring (4 unique vertices + closing)
+ * is the minimum for Catmull-Rom on a ring (needs npoints-1 >= 4).
+ */
+static void
+do_test_catmull_rom_polygons(void)
+{
+ LWGEOM *geom_in, *geom_out;
+ LWPOLY *opoly;
+
+ /* Polygon with 5 ring points (4 unique + close): smooth with n_segments=2.
+ * Each ring span produces 2 sub-segments, 4 spans total → 4*2+1 = 9 ring pts. */
+ geom_in = lwgeom_from_wkt(
+ "POLYGON((0 0, 4 0, 4 4, 0 4, 0 0))",
+ LW_PARSER_CHECK_NONE);
+ geom_out = lwgeom_catmull_rom(geom_in, 2);
+
+ opoly = (LWPOLY *)geom_out;
+ CU_ASSERT_EQUAL(opoly->nrings, 1);
+ CU_ASSERT_EQUAL(opoly->rings[0]->npoints, 9);
+
+ /* First and last ring points must be identical (closed) */
+ POINT4D first = getPoint4d(opoly->rings[0], 0);
+ POINT4D last = getPoint4d(opoly->rings[0], opoly->rings[0]->npoints - 1);
+ CU_ASSERT_DOUBLE_EQUAL(first.x, last.x, 1e-10);
+ CU_ASSERT_DOUBLE_EQUAL(first.y, last.y, 1e-10);
+
+ lwgeom_free(geom_in);
+ lwgeom_free(geom_out);
+
+ /* Triangle ring (npoints-1 = 3 < 4): return original ring unchanged */
+ geom_in = lwgeom_from_wkt(
+ "POLYGON((0 0, 4 0, 2 4, 0 0))",
+ LW_PARSER_CHECK_NONE);
+ geom_out = lwgeom_catmull_rom(geom_in, 5);
+ opoly = (LWPOLY *)geom_out;
+ CU_ASSERT_EQUAL(opoly->rings[0]->npoints, 4); /* unchanged */
+ lwgeom_free(geom_in);
+ lwgeom_free(geom_out);
+}
+
+/*
+ * Z and M dimensions must be interpolated by the algorithm.
+ */
+static void
+do_test_catmull_rom_zm(void)
+{
+ LWGEOM *geom_in, *geom_out;
+ LWLINE *oline;
+ POINT4D pt;
+
+ /* Collinear in XYZ: Z should be interpolated linearly on a straight line */
+ geom_in = lwgeom_from_wkt(
+ "LINESTRING Z (0 0 0, 1 0 10, 2 0 20, 3 0 30)",
+ LW_PARSER_CHECK_NONE);
+ geom_out = lwgeom_catmull_rom(geom_in, 2);
+ oline = (LWLINE *)geom_out;
+
+ /* Z at midpoint of first span (index 1) should be ~5 */
+ pt = getPoint4d(oline->points, 1);
+ CU_ASSERT_DOUBLE_EQUAL(pt.z, 5.0, 1e-10);
+
+ /* Original Z values preserved at span boundaries */
+ pt = getPoint4d(oline->points, 2);
+ CU_ASSERT_DOUBLE_EQUAL(pt.z, 10.0, 1e-10);
+
+ lwgeom_free(geom_in);
+ lwgeom_free(geom_out);
+}
+
+/*
+ * Points pass through unchanged.
+ */
+static void
+do_test_catmull_rom_points(void)
+{
+ do_test_catmull_rom("POINT(1 2)", "POINT(1 2)", 5);
+ do_test_catmull_rom("MULTIPOINT((1 2),(3 4))", "MULTIPOINT(1 2,3 4)", 5);
+}
+
+
+void catmullrom_suite_setup(void);
+void catmullrom_suite_setup(void)
+{
+ CU_pSuite suite = CU_add_suite("catmullrom", NULL, NULL);
+ PG_ADD_TEST(suite, do_test_catmull_rom_lines);
+ PG_ADD_TEST(suite, do_test_catmull_rom_vertex_preservation);
+ PG_ADD_TEST(suite, do_test_catmull_rom_polygons);
+ PG_ADD_TEST(suite, do_test_catmull_rom_zm);
+ PG_ADD_TEST(suite, do_test_catmull_rom_points);
+}
diff --git a/liblwgeom/cunit/cu_tester.c b/liblwgeom/cunit/cu_tester.c
index 708b8070c..a32c3f41a 100644
--- a/liblwgeom/cunit/cu_tester.c
+++ b/liblwgeom/cunit/cu_tester.c
@@ -53,6 +53,7 @@ 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 catmullrom_suite_setup(void);
extern void filterm_suite_setup(void);
extern void minimum_bounding_circle_suite_setup(void);
extern void misc_suite_setup(void);
@@ -107,6 +108,7 @@ PG_SuiteSetup setupfuncs[] = {algorithms_suite_setup,
measures_suite_setup,
effectivearea_suite_setup,
chaikin_suite_setup,
+ catmullrom_suite_setup,
filterm_suite_setup,
minimum_bounding_circle_suite_setup,
misc_suite_setup,
diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in
index 9ff00e1a4..401ac78cf 100644
--- a/liblwgeom/liblwgeom.h.in
+++ b/liblwgeom/liblwgeom.h.in
@@ -1132,6 +1132,7 @@ extern LWGEOM* lwgeom_force_4d(const LWGEOM *geom, double zval, double mval);
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);
+extern LWGEOM* lwgeom_catmull_rom(const LWGEOM *igeom, int n_segments);
extern LWGEOM* lwgeom_filter_m(LWGEOM *geom, double min, double max, int returnm);
/*
diff --git a/liblwgeom/lwchaikins.c b/liblwgeom/lwchaikins.c
deleted file mode 100644
index 4312a8631..000000000
--- a/liblwgeom/lwchaikins.c
+++ /dev/null
@@ -1,202 +0,0 @@
-
-/**********************************************************************
- *
- * 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;
- uint32_t 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;
- LWLINE *oline;
-
- if( lwline_is_empty(iline))
- return lwline_clone(iline);
-
- 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;
- int 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;
-}
diff --git a/liblwgeom/lwsmoothing.c b/liblwgeom/lwsmoothing.c
new file mode 100644
index 000000000..8d7e256f5
--- /dev/null
+++ b/liblwgeom/lwsmoothing.c
@@ -0,0 +1,432 @@
+
+/**********************************************************************
+ *
+ * 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/>.
+ *
+ **********************************************************************
+ *
+ * Chaikin smoothing: Copyright 2018 Nicklas Avén
+ * Catmull-Rom smoothing: Copyright 2026 PostGIS contributors
+ *
+ **********************************************************************/
+
+
+#include "liblwgeom_internal.h"
+
+
+/* ========================================================================
+ * Chaikin smoothing
+ * ======================================================================== */
+
+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;
+ uint32_t 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;
+ LWLINE *oline;
+
+ if( lwline_is_empty(iline))
+ return lwline_clone(iline);
+
+ 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;
+ int 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;
+}
+
+
+/* Forward declaration needed by lwcollection_catmull_rom */
+LWGEOM *lwgeom_catmull_rom(const LWGEOM *igeom, int n_segments);
+
+/* ========================================================================
+ * Catmull-Rom spline smoothing
+ *
+ * An interpolating spline: the output curve passes through every original
+ * vertex. Between each pair of original vertices, (n_segments - 1) new
+ * intermediate points are inserted.
+ *
+ * Standard uniform Catmull-Rom formula for t in [0,1] between P1 and P2:
+ * q(t) = 0.5 * [ 2*P1
+ * + (-P0 + P2)*t
+ * + (2*P0 - 5*P1 + 4*P2 - P3)*t^2
+ * + (-P0 + 3*P1 - 3*P2 + P3)*t^3 ]
+ *
+ * Applied per dimension (x, y, and optionally z and m).
+ *
+ * Boundary handling for open lines: phantom endpoints are constructed by
+ * reflecting the adjacent point through the endpoint.
+ * P_{-1} = 2*P0 - P1
+ * P_N = 2*P_{N-1} - P_{N-2}
+ *
+ * Boundary handling for closed rings: wrap-around indexing.
+ * ======================================================================== */
+
+/*
+ * Core Catmull-Rom algorithm on a POINTARRAY.
+ *
+ * isclosed: LW_TRUE for polygon rings (first == last point).
+ * n_segments: number of sub-segments per original span (>= 2).
+ *
+ * Returns a newly allocated POINTARRAY; caller must free.
+ */
+static POINTARRAY *
+ptarray_catmull_rom(const POINTARRAY *inpts, int n_segments, int isclosed)
+{
+ uint32_t ndims = 2 + ptarray_has_z(inpts) + ptarray_has_m(inpts);
+ uint32_t npoints = inpts->npoints;
+
+ /* For closed rings the last point duplicates the first; work with
+ * N = npoints-1 unique points so we don't double-count. */
+ uint32_t N = isclosed ? npoints - 1 : npoints;
+
+ /* Too few unique points: return a clone unchanged */
+ if (N < 4)
+ return ptarray_clone(inpts);
+
+ /* Number of spans between consecutive original vertices */
+ uint32_t nspans = N - (isclosed ? 0 : 1);
+
+ /* Output size: nspans * n_segments + 1 (we always close with the last orig point) */
+ uint32_t out_cap = nspans * (uint32_t)n_segments + 1;
+ POINTARRAY *opts = ptarray_construct_empty(
+ FLAGS_GET_Z(inpts->flags), FLAGS_GET_M(inpts->flags), out_cap);
+
+ /* Helper to get a point by index with phantom-endpoint / wrap-around logic.
+ * idx may be -1 or >= N; we handle accordingly. */
+#define GET_COORD(idx, dim) \
+ (((int32_t)(idx) < 0) \
+ ? (isclosed \
+ ? ((double *)getPoint_internal(inpts, \
+ (uint32_t)((int32_t)N + (int32_t)(idx))))[dim] \
+ : 2.0 * ((double *)getPoint_internal(inpts, 0))[dim] \
+ - ((double *)getPoint_internal(inpts, 1))[dim]) \
+ : ((uint32_t)(idx) >= N) \
+ ? (isclosed \
+ ? ((double *)getPoint_internal(inpts, \
+ (uint32_t)(idx) % N))[dim] \
+ : 2.0 * ((double *)getPoint_internal(inpts, N-1))[dim]\
+ - ((double *)getPoint_internal(inpts, N-2))[dim]) \
+ : ((double *)getPoint_internal(inpts, (uint32_t)(idx)))[dim])
+
+ int first_point = 1; /* track whether we should include t=0 */
+
+ for (uint32_t span = 0; span < nspans; span++)
+ {
+ int32_t i0 = (int32_t)span - 1; /* P0 (leading control) */
+ int32_t i1 = (int32_t)span; /* P1 (start of span) */
+ int32_t i2 = (int32_t)span + 1; /* P2 (end of span) */
+ int32_t i3 = (int32_t)span + 2; /* P3 (trailing control)*/
+
+ int k_start = first_point ? 0 : 1;
+ first_point = 0;
+
+ /* Run k from k_start to n_segments inclusive.
+ * At k=0: t=0.0 (original vertex, only for first span)
+ * At k=n_segments: t=1.0 (next original vertex, closes the span) */
+ for (int k = k_start; k <= n_segments; k++)
+ {
+ double t = k / (double)n_segments;
+ double t2 = t * t;
+ double t3 = t2 * t;
+ POINT4D pt;
+
+ for (uint32_t d = 0; d < ndims; d++)
+ {
+ double c0 = GET_COORD(i0, d);
+ double c1 = GET_COORD(i1, d);
+ double c2 = GET_COORD(i2, d);
+ double c3 = GET_COORD(i3, d);
+
+ ((double *)&pt)[d] = 0.5 * (
+ 2.0 * c1
+ + (-c0 + c2) * t
+ + (2.0*c0 - 5.0*c1 + 4.0*c2 - c3) * t2
+ + (-c0 + 3.0*c1 - 3.0*c2 + c3) * t3
+ );
+ }
+ /* Zero out unused dimensions so POINT4D is clean */
+ if (ndims < 3) pt.z = 0.0;
+ if (ndims < 4) pt.m = 0.0;
+
+ ptarray_append_point(opts, &pt, LW_TRUE);
+ }
+ }
+ /* At k=n_segments (t=1.0), the Catmull-Rom formula evaluates to c2 (the span
+ * endpoint), so all original vertices are naturally included in the output,
+ * and closed rings are automatically closed by the last span's t=1 giving P0. */
+
+#undef GET_COORD
+
+ return opts;
+}
+
+
+static LWLINE *
+lwline_catmull_rom(const LWLINE *iline, int n_segments)
+{
+ if (lwline_is_empty(iline))
+ return lwline_clone(iline);
+
+ if (iline->points->npoints < 4)
+ return lwline_clone(iline);
+
+ POINTARRAY *pa = ptarray_catmull_rom(iline->points, n_segments, LW_FALSE);
+ LWLINE *oline = lwline_construct(iline->srid, NULL, pa);
+ oline->type = iline->type;
+ return oline;
+}
+
+
+static LWPOLY *
+lwpoly_catmull_rom(const LWPOLY *ipoly, int n_segments)
+{
+ 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 (uint32_t i = 0; i < ipoly->nrings; i++)
+ {
+ POINTARRAY *ring = ipoly->rings[i];
+ POINTARRAY *pa_out;
+
+ /* npoints-1 unique points; need at least 4 to smooth */
+ if (ring->npoints - 1 < 4)
+ pa_out = ptarray_clone(ring);
+ else
+ pa_out = ptarray_catmull_rom(ring, n_segments, LW_TRUE);
+
+ if (pa_out->npoints >= 4)
+ {
+ if (lwpoly_add_ring(opoly, pa_out) == LW_FAILURE)
+ return NULL;
+ }
+ else
+ {
+ ptarray_free(pa_out);
+ }
+ }
+
+ opoly->type = ipoly->type;
+
+ if (lwpoly_is_empty(opoly))
+ return NULL;
+
+ return opoly;
+}
+
+
+static LWCOLLECTION *
+lwcollection_catmull_rom(const LWCOLLECTION *igeom, int n_segments)
+{
+ 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;
+
+ for (uint32_t i = 0; i < igeom->ngeoms; i++)
+ {
+ LWGEOM *ngeom = lwgeom_catmull_rom(igeom->geoms[i], n_segments);
+ if (ngeom)
+ out = lwcollection_add_lwgeom(out, ngeom);
+ }
+
+ return out;
+}
+
+
+LWGEOM *
+lwgeom_catmull_rom(const LWGEOM *igeom, int n_segments)
+{
+ switch (igeom->type)
+ {
+ case POINTTYPE:
+ case MULTIPOINTTYPE:
+ return lwgeom_clone(igeom);
+ case LINETYPE:
+ return (LWGEOM *)lwline_catmull_rom((LWLINE *)igeom, n_segments);
+ case POLYGONTYPE:
+ return (LWGEOM *)lwpoly_catmull_rom((LWPOLY *)igeom, n_segments);
+ case MULTILINETYPE:
+ case MULTIPOLYGONTYPE:
+ case COLLECTIONTYPE:
+ return (LWGEOM *)lwcollection_catmull_rom((LWCOLLECTION *)igeom, n_segments);
+ default:
+ lwerror("lwgeom_catmull_rom: 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 13e31cf01..6d6ca9e55 100644
--- a/postgis/lwgeom_functions_analytic.c
+++ b/postgis/lwgeom_functions_analytic.c
@@ -163,6 +163,38 @@ Datum LWGEOM_ChaikinSmoothing(PG_FUNCTION_ARGS)
}
+PG_FUNCTION_INFO_V1(LWGEOM_CatmullRomSmoothing);
+Datum LWGEOM_CatmullRomSmoothing(PG_FUNCTION_ARGS)
+{
+ GSERIALIZED *geom = PG_GETARG_GSERIALIZED_P(0);
+ GSERIALIZED *result;
+ int type = gserialized_get_type(geom);
+ LWGEOM *in, *out;
+ int n_segments = 5;
+
+ if (type == POINTTYPE || type == MULTIPOINTTYPE)
+ PG_RETURN_POINTER(geom);
+
+ if ((PG_NARGS() > 1) && (!PG_ARGISNULL(1)))
+ n_segments = PG_GETARG_INT32(1);
+
+ if (n_segments < 2)
+ elog(ERROR, "nSegments must be >= 2: %s", __func__);
+
+ in = lwgeom_from_gserialized(geom);
+ out = lwgeom_catmull_rom(in, n_segments);
+ 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 at kbt.io;
***********************************************************************/
diff --git a/postgis/postgis.sql.in b/postgis/postgis.sql.in
index 649ed3fc5..73d41a465 100644
--- a/postgis/postgis.sql.in
+++ b/postgis/postgis.sql.in
@@ -3600,6 +3600,13 @@ CREATE OR REPLACE FUNCTION ST_ChaikinSmoothing(geometry, integer default 1, bool
LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
_COST_MEDIUM;
+-- Availability: 3.6.0
+CREATE OR REPLACE FUNCTION ST_CatmullRomSmoothing(geometry, integer default 5)
+ RETURNS geometry
+ AS 'MODULE_PATHNAME', 'LWGEOM_CatmullRomSmoothing'
+ LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
+ _COST_MEDIUM;
+
-- 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/core/catmullrom.sql b/regress/core/catmullrom.sql
new file mode 100644
index 000000000..fe4d9d533
--- /dev/null
+++ b/regress/core/catmullrom.sql
@@ -0,0 +1,39 @@
+-- ST_CatmullRomSmoothing regression tests
+
+-- Basic 4-point collinear line, default nSegments=5
+-- Scale by 5 so interpolated values are integers
+SELECT '1', ST_AsText(ST_CatmullRomSmoothing('LINESTRING(0 0, 5 0, 10 0, 15 0)'));
+
+-- Explicit nSegments=2, outputs at half-integer positions
+SELECT '2', ST_AsText(ST_CatmullRomSmoothing('LINESTRING(0 0, 1 0, 2 0, 3 0)', 2));
+
+-- Fewer than 4 points: return original unchanged
+SELECT '3', ST_AsText(ST_CatmullRomSmoothing('LINESTRING(0 0, 1 1, 2 0)'));
+
+-- 2-point line: return original unchanged
+SELECT '4', ST_AsText(ST_CatmullRomSmoothing('LINESTRING(0 0, 10 10)'));
+
+-- Point: pass through unchanged
+SELECT '5', ST_AsText(ST_CatmullRomSmoothing('POINT(0 0)'));
+
+-- Empty line: return empty
+SELECT '6', ST_AsText(ST_CatmullRomSmoothing('LINESTRING EMPTY'));
+
+-- Polygon with square ring, nSegments=2
+SELECT '7', ST_AsText(ST_CatmullRomSmoothing('POLYGON((0 0, 4 0, 4 4, 0 4, 0 0))', 2));
+
+-- GeometryCollection: point passthrough + line smoothing
+SELECT '8', ST_AsText(ST_CatmullRomSmoothing(
+ 'GEOMETRYCOLLECTION(POINT(1 1), LINESTRING(0 0, 1 0, 2 0, 3 0))', 2));
+
+-- Z dimension is preserved and interpolated
+SELECT '9', ST_AsText(ST_CatmullRomSmoothing('LINESTRING Z (0 0 0, 1 0 10, 2 0 20, 3 0 30)', 2));
+
+-- The geometry bbox is updated
+WITH geom AS (
+ SELECT ST_CatmullRomSmoothing('LINESTRING(0 0, 1 0, 2 0, 3 0)', 2) AS g
+)
+SELECT '10', postgis_getbbox(g) AS box FROM geom;
+
+-- Triangle polygon ring has npoints-1=3 unique pts (<4), return original
+SELECT '11', ST_AsText(ST_CatmullRomSmoothing('POLYGON((0 0, 4 0, 2 4, 0 0))', 5));
diff --git a/regress/core/catmullrom_expected b/regress/core/catmullrom_expected
new file mode 100644
index 000000000..ef9dcc837
--- /dev/null
+++ b/regress/core/catmullrom_expected
@@ -0,0 +1,11 @@
+1|LINESTRING(0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 0,10 0,11 0,12 0,13 0,14 0,15 0)
+2|LINESTRING(0 0,0.5 0,1 0,1.5 0,2 0,2.5 0,3 0)
+3|LINESTRING(0 0,1 1,2 0)
+4|LINESTRING(0 0,10 10)
+5|POINT(0 0)
+6|LINESTRING EMPTY
+7|POLYGON((0 0,2 -0.5,4 0,4.5 2,4 4,2 4.5,0 4,-0.5 2,0 0))
+8|GEOMETRYCOLLECTION(POINT(1 1),LINESTRING(0 0,0.5 0,1 0,1.5 0,2 0,2.5 0,3 0))
+9|LINESTRING Z (0 0 0,0.5 0 5,1 0 10,1.5 0 15,2 0 20,2.5 0 25,3 0 30)
+10|BOX(0 0,3 0)
+11|POLYGON((0 0,4 0,2 4,0 0))
diff --git a/regress/core/tests.mk.in b/regress/core/tests.mk.in
index 5e334ee5e..80d136cbc 100644
--- a/regress/core/tests.mk.in
+++ b/regress/core/tests.mk.in
@@ -33,6 +33,7 @@ TESTS += \
$(top_srcdir)/regress/core/bestsrid \
$(top_srcdir)/regress/core/binary \
$(top_srcdir)/regress/core/boundary \
+ $(top_srcdir)/regress/core/catmullrom \
$(top_srcdir)/regress/core/chaikin \
$(top_srcdir)/regress/core/clean \
$(top_srcdir)/regress/core/clipbybox2d \
-----------------------------------------------------------------------
Summary of changes:
NEWS | 1 +
doc/html/images/Makefile.in | 4 +
doc/html/images/wkt/st_catmullromsmoothing01.wkt | 2 +
doc/html/images/wkt/st_catmullromsmoothing02.wkt | 2 +
doc/html/images/wkt/st_catmullromsmoothing03.wkt | 2 +
doc/html/images/wkt/st_catmullromsmoothing04.wkt | 2 +
doc/reference_processing.xml | 134 +++++++
liblwgeom/Makefile.in | 2 +-
liblwgeom/cunit/Makefile.in | 1 +
liblwgeom/cunit/cu_catmullrom.c | 201 +++++++++++
liblwgeom/cunit/cu_tester.c | 2 +
liblwgeom/liblwgeom.h.in | 1 +
liblwgeom/lwchaikins.c | 202 -----------
liblwgeom/lwsmoothing.c | 432 +++++++++++++++++++++++
postgis/lwgeom_functions_analytic.c | 32 ++
postgis/postgis.sql.in | 7 +
regress/core/catmullrom.sql | 39 ++
regress/core/catmullrom_expected | 11 +
regress/core/tests.mk.in | 1 +
19 files changed, 875 insertions(+), 203 deletions(-)
create mode 100644 doc/html/images/wkt/st_catmullromsmoothing01.wkt
create mode 100644 doc/html/images/wkt/st_catmullromsmoothing02.wkt
create mode 100644 doc/html/images/wkt/st_catmullromsmoothing03.wkt
create mode 100644 doc/html/images/wkt/st_catmullromsmoothing04.wkt
create mode 100644 liblwgeom/cunit/cu_catmullrom.c
delete mode 100644 liblwgeom/lwchaikins.c
create mode 100644 liblwgeom/lwsmoothing.c
create mode 100644 regress/core/catmullrom.sql
create mode 100644 regress/core/catmullrom_expected
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list