三次贝塞尔曲线 - 给定 X 得到 Y - 控制点 X 增加的特殊情况

Cubic bezier curves - get Y for given X - special case where X of control points is increasing

我已阅读 a few discussions regarding finding Y at X for a cubic Bezier curve, and have also read this article 关于此事的内容。

我的情况比一般情况更受限制,我想知道是否有比上述讨论中提到的一般情况更好的解决方案。

我的情况:

有没有考虑到这些约束的有效算法?

如果二分搜索太复杂,仍然有一个 O(1) 方法,但它相当有限。我假设您使用的是 4 个控制点 (p0(x0,y0),p1(x1,y1),p2(x2,y2),p3(x3,y3)) 三次贝塞尔曲线,由区间 [0.0 , 1.0] 中的某些 t 参数化,因此:

t = 0.0 -> x(t) = x0, y(t) = y0;
t = 1.0 -> x(t) = x3, y(t) = y3;

首先让我们暂时忘掉 Beziers 并使用 a catmull-rom curve 代替,这只是表示相同曲线的另一种方法。要在 2 个立方体之间进行转换,请使用这些:

// BEzier to Catmull-Rom
const double m=6.0;
X0 = x3+(x0-x1)*m; Y0 = y3+(y0-y1)*m;
X1 = x0;           Y1 = y0;
X2 = x3;           Y2 = y3;
X3 = x0+(x3-x2)*m; Y3 = y0+(y3-y2)*m;

// Catmull-Rom to Bezier
const double m=1.0/6.0;
x0 = X1;           y0 = Y1;
x1 = X1-(X0-X2)*m; y1 = Y1-(Y0-Y2)*m;
x2 = X2+(X1-X3)*m; y2 = Y2+(Y1-Y3)*m;
x3 = X2;           y3 = Y2;

其中 (xi,yi) 是 Bezier 控制点,(Xi,Yi) 是 Catmull-Rom 点。现在如果所有控制点之间的 X 距离具有相同的距离:

(X3-X2) == (X2-X1) == (X1-X0)

X坐标与t成线性关系。这意味着我们可以直接从 X:

计算 t
t = (X-X1)/(X2-X1);

现在我们可以直接计算任何 XY。因此,如果您可以选择控制点,则选择它们以满足 X 距离条件。

如果不满足条件,您可以尝试更改控制点使其满足(通过二进制搜索,通过将立方体细分为更多块等...)但要注意更改控制点会改变形状一不小心就会出现曲线。

首先:这个答案之所以有效,是因为您的控制点约束意味着我们总是在处理与普通函数等效的参数。在这种情况下,这显然是您想要的,但是将来找到此答案的任何人都应该知道,此答案基于 假设 任何给定的只有一个 y 值x 值...

这对于一般的贝塞尔曲线来说绝对不是这样

话虽如此,我们知道——即使我们将这条曲线表示为二维参数曲线——我们正在处理一条曲线,就所有意图和目的而言,它必须具有某种形式的未知函数y = f(x)。我们还知道,只要我们知道唯一属于特定 x 的“t”值(这只是因为您严格单调递增的系数 属性),我们就可以将 y 计算为 y = By(t),所以问题是:给定一些已知的 x 值,我们能否计算出需要插入 By(t)t 值?

答案是:是的,我们可以。

首先,我们开始的任何 x 值都可以说来自 x = Bx(t),所以假设我们知道 x,我们应该能够找到相应的值t 的 (s) 导致 x.

让我们看看 x(t) 的函数:

x(t) = a(1-t)³ + 3b(1-t)²t + 3c(1-t)t² + dt³

我们可以将其重写为简单的多项式形式:

x(t) = (-a + 3b- 3c + d)t³ + (3a - 6b + 3c)t² + (-3a + 3b)t + a

这是一个标准的三次多项式,只有已知常数作为系数,我们可以简单地将其重写为:

(-a + 3b- 3c + d)t³ + (3a - 6b + 3c)t² + (-3a + 3b)t + (a-x) = 0

您可能想知道“其他值 a、b、c 和 d 的所有 -x 去哪儿了?”答案是它们都抵消了,所以我们 实际上 最终需要减去的唯一一个是最后一个。好用!

所以现在我们只是...解决这个方程:我们知道除了 t 之外的一切,我们只需要一些数学洞察力来告诉我们如何做到这一点。

...当然,“just”在这里不是正确的限定词,求三次函数的根没有什么“just”,但值得庆幸的是,Gerolano Cardano 为确定追溯到 16th 世纪,使用复数。在任何人甚至发明复数之前。相当了不起!但这是一个编程答案,而不是历史课,所以让我们开始实施吧:

