通过点坐标, 反求贝塞尔曲线对应的时间参数

2458 字

12 分钟

通过点坐标, 反求贝塞尔曲线对应的时间参数

2024-10-03

笔记

.NET

/

CSharp

/

算法

贝塞尔曲线在二维上的算法, 可以分解成两个相同的一维上的算法. 而任意阶的贝塞尔曲线都可以写作一个方程. 有方程, 就能尝试求解.

本文所有代码均已上传至 GitHub 仓库: SlimeNull/SolveBezierCurve

算法概述#

bezier2

以二阶贝塞尔曲线开始, 它有三个控制点, 构成了两个线段. 计算指定时间 t 上的点, 步骤如下:

在两个线段上, 取得进度为 t 的点, 得到两个点, 即一个线段
在一个线段上, 取得进度为 t 的点, 得到一个点, 即结果.

这里说的 “进度为 t 的点”, 则是对两点进行插值, t 为参数. 当 t 越接近 0 时, 结果越接近起始点. 当 t 越接近 1 时, 结果越接近结束点.

写成 C# 代码如下:

1
// 定义表示 Point 的结构体
2
public record struct Point(double X, double Y);
3

4
public class BezierUtilities
5
{
6
    // 对 double 进行插值的函数
7
    public static double Lerp(double a, double b, double t)
8
    {
9
        return a * (1 - t) + b * t;
10
    }
11

12
    // 对 Point 进行插值的函数
13
    public static Point Lerp(Point a, Point b, double t)
14
    {
15
        return new Point(Lerp(a.X, b.X, t), Lerp(a.Y, b.Y, t));
16
    }
17
}

同理, 三阶贝塞尔曲线由四个控制点, 也就是三个线段构成, 其计算过程如下:

在三个线段上, 取得进度为 t 的点, 得到三个点, 即两个线段
在两个线段上, 取得进度为 t 的点, 得到两个点, 即一个线段
在一个线段上, 取得进度为 t 的点, 得到一个点, 即结果.

可以看到, 不论多少阶的贝塞尔曲线, 其算法都是单调的

逻辑拆解#

通过前面的插值算法可以知道, 贝塞尔曲线不同维度上的计算结果互不影响, 也就是说, 我们完全可以将二维上的二阶贝塞尔曲线, 降维到一维. 这样, 原本由两个值构成的 “控制点”, 现在变成了纯数字.

二阶贝塞尔曲线#

那么接下来, 设一维二阶贝塞尔曲线三个控制点的值分别为 a, b, c , 时间参数为 t, 计算逻辑如下:

计算线段 a-b 之间的插值点 $d = a \times (1 - t) + b \times t$
计算线段 b-c 之间的插值点 $e = b \times (1 - t) + c \times t$
计算线段 d-e 之间的插值点, 即最终的插值点: $r = d \times (1 - t) + e \times t$

将它们展开, 最终的计算公式是:

$r = (a \times (1 - t) + b \times t) \times (1 - t) + (b \times (1 - t) + c \times t) \times t$

即:

$r = a t^2 - 2 a t + a - 2 b t^2 + 2 b t + c t^2$

对于编程中的 “动画” 所使用的贝塞尔曲线, a 为 0, c 为 1, 此时公式是:

$r = - 2 b t^2 + 2 b t + t^2$

当 b 值确定时, 这就是一个简单的一元二次方程. 将其化为标准形式:

$(- 2b + 1) t^2 + 2 b t - r = 0$

随便扔给 WolframAlpha 就能得出:

当 b 为 0.5 时, 此时解为: $t = r$
当 b 不为 0.5 时, 解: $r = \frac{b \pm \sqrt{b^2 - 2br + r}}{2b - 1}$

三阶贝塞尔曲线#

对于三阶贝塞尔曲线,

那么接下来, 设一维二阶贝塞尔曲线三个控制点的值分别为 a, b, c , 时间参数为 t, 计算逻辑如下:

计算线段 a-b 之间的插值点 $e = a \times (1 - t) + b \times t$
计算线段 b-c 之间的插值点 $f = b \times (1 - t) + c \times t$
计算线段 c-d 之间的插值点 $g = c \times (1 - t) + d \times t$
计算线段 e-f 之间的插值点 $h = e \times (1 - t) + f \times t$
计算线段 f-g 之间的插值点 $i = f \times (1 - t) + g \times t$
计算线段 h-i 之间的插值点, 即最终的插值点: $r = h \times (1 - t) + i \times t$

将它们展开, 最终的计算公式是:

$r = ((a \times (1 - t) + b \times t) \times (1 - t) + (b \times (1 - t) + c \times t) \times t) \times (1 - t) + ((b \times (1 - t) + c \times t) \times (1 - t) + (c \times (1 - t) + d \times t) \times t) \times t$

即:

$r = -a t^3 + 3 a t^2 - 3 a t + a + 3 b t^3 - 6 b t^2 + 3 b t - 3 c t^3 + 3 c t^2 + d t^3$

对于编程中的 “动画” 所使用的贝塞尔曲线, a 为 0, d 为 1, 此时公式是:

$r = 3 b t^3 - 6 b t^2 + 3 b t - 3 c t^3 + 3 c t^2 + t^3$

当 b 值确定时, 这就是一个简单的一元二次方程. 将其化为标准形式:

$(3b - 3c + 1)t^3 + (-6b + 3c)t^2 + 3bt - r = 0$

当然, 三次方程的求根公式就比二次方程的复杂多了, 在这就不写出来了. 在下一部分中给出的逻辑代码, 使用的是盛金法求三次方程.

选取合适的根#

对于贝塞尔曲线, t 的值应该总是大于等于 0, 小于等于 1 的. 所以在解方程之后, 不在此范围的根应该抛弃掉.

代码实现#

简单抽象一个 “点”:

1
public record struct Point(double X, double Y)
2
{
3
    public bool IsEmpty
4
    {
5
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
6
        get => Is(0, 0);
7
    }
8

9
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
10
    public bool Is(double x, double y) => X == x && Y == y;
11

12
    public static Point operator +(Point p1, Point p2)
13
    {
14
        return new Point(p1.X + p2.X, p1.Y + p2.Y);
15
    }
16

17
    public static Point operator -(Point p1, Point p2)
18
    {
19
        return new Point(p1.X - p2.X, p1.Y - p2.Y);
20
    }
21

22
    public static Point operator *(Point point, double factor)
23
    {
24
        return new Point(point.X * factor, point.Y * factor);
25
    }
26

27
    public static Point operator /(Point point, double factor)
28
    {
29
        return new Point(point.X / factor, point.Y / factor);
30
    }
31

32
    public static Point Lerp(Point a, Point b, double t)
33
    {
34
        return a * (1 - t) + b * t;
35
    }
36
}

一次方程, 二次方程, 三次方程的求解逻辑:

