计算三次贝塞尔曲线的拐点?
Calculating the inflection point of a cubic bezier curve?
我有四个点构成三次贝塞尔曲线:
P1 = (10, 5)
P2 = (9, 12)
P3 = (24, -2)
P4 = (25, 3)
现在我想找到这条曲线的拐点。我用谷歌搜索 and/but 每个人都指的是一个网站:http://www.caffeineowl.com/graphics/2d/vectorial/cubic-inflexion.html
不幸的是,这些小程序无法正常工作,我根本无法将它们组合在一起。有人可以教我如何计算曲线的拐点吗?
对于贝塞尔曲线,您有两个参数方程 X(t) 和 Y(t)。
要确定参数曲线的拐点,您需要找到曲线曲率(wiki)改变符号的位置。所以你需要找到上述函数的一阶和二阶导数并求解方程:
C(t) = X' * Y'' - X'' * Y' = 0
一阶导数是二次导数,二阶导数是线性导数,所以t
的方程是三次导数,最多可以有3个解。
编辑:已阅读链接文章,方程已简化为二次方程,最多可有 2 个解。
如果 t 范围 0..1 中存在解,您还必须检查它是否是真正的拐点 - 检查 C'(t) <> 0
在此 t 值。
示例:蓝色圆圈是拐点(抓到两个)
实码(Delphi)
给定:
P是控制点数组
Cf 是贝塞尔曲线的幂基系数
P, Cf: array[0..3] of TPoint
//calculate Bezier coefficients
Cf[3].x := p[3].x - 3 * p[2].x + 3 * p[1].x - p[0].x;
Cf[2].x := 3 * (p[0].x - 2 * p[1].x + p[2].x);
Cf[1].x := 3 * (p[1].x - p[0].x);
Cf[0].x := p[0].x;
//the same for Y
//find parameters of quadratic equation
// a*t^2 + b*t + c = 0
a := 3 * (cf[2].X *cf[3].Y - cf[2].Y *cf[3].X);
b := 3 * (cf[1].X *cf[3].Y - cf[1].Y *cf[3].X);
c := cf[1].X *cf[2].Y - cf[1].Y *cf[2].X;
//here solve quadratic equations, find t parameters
//don't forget a lot of special cases like a=0, D<0, D=0, t outside 0..1 range
Discriminant := b * b - 4 * a * c;
....
为了详细说明 MBo 的答案,C(t) = X' * Y'' - X'' * Y' = 0
几乎就是您想要的伪代码,因为一阶和二阶导数非常容易计算。
根据 this explanation of Bezier derivatives 和坐标为 (x1,y1)...(x4,y4) 的通用贝塞尔函数,我们得到:
fx(t) = x1 (1-t)³ + 3·x2·(1-t)²·t + 3·x3·(1-t)·t² + x4·t³
fx'(t) = a·(1-t)² + 2·b·(1-t)·t + c·t²
其中 a = 3(x2-x1),b = 3(x3-x2),c = 3(x4-x3),并且:
fx''(t) = u·(1-t) + v·t
其中 u = 2(b-a) 且 v = 2(c-b)。当然,y 分量也是如此:
fy(t) = y1 (1-t)³ + 3·y2·(1-t)²·t + 3·y3·(1-t)·t² + y4·t³
fy'(t) = a'·(1-t)² + 2·b'·(1-t)·t + c'·t²
fy''(t) = u'·(1-t) + v'·t
其中 a'
与 a
相同,但具有 y
值等。
计算 C(t) = fx'(t)·fy''(t) - fx''(t)·fy'(t)
的数学很烦人,但这就是我们拥有计算机的原因。如果您拥有 raspberry pi,则您拥有 Mathematica 的许可证,所以让我们使用它:
这是一个庞大的公式,但是找到 "arbitrary" 曲线的拐点有点傻,因为贝塞尔曲线对于线性仿射变换是不变的,所以拐点的 t
值无论我们检查 "the real curve" 还是 rotate/translate/scale 曲线都保持不变,以便它具有更方便的坐标。就像翻译它使 (x1,y1) 最终成为 (0,0),并且 (x4,y4) 位于 x 轴上,因此 y4 为零。
如果我们这样做,那么我们会得到一个简单得多的公式:
这有多简单?嗯:
18 times:
- x3 * y2
+ 3 * x3 * y2 * t
- 3 * x3 * y2 * t^2
- x4 * y2 * t
+ 2 * x4 * y2 * t^2
+ x2 * y3
- 3 * x2 * y3 * t
+ 3 * x2 * y3 * t^2
- x4 * y3 * t^2
由于我们正在编程,因此其中包含大量可缓存的值。服用:
a = x3 * y2
b = x4 * y2
c = x2 * y3
d = x4 * y3
我们可以将C(t)
简化为:
1/18 * C(t) = -a + 3at - 3at^2 - bt + 2bt^2 + c - 3ct + 3ct^2 - dt^2
= -3at^2 + 2bt^2 + 3ct^2 - dt^2 + 3at - bt - 3ct - a + c
= (-3a + 2b + 3c - d)t^2 + (3a - b - 3c)t + (c - a)
把那个因子 18 放回去,这只是一个简单的二次公式,我们可以使用具有更简单值的二次根恒等式求根:
v1 = (-3a + 2b + 3c - d) * 18
v2 = (3a - b - 3c) * 18
v3 = (c - a) * 18
并且,假设 3a +d
不等于 2b+3c
(因为如果是这样的话就没有根),我们得到:
sqr = sqrt(v2^2 - 4 * v1 * v3)
d = 2 * v1
root1 = (sqr - v2) / d
root2 = -(sqr + v2) / d
丢弃不在贝塞尔区间 [0,1] 内的根,剩下的是原始曲线拐点的 t
值。
It's just that I am not able to put things together in a way so that I end up having some nice pseudo-code. That is also why I gave some points and coordinates. I'd like to see someone calculating this with real numbers
尤其是在 Whosebug 上,最好不要一开始就懒惰,而是承诺可能必须学习新事物。
我有四个点构成三次贝塞尔曲线:
P1 = (10, 5)
P2 = (9, 12)
P3 = (24, -2)
P4 = (25, 3)
现在我想找到这条曲线的拐点。我用谷歌搜索 and/but 每个人都指的是一个网站:http://www.caffeineowl.com/graphics/2d/vectorial/cubic-inflexion.html
不幸的是,这些小程序无法正常工作,我根本无法将它们组合在一起。有人可以教我如何计算曲线的拐点吗?
对于贝塞尔曲线,您有两个参数方程 X(t) 和 Y(t)。 要确定参数曲线的拐点,您需要找到曲线曲率(wiki)改变符号的位置。所以你需要找到上述函数的一阶和二阶导数并求解方程:
C(t) = X' * Y'' - X'' * Y' = 0
一阶导数是二次导数,二阶导数是线性导数,所以t
的方程是三次导数,最多可以有3个解。
编辑:已阅读链接文章,方程已简化为二次方程,最多可有 2 个解。
如果 t 范围 0..1 中存在解,您还必须检查它是否是真正的拐点 - 检查 C'(t) <> 0
在此 t 值。
示例:蓝色圆圈是拐点(抓到两个)
实码(Delphi)
给定:
P是控制点数组
Cf 是贝塞尔曲线的幂基系数
P, Cf: array[0..3] of TPoint
//calculate Bezier coefficients
Cf[3].x := p[3].x - 3 * p[2].x + 3 * p[1].x - p[0].x;
Cf[2].x := 3 * (p[0].x - 2 * p[1].x + p[2].x);
Cf[1].x := 3 * (p[1].x - p[0].x);
Cf[0].x := p[0].x;
//the same for Y
//find parameters of quadratic equation
// a*t^2 + b*t + c = 0
a := 3 * (cf[2].X *cf[3].Y - cf[2].Y *cf[3].X);
b := 3 * (cf[1].X *cf[3].Y - cf[1].Y *cf[3].X);
c := cf[1].X *cf[2].Y - cf[1].Y *cf[2].X;
//here solve quadratic equations, find t parameters
//don't forget a lot of special cases like a=0, D<0, D=0, t outside 0..1 range
Discriminant := b * b - 4 * a * c;
....
为了详细说明 MBo 的答案,C(t) = X' * Y'' - X'' * Y' = 0
几乎就是您想要的伪代码,因为一阶和二阶导数非常容易计算。
根据 this explanation of Bezier derivatives 和坐标为 (x1,y1)...(x4,y4) 的通用贝塞尔函数,我们得到:
fx(t) = x1 (1-t)³ + 3·x2·(1-t)²·t + 3·x3·(1-t)·t² + x4·t³
fx'(t) = a·(1-t)² + 2·b·(1-t)·t + c·t²
其中 a = 3(x2-x1),b = 3(x3-x2),c = 3(x4-x3),并且:
fx''(t) = u·(1-t) + v·t
其中 u = 2(b-a) 且 v = 2(c-b)。当然,y 分量也是如此:
fy(t) = y1 (1-t)³ + 3·y2·(1-t)²·t + 3·y3·(1-t)·t² + y4·t³
fy'(t) = a'·(1-t)² + 2·b'·(1-t)·t + c'·t²
fy''(t) = u'·(1-t) + v'·t
其中 a'
与 a
相同,但具有 y
值等。
计算 C(t) = fx'(t)·fy''(t) - fx''(t)·fy'(t)
的数学很烦人,但这就是我们拥有计算机的原因。如果您拥有 raspberry pi,则您拥有 Mathematica 的许可证,所以让我们使用它:
这是一个庞大的公式,但是找到 "arbitrary" 曲线的拐点有点傻,因为贝塞尔曲线对于线性仿射变换是不变的,所以拐点的 t
值无论我们检查 "the real curve" 还是 rotate/translate/scale 曲线都保持不变,以便它具有更方便的坐标。就像翻译它使 (x1,y1) 最终成为 (0,0),并且 (x4,y4) 位于 x 轴上,因此 y4 为零。
如果我们这样做,那么我们会得到一个简单得多的公式:
这有多简单?嗯:
18 times:
- x3 * y2
+ 3 * x3 * y2 * t
- 3 * x3 * y2 * t^2
- x4 * y2 * t
+ 2 * x4 * y2 * t^2
+ x2 * y3
- 3 * x2 * y3 * t
+ 3 * x2 * y3 * t^2
- x4 * y3 * t^2
由于我们正在编程,因此其中包含大量可缓存的值。服用:
a = x3 * y2
b = x4 * y2
c = x2 * y3
d = x4 * y3
我们可以将C(t)
简化为:
1/18 * C(t) = -a + 3at - 3at^2 - bt + 2bt^2 + c - 3ct + 3ct^2 - dt^2
= -3at^2 + 2bt^2 + 3ct^2 - dt^2 + 3at - bt - 3ct - a + c
= (-3a + 2b + 3c - d)t^2 + (3a - b - 3c)t + (c - a)
把那个因子 18 放回去,这只是一个简单的二次公式,我们可以使用具有更简单值的二次根恒等式求根:
v1 = (-3a + 2b + 3c - d) * 18
v2 = (3a - b - 3c) * 18
v3 = (c - a) * 18
并且,假设 3a +d
不等于 2b+3c
(因为如果是这样的话就没有根),我们得到:
sqr = sqrt(v2^2 - 4 * v1 * v3)
d = 2 * v1
root1 = (sqr - v2) / d
root2 = -(sqr + v2) / d
丢弃不在贝塞尔区间 [0,1] 内的根,剩下的是原始曲线拐点的 t
值。
It's just that I am not able to put things together in a way so that I end up having some nice pseudo-code. That is also why I gave some points and coordinates. I'd like to see someone calculating this with real numbers
尤其是在 Whosebug 上,最好不要一开始就懒惰,而是承诺可能必须学习新事物。