登山第三梯：判断点与三角形的关系——不是三角恋

算法介绍：

判断一个点是否在三角形内部，主要有以下四种算法：

1、面积法：若点P在三角形ABC内部，基于四点构建的三个三角形面积之和等于三角形ABC的面积，否则，在三角形外部。

若S△ABC=S△PAB+S△PAC+S△PBC，则点P在三角形内部，否则，在外部。

其中，计算三角形面积可以使用海伦公式和向量法。

（1）海伦公式（以△ABC为例）：

首先，计算出三条边的边长a,b,c；接着使用海伦公式计算其面积：

S△ABC=sqrt(l*(l-a)*(l-b)*(l-c))，其中,l=(a+b+c)*0.5。

（2）向量法（以△ABC为例）：

利用相邻两条边的向量的叉积的模，可计算得到三角形面积

S△ABC=0.5*||ABxAC||

2、内角和法：若点P在三角形ABC内部，基于点P和三个顶点构建三条线段PA、PB、PC，该三条线段与三角形ABC的边的夹角之和等于180°，否则，在三角形外部。

若∠PAB+∠PBA+∠PBC+∠PCB+∠PCA+∠PAC=180°，则点P在三角形内部，否则，在外部。

其中，角度可通过向量的点积进行计算，以∠PAB为例：

cos(∠PAB)=AB*AP/(||AB||*||AP||)，∠PAB=arcos(∠PAB)*180/π。

3、同侧法：若点P在三角形ABC内部，则点P均在三角形边的同一侧。

如上图所示，P点均在向量AB、BC、AC的右侧，则P点在三角形内部；

那如何判断P点在向量的哪一侧？通过观察，发现点P和点C均在向量AB的同一侧。

那么，这个问题就可以转换为判断点P和点C是否在AB的同一侧。

该问题最终可以通过判断两组向量的叉积是否方向一致判断两点是否在同一侧，而叉积结果方向一致可以通过二者的点积来判断。

以AB为向量，然后构建PA和CA，若(ABxPA)*(ABxCA)>0，则两点均在同一侧。

同样的判断在以另外两条边作为参考边也满足，即点P在三角形内部，否则，点P在外部。

4、重心法：若点P在三角形ABC内部，以任意两条边的向量作为基向量，分别乘上系数u和v可表示点P，其中，u≥0，v≥0，u+v≤1

假设以点A为原点，则三角形内部任何一点P均可表示为：

AP=u*AB+v*AC，其中u≥0，v≥0，u+v≤1。

如上图所示，u=0,v=0，点P在A点，u=1，v=0，点P在B点，u=0,v=1，点P在C点

两个未知数需要两个方程，则可以两边分别乘上AC和AB构建两个方程：

（1）AP*AB=u*AB*AB+v*AC*AB

（2）AP*AC=u*AB*AC+v*AC*AC

最后解的：

u=((AC*AC)(AP*AB)-(AC*AB)(AP*AC))/((AB*AB)(AC*AC)-(AB*AC)(AC*AB))

v=((AB*AB)(AP*AC)-(AB*AC)(AP*AB))/((AB*AB)(AC*AC)-(AB*AC)(AC*AB))

效率对比：

100w点判断耗时统计表
算法	面积法	内角和法	同侧法	重心法
耗时(ms)	31	215	19	11

代码提供：

面积法：

/**
* @brief:基于海伦公式计算三角形面积.
* @param[in]:a
* @param[in]:b
* @param[in]:c   输入的为边长
* @return: 返回面积
*/
float computeTriangleAreaWithHeiLen(const float a, const float b, const float c)
{
       float p = 0.5 * (a + b + c);
       float area = sqrt(p * (p - a) * (p - b) * (p - c));
       return area;
}
/**
* @brief:利用面积法判断点是否在三角形内部.
* @param[in]:tri
* @param[in]:pt
* @return: 返回值为-1，则在外部，为0在角点，为1在内部
*/
int judgePtInTriangleWithArea(std::vector<HType::Point3f>& tri, HType::Point3f&  pt)
{
       if (3 != tri.size())
       {
              return -2;
       }
       float zeroEps(0.0000001f);
       HType::Point3f PA = tri[0] - pt;
       float pa = PA.norm();
       HType::Point3f PB = tri[1] - pt;
       float pb = PB.norm();
       HType::Point3f PC = tri[2] - pt;
       float pc = PC.norm();
       if (pa < zeroEps || pb < zeroEps || pc < zeroEps)
       {
              //与角点重合
              return 0;
       }
       HType::Point3f AB = tri[1] - tri[0];
       float a = AB.norm();
       HType::Point3f BC = tri[2] - tri[1];
       float b = BC.norm();
       HType::Point3f CA = tri[0] - tri[2];
       float c = CA.norm();
       float area_ABC = computeTriangleAreaWithHeiLen(a,b,c);
       float area_PAB = computeTriangleAreaWithHeiLen(a, pa, pb);
       float area_PBC = computeTriangleAreaWithHeiLen(b, pb, pc);
       float area_PCA = computeTriangleAreaWithHeiLen(c, pa, pc);
       if((area_PAB+ area_PBC+ area_PCA)-area_ABC>zeroEps)
       {
              //在外部
              return -1;
       }
       return 1;
}

