根据贝塞尔曲线上的点反算t值

这是一个项目中遇到的实际需求。
场景是一个智能仓库管理系统，场景里面有直线和曲线构成的环穿轨道。环穿轨道上面会有小车运动，后台推动小车的两个点位A和B，其中A和B都会在轨道上面，前端需要根据这两个推送点，自动播放小车从A点沿轨道到B点的动画。下面是项目截图：

项目中使用的是二次贝塞尔曲线，所以本文也主要以二次贝塞尔曲线为讲解重点。
要实现上述动画，需要首先确定A点和B点在曲线上面的比例值t_a和t_b

最终的需求变成：“根据贝塞尔曲线上的点反算t值”。 大概有以下几种方法。现假设贝塞尔曲线上的点为点P（后续会用到该点）。

分片迭代

分片迭代是一种近似的方法。我们知道，二次贝塞尔曲线的公式如下：
B(t) = (1-t)2 P0 + 2t(1-t) P1 + t2 * P2
其中： $t \in $[0,1]，P0为二次贝塞尔曲线的起始点，P1为控制点，P2为终止点。

如果你对于上面的知识点不是很熟悉，建议学习贝塞尔曲线相关知识。推荐学习本人的专栏Canvas高级进阶, 里面有专门的章节对贝塞尔曲线进行了全面详细的讲解。本文也是从该专栏的文章中摘录并适当改编而成的。

从以上公式，我们可以得到，对于任意给定的比例值t，可以求出对应该比例值的点B(t)。分片迭代思路是：现在加设把范围[0,1]平均分成N（比如100）等份，形成一系列的比例值t，对于每一个t值，求取对应的点B(t) ，然后让点B(t)和已知在贝塞尔曲线上的点P进行比较，如果点B(t)和点P之间的直线距离在一定的误差范围之内，则认为B(t)等于P，而此时的t值，就是我们要求的t值。
以下是主要代码：

function computeT(p0,p1,p2,p) {
  var t = 0;
  for(var i = 0;i < 1000;i ++){
      var point = getPointOnQuadraticCurve(p0,p1,p2,t);//根据二次贝塞尔曲线公式求B(t)，其中point = B(t)
      if(distance(point,p) < 0.01){ // 判断point和p点的距离是否在特定误差之内
        return t;
      }
      t+= 0.001;
  }
  return null;
}

上述分片迭代的方法，思路最简单，最直观。在精度要求不高的情况下是可以满足的。而在精度要求高的时候，即代码中的“特定误差”值要很小，可能会出现函数返回值为null的情况，在精度要求高的时候要能够计算出值，就要增加迭代次数，此时会极大增加性能消耗。比如上面代码的迭代次数可能会变成10000甚至10000。

迭代方法同样适用于三次贝塞尔曲线和更加高阶的贝塞尔曲线。

分片迭代优化版本

上面提到在精度要求高的情况下，要得到正确结果，要极大的增加迭代次数，造成性能的极大消耗。 有没有办法既提高精度，又不大量增加迭代次数呢？ 经过笔者的思考，发现是可以的。想想假设要求的t值在0.5附近，那么我们只需要在0.5附近加大分片的数量，而不需要在其他地方（0.1~0.4，0.6~1.0）增加分片的数量。 应此升级版本的思路就是，先用比较粗的分片初步确定t值的一个大致范围，再在该范围之类，比较细的分片确定t值。注意这是个递归的过程，如果在第二次比较细的分片情况下，仍然不能确定t值，那么就确定一个t值的更小分范围；重复上面过程，直到找到t值为止。
大致步骤如下：

• 首先，通过一个小的迭代次数进行分片迭代；
• 在迭代的过程中如果找到了符合的比例值t，直接返回；
• 在迭代的过程中同时记录离目标点P最近的t值，如果上一步未找到符合的t值，则进行下一步操作。
• 上一步找到了离目标点P最近的t值，在t值的附近(t - step,t + step)(其中step为上一次分片的步进值)进行分片迭代查找，在迭代的过程中如果找到了符合的比例值t，直接返回。
• 如果没找到，重复上面的不断缩小范围并加大分片精度的过程。 直到找到t值为止。

下面是示例代码：

