沿着等于给定切向量的二维立方贝塞尔曲线计算时间 t

Calculate time t along a 2D cubic bezier equal to a given tangent vector

我有一个由四个点定义的立方贝塞尔曲线。我需要沿着立方贝塞尔曲线找到时间 t,其中切线等于给定向量。这个问题并不像乍看起来那么简单。我将首先解释我是如何处理它的基础数学,这样您就可以找到缺陷并​​可能找到更好的解决方案。

二维立方贝塞尔曲线及其切线可以由 these equations 定义。具体切线:

T(t) = -3(1-t)^2 * P0 + 3(1-t)^2 * P1 - 6t(1-t) * P1 - 3t^2 * P2 + 6t(1-t) * P2 + 3t^2 * P3

并扩展为二维向量:

T_x(t) = -3(1-t)^2 * x0 + 3(1-t)^2 * x1 - 6t(1-t) * x1 - 3t^2 * x2 + 6t(1-t) * x2 + 3t^2 * x3
T_y(t) = -3(1-t)^2 * y0 + 3(1-t)^2 * y1 - 6t(1-t) * y1 - 3t^2 * y2 + 6t(1-t) * y2 + 3t^2 * y3

然后我们还有一个向量 (x, y) 表示我们想要找到时间 t 的切线。

这些是简单的二次方程,所以我们只需要一个方程来求解。我们可以取两者之间的叉积 (vx0 * vy1 - vy0 * vx1) 并求解 0。这会发现三次贝塞尔曲线的切线何时等于我们给定的切线向量,我们将求解 t。 (我不关心向量是否与切线相反,所以如果我们的向量是 (1, 0) 那么它也会寻找 (-1, 0))。在 Mathematica 中,使用这种叉积方法求解 t 如下所示:

Solve[(-3(1-t)^2*x0+3(1-t)^2*x1-6t(1-t)*x1-3t^2*x2+6t(1-t)*x2+3t^2*x3)*y-(-3(1-t)^2*y0+3(1-t)^2*y1-6t(1-t)*y1-3t^2*y2+6t(1-t)*y2+3t^2*y3)*x==0,t,Reals]

Mathematica 将输出:

