计算三次贝塞尔曲线的拐点？

Question

我有四个点可以组成一个三次贝塞尔曲线：

P1 = (10, 5)
P2 = (9, 12)
P3 = (24, -2)
P4 = (25, 3)

现在我想找出这条曲线的拐点。我谷歌搜索了一下，但是每个人都只提到一个网站：http://www.caffeineowl.com/graphics/2d/vectorial/cubic-inflexion.html。不幸的是，这些小应用程序无法正常工作，我也无法将它们整合到一起。请问是否有人能够友情提示如何计算我的曲线的拐点呢？

- user1602492

Answer 1

对于贝塞尔曲线，你有两个参数方程X(t)和Y(t)。 要确定参数曲线的拐点，需要找到曲率（wiki）改变符号的位置。因此，需要找到上述函数的一阶和二阶导数，并解方程式：

C(t) = X' * Y'' - X'' * Y' = 0

第一阶导数是二次的，第二阶导数是线性的，因此方程对于 t 是三次的，可能有最多三个解。

编辑：已经阅读了链接的文章，方程被简化为二次方程，可能有最多两个解。

如果在 t 范围 0..1 内存在解，则还必须检查它是否为实际的拐点-检查在该 t 值处 C'(t) <> 0。

例如：蓝色圆圈是拐点（捕捉到两个）。

真实的代码片段（Delphi）
给定：
P是控制点数组
Cf是Bezier在幂基础上的系数

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

---

你是否熟悉参数曲线表示法？你能否将贝塞尔曲线从伯恩斯坦基转换为多项式基表示法？ - MBo

Answer 2

为了进一步解释MBo的答案，C(t) = X' * Y'' - X'' * Y' = 0 实际上已经是你需要的伪代码，因为计算一阶和二阶导数非常容易。

根据Bezier导数的说明，对于具有坐标(x1,y1)...(x4,y4)的通用Bezier函数，我们得到如下结果:

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)的数学公式很麻烦，但这就是为什么我们拥有电脑。如果你拥有树莓派，那么你也有Mathematica的许可证，让我们把它用起来：

这是一个巨大的公式，但是对于“任意”曲线找到拐点有些愚蠢，因为贝塞尔曲线在线性仿射变换下不变，所以拐点的t值无论我们是检查“真实曲线”还是旋转/平移/缩放曲线，使其具有更方便的坐标，都保持不变。例如将其平移到(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值。

问题在于我无法将事物组合在一起，以便最终得到一些漂亮的伪代码。这也是为什么我给出了一些点和坐标。我想看到有人用实际数字计算这个。

在Stackoverflow上，最好不要懒惰，而是承诺可能需要学习新事物。

- Mike 'Pomax' Kamermans

---

我已经用Python编写了代码，但是结果不正确。这是我还是漏掉了什么，还是有打字错误？？http://hastebin.com/osujiwagix.py - user1602492

是的，显然这些输入将失败。请重新阅读上文中的段落：“如果我们这样做，那么我们就可以得到一个简单得多的公式”，因为您构建代码的分析需要您先轴向对齐曲线。 - Mike 'Pomax' Kamermans

现在我明白了，事情变得很明显。谢谢你，Mike！ - user1602492

贝塞尔曲线如果有三个维度怎么办？如何计算三维立方贝塞尔曲线的拐点？ - Ctrl_Alt_Del

为什么有这么多问号？一个就足够了。话虽如此：你需要找到所有平面曲率变化符号的点，这很容易计算“对于任何给定的点”，但不是你可以在几步内符号地解决整个曲线的问题，所以你只能用数值分析（在编程中意味着：编写一个for循环来采样曲线，找到转换对，并进行细化）。 - Mike 'Pomax' Kamermans