雁过无痕

  C++博客 :: 首页 :: 新随笔 :: 联系 :: 聚合  :: 管理 ::
问题:判断点P是否在三角形ABC


判断一个点是否在在三角形内,最常用的两种方法:面积法、向量同向法。算法虽然很简单,但要做到高效却不容易,要考虑到二维、三维的区别,还要考虑到坐标是用浮点数还是用整数来表示。

 

在二维平面上,问题相对简单,一般只需6次乘法计算。但在三维平面时问题要复杂很多,在网上看到的算法,一般都需要30次乘法计算(如果已知点P在平面ABC上,则需21次)。实际上,在三维坐标系下,可以做到增加1次比较,将乘法计算降到13次(如果点P在平面ABC上,则最多只要8次乘法计算)。

 

最常用的两种方法:面积法和向量同向法本质上是等价的。

向量同向法:若点P在三角形内,则三个向量:ab × apap × acpb × pc平行同向(它们也与向量ab × ac平行同向),由于这三个向量均有可能为0,直接判断它们平行同向相当麻烦,但考虑到ab × ac不可能为0,直接判断“向量:ab × apap × acpb × 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)|

 

对任意向量abc   |a + b + c| = |a| + |b| + |c| <==> 向量abc 平行同向

   因而,面积法和向量同向法本质上是等价的。

 

 

下面先讨论二维坐标系(每个点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_triangle(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_triangle(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

 (编译器可以利用cdqsar等指令直接由x计算出y值)

  abs(x)  =  (x xor y) – y

      或:   = (x + y) xor y

      或:   = x – (2 * x & y)

 

  对整数abc 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_triangle(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;

}

 


作者: flyinghearts
出处: http://www.cnblogs.com/flyinghearts/
本文采用知识共享署名-非商业性使用-相同方式共享 2.5 中国大陆许可协议进行许可,欢迎转载,但未经作者同意必须保留此段声明,且在文章页面明显位置给出原文连接,否则保留追究法律责任的权利。
posted on 2011-07-07 23:16 flyinghearts 阅读(3460) 评论(9)  编辑 收藏 引用 所属分类: 算法编程之美

评论

# re: 判断一个点是否在指定三角形内(1) 2011-07-08 10:39 zober
你这篇文章,我在某个ios blog看过 ~~  回复  更多评论
  

# re: 判断一个点是否在指定三角形内(1) 2011-07-08 11:35 千暮(zblc)
我曾经写过判断一个点是否在任意凹、凸非自交n多边形里面....3个点的话,也不用这么麻烦的- -....直接用射线奇偶法,一个return 就解决了  回复  更多评论
  

# re: 判断一个点是否在指定三角形内(1) 2011-07-08 15:34 陈梓瀚(vczh)
@千暮(zblc)
不行的,三维里面的一个点,有可能在三角形所在的平面外面,你那个方法不管用。  回复  更多评论
  

# re: 判断一个点是否在指定三角形内(1) 2011-07-08 20:13 千暮(zblc)
@陈梓瀚(vczh)
诶诶,人家是2D,不是3D。  回复  更多评论
  

# re: 判断一个点是否在指定三角形内(1) 2011-07-08 22:55 flyinghearts
@zober
原创文章。若你看过相似的,请给出链接。
  回复  更多评论
  

# re: 判断一个点是否在指定三角形内(1) 2011-07-08 22:56 flyinghearts
@千暮(zblc)
请贴出你的代码,谢谢。
  回复  更多评论
  

# re: 判断一个点是否在指定三角形内(1) 2011-07-09 01:23 千暮(zblc)
@flyinghearts
http://www.cppblog.com/zblc/archive/2009/12/23/103822.html
不是独立的代码,而是被封装进了一个09年我写的数学子库里,包含任意凹凸非自交n多边形检测(当然里面公式没完全压缩化简)。
不好意思哈,太晚了,关于一行return我就不鸟你了嘿嘿(算我当时没考虑美学观点:)。
  回复  更多评论
  

# re: 判断一个点是否在指定三角形内(1) 2011-07-11 22:49 flyinghearts
@千暮(zblc)
你的代码有点复杂呀。
 对凸多边形,射线法 比 矢量同向法,要复杂很多。
 参考我的实现代码: http://www.cppblog.com/flyinghearts/archive/2011/07/11/150716.html
  回复  更多评论
  

# re: 判断一个点是否在指定三角形内(1) 2011-07-11 23:48 千暮(zblc)
@flyinghearts
是的,所以你的方法如果化简下代码,用于凸边形更实用,因为我封装进库的时候,不只是为了实现多边形这一种功能,而是相关数学子库设计:顺便把这个功能作为Demo而已。注意Demo的代码就那么点  回复  更多评论
  


只有注册用户登录后才能发表评论。
【推荐】超50万行VC++源码: 大型组态工控、电力仿真CAD与GIS源码库
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理