diff --git a/osu.Framework.Benchmarks/BenchmarkPiecewiseLinearToBezier.cs b/osu.Framework.Benchmarks/BenchmarkPiecewiseLinearToBezier.cs new file mode 100644 index 0000000000..918f537d39 --- /dev/null +++ b/osu.Framework.Benchmarks/BenchmarkPiecewiseLinearToBezier.cs @@ -0,0 +1,43 @@ +// Copyright (c) ppy Pty Ltd . Licensed under the MIT Licence. +// See the LICENCE file in the repository root for full licence text. + +using System.Collections.Generic; +using BenchmarkDotNet.Attributes; +using osu.Framework.Utils; +using osuTK; + +namespace osu.Framework.Benchmarks +{ + public class BenchmarkPiecewiseLinearToBezier : BenchmarkTest + { + private Vector2[] inputPath = null!; + + [Params(5, 25)] + public int NumControlPoints; + + [Params(5, 200)] + public int NumTestPoints; + + [Params(0, 100, 200)] + public int MaxIterations; + + public override void SetUp() + { + base.SetUp(); + + Vector2[] points = new Vector2[5]; + points[0] = new Vector2(50, 250); + points[1] = new Vector2(150, 230); + points[2] = new Vector2(100, 150); + points[3] = new Vector2(200, 80); + points[4] = new Vector2(250, 50); + inputPath = PathApproximator.LagrangePolynomialToPiecewiseLinear(points).ToArray(); + } + + [Benchmark] + public List PiecewiseLinearToBezier() + { + return PathApproximator.PiecewiseLinearToBezier(inputPath, NumControlPoints, NumTestPoints, MaxIterations); + } + } +} diff --git a/osu.Framework.Tests/Visual/Drawables/TestScenePiecewiseLinearToBSpline.cs b/osu.Framework.Tests/Visual/Drawables/TestScenePiecewiseLinearToBSpline.cs new file mode 100644 index 0000000000..b6a21c7459 --- /dev/null +++ b/osu.Framework.Tests/Visual/Drawables/TestScenePiecewiseLinearToBSpline.cs @@ -0,0 +1,201 @@ +// Copyright (c) ppy Pty Ltd . Licensed under the MIT Licence. +// See the LICENCE file in the repository root for full licence text. + +using System; +using System.Collections.Generic; +using osu.Framework.Graphics; +using osu.Framework.Graphics.Lines; +using osu.Framework.Graphics.Sprites; +using osu.Framework.Utils; +using osu.Framework.Testing; +using osuTK; +using osuTK.Graphics; + +namespace osu.Framework.Tests.Visual.Drawables +{ + public partial class TestScenePiecewiseLinearToBSpline : GridTestScene + { + private int numControlPoints = 5; + private int degree = 2; + private int numTestPoints = 100; + private int maxIterations = 100; + private float learningRate = 8; + private float b1 = 0.8f; + private float b2 = 0.99f; + + private readonly List doubleApproximatedPathTests = new List(); + + public TestScenePiecewiseLinearToBSpline() + : base(2, 2) + { + doubleApproximatedPathTests.Add(new DoubleApproximatedPathTest(PathApproximator.BezierToPiecewiseLinear)); + Cell(0).AddRange(new[] + { + createLabel(nameof(PathApproximator.BezierToPiecewiseLinear)), + new ApproximatedPathTest(PathApproximator.BezierToPiecewiseLinear), + doubleApproximatedPathTests[^1], + }); + + doubleApproximatedPathTests.Add(new DoubleApproximatedPathTest(PathApproximator.CatmullToPiecewiseLinear)); + Cell(1).AddRange(new[] + { + createLabel(nameof(PathApproximator.CatmullToPiecewiseLinear)), + new ApproximatedPathTest(PathApproximator.CatmullToPiecewiseLinear), + doubleApproximatedPathTests[^1], + }); + + doubleApproximatedPathTests.Add(new DoubleApproximatedPathTest(PathApproximator.CircularArcToPiecewiseLinear)); + Cell(2).AddRange(new[] + { + createLabel(nameof(PathApproximator.CircularArcToPiecewiseLinear)), + new ApproximatedPathTest(PathApproximator.CircularArcToPiecewiseLinear), + doubleApproximatedPathTests[^1], + }); + + doubleApproximatedPathTests.Add(new DoubleApproximatedPathTest(PathApproximator.LagrangePolynomialToPiecewiseLinear)); + Cell(3).AddRange(new[] + { + createLabel(nameof(PathApproximator.LagrangePolynomialToPiecewiseLinear)), + new ApproximatedPathTest(PathApproximator.LagrangePolynomialToPiecewiseLinear), + doubleApproximatedPathTests[^1], + }); + + AddSliderStep($"{nameof(numControlPoints)}", 3, 25, 5, v => + { + numControlPoints = v; + updateTests(); + }); + + AddSliderStep($"{nameof(degree)}", 1, 5, 3, v => + { + degree = v; + updateTests(); + }); + + AddSliderStep($"{nameof(numTestPoints)}", 10, 200, 100, v => + { + numTestPoints = v; + updateTests(); + }); + + AddSliderStep($"{nameof(maxIterations)}", 0, 200, 10, v => + { + maxIterations = v; + updateTests(); + }); + + AddSliderStep($"{nameof(learningRate)}", 0, 10, 8f, v => + { + learningRate = v; + updateTests(); + }); + + AddSliderStep($"{nameof(b1)}", 0, 0.999f, 0.8f, v => + { + b1 = v; + updateTests(); + }); + + AddSliderStep($"{nameof(b2)}", 0, 0.999f, 0.99f, v => + { + b2 = v; + updateTests(); + }); + + AddStep("Enable optimization", () => + { + foreach (var test in doubleApproximatedPathTests) + test.OptimizePath = true; + updateTests(); + }); + } + + private void updateTests() + { + foreach (var test in doubleApproximatedPathTests) + { + test.NumControlPoints = numControlPoints; + test.Degree = degree; + test.NumTestPoints = numTestPoints; + test.MaxIterations = maxIterations; + test.LearningRate = learningRate; + test.B1 = b1; + test.B2 = b2; + test.UpdatePath(); + } + } + + private Drawable createLabel(string text) => new SpriteText + { + Text = text + "ToBSpline", + Font = new FontUsage(size: 20), + Colour = Color4.White, + }; + + public delegate List ApproximatorFunc(ReadOnlySpan controlPoints); + + private partial class ApproximatedPathTest : SmoothPath + { + public ApproximatedPathTest(ApproximatorFunc approximator) + { + Vector2[] points = new Vector2[5]; + points[0] = new Vector2(50, 250); + points[1] = new Vector2(150, 230); + points[2] = new Vector2(100, 150); + points[3] = new Vector2(200, 80); + points[4] = new Vector2(250, 50); + + AutoSizeAxes = Axes.None; + RelativeSizeAxes = Axes.Both; + PathRadius = 2; + Vertices = approximator(points); + Colour = Color4.White; + } + } + + private partial class DoubleApproximatedPathTest : SmoothPath + { + private readonly Vector2[] inputPath; + + public int NumControlPoints { get; set; } + + public int Degree { get; set; } + + public int NumTestPoints { get; set; } + + public int MaxIterations { get; set; } + + public float LearningRate { get; set; } + + public float B1 { get; set; } + + public float B2 { get; set; } + + public bool OptimizePath { get; set; } + + public DoubleApproximatedPathTest(ApproximatorFunc approximator) + { + Vector2[] points = new Vector2[5]; + points[0] = new Vector2(50, 250); + points[1] = new Vector2(150, 230); + points[2] = new Vector2(100, 150); + points[3] = new Vector2(200, 80); + points[4] = new Vector2(250, 50); + + AutoSizeAxes = Axes.None; + RelativeSizeAxes = Axes.Both; + PathRadius = 2; + Colour = Color4.Magenta; + inputPath = approximator(points).ToArray(); + } + + public void UpdatePath() + { + if (!OptimizePath) return; + + var controlPoints = PathApproximator.PiecewiseLinearToBSpline(inputPath, NumControlPoints, Degree, NumTestPoints, MaxIterations, LearningRate, B1, B2); + Vertices = PathApproximator.BSplineToPiecewiseLinear(controlPoints.ToArray(), Degree); + } + } + } +} diff --git a/osu.Framework/Utils/PathApproximator.cs b/osu.Framework/Utils/PathApproximator.cs index d286da0d86..1228f5c75d 100644 --- a/osu.Framework/Utils/PathApproximator.cs +++ b/osu.Framework/Utils/PathApproximator.cs @@ -4,6 +4,7 @@ using System; using System.Collections.Generic; using System.Linq; +using System.Numerics.Tensors; using osu.Framework.Graphics.Primitives; using osuTK; @@ -297,6 +298,471 @@ public static List LagrangePolynomialToPiecewiseLinear(ReadOnlySpan + /// Creates a bezier curve approximation from a piecewise-linear path. + /// + /// The piecewise-linear path to approximate. + /// The number of control points to use in the bezier approximation. + /// The number of points to evaluate the bezier path at for optimization, basically a resolution. + /// The number of optimization steps. + /// The rate of optimization. Larger values converge faster but can be unstable. + /// The B1 parameter for the Adam optimizer. Between 0 and 1. + /// The B2 parameter for the Adam optimizer. Between 0 and 1. + /// The initial bezier control points to use before optimization. The length of this list should be equal to . + /// A List of vectors representing the bezier control points. + public static List PiecewiseLinearToBezier(ReadOnlySpan inputPath, + int numControlPoints, + int numTestPoints = 100, + int maxIterations = 100, + float learningRate = 8f, + float b1 = 0.8f, + float b2 = 0.99f, + List? initialControlPoints = null) + { + return piecewiseLinearToSpline(inputPath, generateBezierWeights(numControlPoints, numTestPoints), maxIterations, learningRate, b1, b2, initialControlPoints); + } + + /// + /// Creates a B-spline approximation from a piecewise-linear path. + /// + /// The piecewise-linear path to approximate. + /// The number of control points to use in the B-spline approximation. + /// The polynomial order. + /// The number of points to evaluate the B-spline path at for optimization, basically a resolution. + /// The number of optimization steps. + /// The rate of optimization. Larger values converge faster but can be unstable. + /// The B1 parameter for the Adam optimizer. Between 0 and 1. + /// The B2 parameter for the Adam optimizer. Between 0 and 1. + /// The initial B-spline control points to use before optimization. The length of this list should be equal to . + /// A List of vectors representing the B-spline control points. + public static List PiecewiseLinearToBSpline(ReadOnlySpan inputPath, + int numControlPoints, + int degree, + int numTestPoints = 100, + int maxIterations = 100, + float learningRate = 8f, + float b1 = 0.8f, + float b2 = 0.99f, + List? initialControlPoints = null) + { + degree = Math.Min(degree, numControlPoints - 1); + return piecewiseLinearToSpline(inputPath, generateBSplineWeights(numControlPoints, numTestPoints, degree), maxIterations, learningRate, b1, b2, initialControlPoints); + } + + /// + /// Creates an arbitrary spline approximation from a piecewise-linear path. + /// Works for any spline type where the interpolation is a linear combination of the control points. + /// + /// The piecewise-linear path to approximate. + /// A 2D matrix that contains the spline basis functions at multiple positions. The length of the first dimension is the number of test points, and the length of the second dimension is the number of control points. + /// The number of optimization steps. + /// The rate of optimization. Larger values converge faster but can be unstable. + /// The B1 parameter for the Adam optimizer. Between 0 and 1. + /// The B2 parameter for the Adam optimizer. Between 0 and 1. + /// The initial control points to use before optimization. The length of this list should be equal to the number of test points. + /// A List of vectors representing the spline control points. + private static List piecewiseLinearToSpline(ReadOnlySpan inputPath, + float[,] weights, + int maxIterations = 100, + float learningRate = 8f, + float b1 = 0.8f, + float b2 = 0.99f, + List? initialControlPoints = null) + { + int numControlPoints = weights.GetLength(1); + int numTestPoints = weights.GetLength(0); + + // Generate transpose weight matrix + float[,] weightsTranspose = new float[numControlPoints, numTestPoints]; + + for (int i = 0; i < numControlPoints; i++) + { + for (int j = 0; j < numTestPoints; j++) + { + weightsTranspose[i, j] = weights[j, i]; + } + } + + // Create efficient interpolation on the input path + var interpolator = new Interpolator(inputPath, numTestPoints); + + // Initialize control points + float[,] labels = new float[2, numTestPoints]; + float[,] controlPoints = new float[2, numControlPoints]; + + if (initialControlPoints is not null) + { + for (int i = 0; i < numControlPoints; i++) + { + controlPoints[0, i] = initialControlPoints[i].X; + controlPoints[1, i] = initialControlPoints[i].Y; + } + } + else + // Create initial control point positions equally spaced along the input path + interpolator.Interpolate(linspace(0, 1, numControlPoints), controlPoints); + + // Initialize Adam optimizer variables + float[,] m = new float[2, numControlPoints]; + float[,] v = new float[2, numControlPoints]; + float[,] learnableMask = new float[2, numControlPoints]; + + for (int i = 1; i < numControlPoints - 1; i++) + { + learnableMask[0, i] = 1; + learnableMask[1, i] = 1; + } + + // Initialize intermediate variables + float[,] points = new float[2, numTestPoints]; + float[,] grad = new float[2, numControlPoints]; + float[] distanceDistribution = new float[numTestPoints]; + + for (int step = 0; step < maxIterations; step++) + { + matmul(controlPoints, weights, points); + + // Update labels to shift the distance distribution between points + if (step % 11 == 0) + { + getDistanceDistribution(points, distanceDistribution); + interpolator.Interpolate(distanceDistribution, labels); + } + + // Calculate the gradient on the control points + matDiff(labels, points, points); + matmul(points, weightsTranspose, grad); + matScale(grad, -1f / numControlPoints, grad); + + // Apply learnable mask to prevent moving the endpoints + matProduct(grad, learnableMask, grad); + + // Update control points with Adam optimizer + matLerp(grad, m, b1, m); + matProduct(grad, grad, grad); + matLerp(grad, v, b2, v); + adamUpdate(controlPoints, m, v, step, learningRate, b1, b2); + } + + // Convert the resulting control points array + var result = new List(numControlPoints); + + for (int i = 0; i < numControlPoints; i++) + { + result.Add(new Vector2(controlPoints[0, i], controlPoints[1, i])); + } + + return result; + } + + private static void adamUpdate(float[,] parameters, float[,] m, float[,] v, int step, float learningRate, float b1, float b2) + { + const float epsilon = 1E-8f; + float mMult = 1 / (1 - MathF.Pow(b1, step + 1)); + float vMult = 1 / (1 - MathF.Pow(b2, step + 1)); + int m0 = m.GetLength(0); + int m1 = m.GetLength(1); + + for (int i = 0; i < m0; i++) + { + for (int j = 0; j < m1; j++) + { + float mCorr = m[i, j] * mMult; + float vCorr = v[i, j] * vMult; + parameters[i, j] -= learningRate * mCorr / (MathF.Sqrt(vCorr) + epsilon); + } + } + } + + private static unsafe void matLerp(float[,] mat1, float[,] mat2, float t, float[,] result) + { + // mat1 can not be the same array as result, or it will not work correctly + if (ReferenceEquals(mat1, result)) + throw new ArgumentException($"{nameof(mat1)} can not be the same array as {nameof(result)}."); + + fixed (float* mat1P = mat1, mat2P = mat2, resultP = result) + { + var span1 = new Span(mat1P, mat1.Length); + var span2 = new Span(mat2P, mat2.Length); + var spanR = new Span(resultP, result.Length); + TensorPrimitives.Multiply(span2, t, spanR); + TensorPrimitives.MultiplyAdd(span1, 1 - t, spanR, spanR); + } + } + + private static unsafe void matProduct(float[,] mat1, float[,] mat2, float[,] result) + { + fixed (float* mat1P = mat1, mat2P = mat2, resultP = result) + { + var span1 = new Span(mat1P, mat1.Length); + var span2 = new Span(mat2P, mat2.Length); + var spanR = new Span(resultP, result.Length); + TensorPrimitives.Multiply(span1, span2, spanR); + } + } + + private static unsafe void matScale(float[,] mat, float scalar, float[,] result) + { + fixed (float* matP = mat, resultP = result) + { + var span1 = new Span(matP, mat.Length); + var spanR = new Span(resultP, result.Length); + TensorPrimitives.Multiply(span1, scalar, spanR); + } + } + + private static unsafe void matDiff(float[,] mat1, float[,] mat2, float[,] result) + { + fixed (float* mat1P = mat1, mat2P = mat2, resultP = result) + { + var span1 = new Span(mat1P, mat1.Length); + var span2 = new Span(mat2P, mat2.Length); + var spanR = new Span(resultP, result.Length); + TensorPrimitives.Subtract(span1, span2, spanR); + } + } + + // This matmul operation is not standard because it computes (m, p) * (n, p) -> (m, n) + // This is because the memory for the reduced dimension must be contiguous + private static unsafe void matmul(float[,] mat1, float[,] mat2, float[,] result) + { + int m = mat1.GetLength(0); + int n = mat2.GetLength(0); + int p = mat1.GetLength(1); + + for (int i = 0; i < m; i++) + { + for (int j = 0; j < n; j++) + { + fixed (float* mat1P = mat1, mat2P = mat2) + { + var span1 = new Span(mat1P + i * p, p); + var span2 = new Span(mat2P + j * p, p); + result[i, j] = TensorPrimitives.Dot(span1, span2); + } + } + } + } + + private static float[] linspace(float start, float end, int count) + { + float[] result = new float[count]; + + for (int i = 0; i < count; i++) + { + result[i] = start + (end - start) * i / (count - 1); + } + + return result; + } + + /// + /// Calculates a normalized cumulative distribution for the Euclidean distance between points on a piecewise-linear path. + /// + /// (2, n) shape array which represents the points of the piecewise-linear path. + /// n-length array to write the result to. + private static void getDistanceDistribution(float[,] points, float[] result) + { + int m = points.GetLength(1); + float accumulator = 0; + result[0] = 0; + + for (int i = 1; i < m; i++) + { + float dist = MathF.Sqrt(MathF.Pow(points[0, i] - points[0, i - 1], 2) + MathF.Pow(points[1, i] - points[1, i - 1], 2)); + accumulator += dist; + result[i] = accumulator; + } + + var spanR = result.AsSpan(); + TensorPrimitives.Divide(spanR, accumulator, spanR); + } + + private class Interpolator + { + private readonly int ny; + private readonly float[,] ys; + + public Interpolator(ReadOnlySpan inputPath, int resolution = 1000) + { + float[,] arr = new float[2, inputPath.Length]; + + for (int i = 0; i < inputPath.Length; i++) + { + arr[0, i] = inputPath[i].X; + arr[1, i] = inputPath[i].Y; + } + + float[] dist = new float[inputPath.Length]; + getDistanceDistribution(arr, dist); + ny = resolution; + ys = new float[resolution, 2]; + int current = 0; + + for (int i = 0; i < resolution; i++) + { + float target = (float)i / (resolution - 1); + + while (dist[current] < target) + current++; + + int prev = Math.Max(0, current - 1); + float currDist = dist[current]; + float prevDist = dist[prev]; + float t = (currDist - target) / (currDist - prevDist); + + if (float.IsNaN(t)) + t = 0; + + ys[i, 0] = t * arr[0, prev] + (1 - t) * arr[0, current]; + ys[i, 1] = t * arr[1, prev] + (1 - t) * arr[1, current]; + } + } + + public void Interpolate(float[] x, float[,] result) + { + int nx = x.Length; + + for (int i = 0; i < nx; i++) + { + float idx = x[i] * (ny - 1); + int idxBelow = (int)idx; + int idxAbove = Math.Min(idxBelow + 1, ny - 1); + idxBelow = Math.Max(idxAbove - 1, 0); + + float t = idx - idxBelow; + + result[0, i] = t * ys[idxAbove, 0] + (1 - t) * ys[idxBelow, 0]; + result[1, i] = t * ys[idxAbove, 1] + (1 - t) * ys[idxBelow, 1]; + } + } + } + + /// + /// Calculate a matrix of B-spline basis function values. + /// + /// The number of control points. + /// The number of points to evaluate the spline at. + /// The order of the B-spline. + /// Matrix array of B-spline basis function values. + private static float[,] generateBSplineWeights(int numControlPoints, int numTestPoints, int degree) + { + if (numControlPoints < 2) + throw new ArgumentOutOfRangeException(nameof(numControlPoints), $"{nameof(numControlPoints)} must be >=2 but was {numControlPoints}."); + + if (numTestPoints < 2) + throw new ArgumentOutOfRangeException(nameof(numTestPoints), $"{nameof(numTestPoints)} must be >=2 but was {numTestPoints}."); + + if (degree < 0 || degree >= numControlPoints) + throw new ArgumentOutOfRangeException(nameof(degree), $"{nameof(degree)} must be >=0 and <{nameof(numControlPoints)} but was {degree}."); + + // Calculate the basis function values using the Cox-de Boor recursion formula + // Generate an open uniform knot vector from 0 to 1 + float[] x = linspace(0, 1, numTestPoints); + float[] knots = new float[numControlPoints + degree + 1]; + + for (int i = 0; i < degree; i++) + { + knots[i] = 0; + knots[numControlPoints + degree - i] = 1; + } + + for (int i = degree; i < numControlPoints + 1; i++) + { + knots[i] = (float)(i - degree) / (numControlPoints - degree); + } + + // Calculate the first order basis + float[,] prevOrder = new float[numTestPoints, numControlPoints]; + + for (int i = 0; i < numTestPoints; i++) + { + prevOrder[i, (int)MathHelper.Clamp(x[i] * (numControlPoints - degree), 0, numControlPoints - degree - 1)] = 1; + } + + // Calculate the higher order basis + for (int q = 1; q < degree + 1; q++) + { + for (int i = 0; i < numTestPoints; i++) + { + // This code multiplies the previous order by equal length arrays of alphas and betas, + // then shifts the alpha array by one index, and adds the results, resulting in one extra length. + // nextOrder = (prevOrder * alphas).shiftRight() + (prevOrder * betas) + float prevAlpha = 0; + + for (int j = 0; j < numControlPoints - degree + q - 1; j++) + { + float alpha = (x[i] - knots[degree - q + 1 + j]) / (knots[degree + 1 + j] - knots[degree - q + 1 + j]); + float alphaVal = alpha * prevOrder[i, j]; + float betaVal = (1 - alpha) * prevOrder[i, j]; + prevOrder[i, j] = prevAlpha + betaVal; + prevAlpha = alphaVal; + } + + prevOrder[i, numControlPoints - degree + q - 1] = prevAlpha; + } + } + + return prevOrder; + } + + private static float[,] generateBezierWeights(int numControlPoints, int numTestPoints) + { + if (numControlPoints < 2) + throw new ArgumentOutOfRangeException(nameof(numControlPoints), $"{nameof(numControlPoints)} must be >=2 but was {numControlPoints}."); + + if (numTestPoints < 2) + throw new ArgumentOutOfRangeException(nameof(numTestPoints), $"{nameof(numTestPoints)} must be >=2 but was {numTestPoints}."); + + long[] coefficients = binomialCoefficients(numControlPoints - 1); + float[,] p = new float[numTestPoints, numControlPoints]; + + for (int i = 0; i < numTestPoints; i++) + { + p[i, 0] = 1; + float t = (float)i / (numTestPoints - 1); + + for (int j = 1; j < numControlPoints; j++) + { + p[i, j] = p[i, j - 1] * t; + } + } + + float[,] result = new float[numTestPoints, numControlPoints]; + + for (int i = 0; i < numTestPoints; i++) + { + for (int j = 0; j < numControlPoints; j++) + { + result[i, j] = coefficients[j] * p[i, j] * p[numTestPoints - i - 1, numControlPoints - j - 1]; + } + } + + return result; + } + + /// + /// Computes an array with all binomial coefficients from 0 to n inclusive. + /// + /// n+1 length array with the binomial coefficients. + private static long[] binomialCoefficients(int n) + { + long[] coefficients = new long[n + 1]; + coefficients[0] = 1; + + for (int i = 1; i < (n + 2) / 2; i++) + { + coefficients[i] = coefficients[i - 1] * (n + 1 - i) / i; + } + + for (int i = n; i > n / 2; i--) + { + coefficients[i] = coefficients[n - i]; + } + + return coefficients; + } + private static Stack bSplineToBezierInternal(ReadOnlySpan controlPoints, ref int degree) { Stack result = new Stack(); diff --git a/osu.Framework/osu.Framework.csproj b/osu.Framework/osu.Framework.csproj index de287ea885..591150e5c3 100644 --- a/osu.Framework/osu.Framework.csproj +++ b/osu.Framework/osu.Framework.csproj @@ -52,5 +52,6 @@ +