给定 x 的一些已知值,以及坐标 a、b、c 和 d 的知识,我们可以按如下方式实现求根:

// Find the roots for a cubic polynomial with bernstein coefficients
// {pa, pb, pc, pd}. The function will first convert those to the
// standard polynomial coefficients, and then run through Cardano's
// formula for finding the roots of a depressed cubic curve.
double[] findRoots(double x, double pa, double pb, double pc, double pd) {
  double
    pa3 = 3 * pa,
    pb3 = 3 * pb,
    pc3 = 3 * pc,
    a = -pa  +   pb3 - pc3 + pd,     
    b =  pa3 - 2*pb3 + pc3, 
    c = -pa3 +   pb3, 
    d =  pa  -     x;

  // Fun fact: any Bezier curve may (accidentally or on purpose)
  // perfectly model any lower order curve, so we want to test 
  // for that: lower order curves are much easier to root-find.
  if (approximately(a, 0)) {
    // this is not a cubic curve.
    if (approximately(b, 0)) {
      // in fact, this is not a quadratic curve either.
      if (approximately(c, 0)) {
        // in fact in fact, there are no solutions.
        return new double[]{};
      }
      // linear solution:
      return new double[]{-d / c};
    }
    // quadratic solution:
    double
      q = sqrt(c * c - 4 * b * d), 
      b2 = 2 * b;
    return new double[]{
      (q - c) / b2, 
      (-c - q) / b2
    };
  }

  // At this point, we know we need a cubic solution,
  // and the above a/b/c/d values were technically
  // a pre-optimized set because a might be zero and
  // that would cause the following divisions to error.

  b /= a;
  c /= a;
  d /= a;

  double
    b3 = b / 3,
    p = (3 * c - b*b) / 3, 
    p3 = p / 3, 
    q = (2 * b*b*b - 9 * b * c + 27 * d) / 27, 
    q2 = q / 2, 
    discriminant = q2*q2 + p3*p3*p3, 
    u1, v1;

  // case 1: three real roots, but finding them involves complex
  // maths. Since we don't have a complex data type, we use trig
  // instead, because complex numbers have nice geometric properties.
  if (discriminant < 0) {
    double
      mp3 = -p/3,
      r = sqrt(mp3*mp3*mp3), 
      t = -q / (2 * r), 
      cosphi = t < -1 ? -1 : t > 1 ? 1 : t, 
      phi = acos(cosphi), 
      crtr = crt(r), 
      t1 = 2 * crtr;
    return new double[]{
      t1 * cos(phi / 3) - b3,
      t1 * cos((phi + TAU) / 3) - b3,
      t1 * cos((phi + 2 * TAU) / 3) - b3
    };
  }

  // case 2: three real roots, but two form a "double root",
  // and so will have the same resultant value. We only need
  // to return two values in this case.
  else if (discriminant == 0) {
    u1 = q2 < 0 ? crt(-q2) : -crt(q2);
    return new double[]{
      2 * u1 - b3,
      -u1 - b3
    };
  }

  // case 3: one real root, 2 complex roots. We don't care about
  // complex results so we just ignore those and directly compute
  // that single real root.
  else {
    double sd = sqrt(discriminant);
    u1 = crt(-q2 + sd);
    v1 = crt(q2 + sd);
    return new double[]{u1 - v1 - b3};
  }
}

好的,这就是一大堆代码,还有很多额外的东西:

  • crt() 是立方根函数。在这种情况下,我们实际上并不关心复数,因此更简单的实现方法是使用 def、宏、三元或任何 shorthand 您选择的语言提供的任何东西:crt(x) = x < 0 ? -pow(-x, 1f/3f) : pow(x, 1f/3f);.
  • tau 就是 2π。做几何编程的时候有它在身边很有用。
  • approximately 是一个函数,它将一个值与目标 周围的非常小的间隔 进行比较,因为 IEEE 浮点数字是 jerks.基本上我们在谈论 approximately(a,b) = return abs(a-b) < 0.000001 什么的。

其余部分应该是不言自明的,如果有点 java 风格(我使用 Processing 来处理这类事情)。

有了这个实现,我们可以编写我们的实现来找到给定 x 的 y。这比调用函数要复杂一些,因为立方根是很复杂的事情。您最多可以返回三个根。但是其中只有一个适用于我们的“t区间” [0,1]:

double x = some value we know!
double[] roots = getTforX(x);
double t;
if (roots.length > 0) {
  for (double _t: roots) {
    if (_t<0 || _t>1) continue;
    t = _t;
    break;
  }
}

就是这样,我们完成了:我们现在有了“t”值,我们可以使用它来获取关联的“y”值。