Add interpolated versions of polyline and polygon smoothing

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

@@ -174,6 +174,102 @@ Vector<Point2> GoostGeometry2D::simplify_polyline(const Vector<Point2> &p_polyli
 	return ret;
 }
 
+// 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 = Math::pow(v1.x - v0.x, 2.0f) + Math::pow(v1.y - v0.y, 2.0f);
+			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...
+		c = 0.5 * (2 * p1 + (-1 * p0 + p2) * p_t + (2 * p0 - 5 * p1 + 4 * p2 - p3) * Math::pow(p_t, 2) + (-1 * p0 + 3 * p1 - 3 * p2 + p3) * Math::pow(p_t, 3));
+	}
+	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.");
+
+	const int point_count = p_polyline.size() * p_density;
+	if (point_count <= p_polyline.size()) {
+		// No need to interpolate.
+		return p_polyline;
+	}
+	Vector<Point2> pts = p_polyline;
+	// Extrapolate first and last points to act as control points.
+	pts.insert(0, pts[1] - 2 * pts[0]);
+	pts.insert(pts.size(), 2 * pts[pts.size() - 1] - pts[pts.size() - 2]);
+
+	const real_t length = polyline_length(p_polyline);
+
+	Vector<Point2> smoothed;
+	for (int i = 0; i < pts.size() - 3; ++i) {
+		// Weighted distribution.
+		const real_t segment_length = pts[i + 1].distance_to(pts[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(
+					pts[i + 0], pts[i + 1], pts[i + 2], pts[i + 3], t, p_alpha));
+		}
+	}
+	smoothed.push_back(pts[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() < 4, Vector<Point2>(),
+			"Cannot smooth polygon: requires at least 4 points for interpolation.");
+
+	const int point_count = p_polygon.size() * p_density;
+	if (point_count <= p_polygon.size()) {
+		// No need to interpolate.
+		return p_polygon;
+	}
+	Vector<Point2> pts = p_polygon;
+	const real_t perimeter = polygon_perimeter(p_polygon);
+
+	auto pt = [&](int i) {
+		return pts[((i % pts.size()) + pts.size()) % pts.size()];
+	};
+	Vector<Point2> smoothed;
+	for (int i = 0; i < pts.size(); ++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 algorithm.
 // https://www.cs.unc.edu/~dm/UNC/COMP258/LECTURES/Chaikins-Algorithm.pdf
 //

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 */
-	static Vector<Point2> simplify_polyline(const Vector<Point2> &p_polyline, real_t p_epsilon);
+	/* Polygon/Polyline smoothing and simplification */
+	static Vector<Point2> smooth_polygon(const Vector<Point2> &p_polygon, float p_density, float p_alpha = 0.5);
+	static Vector<Point2> smooth_polyline(const Vector<Point2> &p_polyline, float p_density, float p_alpha = 0.5);
 	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);
+	static Vector<Point2> simplify_polyline(const Vector<Point2> &p_polyline, real_t p_epsilon);
 
 	/* Polygon/Polyline attributes */
 	static Point2 polygon_centroid(const Vector<Point2> &p_polygon);

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

@@ -106,6 +106,14 @@ Vector<Point2> _GoostGeometry2D::simplify_polyline(const Vector<Point2> &p_polyl
 	return GoostGeometry2D::simplify_polyline(p_polyline, p_epsilon);
 }
 
+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, real_t 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, real_t cut_distance) const {
 	return GoostGeometry2D::smooth_polygon_approx(p_polygon, p_iterations, cut_distance);
 }
@@ -172,6 +180,8 @@ 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", "polygon", "density", "alpha"), &_GoostGeometry2D::smooth_polygon, DEFVAL(0.5));
+	ClassDB::bind_method(D_METHOD("smooth_polyline", "polyline", "density", "alpha"), &_GoostGeometry2D::smooth_polyline, DEFVAL(0.5));
 	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));
 

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

@@ -31,6 +31,8 @@ class _GoostGeometry2D : public Object {
 	Array decompose_polygon(const Vector<Point2> &p_polygon) const;
 
 	Vector<Point2> simplify_polyline(const Vector<Point2> &p_polyline, real_t p_epsilon) const;
+	Vector<Point2> smooth_polygon(const Vector<Point2> &p_polygon, float p_density, float p_alpha = 0.5) const;
+	Vector<Point2> smooth_polyline(const Vector<Point2> &p_polyline, float p_density, float p_alpha = 0.5) 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;
 

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">