diff --git a/core/math/2d/geometry/goost_geometry_2d.cpp b/core/math/2d/geometry/goost_geometry_2d.cpp index 6d5b15e8..c0574c93 100644 --- a/core/math/2d/geometry/goost_geometry_2d.cpp +++ b/core/math/2d/geometry/goost_geometry_2d.cpp @@ -174,14 +174,113 @@ Vector GoostGeometry2D::simplify_polyline(const Vector &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 GoostGeometry2D::smooth_polyline(const Vector &p_polyline, float p_density, float p_alpha) { + ERR_FAIL_COND_V_MSG(p_polyline.size() < 3, Vector(), + "Cannot smooth polyline: requires at least 3 points for interpolation."); + + if (p_density <= 1.0f) { + return p_polyline; // No need to interpolate. + } + Vector 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 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 GoostGeometry2D::smooth_polygon(const Vector &p_polygon, float p_density, float p_alpha) { + ERR_FAIL_COND_V_MSG(p_polygon.size() < 3, Vector(), "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 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 GoostGeometry2D::smooth_polygon_approx(const Vector &p_polygon, int p_iterations, real_t p_cut_distance) { +Vector GoostGeometry2D::smooth_polygon_approx(const Vector &p_polygon, int p_iterations, float p_cut_distance) { ERR_FAIL_COND_V_MSG(p_polygon.size() < 3, Vector(), "Bad polygon!"); Vector 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 smoothed; for (int j = 0; j < subject.size(); ++j) { @@ -195,17 +294,17 @@ Vector GoostGeometry2D::smooth_polygon_approx(const Vector &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 GoostGeometry2D::smooth_polyline_approx(const Vector &p_polyline, int p_iterations, real_t p_cut_distance) { +Vector GoostGeometry2D::smooth_polyline_approx(const Vector &p_polyline, int p_iterations, float p_cut_distance) { if (p_polyline.size() <= 2) { return p_polyline; } Vector 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 smoothed; smoothed.push_back(subject[0]); // Always add first point. diff --git a/core/math/2d/geometry/goost_geometry_2d.h b/core/math/2d/geometry/goost_geometry_2d.h index 4f02d01f..b5fe71de 100644 --- a/core/math/2d/geometry/goost_geometry_2d.h +++ b/core/math/2d/geometry/goost_geometry_2d.h @@ -22,10 +22,12 @@ class GoostGeometry2D { static Vector> triangulate_polygon(const Vector &p_polygon); static Vector> decompose_polygon(const Vector &p_polygon); - /* Polygon/Polyline simplification and smoothing */ + /* Polygon/Polyline smoothing and simplification */ + static Vector smooth_polygon(const Vector &p_polygon, float p_density, float p_alpha = 0.5f); + static Vector smooth_polyline(const Vector &p_polyline, float p_density, float p_alpha = 0.5f); + static Vector smooth_polygon_approx(const Vector &p_polygon, int p_iterations = 1, float p_cut_distance = 0.25f); + static Vector smooth_polyline_approx(const Vector &p_polyline, int p_iterations = 1, float p_cut_distance = 0.25f); static Vector simplify_polyline(const Vector &p_polyline, real_t p_epsilon); - static Vector smooth_polygon_approx(const Vector &p_polygon, int p_iterations = 1, real_t p_cut_distance = 0.25); - static Vector smooth_polyline_approx(const Vector &p_polyline, int p_iterations = 1, real_t p_cut_distance = 0.25); /* Polygon/Polyline attributes */ static Point2 polygon_centroid(const Vector &p_polygon); diff --git a/core/math/2d/geometry/goost_geometry_2d_bind.cpp b/core/math/2d/geometry/goost_geometry_2d_bind.cpp index f3465714..d2909f5c 100644 --- a/core/math/2d/geometry/goost_geometry_2d_bind.cpp +++ b/core/math/2d/geometry/goost_geometry_2d_bind.cpp @@ -106,11 +106,19 @@ Vector _GoostGeometry2D::simplify_polyline(const Vector &p_polyl return GoostGeometry2D::simplify_polyline(p_polyline, p_epsilon); } -Vector _GoostGeometry2D::smooth_polygon_approx(const Vector &p_polygon, int p_iterations, real_t cut_distance) const { +Vector _GoostGeometry2D::smooth_polygon(const Vector &p_polygon, float p_density, float p_alpha) const { + return GoostGeometry2D::smooth_polygon(p_polygon, p_density, p_alpha); +} + +Vector _GoostGeometry2D::smooth_polyline(const Vector &p_polyline, float p_density, float p_alpha) const { + return GoostGeometry2D::smooth_polyline(p_polyline, p_density, p_alpha); +} + +Vector _GoostGeometry2D::smooth_polygon_approx(const Vector &p_polygon, int p_iterations, float cut_distance) const { return GoostGeometry2D::smooth_polygon_approx(p_polygon, p_iterations, cut_distance); } -Vector _GoostGeometry2D::smooth_polyline_approx(const Vector &p_polyline, int p_iterations, real_t cut_distance) const { +Vector _GoostGeometry2D::smooth_polyline_approx(const Vector &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); diff --git a/core/math/2d/geometry/goost_geometry_2d_bind.h b/core/math/2d/geometry/goost_geometry_2d_bind.h index 134acde2..eb8dded4 100644 --- a/core/math/2d/geometry/goost_geometry_2d_bind.h +++ b/core/math/2d/geometry/goost_geometry_2d_bind.h @@ -30,9 +30,11 @@ class _GoostGeometry2D : public Object { Array triangulate_polygon(const Vector &p_polygon) const; Array decompose_polygon(const Vector &p_polygon) const; + Vector smooth_polygon(const Vector &p_polygon, float p_density, float p_alpha = 0.5f) const; + Vector smooth_polyline(const Vector &p_polyline, float p_density, float p_alpha = 0.5f) const; + Vector smooth_polygon_approx(const Vector &p_polygon, int p_iterations = 1, float cut_distance = 0.25f) const; + Vector smooth_polyline_approx(const Vector &p_polyline, int p_iterations = 1, float cut_distance = 0.25f) const; Vector simplify_polyline(const Vector &p_polyline, real_t p_epsilon) const; - Vector smooth_polygon_approx(const Vector &p_polygon, int p_iterations = 1, real_t cut_distance = 0.25) const; - Vector smooth_polyline_approx(const Vector &p_polyline, int p_iterations = 1, real_t cut_distance = 0.25) const; Vector2 polygon_centroid(const Vector &p_polygon) const; real_t polygon_area(const Vector &p_polygon) const; diff --git a/doc/GoostGeometry2D.xml b/doc/GoostGeometry2D.xml index 84f6fd44..65aff192 100644 --- a/doc/GoostGeometry2D.xml +++ b/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. + + + + + + + + + + + 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 < 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]. + + @@ -232,6 +248,23 @@ 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]. + + + + + + + + + + + + + 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 < 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]. @@ -245,7 +278,8 @@ 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]. diff --git a/tests/project/goost/core/math/2d/geometry/test_geometry_2d.gd b/tests/project/goost/core/math/2d/geometry/test_geometry_2d.gd index 0ccce76f..29f82f6e 100644 --- a/tests/project/goost/core/math/2d/geometry/test_geometry_2d.gd +++ b/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)