判断一个点是否在指定三角形内（1）


问题：判断点P是否在三角形ABC内 
 
判断一个点是否在在三角形内，最常用的两种方法：面积法、向量同向法。算法虽然很简单，但要做到高效却不容易，要考虑到二维、三维的区别，还要考虑到坐标是用浮点数还是用整数来表示。
 
在二维平面上，问题相对简单，一般只需6次乘法计算。但在三维平面时问题要复杂很多，在网上看到的算法，一般都需要30次乘法计算（如果已知点P在平面ABC上，则需21次）。实际上，在三维坐标系下，可以做到增加1次比较，将乘法计算降到13次（如果点P在平面ABC上，则最多只要8次乘法计算）。
 
最常用的两种方法：面积法和向量同向法本质上是等价的。
向量同向法：若点P在三角形内，则三个向量：ab × ap、ap × ac、pb × pc平行同向（它们也与向量ab × ac平行同向），由于这三个向量均有可能为0，直接判断它们平行同向相当麻烦，但考虑到ab × ac不可能为0，直接判断“向量：ab × ap、ap × ac、pb × pc均与ab × ac平行同向”反而更简单。
 
面积法：当点p在三角形abc内时，4个三角形的面积满足： abc = abp + apc + pbc
对面积的计算，可以通过向量的向量积计算得到： 面积 abc = |ab × ac| / 2
表面上，要计算4个三角形的面积，但根据下面的公式：
           ap × ap = 0, pb × pc = (ab - ap) × (ac - ap) = ab × ac - ab × ap - ap × ac
可以少算一次矢量积。
公式： |ab × ac| = |ab × ap| + |ap × ac| + |(ab × ac - ab × ap - ap × ac)| 
 
对任意向量a、b、c：   |a + b + c| = |a| + |b| + |c| <==> 向量a、b、c 平行同向
   因而，面积法和向量同向法本质上是等价的。
 
 
下面先讨论二维坐标系（每个点X，都看作是原点O到该点X的二维向量OX）。
先定义一个二维向量模板：
 
template<typename T> class Vec2 {
 T x, y;
public:
 typedef T value_type;
 Vec2(T xx = 0, T yy = 0) : x(xx), y(yy) {};
 T cross(const Vec2& v) const { return x * v.y - y * v.x;} // 矢量积
 Vec2 operator-(const Vec2& v) const { return Vec2(x - v.x, y - v.y); }
};
 
 
如果坐标采用浮点数，考虑到浮点数取绝对值方便（有专门的浮点指令），但彼此间比较大小存在误差，采用面积法比较方便：
 
typedef Vec2<double> Vd2;
 
bool is_in_易做图(const Vd2& a, const Vd2& b, const Vd2& c, const Vd2& p)
{
 Vd2 ab(b -a), ac(c - a), ap(p - a);
 
 //用矢量积计算面积，下面4个值的绝对值，是对应的三角形的面积的两倍，
 double abc = ab.cross(ac); 
 double abp = ab.cross(ap);
 double apc = ap.cross(ac); 
 double pbc = abc - abp - apc;   //等于pb.cross(pc)
 //面积法：4个三角形的面积差 等于 0
 double delta = fabs(abc) - fabs(abp) - fabs(apc) - fabs(pbc); 
 return fabs(delta) < DBL_EPSILON;         
}
 
 
如果坐标采用整数表示，代码相对麻烦点：
 
typedef Vec2<int> Vi2;
 
bool is_in_易做图(const Vi2& a, const Vi2& b, const Vi2& c, const Vi2& p)
{
 Vi2 ab(b -a), ac(c - a), ap(p - a);
 
 //用矢量积计算面积，下面4个值的绝对值，是对应的三角形的面积的两倍，
 int abc = ab.cross(ac); 
 int abp = ab.cross(ap);
 int apc = ap.cross(ac); 
 int pbc = abc - abp - apc;   //等于pb.cross(pc)
 
 //方法1： 面积法：4个三角形的面积差 等于 0
 return abs(abc) == abs(abp) + abs(apc) + abs(pbc)
 
 //方法2： 矢量同向法： abp apc pbc 均与 abc 同向：
 if (abc < 0) { abp = -abp; apc = -apc; pbc = -pbc; }
 return (abp >= 0) & (apc >= 0) & (pbc >= 0); 
}
 
方法1：要计算4次绝对值，看似需要4次条件跳转，但主流的编译器，都能采用位运算直接计算绝对值（注意：GCC需要加额外的参数），不需要任何条件跳转。
方法2：比方法1指令少，但多1次条件跳转。
哪种方法效率较高，与编译器生成的具体代码有关。
 
上面代码中，可采用的两种优化方法：
①　对整数x取绝对值，可以利用位运算：
    设 y = 0 （当x >= 0）
         = -1 （当x < 0） 
 （编译器可以利用cdq或sar等指令直接由x计算出y值）
 则 abs(x)  =  (x xor y) – y
      或：   = (x + y) xor y
      或：   = x – (2 * x & y)
 
 ② 对整数a、b、c， a >= 0 && b >= 0 && c >= 0 等价于 
(a >= 0) & (b >= 0) & (c >= 0) 等价于：
(a | b | c) >= 0
 
 为避免编译器没有进行相关优化，直接手动优化，可得：
 
inline int chg_sign(int x, int sign) //sign只能取0或-1，函数分别返回x、-x
{
 return (x + sign) ^ sign; 
 //return (x ^ sign) - sign; 
}
 
bool is_in_易做图(const Vi2& a, const Vi2& b, const Vi2& c, const Vi2& p)
{
 Vi2 ab(b -a), ac(c - a), ap(p - a);
 
 //用矢量积计算面积，下面4个值的绝对值，是对应的三角形的面积的两倍，
 int abc = ab.cross(ac); 
 int abp = ab.cross(ap);
 int apc = ap.cross(ac); 
 int pbc = abc - abp - apc;   //等于pb.cross(pc)
 
 //方法3： 矢量同向法（优化版）
 const int sign = (abc >= 0) - 1; 
 //const int sign = abc >> (sizeof(abc) * CHAR_BIT - 1);
 return (chg_sign(abp, sign) | chg_sign(apc, sign) | chg_sign(pbc, sign)) >= 0;
}

补充：软件开发 , C语言 ,