1
internal static class MathUtilities
2
{
3
    /// <summary>
4
    /// 解一元一次方程 (标准形式)
5
    /// </summary>
6
    /// <param name="a"></param>
7
    /// <param name="b"></param>
8
    /// <returns></returns>
9
    public static IEnumerable<double> SolveLinearFunction(double a, double b)
10
    {
11
        yield return -b / a;
12
    }
13

14
    /// <summary>
15
    /// 解一元二次方程 (标准形式)
16
    /// </summary>
17
    /// <param name="a">二次项系数</param>
18
    /// <param name="b">一次项系数</param>
19
    /// <param name="c">常数项</param>
20
    /// <returns></returns>
21
    public static IEnumerable<double> SolveQuadraticFunction(double a, double b, double c)
22
    {
23
        if (a == 0)
24
        {
25
            foreach (var root in SolveLinearFunction(b, c))
26
            {
27
                yield return root;
28
            }
29

30
            yield break;
31
        }
32

33
        var delta = b * b - 4 * a * c;
34
        if (delta < 0)
35
        {
36
            // 没有实数根
37
            yield break;
38
        }
39

40
        var x1 = (-b + Math.Sqrt(delta)) / (2 * a);
41
        yield return x1;
42

43
        if (delta == 0)
44
        {
45
            yield break;
46
        }
47

48
        var x2 = (-b - Math.Sqrt(delta)) / (2 * a);
49
        yield return x2;
50
    }
51

52
    /// <summary>
53
    /// 解一元三次方程 (标准形式)
54
    /// </summary>
55
    /// <param name="a">三次项系数</param>
56
    /// <param name="b">二次项系数</param>
57
    /// <param name="c">一次项系数</param>
58
    /// <param name="d">常数项</param>
59
    /// <returns></returns>
60
    public static IEnumerable<double> SolveCubicFunction(double a, double b, double c, double d)
61
    {
62
        if (a == 0)
63
        {
64
            foreach (var root in SolveQuadraticFunction(b, c, d))
65
            {
66
                yield return root;
67
            }
68

69
            yield break;
70
        }
71

72
        // 盛金公式法
73

74
        double bigA = b * b - 3 * a * c;
75
        double bigB = b * c - 9 * a * d;
76
        double bigC = c * c - 3 * b * d;
77

78
        double delta = bigB * bigB - 4 * bigA * bigC;
79

80
        if (bigA == bigB && bigA == 0)
81
        {
82
            // 三个相同的实数根
83
            yield return -c / b;
84
            yield break;
85
        }
86
        else if (delta > 0)
87
        {
88
            var bigY1 = bigA * b + 3 * a * ((-bigB + Math.Sqrt(delta)) / 2);
89
            var bigY2 = bigA * b + 3 * a * ((-bigB - Math.Sqrt(delta)) / 2);
90

91
            // 只有一个实数根, 剩下两个是虚的
92
            yield return (-b - (Math.Pow(bigY1, 1.0 / 3.0) + Math.Pow(bigY2, 1.0 / 3.0))) / (3 * a);
93
        }
94
        else if (delta == 0)
95
        {
96
            var bigK = bigB / bigA;
97

98
            yield return -b / a + bigK;
99

100
            // 两个相同的实数根
101
            yield return -bigK / 2;
102
        }
103
        else
104
        {
105
            var bigT = (2 * bigA * b - 3 * a * bigB) / (2 * Math.Pow(bigA, 3.0 / 2.0));
106
            var theta = Math.Acos(bigT);
107

108
            // 三个根
109
            yield return (-b - 2 * Math.Sqrt(bigA) * Math.Cos(theta / 3)) / (3 * a);
110
            yield return (-b + Math.Sqrt(bigA) * (Math.Cos(theta / 3) + Math.Sqrt(3) * Math.Sin(theta / 3))) / (3 * a);
111
            yield return (-b + Math.Sqrt(bigA) * (Math.Cos(theta / 3) - Math.Sqrt(3) * Math.Sin(theta / 3))) / (3 * a);
112
        }
113
    }
114
}

贝塞尔曲线抽象:

