5

我已经阅读了 一些 关于在 X 处为三次贝塞尔曲线找到 Y 的讨论，并且还阅读了有关此事的这篇文章。

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

我的情况：

  • X不同控制点的值在增加。即： X3 > X2 > X1 > X0。
  • 此外，由于上述原因，X(t)也是严格单调递增的。

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

4

2 回答 2

12

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

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

话虽如此，我们知道——即使我们已经将这条曲线表示为二维的参数曲线——我们正在处理的曲线对于所有意图和目的都必须具有某种形式的未知函数y = f(x)。我们还知道，只要我们知道唯一属于特定 x 的“t”值（由于您严格单调递增的系数属性，这只是这种情况），我们可以将 y 计算为y = By(t)，所以问题是：我们可以给定一些已知值，计算t我们需要插入的值？By(t)x

答案是：是的，我们可以。

首先，x我们开始的任何值都可以说来自x = Bx(t)，因此，如果我们知道x，我们应该能够找到导致 的对应t值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，我们只需要一些数学洞察力来告诉我们如何做到这一点。

...当然“公正”在这里不是正确的限定词，找到三次函数的根并不是“公正”的，但幸运的是，杰罗拉诺·卡尔达诺（ Gerolano Cardano）早在 16世纪就为确定根奠定了基础，使用复数。在任何人甚至发明复数之前。相当的壮举！但这是一个编程答案，而不是历史课，所以让我们开始实施：

给定 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、宏、三进制或您选择的语言提供的任何简写：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 风格（我正在使用处理这些事情）。

通过这个实现，我们可以编写我们的实现来找到给定 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”值。

于 2018-08-16T18:32:48.240 回答
1

如果二分搜索太复杂，仍然有一种O(1)方法，但它相当有限。p0(x0,y0),p1(x1,y1),p2(x2,y2),p3(x3,y3)我假设您使用的是 4 个控制t点[0.0 , 1.0]（

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

首先让我们暂时忘掉贝塞尔曲线，改用catmull-rom 曲线，这只是表示同一曲线的另一种方式。要在 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。这意味着我们可以t直接计算X：

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

现在我们可以直接计算Y任意值X。因此，如果您可以选择控制点，请选择它们，使它们满足 X 距离条件。

如果不满足条件，您可以尝试更改控制点以使其满足（通过二分搜索，通过将三次细分为更多块等），但要注意如果不小心更改控制点可能会改变结果曲线的形状.

于 2018-08-16T16:33:40.270 回答