{{t->ConditionalExpression[(x0 y-2 x1 y+x2 y-x y0+2 x y1-x y2)/(x0 y-3 x1 y+3 x2 y-x3 y-x y0+3 x y1-3 x y2+x y3)-\[Sqrt]((x1^2 y^2-x0 x2 y^2-x1 x2 y^2+x2^2 y^2+x0 x3 y^2-x1 x3 y^2+x x2 y y0-x x3 y y0-2 x x1 y y1+x x2 y y1+x x3 y y1+x^2 y1^2+x x0 y y2+x x1 y y2-2 x x2 y y2-x^2 y0 y2-x^2 y1 y2+x^2 y2^2-x x0 y y3+x x1 y y3+x^2 y0 y3-x^2 y1 y3)/(x0 y-3 x1 y+3 x2 y-x3 y-x y0+3 x y1-3 x y2+x y3)^2),(x>(x2 y-x3 y)/(y2-y3)&&x0>(x1^2 y^2-x1 x2 y^2+x2^2 y^2-x1 x3 y^2+x x2 y y0-x x3 y y0-2 x x1 y y1+x x2 y y1+x x3 y y1+x^2 y1^2+x x1 y y2-2 x x2 y y2-x^2 y0 y2-x^2 y1 y2+x^2 y2^2+x x1 y y3+x^2 y0 y3-x^2 y1 y3)/(x2 y^2-x3 y^2-x y y2+x y y3)&&y>0&&y2>y3)||(x<(x2 y-x3 y)/(y2-y3)&&x0<(x1^2 y^2-x1 x2 y^2+x2^2 y^2-x1 x3 y^2+x x2 y y0-x x3 y y0-2 x x1 y y1+x x2 y y1+x x3 y y1+x^2 y1^2+x x1 y y2-2 x x2 y y2-x^2 y0 y2-x^2 y1 y2+x^2 y2^2+x x1 y y3+x^2 y0 y3-x^2 y1 y3)/(x2 y^2-x3 y^2-x y y2+x y y3)&&y>0&&y2>y3)||(x<(x2 y-x3 y)/(y2-y3)&&x0<(x1^2 y^2-x1 x2 y^2+x2^2 y^2-x1 x3 y^2+x x2 y y0-x x3 y y0-2 x x1 y y1+x x2 y y1+x x3 y y1+x^2 y1^2+x x1 y y2-2 x x2 y y2-x^2 y0 y2-x^2 y1 y2+x^2 y2^2+x x1 y y3+x^2 y0 y3-x^2 y1 y3)/(x2 y^2-x3 y^2-x y y2+x y y3)&&y<0&&y2<y3)||(x<(x2 y-x3 y)/(y2-y3)&&y<0&&x0>(x1^2 y^2-x1 x2 y^2+x2^2 y^2-x1 x3 y^2+x x2 y y0-x x3 y y0-2 x x1 y y1+x x2 y y1+x x3 y y1+x^2 y1^2+x x1 y y2-2 x x2 y y2-x^2 y0 y2-x^2 y1 y2+x^2 y2^2+x x1 y y3+x^2 y0 y3-x^2 y1 y3)/(x2 y^2-x3 y^2-x y y2+x y y3)&&y2>y3)||(x<(x2 y-x3 y)/(y2-y3)&&y2<y3&&x0>(x1^2 y^2-x1 x2 y^2+x2^2 y^2-x1 x3 y^2+x x2 y y0-x x3 y y0-2 x x1 y y1+x x2 y y1+x x3 y y1+x^2 y1^2+x x1 y y2-2 x x2 y y2-x^2 y0 y2-x^2 y1 y2+x^2 y2^2+x x1 y y3+x^2 y0 y3-x^2 y1 y3)/(x2 y^2-x3 y^2-x y y2+x y y3)&&y>0)||(x0<(x1^2 y^2-x1 x2 y^2+x2^2 y^2-x1 x3 y^2+x x2 y y0-x x3 y y0-2 x x1 y y1+x x2 y y1+x x3 y y1+x^2 y1^2+x x1 y y2-2 x x2 y y2-x^2 y0 y2-x^2 y1 y2+x^2 y2^2+x x1 y y3+x^2 y0 y3-x^2 y1 y3)/(x2 y^2-x3 y^2-x y y2+x y y3)&&y<0&&x>(x2 y-x3 y)/(y2-y3)&&y2>y3)||(x0<(x1^2 y^2-x1 x2 y^2+x2^2 y^2-x1 x3 y^2+x x2 y y0-x x3 y y0-2 x x1 y y1+x x2 y y1+x x3 y y1+x^2 y1^2+x x1 y y2-2 x x2 y y2-x^2 y0 y2-x^2 y1 y2+x^2 y2^2+x x1 y y3+x^2 y0 y3-x^2 y1 y3)/(x2 y^2-x3 y^2-x y y2+x y y3)&&y2<y3&&x>(x2 y-x3 y)/(y2-y3)&&y>0)||(y<0&&y2<y3&&x>(x2 y-x3 y)/(y2-y3)&&x0>(x1^2 y^2-x1 x2 y^2+x2^2 y^2-x1 x3 y^2+x x2 y y0-x x3 y y0-2 x x1 y y1+x x2 y y1+x x3 y y1+x^2 y1^2+x x1 y y2-2 x x2 y y2-x^2 y0 y2-x^2 y1 y2+x^2 y2^2+x x1 y y3+x^2 y0 y3-x^2 y1 y3)/(x2 y^2-x3 y^2-x y y2+x y y3))]},

{t->ConditionalExpression[(x0 y-2 x1 y+x2 y-x y0+2 x y1-x y2)/(x0 y-3 x1 y+3 x2 y-x3 y-x y0+3 x y1-3 x y2+x y3)+\[Sqrt]((x1^2 y^2-x0 x2 y^2-x1 x2 y^2+x2^2 y^2+x0 x3 y^2-x1 x3 y^2+x x2 y y0-x x3 y y0-2 x x1 y y1+x x2 y y1+x x3 y y1+x^2 y1^2+x x0 y y2+x x1 y y2-2 x x2 y y2-x^2 y0 y2-x^2 y1 y2+x^2 y2^2-x x0 y y3+x x1 y y3+x^2 y0 y3-x^2 y1 y3)/(x0 y-3 x1 y+3 x2 y-x3 y-x y0+3 x y1-3 x y2+x y3)^2),(x>(x2 y-x3 y)/(y2-y3)&&x0>(x1^2 y^2-x1 x2 y^2+x2^2 y^2-x1 x3 y^2+x x2 y y0-x x3 y y0-2 x x1 y y1+x x2 y y1+x x3 y y1+x^2 y1^2+x x1 y y2-2 x x2 y y2-x^2 y0 y2-x^2 y1 y2+x^2 y2^2+x x1 y y3+x^2 y0 y3-x^2 y1 y3)/(x2 y^2-x3 y^2-x y y2+x y y3)&&y>0&&y2>y3)||(x<(x2 y-x3 y)/(y2-y3)&&x0<(x1^2 y^2-x1 x2 y^2+x2^2 y^2-x1 x3 y^2+x x2 y y0-x x3 y y0-2 x x1 y y1+x x2 y y1+x x3 y y1+x^2 y1^2+x x1 y y2-2 x x2 y y2-x^2 y0 y2-x^2 y1 y2+x^2 y2^2+x x1 y y3+x^2 y0 y3-x^2 y1 y3)/(x2 y^2-x3 y^2-x y y2+x y y3)&&y>0&&y2>y3)||(x<(x2 y-x3 y)/(y2-y3)&&x0<(x1^2 y^2-x1 x2 y^2+x2^2 y^2-x1 x3 y^2+x x2 y y0-x x3 y y0-2 x x1 y y1+x x2 y y1+x x3 y y1+x^2 y1^2+x x1 y y2-2 x x2 y y2-x^2 y0 y2-x^2 y1 y2+x^2 y2^2+x x1 y y3+x^2 y0 y3-x^2 y1 y3)/(x2 y^2-x3 y^2-x y y2+x y y3)&&y<0&&y2<y3)||(x<(x2 y-x3 y)/(y2-y3)&&y<0&&x0>(x1^2 y^2-x1 x2 y^2+x2^2 y^2-x1 x3 y^2+x x2 y y0-x x3 y y0-2 x x1 y y1+x x2 y y1+x x3 y y1+x^2 y1^2+x x1 y y2-2 x x2 y y2-x^2 y0 y2-x^2 y1 y2+x^2 y2^2+x x1 y y3+x^2 y0 y3-x^2 y1 y3)/(x2 y^2-x3 y^2-x y y2+x y y3)&&y2>y3)||(x<(x2 y-x3 y)/(y2-y3)&&y2<y3&&x0>(x1^2 y^2-x1 x2 y^2+x2^2 y^2-x1 x3 y^2+x x2 y y0-x x3 y y0-2 x x1 y y1+x x2 y y1+x x3 y y1+x^2 y1^2+x x1 y y2-2 x x2 y y2-x^2 y0 y2-x^2 y1 y2+x^2 y2^2+x x1 y y3+x^2 y0 y3-x^2 y1 y3)/(x2 y^2-x3 y^2-x y y2+x y y3)&&y>0)||(x0<(x1^2 y^2-x1 x2 y^2+x2^2 y^2-x1 x3 y^2+x x2 y y0-x x3 y y0-2 x x1 y y1+x x2 y y1+x x3 y y1+x^2 y1^2+x x1 y y2-2 x x2 y y2-x^2 y0 y2-x^2 y1 y2+x^2 y2^2+x x1 y y3+x^2 y0 y3-x^2 y1 y3)/(x2 y^2-x3 y^2-x y y2+x y y3)&&y<0&&x>(x2 y-x3 y)/(y2-y3)&&y2>y3)||(x0<(x1^2 y^2-x1 x2 y^2+x2^2 y^2-x1 x3 y^2+x x2 y y0-x x3 y y0-2 x x1 y y1+x x2 y y1+x x3 y y1+x^2 y1^2+x x1 y y2-2 x x2 y y2-x^2 y0 y2-x^2 y1 y2+x^2 y2^2+x x1 y y3+x^2 y0 y3-x^2 y1 y3)/(x2 y^2-x3 y^2-x y y2+x y y3)&&y2<y3&&x>(x2 y-x3 y)/(y2-y3)&&y>0)||(y<0&&y2<y3&&x>(x2 y-x3 y)/(y2-y3)&&x0>(x1^2 y^2-x1 x2 y^2+x2^2 y^2-x1 x3 y^2+x x2 y y0-x x3 y y0-2 x x1 y y1+x x2 y y1+x x3 y y1+x^2 y1^2+x x1 y y2-2 x x2 y y2-x^2 y0 y2-x^2 y1 y2+x^2 y2^2+x x1 y y3+x^2 y0 y3-x^2 y1 y3)/(x2 y^2-x3 y^2-x y y2+x y y3))]}}