1
public class BezierCurve
2
{
3
    private readonly Point[] _points;
4

5
    public BezierCurve(IEnumerable<Point> points)
6
    {
7
        ArgumentNullException.ThrowIfNull(points);
8

9
        _points = points.ToArray();
10

11
        if (_points.Length < 1)
12
        {
13
            throw new ArgumentException("Point count can not be less than 1");
14
        }
15
    }
16

17
    private double SolveTimeForPoint(double controlPoint1, double controlPoint2, double point)
18
    {
19
        double a = -controlPoint1 + controlPoint2;
20
        double b = controlPoint1  - point;
21

22

23
        foreach (var root in MathUtilities.SolveLinearFunction(a, b))
24
        {
25
            if (root < 0 || root > 1)
26
            {
27
                continue;
28
            }
29

30
            return root;
31
        }
32

33
        return double.NaN;
34
    }
35

36
    private double SolveTimeForPoint(double controlPoint1, double controlPoint2, double controlPoint3, double point)
37
    {
38
        double a = controlPoint1 - 2 * controlPoint2 + controlPoint3;
39
        double b = -2 * controlPoint1 + 2 * controlPoint2;
40
        double c = controlPoint1 - point;
41

42
        foreach (var root in MathUtilities.SolveQuadraticFunction(a, b, c))
43
        {
44
            if (root < 0 || root > 1)
45
            {
46
                continue;
47
            }
48

49
            return root;
50
        }
51

52
        return double.NaN;
53
    }
54

55
    private double SolveTimeForPoint(double controlPoint1, double controlPoint2, double controlPoint3, double controlPoint4, double point)
56
    {
57
        double a = -controlPoint1 + 3 * controlPoint2 - 3 * controlPoint3 + controlPoint4;
58
        double b = 3 * controlPoint1 - 6 * controlPoint2 + 3 * controlPoint3;
59
        double c = -3 * controlPoint1 + 3 * controlPoint2;
60
        double d = controlPoint1 - point;
61

62
        foreach (var root in MathUtilities.SolveCubicFunction(a, b, c, d))
63
        {
64
            if (root < 0 || root > 1)
65
            {
66
                continue;
67
            }
68

69
            return root;
70
        }
71

72
        return double.NaN;
73
    }
74

75
    /// <summary>
76
    /// 传入进度 t 值, 求曲线上点坐标
77
    /// </summary>
78
    /// <param name="t"></param>
79
    /// <returns></returns>
80
    public Point Solve(double t)
81
    {
82
        Point[] buffer = new Point[_points.Length];
83
        _points.CopyTo(buffer, 0);
84

85
        int pointCount = _points.Length;
86
        while (pointCount > 1)
87
        {
88
            for (int i = 0; i < pointCount - 1; i++)
89
            {
90
                buffer[i] = Point.Lerp(buffer[i], buffer[i + 1], t);
91
            }
92

93
            pointCount--;
94
        }
95

96
        return buffer[0];
97
    }
98

99
    /// <summary>
100
    /// 根据曲线上某点的 X 值求解 t 参数
101
    /// </summary>
102
    /// <param name="pointX"></param>
103
    /// <returns></returns>
104
    /// <exception cref="NotSupportedException"></exception>
105
    public double SolveTimeForPointX(double pointX)
106
    {
107
        if (_points.Length == 2)
108
        {
109
            return SolveTimeForPoint(_points[0].X, _points[1].X, pointX);
110
        }
111
        else if (_points.Length == 3)
112
        {
113
            return SolveTimeForPoint(_points[0].X, _points[1].X, _points[2].X, pointX);
114
        }
115
        else if (_points.Length == 4)
116
        {
117
            return SolveTimeForPoint(_points[0].X, _points[1].X, _points[2].X, _points[3].X, pointX);
118
        }
119
        else
120
        {
121
            throw new NotSupportedException();
122
        }
123
    }
124

125
    /// <summary>
126
    /// 根据曲线上某点的 Y 值求解 t 参数
127
    /// </summary>
128
    /// <param name="pointY"></param>
129
    /// <returns></returns>
130
    /// <exception cref="NotSupportedException"></exception>
131
    public double SolveTimeForPointY(double pointY)
132
    {
133
        if (_points.Length == 2)
134
        {
135
            return SolveTimeForPoint(_points[0].Y, _points[1].Y, pointY);
136
        }
137
        else if (_points.Length == 3)
138
        {
139
            return SolveTimeForPoint(_points[0].Y, _points[1].Y, _points[2].Y, pointY);
140
        }
141
        else if (_points.Length == 4)
142
        {
143
            return SolveTimeForPoint(_points[0].Y, _points[1].Y, _points[2].Y, _points[3].Y, pointY);
144
        }
145
        else
146
        {
147
            throw new NotSupportedException();
148
        }
149
    }
150
}