检测一个点是否在三角形内

880 阅读1分钟

检测一个点是否在三角形内

1. 夹角求和

1.png

ABC\triangle_{ABC}中有一点PP,很容易得到如下关系:

APB+BPC+CPA=2π\angle_{APB} + \angle_{BPC} + \angle_{CPA} = 2 \pi

但是,PP点在三角形外时:

2.png

{APB+BPC+CPA=2APBAPB<πAPB+BPC+CPA<2π\left\{ \begin{aligned} \angle_{APB} + \angle_{BPC} + \angle_{CPA} &= 2 \angle_{APB} \\ \angle_{APB} &< \pi \end{aligned} \right. \Rightarrow \angle_{APB} + \angle_{BPC} + \angle_{CPA} < 2 \pi

实现起来很简单:

function sumCheck(triangle, p) {
  // 求出三个向量
  const pA = new THREE.Vector3().subVectors(triangle.a, p);
  const pB = new THREE.Vector3().subVectors(triangle.b, p);
  const pC = new THREE.Vector3().subVectors(triangle.c, p);

  // 计算三个角度
  const angleAPB = pA.angleTo(pB);
  const angleBPC = pB.angleTo(pC);
  const angleCPA = pC.angleTo(pA);

  // 由于浮点计算误差 我们将差距小于1e-10的认为是相等
  return Math.abs(angleAPB + angleBPC + angleCPA - Math.PI * 2) < 1e-10
}

关于向量之间的夹角,我们这里直接使用了Three.js定义好的方法,其实就是通过向量“点积”求出余弦值,然后通过反三角函数求出角度。这里就不再赘述了。

这个方法的好处就是算法简单,便于理解。但是,在求夹角的时候使用了反三角函数,这是比较耗时,接下来介绍一个更快的算法。

2. 同侧检测

1.png

上图中,点PPABC\triangle_{ABC}中的充要条件是:

  1. PPABAB的下侧
  2. PPBCBC的左侧
  3. ppACAC的右侧

但是,上下左右这样的描述并不准确,也不通用;我们将它转换为更通用的描述方式:

  1. AP\overrightarrow{AP}AB\overrightarrow{AB}的顺时针方向
  2. BP\overrightarrow{BP}BC\overrightarrow{BC}的顺时针方向
  3. CP\overrightarrow{CP}CA\overrightarrow{CA}的顺时针方向

好像看起来通用了,但是,如果是一个空间中的三角形,从一个方向我们上面的说法没有问题;从另一个方向看,如果要保证正确,我们就要把上面的说法全改成逆时针。

怎么知道它在那个方向? 当我们比较PPABAB时,我们可以利用CC,如果PP在三角形内,PPCC一定是同侧的!!! 所以,我们再次修改我们上面的算法:

  1. AP\overrightarrow{AP}AC\overrightarrow{AC}AB\overrightarrow{AB}的同侧
  2. BP\overrightarrow{BP}BA\overrightarrow{BA}BC\overrightarrow{BC}的同侧
  3. CP\overrightarrow{CP}CB\overrightarrow{CB}CA\overrightarrow{CA}的同侧

那用数学怎么判断同侧呢?向量叉乘可以判断方向 要判断AP\overrightarrow{AP}AC\overrightarrow{AC}是不是在AB\overrightarrow{AB}的同一侧,只需要将AP\overrightarrow{AP}AC\overrightarrow{AC}分别与AB\overrightarrow{AB}做向量积,如果它们结果的方向相同,那就是在同一侧。即,direction(AP×AB)direction(\overrightarrow{AP} \times \overrightarrow{AB})是否等于direction(AP×AB)direction(\overrightarrow{AP} \times \overrightarrow{AB})

下面我们用代码实现它:


// 判断向量v1,v2是不是在向量base的同一侧
function sameSide(v1, v2, base) {
  const direction1 = new THREE.Vector3().crossVectors(v1, base);
  const direction2 = new THREE.Vector3().crossVectors(v2, base);

  // 如果方向相同 点积后的结果一定大于零
  return direction1.dot(direction2) > 0;
}

function sameSideCheck(triangle, p) {

  const _AB = new THREE.Vector3().subVectors(triangle.b, triangle.a);
  const _BC = new THREE.Vector3().subVectors(triangle.c, triangle.b);
  const _CA = new THREE.Vector3().subVectors(triangle.a, triangle.c);
  const _BA = _AB.clone().negate();
  const _AC = _CA.clone().negate();
  const _CB = _BC.clone().negate();
  
  const _AP = new THREE.Vector3().subVectors(p, triangle.a);
  const _BP = new THREE.Vector3().subVectors(p, triangle.b);
  const _CP = new THREE.Vector3().subVectors(p, triangle.c);

  return sameSide(_AP, _AC, _AB) && sameSide(_BP, _BA, _BC) && sameSide(_CP, _CB, _CA);
}

相比于之前的夹角求和的方式,我们这里只是用了向量的内积和外积,运算速度更快。

3. 重心坐标

ABCA,B,C三个点可以确定一个平面,平面上一点PP可以通过A,B,CA,B,C的线性组合表示:

设,AP=uAB+vAC则,PA=u(BA)+v(CA)化简得,P=(1uv)A+uB+vC\begin{aligned} 设,& \overrightarrow{AP} = u \overrightarrow{AB} + v \overrightarrow{AC} \\ 则,& P - A = u(B - A) + v(C - A) \\ 化简得,& P = (1 - u - v)A + uB + vC \end{aligned}

即,P=αA+βB+γCα+β+γ=1P = \alpha A + \beta B + \gamma C,\alpha + \beta + \gamma = 1

坐标(α,β,γ)(\alpha, \beta, \gamma)PP关于A,B,CA,B,C的重心坐标。

那为什么叫重心坐标呢?我们看一下计算重心的方式:

设,有n个质点,质点的坐标位置和质量分别为PiP_imi(i=1...n)m_i(i=1...n);则,这些质点的重心的坐标为:

P=i=1nmiPii=1nmiP = \frac{\sum\limits_{i=1}^n m_i \cdot P_i}{\sum\limits_{i=1}^n m_i}

从形式上看,就是将质量当作权值的一个加权平均数:

P=m1P1i=1nmi+m1P1i=1nmi+...+mnPni=1nmi记,w1=m1i=1nmi,w2=m2i=1nmi,...,wn=m1i=1nmi则,w1+w2+...+wn=1\begin{aligned} &P = \frac{m_1 \cdot P_1}{\sum\limits_{i=1}^n m_i} + \frac{m_1 \cdot P_1}{\sum\limits_{i=1}^n m_i} + ... + \frac{m_n \cdot P_n}{\sum\limits_{i=1}^n m_i} \\ &记, w_1=\frac{m_1}{\sum\limits_{i=1}^n m_i}, w_2=\frac{m_2}{\sum\limits_{i=1}^n m_i}, ..., w_n=\frac{m_1}{\sum\limits_{i=1}^n m_i} \\ &则, w_1+w_2+...+w_n = 1 \end{aligned}

n=3n=3时,正是我们上面推导的公式:

P=αA+βB+γCα+β+γ=1P = \alpha A + \beta B + \gamma C,\alpha + \beta + \gamma = 1

(w1,w2,w3)(w_1, w_2,w_3)就是重心坐标(α,β,γ)(\alpha,\beta,\gamma)

重心坐标在三角形内的判断条件

我们把上面的公式变一下形式:P=A+uAB+vACP = A + u \overrightarrow{AB} + v \overrightarrow{AC}

从几何上解释这个公式就是:

  1. AA为起点
  2. 沿着AB\overrightarrow{AB}的方向走了uu
  3. 沿着AC\overrightarrow{AC}的方向走了vv

3.png

u,vu,v分别在什么范围内可以保证PP不“走出”三角形呢?

如果u,vu,v有一个小于00,那一定落在三角形外面。

5.png

如果u,vu,v有一个大于11,那一定落在三角形外面

6.png

现在,我们知道了u,vu,v都在[0,1][0, 1]区间内,接下来讨论一下00 如果它们其中有一个为00PP点就会落在ABABACAC边所在的直线上,另一个必须是[0,1][0, 1]区间内的一个值,才能落在三角形的边上。

4.png 当它们都为00时,PP就落在了AA点。

关于在边上是不是属于三角形内,这个怎么认为都可以,我们这里就认为边上的点在三角形内。

因此,我们得到了u,vu,v都必须在[0,1][0, 1]区间内。

现在我们还差最后一个坐标的范围,它很简单,就是1(u+v)1-(u + v),它的取值范围是什么呢?不防先来看下u+vu+v

u+v=1u + v = 1时,就表示P=(1uv)A+uB+vC=uB+(1u)CP=(1-u-v)A+uB+vC=uB+(1-u)C

P=uB+(1u)CP=uB+(1-u)C这是什么?正是BCBC边,也就是说,当u+v=1u + v=1时,PP点必定落在BCBC边上。一旦u+v>1u + v > 1,那就“走出”了三角形。

因此,我们可以确定0<u+v<10<1(u+v)<10< u + v <1 \Rightarrow 0 < 1 - (u + v) < 1

综上,我们可以得出,当PP的重心坐标的三个值都属于[0,1][0, 1]区间,则PP点在三角形内。

所以,一旦我们知道了重心坐标,判断起来会相当简单。

function baryCoordCheck(triangle, p) {
  const baryCoord = getBaryCoordCheck(triangle, p);

  return baryCoord.x >= 0 && baryCoord.y >= 0 && baryCoord.z <= 1
}

现在的问题是,怎么求重心坐标?

重心坐标的求法

PA=u(CA)+v(BA)令,v0=CAv1=BAv2=PA得,v2=uv0+vv1\begin{aligned} P - A &= u(C - A) + v(B - A) \\ 令, \overrightarrow{v_0} &= C-A \\ \overrightarrow{v_1} &= B-A \\ \overrightarrow{v_2} &= P-A \\ 得, \overrightarrow{v_2}&=u\overrightarrow{v_0} + v\overrightarrow{v_1} \end{aligned}
{v2v0=(uv0+vv1)v0v2v1=(uv0+vv1)v1{v2v0=u(v0v0)+v(v0v1)v2v1=u(v0v1)+v(v1v1)\begin{aligned} &\left\{ \begin{aligned} \overrightarrow{v_2} \cdot \overrightarrow{v_0} &= (u\overrightarrow{v_0} + v\overrightarrow{v_1}) \cdot \overrightarrow{v_0} \\ \overrightarrow{v_2} \cdot \overrightarrow{v_1} &= (u\overrightarrow{v_0} + v\overrightarrow{v_1}) \cdot \overrightarrow{v_1} \end{aligned} \right. \\ \Rightarrow &\left\{ \begin{aligned} \overrightarrow{v_2} \cdot \overrightarrow{v_0} &= u(\overrightarrow{v_0} \cdot \overrightarrow{v_0}) + v(\overrightarrow{v_0} \cdot \overrightarrow{v_1}) \\ \overrightarrow{v_2} \cdot \overrightarrow{v_1} &= u(\overrightarrow{v_0} \cdot \overrightarrow{v_1}) + v(\overrightarrow{v_1} \cdot \overrightarrow{v_1}) \\ \end{aligned} \right. \end{aligned}

两个未知数u,vu,v可以通过求解上述方程组得到:

u=(v1v1)(v2v0)(v1v0)(v2v1)(v0v0)(v1v1)(v0v1)(v1v0)v=(v0v0)(v2v1)(v0v1)(v2v0)(v0v0)(v1v1)(v0v1)(v1v0)u = \frac{(v_1 \cdot v_1)(v_2 \cdot v_0)-(v_1 \cdot v_0)(v_2 \cdot v_1)}{(v_0 \cdot v_0)(v_1 \cdot v_1) - (v_0 \cdot v_1)(v_1 \cdot v_0)} \\ v = \frac{(v_0 \cdot v_0)(v_2 \cdot v_1)-(v_0 \cdot v_1)(v_2 \cdot v_0)}{(v_0 \cdot v_0)(v_1 \cdot v_1) - (v_0 \cdot v_1)(v_1 \cdot v_0)}

Three.js使用的就是这样的算法,我们直接看源码:

// src/math/Triangle.js

static getBarycoord( point, a, b, c, target ) {
  // point P点坐标
  // a, b, c为三角形的三个顶点
  // 将结果保存到target中

  _v0.subVectors( c, a );
  _v1.subVectors( b, a );
  _v2.subVectors( point, a );

  const dot00 = _v0.dot( _v0 );
  const dot01 = _v0.dot( _v1 );
  const dot02 = _v0.dot( _v2 );
  const dot11 = _v1.dot( _v1 );
  const dot12 = _v1.dot( _v2 );

  const denom = ( dot00 * dot11 - dot01 * dot01 );

  // collinear or singular triangle
  if ( denom === 0 ) {

    // arbitrary location outside of triangle?
    // not sure if this is the best idea, maybe should be returning undefined
    return target.set( - 2, - 1, - 1 );

  }

  const invDenom = 1 / denom;
  const u = ( dot11 * dot02 - dot01 * dot12 ) * invDenom;
  const v = ( dot00 * dot12 - dot01 * dot02 ) * invDenom;

  // barycentric coordinates must always sum to 1
  return target.set( 1 - u - v, v, u );

}

同样,检查一个点是否在三角形内,也有对应得实现:

// src/math/Triangle.js
static containsPoint( point, a, b, c ) {

  this.getBarycoord( point, a, b, c, _v3 );

  return ( _v3.x >= 0 ) && ( _v3.y >= 0 ) && ( ( _v3.x + _v3.y ) <= 1 );
}

验证代码

我们这里借助canvas写一段模拟代码,用来验证我们的函数。 首先,我们需要定义我们的三角形。我们假定视窗大小为8,然后用canvas将这个尺寸为8的视窗画出来。

为了坐标的统一,需要将坐标除以视窗大小,规范化到[0,1][0, 1]

参考文章

Point in triangle test (blackpawn.com)