function computeT(p0, p1, p2, p,startT = 0,endT = 1) {
  var t = startT;
  var minDistance = Infinity,
      minDistanceT = null;
  var step = (endT - startT) / 100;
  for (var i = 0; i < 100; i++) {
    var point = getPointOnQuadraticCurve(p0, p1, p2, t);
    var dst = distance(point,p);
    if (dst < minDistance) {
      minDistance = dst;
      minDistanceT = t;
    }
    if (dst < 0.0001) {
      return t;
    }
    t += step;
  }
  return computeT(p0, p1, p2, p, minDistanceT - step,minDistanceT + step);
}

以上过程虽然增加了一定的迭代次数，但是是常量级别的增加，而非数量级别的增加，所以会极大提高性能。 比如目标t值在0.5附近，第一次通过100次迭代可以确定t值的范围在0.4 ~ 0.6之间；然后进行第二次迭代，第二次迭代此次数仍然为100次，假设确定t值的范围在0.51 ~ 0.53之间；然后进行第三次迭代，第三次迭代此次数仍然为100次，此时可以获取t值为0.516，可以看出最多值迭代了300次。 假设总共经过第N次迭代，每次迭代次数为M，才找到t值，那么总共的迭代次数是N * M。

该迭代方法同样适用于三次贝塞尔曲线和更加高阶的贝塞尔曲线。而且相对于未优化的版本，该方法的性能好了很多。是适合所有贝塞尔曲线的比较好的反算t值的方法。

二分法

二分法的思路是：

• 首先确定一个起始t值和结束t值t₀和t₁，初始值t₀ = 0，t₁ = 1。
• 取t₀和t₁的中间值t_m = (t₀+t₁)/2
• 通过t_m计算出点P_m，如果P_m和目标点P之间的距离在误差值范围之内，则t_m为需要计算的目标t值。
• 如果上一步P_m和目标点P之间的距离不在误差值范围之内，则判断P_m和目标点P的前后顺序，如果P_m在目标点P的前面，则把t_m赋值给t₁；否则把t_m赋值给t₀。
• 重复以上步骤直到找到合适的t_m值。

上述步骤有一个难点： 如何判断P_m和目标点P的前后顺序？
对于二次贝塞尔曲线，如下图所示：

其中，P0为起始点，P2为终止点，P1为控制点。 二次贝塞尔曲线有如下特点：
线段(P1，P0)、(P1，P2)和曲线相切，这也就意味着曲线一定在三角形(P0，P1，P2)之内，而且二次贝塞尔曲线本身不会自身相交，所有我们可以有如下结论，

对于曲线上面的点A，直线(P1,A)和线段(P0,P1)相交于点a；对于曲线上面的点B，直线(P1,B)和线段(P0,P1)相交于点b。点A和点B的先后顺序与点a和点b的先后顺序是一致的，而直线上面的点（a和b）的前后顺序是容易判断的。 也就是说如果点a在点b的前面，则点A也在点B的前面，反之亦然。如下图所示：

有了以上的结论，我们就找到了判断P_m和目标点P的前后顺序的方法。

如果你对上述结论不熟悉，建议学习贝塞尔曲线的相关知识，推荐学习本人的专栏Canvas高级进阶, 里面有专门的章节对贝塞尔曲线进行了全面详细的讲解。本文也是从该专栏的文章中摘录并适当改编而成的。

有了这个方法，加上前面描述的二分查找的步骤，可以得出示例代码如下：

function computeT2(p0,p1,p2,p,startT = 0,endT = 1) {
   var halfT  = (startT + endT) / 2;
   var halfPoint = getPointOnQuadraticCurve(p0,p1,p2,halfT);
    if(distance(halfPoint,p) < 0.0001){
      return halfT;
   }
   //求交点：
   var inter1 = segmentsIntr(p0,p2,p1,p);
   var inter2 = segmentsIntr(p0,p2,p1,halfPoint);
  
   var r1 = interpolationRate(p0,inter1,p2),
       r2 = interpolationRate(p0,inter2,p2);
   
   if(r1 > r2){
       startT = halfT;  
   }else {
     endT = halfT;
   }
   return computeT2(p0,p1,p2,p,startT,endT);
}

