Add interpolated versions of polyline and polygon smoothing

core/math/2d/geometry/goost_geometry_2d.cpp

@@ -174,14 +174,113 @@ Vector<Point2> GoostGeometry2D::simplify_polyline(const Vector<Point2> &p_polyli
 	return ret;
 }
 
-// Approximate polygon smoothing using Chaikin algorithm.
+// Catmull-Rom interpolation. See also:
+//
+// "On the Parameterization of Catmull-Rom Curves" by Cem Yuksel, Scott Schaefer, John Keyser.
+// https://people.engr.tamu.edu/schaefer/research/catmull_rom.pdf
+//
+static Vector2 catmull_rom(const Vector2 &p0, const Vector2 &p1, const Vector2 &p2, const Vector2 &p3, float p_t, float p_alpha) {
+	Vector2 c;
+	if (p_alpha > 0.0f) {
+		// Centripetal (alpha == 0.5) or chordal (alpha > 0.5).
+#ifdef DEBUG_ENABLED
+		// Division by zero...
+		ERR_FAIL_COND_V_MSG(p0 == p1 || p1 == p2 || p2 == p3, Vector2(), "Duplicate points detected, cannot interpolate.");
+#endif
+		auto compute_t = [&](float t, float alpha, const Vector2 &v0, const Vector2 &v1) {
+			real_t a = (v1.x - v0.x) * (v1.x - v0.x) + (v1.y - v0.y) * (v1.y - v0.y);
+			real_t b = Math::pow(a, alpha * 0.5f);
+			return b + t;
+		};
+		real_t t0 = 0.0;
+		real_t t1 = compute_t(t0, p_alpha, p0, p1);
+		real_t t2 = compute_t(t1, p_alpha, p1, p2);
+		real_t t3 = compute_t(t2, p_alpha, p2, p3);
+		real_t t = Math::lerp(t1, t2, p_t);
+		Vector2 a1 = (t1 - t) / (t1 - t0) * p0 + (t - t0) / (t1 - t0) * p1;
+		Vector2 a2 = (t2 - t) / (t2 - t1) * p1 + (t - t1) / (t2 - t1) * p2;
+		Vector2 a3 = (t3 - t) / (t3 - t2) * p2 + (t - t2) / (t3 - t2) * p3;
+		Vector2 b1 = (t2 - t) / (t2 - t0) * a1 + (t - t0) / (t2 - t0) * a2;
+		Vector2 b2 = (t3 - t) / (t3 - t1) * a2 + (t - t1) / (t3 - t1) * a3;
+		c = (t2 - t) / (t2 - t1) * b1 + (t - t1) / (t2 - t1) * b2;
+	} else {
+		// Uniform, faster to compute, duplicate points allowed but not recommended...
+		real_t t2 = p_t * p_t;
+		real_t t3 = t2 * p_t;
+		c = 0.5 * (2 * p1 + (-1 * p0 + p2) * p_t + (2 * p0 - 5 * p1 + 4 * p2 - p3) * t2 + (-1 * p0 + 3 * p1 - 3 * p2 + p3) * t3);
+	}
+	return c;
+}
+
+Vector<Point2> GoostGeometry2D::smooth_polyline(const Vector<Point2> &p_polyline, float p_density, float p_alpha) {
+	ERR_FAIL_COND_V_MSG(p_polyline.size() < 3, Vector<Point2>(),
+			"Cannot smooth polyline: requires at least 3 points for interpolation.");
+
+	if (p_density <= 1.0f) {
+		return p_polyline; // No need to interpolate.
+	}
+	Vector<Point2> pts = p_polyline;
+	// Extrapolate first and last points to act as control points.
+	const Vector2 &d1 = pts[0] - pts[1];
+	pts.insert(0, pts[0] + d1);
+	const Vector2 &d2 = pts[pts.size() - 1] - pts[pts.size() - 2];
+	pts.insert(pts.size(), pts[pts.size() - 1] + d2);
+
+	const int point_count = p_polyline.size() * p_density;
+	const real_t length = polyline_length(p_polyline);
+
+	const Point2 *p = pts.ptr();
+	Vector<Point2> smoothed;
+	for (int i = 0; i < pts.size() - 3; ++i) {
+		// Weighted distribution.
+		const real_t segment_length = p[i + 1].distance_to(p[i + 2]);
+		const int pc = Math::ceil(point_count * segment_length / length);
+		for (int j = 0; j < pc; ++j) {
+			real_t t = 1.0 / pc * j;
+			smoothed.push_back(catmull_rom(
+					p[i + 0], p[i + 1], p[i + 2], p[i + 3], t, p_alpha));
+		}
+	}
+	smoothed.push_back(p[pts.size() - 2]);
+	return smoothed;
+}
+
+Vector<Point2> GoostGeometry2D::smooth_polygon(const Vector<Point2> &p_polygon, float p_density, float p_alpha) {
+	ERR_FAIL_COND_V_MSG(p_polygon.size() < 3, Vector<Point2>(), "Bad polygon!");
+
+	if (p_density <= 1.0f) {
+		return p_polygon; // No need to interpolate.
+	}
+	const int point_count = p_polygon.size() * p_density;
+	const real_t perimeter = polygon_perimeter(p_polygon);
+
+	const int s = p_polygon.size();
+	const Point2 *p = p_polygon.ptr();
+	auto pt = [&](int idx) {
+		return p[(idx % s + s) % s];
+	};
+	Vector<Point2> smoothed;
+	for (int i = 0; i < s; ++i) {
+		// Weighted distribution.
+		const real_t segment_length = pt(i + 0).distance_to(pt(i + 1));
+		const int pc = Math::ceil(point_count * segment_length / perimeter);
+		for (int j = 0; j < pc; ++j) {
+			real_t t = 1.0 / pc * j;
+			smoothed.push_back(catmull_rom(
+					pt(i - 1), pt(i + 0), pt(i + 1), pt(i + 2), t, p_alpha));
+		}
+	}
+	return smoothed;
+}
+
+// Approximate polygon smoothing using Chaikin's corner-cutting algorithm.
 // https://www.cs.unc.edu/~dm/UNC/COMP258/LECTURES/Chaikins-Algorithm.pdf
 //