Here's an image that's easier to see. 也就是说,大多数情况下都有重复的变量,所以它比看起来简单得多。 (两种条件情况相同,并且方程中的解是正例或负例,因为它求解的是二次方程)。在代码形式中,这很容易看出:

var temp1 = (tx2 - tx3) / (y2 - y3);
var temp2 = (tx1 * tx1 + tx2 * tx2 + tx2 * (ty0 + ty1 - 2 * ty2) + tx1 * (-tx2 - tx3 - 2 * ty1 + ty2 + ty3) + tx3 * (ty1 - ty0) + ty1 * ty1 - ty0 * ty2 + ty2 * ty2 + ty0 * ty3 - ty1 * (ty2 + ty3)) / (tangent.y * (tx2 - tx3 - ty2 + ty3));
console.log ('Temp1: ', temp1, ' Temp2: ', temp2);
if
(
    tangent.x < temp1 &&
    (
        tangent.y < 0 && 
        (
            x0 < temp2 && y2 < y3 ||
            x0 > temp2 && y2 > y3
        ) ||
        tangent.y > 0 &&
        (
            x0 < temp2 && y2 > y3 ||
            x0 > temp2 && y2 < y3
        )
    ) ||
    tangent.x > temp1 &&
    (
        tangent.y < 0 &&
        (
            x0 < temp2 && y2 > y3 ||
            x0 > temp2 && y2 < y3
        ) ||
        tangent.y > 0 &&
        (
            x0 < temp2 && y2 < y3 ||
            x0 > temp2 && y2 > y3
        )
    )
)
{
    var tx0ty0 = tx0 - ty0;
    var ty1tx1 = ty1 - tx1;
    var tx2ty2 = tx2 - ty2;

    var temp6 = 2 * (tx0ty0 + tx2ty2) + 4 * ty1tx1;
    var temp5 = tx0ty0 + 3 * (tx2ty2 + ty1tx1) + ty3 - tx3;
    var temp7 = temp6 * temp6 - 4 * (tx0ty0 + ty1tx1) * temp5;
    var temp3 = Math.sqrt(temp7);
    var temp4 = 2 * temp5;
    var t1 = (temp6 - temp3) / temp4;
    var t2 = (temp6 + temp3) / temp4;
}