解方程

前面说过，贝塞尔曲线的公式如下：
B(t) = (1-t)2 P0 + 2t(1-t) P1 + t2 * P2
其中： $t \in $[0,1]，P0为二次贝塞尔曲线的起始点，P1为控制点，P2为终止点。
分别表示成x和y的方程，则可以表示如下：

• x_P = (1-t)² x_P0 + 2t(1-t) x_P1 + t² * x_P2
• y_P = (1-t)² y_P0 + 2t(1-t) y_P1 + t² * y_P2

实际上就是两个变量t的二次元方程，取上面任意一个方程，带入相关的值解方程，方程的解即为我们要求的目标t值。

整理方程: x_P = (1-t)² x_P0 + 2t(1-t) x_P1 + t² * x_P2，可以得出二次方程如下：
(x_P2 + x_P0 - 2 x_P1 ) t² + 2(x_P1 - x_P0) t + (x_P0 - x_P) = 0。
我们已知二次方程的： at² + b t + c = 0的解为：

• 如果a = 0，则解为 -c/b
• 如果a != 0，解如下图所示：

应此令：

• a = (x_P2 + x_P0 - 2 * x_P1 )
• b = 2*(x_P1 - x_P0)
• c = (x_P0 - x_P)

可以方便求出方程的解。

需要注意的是，二次方程的解可能会有两个。如果求出的解有两个怎么办呢。 首先我们知道贝塞尔曲线的t值的范围是$t \in $[0,1]，所以如果有两个解：

• 其中一个不再[0,1]的范围之内，则另外一个解就是目标t值。（注意不可能两个都不在[0,1]范围之内，因为我们知道，目标点P在贝塞尔曲线上面）。
• 如果两个解都在[0,1]范围之内，那就把两个解再带入贝塞尔曲线的公式，分别求出两个B(t)点，那个离目标点P近，就取那个解。

下面是示例代码,其中函数equation2用于解曲线的方程：

function computeT(p0,p1,p2,p) {
    let interpolationx =  (p1.x - p0.x) / (p2.x - p0.x);
    let tt;
    if(interpolationx >= 0 && interpolationx <= 1){
      let ty = equation2(p0.y,p1.y,p2.y,p.y);
      return ty;
    }else{
      tt = equation2(p0.x,p1.x,p2.x,p.x);
      if(tt.tt1){
        var pointTest = getPointOnQuadraticCurve(p0,p1,p2,tt.tt1);
        if(distance(pointTest,p) < 0.01){
          return tt.tt1;
        }else{
          return tt.tt2;
        }
      }else{
        return tt;
      }
    }
}

function equation2(z0,z1,z2,zp){ // z0、z1，z2代表P0、P1、P2的x坐标值或者y坐标值，zp代表目标点P的x坐标值或者y坐标值
    var a = z0 - z1 * 2 + z2,
      b = 2*(z1 - z0),
      c = z0 - zp;
  var tt = null;
  if(a == 0 && b != 0){
    tt = - c / b;
  } else {
    var sq = Math.sqrt( b * b - 4 * a * c );
    var tt1 = (sq - b)/ (2 * a),
        tt2 = (-sq - b) / (2 * a);
    // console.log("tt1,tt2:",tt1,tt2); 
    if((tt1 <= 1 && tt1>= 0) && (tt2 <= 1 && tt2>= 0)){
      return {tt1,tt2};
    }else if(tt1 <= 1 && tt1>= 0){
      tt = tt1;
    }else {
      tt = tt2;
    }
  }
  return tt;
}

几种方法的比较

从性能方面来说：

• 解方程的方式是最快的
• 二分法和分片迭代的优化版次之
• 原始的分片迭代方法最差

从通用性来说，分片迭代的方式是适合任意阶的贝塞尔曲线。但是考虑到性能问题所以分片迭代的优化版是通用性最好的求解方法。

欢迎关注公众号“ITman彪叔”。彪叔，拥有10多年开发经验，现任公司系统架构师、技术总监、技术培训师、职业规划师。在计算机图形学、WebGL、前端可视化方面有深入研究。对程序员思维能力训练和培训、程序员职业规划有浓厚兴趣。