二維平面上判斷點是否在三角形內


最近在項目中碰到的這個問題,在此記錄一下。已知三角形的三個頂點坐標,判斷某個點是否在三角形中(在三角形的邊上,我們也視作在三角形中),本文給出了三種方法。

image

 

算法1

利用面積法,如上圖所示,如果點P在三角形ABC的內部,則三個小三角形PAB, PBC, PAC的面積之和 = ABC的面積,反之則不相等。

已知三角形的三個頂點坐標求其面積,可以根據向量的叉乘,參考here

該算法詳見后面代碼中的函數:IsPointInTriangle1

 

算法2

首先看一下這個問題,如何判斷某兩個點在某條直線的同一側(代碼中函數:IsPointAtSameSideOfLine)?

image

根據向量的叉乘以及右手螺旋定則,AB^AM (^表示叉乘,這里向量省略了字母上面的箭頭符號)的方向為向外指出屏幕,AB^AN也是向外指出屏幕,但AB^AO的方向是向內指向屏幕,因此M,N在直線AB的同側,M ,O在直線AB的兩側。實際計算時,只需要考慮叉積的數值正負

假設以上各點坐標為A(0,0), B(4,0), M(1,2), N(3,4), O(3,-4), 則:

AB^AM = (4,0)^(1,2) = 4*2 - 0*1 = 8

AB^AN = (4,0)^(3,4) = 4*4 – 0*3 = 16

AB^AO = (4,0)^(3,-4) = 4*-4 – 0*3 = –16

由上面的數值可知,可以根據數值的正負判斷叉乘后向量的方向。即,如果叉積AB^AM 和 AB^AN的結果同號,那么M,N兩點就在直線的同側,否則不在同一側。特殊地,如果點M在直線AB上,則AB^AM的值為0。(如果是在三維坐標系中,求出的叉積是一個向量,可以根據兩個向量的點積結果正負來判斷兩個向量的是否指向同一側)                                                         本文地址

 

以上的問題解決了,就很容易的用來判斷某個點是否在三角形內,如果P在三角形ABC內部,則滿足以下三個條件:P,A在BC的同側、P,B在AC的同側、PC在AB的同側。某一個不滿足則表示P不在三角形內部。

該算法詳見后面代碼中的函數:IsPointInTriangle2

 

算法3

該方法也用到了向量。對於三角形ABC和一點P,可以有如下的向量表示:

image

p點在三角形內部的充分必要條件是:1 >= u >= 0,   1 >= v >= 0,   u+v <= 1。

 

已知A,B,C,P四個點的坐標,可以求出u,v,把上面的式子分別點乘向量AC和向量AB

image

解方程得到:

image

解出u,v后只需要看他們是否滿足“1 >= u >= 0, 1 >= v >= 0, u+v <= 1”,如滿足,則,p 在三角形內。

(u = 0時,p在AB上, v = 0時,p在AC上,兩者均為0時,p和A重合)

該算法詳見后面代碼中的函數:IsPointInTriangle3

 

 

算法4

該算法和算法2類似,可以看作是對算法2的簡化,也是用到向量的叉乘。假設三角形的三個點按照順時針(或者逆時針)順序是A,B,C。對於某一點P,求出三個向量PA,PB,PC, 然后計算以下三個叉乘(^表示叉乘符號):

t1 = PA^PB,

t2 = PB^PC,

t3 = PC^PA,

如果t1,t2,t3同號(同正或同負),那么P在三角形內部,否則在外部。

該算法詳見后面代碼中的函數:IsPointInTriangle4

 


經過測試,算法4最快,算法3次之,接着算法2,算法1最慢。直觀的從計算量上來看,也是算法4的計算量最少。

以下是代碼,定義了兩個類:二維向量類 和 三角形類

 #include<string>
#include<sstream>
#include<fstream>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<vector>
#include<stack>
#include<queue>
#include<list>
#include<algorithm>
#include<ctime>
#include<unordered_map>
#include<map>
#include<typeinfo>
#include<cmath>
#include<ctime>
#include<climits>
using namespace std;

//類定義:二維向量
class Vector2d
{
public:
double x_;
double y_;

public:
Vector2d(double x, double y):x_(x), y_(y){}
Vector2d():x_(0), y_(0){}

//二維向量叉乘, 叉乘的結果其實是向量,方向垂直於兩個向量組成的平面,這里我們只需要其大小和方向
double CrossProduct(const Vector2d vec)
{
return x_*vec.y_ - y_*vec.x_;
}

//二維向量點積
double DotProduct(const Vector2d vec)
{
return x_ * vec.x_ + y_ * vec.y_;
}

//二維向量減法
Vector2d Minus(const Vector2d vec) const
{
return Vector2d(x_ - vec.x_, y_ - vec.y_);
}

//判斷點M,N是否在直線AB的同一側
static bool IsPointAtSameSideOfLine(const Vector2d &pointM, const Vector2d &pointN,
const Vector2d &pointA, const Vector2d &pointB)
{
Vector2d AB = pointB.Minus(pointA);
Vector2d AM = pointM.Minus(pointA);
Vector2d AN = pointN.Minus(pointA);

//等於0時表示某個點在直線上
return AB.CrossProduct(AM) * AB.CrossProduct(AN) >= 0;
}
};