所以我们有两个可能的时间,因为问题是二次的。 Here's an interactive example in JS. 该示例使用 (0.707, 0.707) 的硬编码切向量。 (因此在该坐标系中指向下方和右侧的矢量)。

上面的代码有问题。即使纠正不等式和平方根计算中的浮点错误,也有一些情况没有明确定义。就像当 y2 - y3 为 0 时导致除以零的情况。这也有一些微妙之处,比如在某些情况下,temp4 将具有非常接近于零的有效结果,要么产生正确的结果,要么由于浮点问题产生比预期大得多的 t1 和 t2 值。在 t1 或 t2 为 0.5 的情况下,我特别注意到了这一点。我认为将它翻转到对角线上并再次求解可能会解决一些边缘情况,但我对这种方法没有信心。

我想要的是一种久经考验的方法,可能带有代码示例,或者另一种没有奇怪边缘情况的解决方法。

这个问题可能有几种特殊情况...例如,贝塞尔弧可以形成一个尖点,甚至可以是一条直线或单个点(只需考虑所有相同的控制点)。没有特殊情况的单一直接公式是不可能的。

将问题设置为 tx*y'(t) - ty*x'(t) = 0,其中 x'(t) 被定义为 a_x*t^2 + b_x*t + c_xy 也类似)并使用 Maxima 求解,我得到了

作为一般情况的两个解决方案。

当然只有当分母不为零时它们才有效,在这种情况下方程的解简化为:

此计算的 Javascript 实现是:

tvals = []; // Array of solutions
var den = 2*ax*ty - 2*ay*tx;
if (Math.abs(den) < 1E-10) {
    var num = ax*cy - ay*cx;
    var den = ax*by - ay*bx;
    if (den != 0) {
        var t = -num / den;
        if (t >= 0 && t <= 1) tvals.push(t);
    }
} else {
    var delta = (bx*bx - 4*ax*cx)*ty*ty + (-2*bx*by + 4*ay*cx + 4*ax*cy)*tx*ty + (by*by - 4*ay*cy)*tx*tx;
    var k = bx*ty - by*tx;
    tvals = [];
    if (delta >= 0 && den != 0) {
        var d = Math.sqrt(delta);
        var t0 = -(k + d) / den;
        var t1 = (-k + d) / den;
        if (t0 >= 0 && t0 < 1) tvals.push(t0);
        if (t1 >= 0 && t1 < 1) tvals.push(t1);
    }
}

您可以在 http://raksy.dyndns.org/beztan.html or as a video in https://youtu.be/5PKQUtytrlQ

中查看有效的交互式示例