-Vector<Point2> GoostGeometry2D::smooth_polygon_approx(const Vector<Point2> &p_polygon, int p_iterations, real_t p_cut_distance) {
+Vector<Point2> GoostGeometry2D::smooth_polygon_approx(const Vector<Point2> &p_polygon, int p_iterations, float p_cut_distance) {
 	ERR_FAIL_COND_V_MSG(p_polygon.size() < 3, Vector<Point2>(), "Bad polygon!");
 
 	Vector<Point2> subject = p_polygon;
-	const real_t cd = CLAMP(p_cut_distance, 0.0, 0.5);
+	const float cd = CLAMP(p_cut_distance, 0.0f, 0.5f);
 	for (int i = 0; i < p_iterations; ++i) {
 		Vector<Point2> smoothed;
 		for (int j = 0; j < subject.size(); ++j) {
@@ -195,17 +294,17 @@ Vector<Point2> GoostGeometry2D::smooth_polygon_approx(const Vector<Point2> &p_po
 	return subject;
 }
 
-// Approximate polyline smoothing using Chaikin algorithm:
+// Approximate polyline smoothing using Chaikin's corner-cutting algorithm:
 // https://www.cs.unc.edu/~dm/UNC/COMP258/LECTURES/Chaikins-Algorithm.pdf
 //
 // Unlike polygon version, the endpoints are always retained.
 //
-Vector<Point2> GoostGeometry2D::smooth_polyline_approx(const Vector<Point2> &p_polyline, int p_iterations, real_t p_cut_distance) {
+Vector<Point2> GoostGeometry2D::smooth_polyline_approx(const Vector<Point2> &p_polyline, int p_iterations, float p_cut_distance) {
 	if (p_polyline.size() <= 2) {
 		return p_polyline;
 	}
 	Vector<Point2> subject = p_polyline;
-	const real_t cd = CLAMP(p_cut_distance, 0.0, 0.5);
+	const float cd = CLAMP(p_cut_distance, 0.0f, 0.5f);
 	for (int i = 0; i < p_iterations; ++i) {
 		Vector<Point2> smoothed;
 		smoothed.push_back(subject[0]); // Always add first point.

core/math/2d/geometry/goost_geometry_2d.h

@@ -22,10 +22,12 @@ class GoostGeometry2D {
 	static Vector<Vector<Point2>> triangulate_polygon(const Vector<Point2> &p_polygon);
 	static Vector<Vector<Point2>> decompose_polygon(const Vector<Point2> &p_polygon);
 
-	/* Polygon/Polyline simplification and smoothing */
+	/* Polygon/Polyline smoothing and simplification */
+	static Vector<Point2> smooth_polygon(const Vector<Point2> &p_polygon, float p_density, float p_alpha = 0.5f);
+	static Vector<Point2> smooth_polyline(const Vector<Point2> &p_polyline, float p_density, float p_alpha = 0.5f);
+	static Vector<Point2> smooth_polygon_approx(const Vector<Point2> &p_polygon, int p_iterations = 1, float p_cut_distance = 0.25f);
+	static Vector<Point2> smooth_polyline_approx(const Vector<Point2> &p_polyline, int p_iterations = 1, float p_cut_distance = 0.25f);
 	static Vector<Point2> simplify_polyline(const Vector<Point2> &p_polyline, real_t p_epsilon);
-	static Vector<Point2> smooth_polygon_approx(const Vector<Point2> &p_polygon, int p_iterations = 1, real_t p_cut_distance = 0.25);
-	static Vector<Point2> smooth_polyline_approx(const Vector<Point2> &p_polyline, int p_iterations = 1, real_t p_cut_distance = 0.25);
 
 	/* Polygon/Polyline attributes */
 	static Point2 polygon_centroid(const Vector<Point2> &p_polygon);

core/math/2d/geometry/goost_geometry_2d_bind.cpp

@@ -106,11 +106,19 @@ Vector<Point2> _GoostGeometry2D::simplify_polyline(const Vector<Point2> &p_polyl
 	return GoostGeometry2D::simplify_polyline(p_polyline, p_epsilon);
 }
 
-Vector<Point2> _GoostGeometry2D::smooth_polygon_approx(const Vector<Point2> &p_polygon, int p_iterations, real_t cut_distance) const {
+Vector<Point2> _GoostGeometry2D::smooth_polygon(const Vector<Point2> &p_polygon, float p_density, float p_alpha) const {
+	return GoostGeometry2D::smooth_polygon(p_polygon, p_density, p_alpha);
+}
+
+Vector<Point2> _GoostGeometry2D::smooth_polyline(const Vector<Point2> &p_polyline, float p_density, float p_alpha) const {
+	return GoostGeometry2D::smooth_polyline(p_polyline, p_density, p_alpha);
+}
+
+Vector<Point2> _GoostGeometry2D::smooth_polygon_approx(const Vector<Point2> &p_polygon, int p_iterations, float cut_distance) const {
 	return GoostGeometry2D::smooth_polygon_approx(p_polygon, p_iterations, cut_distance);
 }
 
-Vector<Point2> _GoostGeometry2D::smooth_polyline_approx(const Vector<Point2> &p_polyline, int p_iterations, real_t cut_distance) const {
+Vector<Point2> _GoostGeometry2D::smooth_polyline_approx(const Vector<Point2> &p_polyline, int p_iterations, float cut_distance) const {
 	return GoostGeometry2D::smooth_polyline_approx(p_polyline, p_iterations, cut_distance);
 }
 
@@ -172,8 +180,10 @@ void _GoostGeometry2D::_bind_methods() {
 	ClassDB::bind_method(D_METHOD("decompose_polygon", "polygon"), &_GoostGeometry2D::decompose_polygon);
 
 	ClassDB::bind_method(D_METHOD("simplify_polyline", "polyline", "epsilon"), &_GoostGeometry2D::simplify_polyline);
-	ClassDB::bind_method(D_METHOD("smooth_polygon_approx", "polygon", "iterations", "cut_distance"), &_GoostGeometry2D::smooth_polygon_approx, DEFVAL(1), DEFVAL(0.25));
-	ClassDB::bind_method(D_METHOD("smooth_polyline_approx", "polyline", "iterations", "cut_distance"), &_GoostGeometry2D::smooth_polyline_approx, DEFVAL(1), DEFVAL(0.25));
+	ClassDB::bind_method(D_METHOD("smooth_polygon", "polygon", "density", "alpha"), &_GoostGeometry2D::smooth_polygon, DEFVAL(0.5f));
+	ClassDB::bind_method(D_METHOD("smooth_polyline", "polyline", "density", "alpha"), &_GoostGeometry2D::smooth_polyline, DEFVAL(0.5f));
+	ClassDB::bind_method(D_METHOD("smooth_polygon_approx", "polygon", "iterations", "cut_distance"), &_GoostGeometry2D::smooth_polygon_approx, DEFVAL(1), DEFVAL(0.25f));
+	ClassDB::bind_method(D_METHOD("smooth_polyline_approx", "polyline", "iterations", "cut_distance"), &_GoostGeometry2D::smooth_polyline_approx, DEFVAL(1), DEFVAL(0.25f));
 
 	ClassDB::bind_method(D_METHOD("polygon_centroid", "polygon"), &_GoostGeometry2D::polygon_centroid);
 	ClassDB::bind_method(D_METHOD("polygon_area", "polygon"), &_GoostGeometry2D::polygon_area);

core/math/2d/geometry/goost_geometry_2d_bind.h

@@ -30,9 +30,11 @@ class _GoostGeometry2D : public Object {
 	Array triangulate_polygon(const Vector<Point2> &p_polygon) const;
 	Array decompose_polygon(const Vector<Point2> &p_polygon) const;
 
+	Vector<Point2> smooth_polygon(const Vector<Point2> &p_polygon, float p_density, float p_alpha = 0.5f) const;
+	Vector<Point2> smooth_polyline(const Vector<Point2> &p_polyline, float p_density, float p_alpha = 0.5f) const;
+	Vector<Point2> smooth_polygon_approx(const Vector<Point2> &p_polygon, int p_iterations = 1, float cut_distance = 0.25f) const;
+	Vector<Point2> smooth_polyline_approx(const Vector<Point2> &p_polyline, int p_iterations = 1, float cut_distance = 0.25f) const;
 	Vector<Point2> simplify_polyline(const Vector<Point2> &p_polyline, real_t p_epsilon) const;
-	Vector<Point2> smooth_polygon_approx(const Vector<Point2> &p_polygon, int p_iterations = 1, real_t cut_distance = 0.25) const;
-	Vector<Point2> smooth_polyline_approx(const Vector<Point2> &p_polyline, int p_iterations = 1, real_t cut_distance = 0.25) const;
 
 	Vector2 polygon_centroid(const Vector<Vector2> &p_polygon) const;
 	real_t polygon_area(const Vector<Vector2> &p_polygon) const;

doc/GoostGeometry2D.xml

@@ -221,6 +221,22 @@
 				Simplifies a polyline by reducing the number of points using the Ramer-Douglas-Peucker (RDP) algorithm. Higher [code]epsilon[/code] values result in fewer points retained.
 			</description>
 		</method>
+		<method name="smooth_polygon" qualifiers="const">
+			<return type="PoolVector2Array">
+			</return>
+			<argument index="0" name="polygon" type="PoolVector2Array">
+			</argument>
+			<argument index="1" name="density" type="float">
+			</argument>
+			<argument index="2" name="alpha" type="float" default="0.5">
+			</argument>
+			<description>
+				Smoothers the polygon using the Catmull-Rom's interpolating spline, resulting in larger number of vertices.
+				The [code]density[/code] parameter configures the desired number of vertices in the output polygon: [code]n = polygon.size() * density[/code], where [code]n[/code] is the point count computed. If [code]density &lt; 1.0[/code], returns original [code]polygon[/code]. The number of vertices is weighted per segment according to the [method polygon_perimeter].
+				The [code]alpha[/code] parameter determines the type of the Catmull-Rom's spline: uniform - [code]alpha == 0[/code], centripetal - [code]alpha == 0.5[/code], chordal - [code]alpha > 0.5[/code]. The default value of [code]0.5[/code] is recommended for eliminating self-intersections and cusps.
+				For faster, approximate smoothing method, see [method smooth_polygon_approx].
+			</description>
+		</method>
 		<method name="smooth_polygon_approx" qualifiers="const">
 			<return type="PoolVector2Array">
 			</return>
@@ -232,6 +248,23 @@
 			</argument>
 			<description>
 				Approximately smoothers the polygon using the Chaikin's algorithm resulting in larger number of vertices. Number of [code]iterations[/code] can be specified to produce smoother polygons. The [code]cut_distance[/code] determines at what distance new control points are selected from segments.
+				Unlike [method smooth_polygon], the resulting curve does not go through input vertices, but instead touches the segments of the original [code]polygon[/code].
+			</description>
+		</method>
+		<method name="smooth_polyline" qualifiers="const">
+			<return type="PoolVector2Array">
+			</return>
+			<argument index="0" name="polyline" type="PoolVector2Array">
+			</argument>
+			<argument index="1" name="density" type="float">
+			</argument>
+			<argument index="2" name="alpha" type="float" default="0.5">
+			</argument>
+			<description>
+				Smoothers the polyline using the Catmull-Rom's interpolating spline, resulting in larger number of vertices.
+				The [code]density[/code] parameter configures the desired number of vertices in the output polyline: [code]n = polyline.size() * density[/code], where [code]n[/code] is the point count computed. If [code]density &lt; 1.0[/code], returns original [code]polyline[/code]. The number of vertices is weighted per segment according to the [method polyline_length].
+				The [code]alpha[/code] parameter determines the type of the Catmull-Rom's spline: uniform - [code]alpha == 0[/code], centripetal - [code]alpha == 0.5[/code], chordal - [code]alpha > 0.5[/code]. The default value of [code]0.5[/code] is recommended for eliminating self-intersections and cusps.
+				For faster, approximate smoothing method, see [method smooth_polyline_approx].
 			</description>
 		</method>
 		<method name="smooth_polyline_approx" qualifiers="const">
@@ -245,7 +278,8 @@
 			</argument>
 			<description>
 				Approximately smoothers the polyline using the Chaikin's algorithm resulting in larger number of vertices. Number of [code]iterations[/code] can be specified to produce smoother polylines. The [code]cut_distance[/code] determines at what distance new control points are selected from segments.
-				Unlike [method smooth_polygon_approx], this method always retains start and end points from the original polyline.
+				Unlike [method smooth_polyline], the resulting curve does not go through input vertices, but instead touches the segments of the original polyline.
+				Unlike [method smooth_polygon_approx], this method always retains start and end points from the original [code]polyline[/code].
 			</description>
 		</method>
 		<method name="triangulate_polygon" qualifiers="const">

tests/project/goost/core/math/2d/geometry/test_geometry_2d.gd

@@ -199,6 +199,34 @@ func test_simplify_polyline():
 		assert_eq(simplified[i], control[i])
 
 
+func test_smooth_polygon():
+	var input = [Vector2(26, 20), Vector2(73, 23), Vector2(72, 62), Vector2(29, 57)]
+	var control = [Vector2(26, 20), Vector2(49.311768, 15.934073), Vector2(73, 23),
+			Vector2(77.531448, 43.002853), Vector2(72, 62), Vector2(50.609352, 64.723686),
+			Vector2(29, 57), Vector2(22.657887, 38.113762)]
+	var smoothed = GoostGeometry2D.smooth_polygon(input, 1.6)
+
+	if smoothed.size() != control.size():
+		assert_eq(smoothed.size(), control.size(), "Point count mismatch")
+		return
+	for i in smoothed.size():
+		assert_eq(smoothed[i], control[i])
+
+
+func test_smooth_polyline():
+	var input = [Vector2(26, 20), Vector2(73, 23), Vector2(72, 62), Vector2(29, 57)]
+	var control = [Vector2(26, 20), Vector2(43.531898, 19.454647), Vector2(61.063797, 18.909294),
+			Vector2(73, 23), Vector2(77.531448, 43.002853), Vector2(72, 62), Vector2(60.854603, 63.835579),
+			Vector2(44.927299, 60.417786), Vector2(29, 57)]
+	var smoothed = GoostGeometry2D.smooth_polyline(input, 1.6)
+
+	if smoothed.size() != control.size():
+		assert_eq(smoothed.size(), control.size(), "Point count mismatch")
+		return
+	for i in smoothed.size():
+		assert_eq(smoothed[i], control[i])
+
+
 func test_smooth_polygon_approx():
 	var input = [Vector2(25, 83), Vector2(49, 16), Vector2(66, 79)]
 	var control = [Vector2(31, 66.25), Vector2(43, 32.75), Vector2(53.25, 31.75), Vector2(61.75, 63.25), Vector2(55.75, 80), Vector2(35.25, 82)]
@@ -234,9 +262,21 @@ class Stress extends "res://addons/gut/test.gd":
 		var input = GoostGeometry2D.circle(1024)
 		for i in input.size():
 			input[i] += Random2D.point_in_circle(100)
-		for i in 100000:
+		for i in 10000:
+			var t1 = OS.get_ticks_msec()
+			var _out = GoostGeometry2D.simplify_polyline(input, 100.0)
+			var t2 = OS.get_ticks_msec()
+			time += t2 - t1
+		gut.p(time / 10000.0)
+
+	func test_smooth_polyline():
+		var time = 0
+		var input = GoostGeometry2D.regular_polygon(1024, 6)
+		for i in input.size():
+			input[i] += Random2D.point_in_circle(100)
+		for i in 100:
 			var t1 = OS.get_ticks_msec()
-			var _simplified = GoostGeometry2D.simplify_polyline(input, 100.0)
+			var _out = GoostGeometry2D.smooth_polyline(input, 20)
 			var t2 = OS.get_ticks_msec()
 			time += t2 - t1
-		gut.p(time / 100000.0)
+		gut.p(time / 100.0)