我已经阅读了a few discussions关于在X处找到Y为立方Bezier曲线,并且还读了关于此问题的article。
我的情况比一般情况更受限制,我想知道是否有比上述讨论中提到的一般解决方案更好的解决方案。
我的情况:
X
值正在增加。即:X3 > X2 > X1 > X0
。X(t)
也严格单调增加。是否有任何有效的算法将这些约束考虑在内?
首先:这个答案才有效,因为你的控制点约束意味着我们总是在处理正常函数的参数等价物。在这种情况下,这显然是你想要的,但是将来找到这个答案的任何人都应该意识到这个答案是基于假设任何给定的x值只有一个y值。当然,对于Bezier曲线来说,绝对不是这样。
话虽如此,我们知道 - 即使我们已经将这条曲线表示为两个维度的参数曲线 - 我们处理的曲线对于所有意图和目的都必须具有y = f(x)
形式的一些未知函数。我们也知道,只要我们知道唯一属于特定x的“t”值(由于你的严格单调增加的系数属性就是这种情况),我们可以将y计算为y = By(t)
,所以问题是:can我们计算了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
您可能想知道“-x对于其他值a,b,c和d的位置是什么?”并且答案是它们都取消了,所以我们实际上最终需要减去的唯一一个是最后一个。便利!
所以现在我们只是...解决这个等式:我们知道除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()
是cuberoot功能。在这种情况下,我们实际上并不关心复数,因此更容易实现这一点的方法是使用def,或宏,或三元,或者您选择的语言提供的简写:crt(x) = x < 0 ? -pow(-x, 1f/3f) : pow(x, 1f/3f);
。tau
只有2π。在进行几何编程时,有用它。approximately
是一个函数,它将值与目标周围的非常小的间隔进行比较,因为IEEE浮点数字是混蛋。基本上我们谈论的是approximately(a,b) = return abs(a-b) < 0.000001
或其他什么。其余的应该是相当不言自明的,如果有一点java风格(我正在使用Processing这些东西)。
通过这个实现,我们可以编写我们的实现来找到y,给定x。这比调用函数更复杂,因为立方根是复杂的东西。你最多可以获得三个回归。但只有其中一个适用于我们的“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”值。
如果二进制搜索太复杂,仍然有一个O(1)
方法,但它相当有限。我假设你正在使用一个4控制点(p0(x0,y0),p1(x1,y1),p2(x2,y2),p3(x3,y3)
)立方贝塞尔参数化t
区间的一些[0.0 , 1.0]
所以:
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
是线性的。这意味着我们可以直接从t
计算X
:
t = (X-X1)/(X2-X1);
现在我们可以直接为任何Y
计算X
。因此,如果您可以选择控制点,则选择它们以使它们符合X距离条件。
如果不满足条件,您可以尝试更改控制点以使其满意(通过二分搜索,通过将立方体细分为更多补丁等...)但要注意,如果不小心,更改控制点可以更改结果曲线的形状。