根据贝塞尔曲线上的点反算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、P一、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、前端可视化方面有深刻研究。对程序员思惟能力训练和培训、程序员职业规划有浓厚兴趣。