二阶贝塞尔曲线的点击判定

2021-05-31codingbezier

引言

最近在做一个可视化项目时,有两点之间绘制链路的需求。考虑到可能有多条链路同时存在,因此选用了二阶贝塞尔曲线。但也有点击此曲线进行交互的功能,因此需要实现二阶贝塞尔曲线的点击判定。

贝塞尔曲线

贝塞尔曲线是一种曲线,用于计算机图形绘制,CSS动画等方面。贝塞尔曲线形状完全由控制点决定,有 nn 个控制点对应 n1n-1 阶贝塞尔曲线。

一阶贝塞尔曲线

在给定两个控制点P0P_0P1P_1的情况下,是两点之间直接连线构成的直线,由下列公式给出:

B(t)=P0+(P1P0)t=(1t)P0+tP1,t[0,1]B(t)=P_0+(P_1-P_0 )t=(1-t) P_0+tP_1 , t∈[0,1]

二阶贝塞尔曲线

二阶贝塞尔曲线,有首末点分别为P0P_0P2P_2的情况下,控制点为P1P_1,由下列公式给出:

B(t)=(1t)2P0+2t(1t)P1+t2P2,t[0,1]B(t)=(1-t)^2 P_0+2t(1-t) P_1+t^2 P_2 , t∈[0,1]

二阶贝塞尔曲线展示图

二阶贝塞尔曲线展示图

点到二阶贝塞尔曲线的准确距离计算

假设有一个二阶贝塞尔曲线,你在这个曲线上旁边点了一点,怎么获得点到这个曲线的最短距离呢? 图片 可以将你点击的位置视作圆心,半径向外不断扩大这个圆,直到这个圆与曲线只有一个交点时,此时这个圆与曲线相切,见上方图片位置。

设二阶贝塞尔曲线三个控制点分别为P0P_0P1P_1P2P_2,点击点位置为PiP_i,与圆相切位置为PtP_t,相切点对应tt值为t0t_0,贝塞尔曲线x轴y轴方向关于t的导数分别为X(t)X'(t)Y(t)Y'(t)。因为相切,向量乘积为0,则应满足以下方程:

X(t)(xixt)+Y(t)(yiyt)=0(1)X'(t)*(x_i-x_t) + Y'(t)*(y_i - y_t) = 0 \qquad (1)

由二阶贝塞尔曲线满足:

B(t)=(P02P1+P2)t2+2(P1P0)t+P0B(t)=(P_0-2P_1+P_2)t^2 + 2*(P_1-P_0)t+P_0
B(t)=2(P02P1+P2)t+2(P1P0)B'(t)=2*(P_0-2P_1+P_2)t+ 2*(P_1-P_0)

易知以上方程可以化为At3+Bt2+Ct+D=0At^3+Bt^2+Ct+D=0的形式,是一个一元三次方程。
P0P_0点的x轴坐标为x0x_0,y轴坐标为y0y_0,其余点同理。为了方便计算,定义如下变量:

