程序師世界是廣大編程愛好者互助、分享、學習的平台,程序師世界有你更精彩!
首頁
編程語言
C語言|JAVA編程
Python編程
網頁編程
ASP編程|PHP編程
JSP編程
數據庫知識
MYSQL數據庫|SqlServer數據庫
Oracle數據庫|DB2數據庫
 程式師世界 >> 編程語言 >> C語言 >> 關於C語言 >> 判斷一個點是否在指定三角形內(1)

判斷一個點是否在指定三角形內(1)

編輯:關於C語言

問題:判斷點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_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)

 (編譯器可以利用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_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;

}

 

 

  1. 上一頁:
  2. 下一頁:
Copyright © 程式師世界 All Rights Reserved