//三角形類
class Triangle
{
private:
Vector2d pointA_, pointB_, pointC_;

public:
Triangle(Vector2d point1, Vector2d point2, Vector2d point3)
:pointA_(point1), pointB_(point2), pointC_(point3)
{
//todo 判斷三點是否共線
}

//計算三角形面積
double ComputeTriangleArea()
{
//依據兩個向量的叉乘來計算,可參考http://blog.csdn.net/zxj1988/article/details/6260576
Vector2d AB = pointB_.Minus(pointA_);
Vector2d BC = pointC_.Minus(pointB_);
return fabs(AB.CrossProduct(BC) / 2.0);
}

bool IsPointInTriangle1(const Vector2d pointP)
{
double area_ABC = ComputeTriangleArea();
double area_PAB = Triangle(pointP, pointA_, pointB_).ComputeTriangleArea();
double area_PAC = Triangle(pointP, pointA_, pointC_).ComputeTriangleArea();
double area_PBC = Triangle(pointP, pointB_, pointC_).ComputeTriangleArea();

if(fabs(area_PAB + area_PBC + area_PAC - area_ABC) < 0.000001)
return true;
else return false;
}

bool IsPointInTriangle2(const Vector2d pointP)
{
return Vector2d::IsPointAtSameSideOfLine(pointP, pointA_, pointB_, pointC_) &&
Vector2d::IsPointAtSameSideOfLine(pointP, pointB_, pointA_, pointC_) &&
Vector2d::IsPointAtSameSideOfLine(pointP, pointC_, pointA_, pointB_);
}

bool IsPointInTriangle3(const Vector2d pointP)
{
Vector2d AB = pointB_.Minus(pointA_);
Vector2d AC = pointC_.Minus(pointA_);
Vector2d AP = pointP.Minus(pointA_);
double dot_ac_ac = AC.DotProduct(AC);
double dot_ac_ab = AC.DotProduct(AB);
double dot_ac_ap = AC.DotProduct(AP);
double dot_ab_ab = AB.DotProduct(AB);
double dot_ab_ap = AB.DotProduct(AP);

double tmp = 1.0 / (dot_ac_ac * dot_ab_ab - dot_ac_ab * dot_ac_ab);

double u = (dot_ab_ab * dot_ac_ap - dot_ac_ab * dot_ab_ap) * tmp;
if(u < 0 || u > 1)
return false;
double v = (dot_ac_ac * dot_ab_ap - dot_ac_ab * dot_ac_ap) * tmp;
if(v < 0 || v > 1)
return false;

return u + v <= 1;
}

bool IsPointInTriangle4(const Vector2d pointP)
{
Vector2d PA = pointA_.Minus(pointP);
Vector2d PB = pointB_.Minus(pointP);
Vector2d PC = pointC_.Minus(pointP);
double t1 = PA.CrossProduct(PB);
double t2 = PB.CrossProduct(PC);
double t3 = PC.CrossProduct(PA);
return t1*t2 >= 0 && t1*t3 >= 0;
}
};


int main()
{
Triangle tri(Vector2d(0,0), Vector2d(6,6), Vector2d(12,0));
srand(time(0));
for(int i = 0; i < 20; ++i)
{
Vector2d point((rand()*1.0 / RAND_MAX) * 12, (rand()*1.0 / RAND_MAX) * 6);
cout<<point.x_<<" "<<point.y_<<": ";
cout<<tri.IsPointInTriangle1(point)<<" ";
cout<<tri.IsPointInTriangle2(point)<<" ";
cout<<tri.IsPointInTriangle3(point)<<" ";
cout<<tri.IsPointInTriangle4(point)<<endl;

}
}

 

參考資料:

http://www.cnblogs.com/graphics/archive/2010/08/05/1793393.html

http://blog.csdn.net/fox64194167/article/details/8147460

http://blog.csdn.net/dracularking/article/details/2217180

 

【版權聲明】轉載請注明出處:http://www.cnblogs.com/TenosDoIt/p/4024413.html


注意!

本站转载的文章为个人学习借鉴使用,本站对版权不负任何法律责任。如果侵犯了您的隐私权益,请联系我们删除。



 
粤ICP备14056181号  © 2014-2020 ITdaan.com