{(a0,a1)=P1P0(b0,b1)=P02P1+P2\left\{ \begin{align*} & (a_0, a_1)=P_1-P_0 \\ & (b_0, b_1)=P_0-2P_1+P_2 \\ \end{align*} \right.

则公式 (1) 通过相关软件可以简单化简为以下公式:

(b02+b12)t3+(3a0b0+3a1b1)t2+[2a02+2a12+b0(x0xi)+b1(y0yi)]t+a0(x0xi)+a1(y0yi)=0(b_0^2 + b_1^2)*t^3 + (3*a_0*b_0 + 3*a_1*b_1)*t^2 + [2*a_0^2 + 2*a_1^2 + b_0*(x_0 - x_i) + b_1*(y_0 - y_i)]*t + a_0*(x_0 - x_i) + a_1*(y_0 - y_i) = 0

接着根据一元三次方程求根公式,可以进一步计算得出对应的tt的值,进而计算出准确距离。

代码实现

搜索可以得到二阶贝塞尔曲线距离计算相关实现实例,其代码中已进行了一些性能优化。因此将其”翻译”到js语言,源代码如下:

// The MIT License
// Copyright © 2018 Inigo Quilez
// Distance to a quadratic bezier segment, which can be solved analyically with a cubic.
// List of some other 2D distances: https://www.shadertoy.com/playlist/MXdSRf
// and www.iquilezles.org/www/articles/distfunctions2d/distfunctions2d.htm
float dot2( in vec2 v ) { return dot(v,v); }
float cross2( in vec2 a, in vec2 b ) { return a.x*b.y - a.y*b.x; }
    
// unsigned distance to a quadratic bezier
float udBezier( in vec2 pos, in vec2 A, in vec2 B, in vec2 C )
{    
    vec2 a = B - A;
    vec2 b = A - 2.0*B + C;
    vec2 c = a * 2.0;
    vec2 d = A - pos;

    float kk = 1.0/dot(b,b);
    float kx = kk * dot(a,b);
    float ky = kk * (2.0*dot(a,a)+dot(d,b))/3.0;
    float kz = kk * dot(d,a);      

    float res = 0.0;

    float p = ky - kx*kx;
    float p3 = p*p*p;
    float q = kx*(2.0*kx*kx - 3.0*ky) + kz;
    float h = q*q + 4.0*p3;

    if( h>=0.0 ) 
    {   // 1 root
        h = sqrt(h);
        vec2 x = (vec2(h,-h)-q)/2.0;
        vec2 uv = sign(x)*pow(abs(x), vec2(1.0/3.0));
        float t = clamp( uv.x+uv.y-kx, 0.0, 1.0 );
        res = dot2(d+(c+b*t)*t);
    }
    else 
    {   // 3 roots
        float z = sqrt(-p);
        float v = acos(q/(p*z*2.0))/3.0;
        float m = cos(v);
        float n = sin(v)*1.732050808;
        vec3  t = clamp( vec3(m+m,-n-m,n-m)*z-kx, 0.0, 1.0 );
        res = min( dot2(d+(c+b*t.x)*t.x),
                   dot2(d+(c+b*t.y)*t.y) );
        // the third root cannot be the closest. See https://www.shadertoy.com/view/4dsfRS
        // res = min(res,dot2(d+(c+b*t.z)*t.z));
    }
    
    return sqrt( res );
}

粗判定

考虑到一些应用场景,比如判断是否点中曲线,可以通过计算点到曲线的最短距离是否小于判定距离。但在点击点到曲线距离较远的情况下,可以明显得出是否点中曲线的判断结果,在此场景下,出于性能优化考虑,对点击点先进行粗判断,在粗判断成功后才开始计算点到曲线的最短距离。

image-20210612231511206

粗判断范围展示图

以二阶贝塞尔的首末节点连线与中垂线建立坐标系,可以得到一个正好将此贝塞尔曲线包含在内的矩形(上图内部蓝色矩形),在此矩形基础上各个边向外扩大一个判定距离,可以得到一个更大的矩形(上图外部红色矩形)。易知红色矩形外的点到此曲线的距离一定大于判定距离,判断点是否在此矩形外即可。

粗判断的部分代码实现如下:

// 根据坐标系的偏移和旋转生成坐标映射函数
function transformCoord(rotation: number, dx: number, dy: number) {
  const cos = Math.cos(rotation);
  const sin = Math.sin(rotation);
  return function transformPoint(x, y) {
    return {
      x: cos * (x - dx) + sin * (y - dy),
      y: -sin * (x - dx) + cos * (y - dy),
    };
  };
}
// 根据dx和dy计算角度
function calculateRotation(dx: number, dy: number) {
  return Math.atan2(dy, dx);
}

const sorter = (a: number, b: number) => a - b;
// 根据贝塞尔曲线创建粗判断函数
function createIsRoughHit(fromX, fromY, cpX, cpY, toX, toY) {
  const transform = transformCoord(
    calculateRotation(toX - fromX, toY - fromY),
    (fromX + toX) / 2,
    (fromY + toY) / 2
  );
  const t_from = transform(fromX, fromY);
  const t_cp = transform(cpX, cpY);
  const t_to = transform(toX, toY);
  const ratio = Math.abs(t_from.x);
  const [bottom, top] = [t_cp.y / 2, 0].sort(sorter);
  let left, right;
  // 左右侧范围可根据数学公式求出
  if (t_cp.x >= t_from.x && t_cp.x <= t_to.y) {
    [left, right] = [t_from.x, t_to.x];
  } else if (t_cp.x > 0) {
    [left, right] = [t_from.x, 0.5 * (t_cp.x + (ratio * ratio) / t_cp.x)];
  } else {
    [left, right] = [0.5 * (t_cp.x + (ratio * ratio) / t_cp.x), t_to.x];
  }
  // 返回粗判断比对函数
  return function isRoughtHit(x: number, y: number, distance: number) {
    const t_point = transform(x, y);
    if (t_point.x < left - distance || t_point.x > right + distance)
      return false;
    if (t_point.y < bottom - distance || t_point.y > top + distance)
      return false;
    return true;
  };
}

具体实现效果

github地址
codepen展示
性能暂未测试,实际使用较为出色,但应还有不少优化空间。

参考

The end