http://www.cnblogs.com/Linkabox/p/3357768.html
三 重心法
解法一:面积法
S=Area(ABC)、S1=Area(ABP)、S2=Area(BCP)、S3=Area(CAP)
若P点在三角形内,S=S1+S2+S3
若P在三角形外,S1+S2+S3>S
计算三角形的面积可以使用海伦公式
而公式里的p为半周长(周长的一半):
但注意了我们可以利用向量积的几何意义以两个向量为邻边的平行四边形的面积(即三角形面积的2倍)
若 |AB×BP| + |BC×CP| + |CA×AP| = |AB×BC|,则得点 P 在 △ABC 内。
double TriArea(point A,point B,point P) { return abs((B.x-A.x)*(P.y-B.y)-(B.y-A.y)*(P.x-B.x)); } bool isInTriangle1(point A,point B,point C,point P) { double sABP,sBCP,sCAP,sABC; sABP=TriArea(A,B,P); sBCP=TriArea(B,C,P); sCAP=TriArea(C,A,P); sABC=TriArea(A,B,C); if (sABC==sABP+sBCP+sCAP) { return true; } return false; }
解法二:叉乘法
沿 △ABC 各有向边按一定方向走(顺时针或逆时针),判断点 P 是否在该边的某侧(右侧或左侧),若点 P 在三条边的同侧,则点 P 在 △ABC 内。
分别计算向量 AB、BC、CA 与向量 AP、BP、CP 的向量积(叉乘),若三个结果均同号(正或负,为零表示 P 在边上),则可得点 P 在 △ABC 内。
其中AB = B - A,AB×AP = AB.x*AP.y - AB.y*AP.x即 x1y2-y1x2。
其实不止三角形对于凸多边形也是可以通过叉乘法来判断的。
补充知识:
所谓凸多边形,就是把一个多边形任意一边向两方无限延长成为一条直线,如果多边形的其他各边均在此直线的同旁,那么这个多边形就叫做凸多边形。
下面是示例代码:注意面积法对输入的点序没要求,叉乘法就要按逆时针的顺序。
inline double VecProduct(point A,point B,point P) { return (B.x-A.x)*(P.y-A.y)-(B.y-A.y)*(P.x-A.x); }
bool isInTriangle2(point A,point B,point C,point P) { if (VecProduct(A,B,P)>=0 && VecProduct(B,C,P)>=0 && VecProduct(C,A,P)>=0) { return true; } return false; }http://www.cnblogs.com/graphics/archive/2010/08/05/1793393.html
假设点P位于三角形内,会有这样一个规律,当我们沿着ABCA的方向在三条边上行走时,你会发现点P始终位于边AB,BC和CA的右侧。我们就利用这一点,但是如何判断一个点在线段的左侧还是右侧呢?我们可以从另一个角度来思考,当选定线段AB时,点C位于AB的右侧,同理选定BC时,点A位于BC的右侧,最后选定CA时,点B位于CA的右侧,所以当选择某一条边时,我们只需验证点P与该边所对的点在同一侧即可。问题又来了,如何判断两个点在某条线段的同一侧呢?可以通过叉积来实现,连接PA,将PA和AB做叉积,再将CA和AB做叉积,如果两个叉积的结果方向一致,那么两个点在同一测。判断两个向量的是否同向可以用点积实现,如果点积大于0,则两向量夹角是锐角,否则是钝角。
三角形的三个点在同一个平面上,如果选中其中一个点,其他两个点不过是相对该点的位移而已,比如选择点A作为起点,那么点B相当于在AB方向移动一段距离得到,而点C相当于在AC方向移动一段距离得到。
所以对于平面内任意一点,都可以由如下方程来表示
P = A + u * (C – A) + v * (B - A) // 方程1
如果系数u或v为负值,那么相当于朝相反的方向移动,即BA或CA方向。那么如果想让P位于三角形ABC内部,u和v必须满足什么条件呢?有如下三个条件
u >= 0
v >= 0
u + v <= 1
几个边界情况,当u = 0且v = 0时,就是点A,当u = 0,v = 1时,就是点B,而当u = 1, v = 0时,就是点C
整理方程1得到P – A = u(C - A) + v(B - A)
令v0 = C – A, v1 = B – A, v2 = P – A,则v2 = u * v0 + v * v1,现在是一个方程,两个未知数,无法解出u和v,将等式两边分别点乘v0和v1的到两个等式
(v2) • v0 = (u * v0 + v * v1) • v0
(v2) • v1 = (u * v0 + v * v1) • v1
注意到这里u和v是数,而v0,v1和v2是向量,所以可以将点积展开得到下面的式子。
v2 • v0 = u * (v0 • v0) + v * (v1 • v0) // 式1
v2 • v1 = u * (v0 • v1) + v * (v1• v1) // 式2
解这个方程得到
u = ((v1•v1)(v2•v0)-(v1•v0)(v2•v1)) / ((v0•v0)(v1•v1) - (v0•v1)(v1•v0))
v = ((v0•v0)(v2•v1)-(v0•v1)(v2•v0)) / ((v0•v0)(v1•v1) - (v0•v1)(v1•v0))
bool PointinTriangle(Vector3 A, Vector3 B, Vector3 C, Vector3 P) { Vector3 v0 = C - A ; Vector3 v1 = B - A ; Vector3 v2 = P - A ; float dot00 = v0.Dot(v0) ; float dot01 = v0.Dot(v1) ; float dot02 = v0.Dot(v2) ; float dot11 = v1.Dot(v1) ; float dot12 = v1.Dot(v2) ; float inverDeno = 1 / (dot00 * dot11 - dot01 * dot01) ; float u = (dot11 * dot02 - dot01 * dot12) * inverDeno ; if (u < 0 || u > 1) // if u out of range, return directly { return false ; } float v = (dot00 * dot12 - dot01 * dot02) * inverDeno ; if (v < 0 || v > 1) // if v out of range, return directly { return false ; } return u + v <= 1 ; }
http://www.cnblogs.com/TenosDoIt/p/4024413.html
算法2
首先看一下这个问题,如何判断某两个点在某条直线的同一侧(代码中函数:IsPointAtSameSideOfLine)?
根据向量的叉乘以及右手螺旋定则,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不在三角形内部。
算法4
该算法和算法2类似,可以看作是对算法2的简化,也是用到向量的叉乘。假设三角形的三个点按照顺时针(或者逆时针)顺序是A,B,C。对于某一点P,求出三个向量PA,PB,PC, 然后计算以下三个叉乘(^表示叉乘符号):
t1 = PA^PB,
t2 = PB^PC,
t3 = PC^PA,
如果t1,t2,t3同号(同正或同负),那么P在三角形内部,否则在外部
//类定义:二维向量
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;
}
};
English