math - 沿着等于给定切线向量的二维三次贝塞尔曲线计算时间 t

Question

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

二维三次贝塞尔曲线及其切线可以由这些方程定义。特别是切线：

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))]}}

这是一个更容易看到的图像。也就是说，大多数情况下都有重复的变量，所以它比看起来要简单得多。（这两种情况是相同的，并且解是方程中的正例或负例，因为它求解了一个二次方程）。在代码形式中，这很容易看到：

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;
}

因此，由于问题是二次的，因此我们所拥有的可能是我们预期的两次。这是 JS 中的一个交互式示例。该示例使用 (0.707, 0.707) 的硬编码切向量。（因此在该坐标系中指向下方和右侧的向量）。

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

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

Answer 1

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

将问题提出为tx*y'(t) - ty*x'(t) = 0wherex'(t)被定义为a_x*t^2 + b_x*t + c_x（和类似的y）并用我得到的 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或https://youtu.be/5PKQUtytrlQ中查看一个有效的交互式示例