内角和法：

/**

* @brief:利用内角和法判断点是否在三角形内部.

* @param[in]:tri

* @param[in]:pt

* @return: 返回值为-1，则在外部，为0在角点，为1在内部

*/

int judgePtInTriangleWithInteriorAngleSum(std::vector<HType::Point3f>& tri,  HType::Point3f& pt)

{

       if (3 != tri.size())

       {

              return -2;

       }

       float zeroEps(0.0000001f);

       HType::Point3f PA = tri[0] - pt;

       float pa = PA.norm();

       HType::Point3f PB = tri[1] - pt;

       float pb = PB.norm();

       HType::Point3f PC = tri[2] - pt;

       float pc = PC.norm();

       if (pa < zeroEps || pb < zeroEps || pc < zeroEps)

       {

              //与角点重合

              return 0;

       }

       HType::Point3f AB = tri[1] - tri[0];

       float a = AB.norm();

       HType::Point3f BC = tri[2] - tri[1];

       float b = BC.norm();

       HType::Point3f CA = tri[0] - tri[2];

       float c = CA.norm();

       float angle_PAB = acos(AB.dot(-PA) / (a * pa)) * 180.0/M_PI;

       float angle_PBA = acos((-AB).dot(-PB)/(a * pb)) * 180.0 / M_PI;

       float angle_PBC = acos(BC.dot(-PB)/(b * pb)) * 180.0 / M_PI;

       float angle_PCB = acos((-BC).dot(-PC)/(b*pc)) * 180.0 / M_PI;

       float angle_PCA = acos(CA.dot(-PC) / (c * pc)) * 180.0 / M_PI;

       float angle_PAC = acos((-CA).dot(-PA)/(c*pa)) * 180.0 / M_PI;

       float angle_all = angle_PAB + angle_PBA + angle_PBC + angle_PCB + angle_PCA  + angle_PAC;

       if (angle_all - 180 > zeroEps)

       {

              return -1;

       }

       return 1;

}

同侧法：

/**
* @brief:利用同侧法判断点是否在三角形内部.
* @param[in]:tri
* @param[in]:pt
* @return: 返回值为-1，则在外部，为0在角点，为1在内部
*/
int judgePtInTriangleWithSameSide(std::vector<HType::Point3f>& tri,  HType::Point3f& pt)
{
       if (3 != tri.size())
       {
              return -2;
       }
       float zeroEps(0.0000001f);
       HType::Point3f PA = tri[0] - pt;
       HType::Point3f PB = tri[1] - pt;
       HType::Point3f PC = tri[2] - pt;
       if (PA.norm() < zeroEps || PB.norm() < zeroEps || PC.norm() < zeroEps)
       {
              //与角点重合
              return 0;
       }
       //以AB为基准
       HType::Point3f AB = tri[1] - tri[0];
       HType::Point3f CA = tri[0] - tri[2];
       float retAB = (AB.cross(CA)).dot(AB.cross(PA));
       if (retAB < 0) return -1;
       //以BC为基准
       HType::Point3f BC = tri[2] - tri[1];
       float retBC = (BC.cross(AB)).dot(BC.cross(PB));
       if (retBC < 0) return -1;
       //以CA为基准
       float retCA = (CA.cross(BC)).dot(CA.cross(PC));
       if (retCA < 0) return -1;
       return 1;
}

 

重心法：

/**

* @brief:利用重心法判断点是否在三角形内部.

* @param[in]:tri

* @param[in]:pt

* @return: 返回值为-1，则在外部，为0在角点，为1在内部

*/

int judgePtInTriangleWithCentroid(std::vector<HType::Point3f>& tri,  HType::Point3f& pt)

{

       if (3 != tri.size())

       {

              return -2;

       }

       float zeroEps(0.0000001f);

       HType::Point3f PA = tri[0] - pt;

       HType::Point3f PB = tri[1] - pt;

       HType::Point3f PC = tri[2] - pt;

       if (PA.norm() < zeroEps || PB.norm() < zeroEps || PC.norm() < zeroEps)

       {

              //与角点重合

              return 0;

       }

       auto t1 = PA[0] * PB[1] - PA[1] * PB[0];

       auto t2 = PB[0] * PC[1] - PB[1] * PC[0];

       auto t3 = PC[0] * PA[1] - PC[1] * PA[0];

       if (t1 * t2 >= 0 && t2 * t3 >= 0 && t3 * t1 >= 0)

       {

              //在三角形边上或里面

              return 1;

       }

       return -1;

}