比如已知三个顶点 A(xa,za) B(xb,zb) C(xc,zc)
以及O点(xo,zo)
判断O点在X-Z平面上是否在三角形ABC内
我想到不少方法,但在找最快最简洁的。
1.面积判别法
如果AOC,BOC,AOB面积之和等于ABC就在内
在此处使用一种常见且简便的方法:如果三角形PAB,PAC和PBC的面积之和与三角形ABC的面积相等,即可判定点P在三角形ABC内(包括在三条边上)。
可知,该方法的关键在于如何计算三角形的面积。幸运地是,当知道三角形顶点(A,B和C)的坐标((Ax, Ay),(Bx, By)和(Cx, Cy))之后,即可计算出其面积:
关键的代码如下,
// Point(java.awt.Point)代表点的坐标。
private static double triangleArea(Point pos1, Point pos2, Point pos3) {
double result = Math.abs((pos1.x * pos2.y + pos2.x * pos3.y + pos3.x * pos1.y
- pos2.x * pos1.y - pos3.x * pos2.y - pos1.x * pos3.y) / 2.0D );
return result;
}
private static boolean inTriangle(Point pos, Point posA, Point posB,
Point posC) {
double triangleArea = triangleArea(posA, posB, posC);
double area = triangleArea(pos, posA, posB);
area += triangleArea(pos, posA, posC);
area += triangleArea(pos, posB, posC);
double epsilon = 0.0001 ; // 由于浮点数的计算存在着误差,故指定一个足够小的数,用于判定两个面积是否(近似)相等。
if (Math.abs(triangleArea - area) < epsilon) {
return true ;
}
return false ;
}
在一个开源脚本了偶然看到这么几行,如下:
function Polygon_contains(point)
{
var i, j, status=false;
for (i=0, j=this.numpoints-1; i<this.numpoints; j=i++) {
if ((((this.points[i].y<=point.y) && (point.y<this.points[j].y)) || ((this.points[j].y<=point.y) && (point.y<this.points[i].y))) && (point.x < (this.points[j].x - this.points[i].x) * (point.y - this.points[i].y) / (this.points[j].y - this.points[i].y) + this.points[i].x))
status = !status;
}
return status;
}
仔细分析了一下,虽然简单,可道理和铅垂线法差不多,但是省去了求交,统计交点个数等过程。
主要逻辑都在那个if条件语句
(((this.points[i].y<=point.y) && (point.y<this.points[j].y)) || ((this.points[j].y<=point.y) && (point.y<this.points[i].y))) 就是判断点point在,过 i 点的水平线和过 j 点的水平线之间(临界条件是一端闭合一端开)
(point.x < (this.points[j].x - this.points[i].x) * (point.y - this.points[i].y) / (this.points[j].y - this.points[i].y) + this.points[i].x) 是判断点点point在线段 ij 的坐标,注意这里的线段 ij 不分方向。我们应该添加少许代码以适应 ij 竖直的情况,而且好像根据斜率的正负,好像要分情况考虑。写一个函数判断点是否在线段左侧:
private bool PointLeft(CPoint p, CPoint i, CPoint j)
{
if (i.X == j.X)
{
return p.X < i.X;
}
double deltY = j.Y - i.Y;
double deltX = j.X - i.X;
if ((deltX >= 0 & deltY > 0) || (deltX <= 0 && deltY < 0))
{
return p.Y > deltY / deltX * (p.X - i.X) + i.Y;
}
else
{
return p.Y < deltY / deltX * (p.X - i.X) + i.Y;
}
}
刚开始假设点在外面,没符合这个条件,里外就反一下,最后status为真就在多边形内。
经过检验,对带环的多边形、逆、顺多边形和凸凹多边形都有效。
添加少许代码就可检验点是否在边上或和某个节点重合。
作者:hyp
0.前言
最近不断遇到类似的几何位置问题,一直没有花时间去总结,本文总结了我常用点跟多边形的位置判断方法以及代码。希望能够对大家有所帮助。
文中所指的多边形均为凸多边形,一些描述可能有误,欢迎指正。
1.测试的多边形
在开始之前,我们需要先构建好测试环境。
我构建了一个比较特殊的多边形,如下。
/ \
| |
|_|
从最上面的顶点顺时针坐标(屏幕坐标系)分别为:(40,10) (60,30)(60,50) (20,50) (20,30)
2.射线判别法
根据对多边形的了解,我们可以得出如下结论:
如果一个点在多边形内部,任意角度做射线肯定会与多边形要么有一个交点,要么有与多边形边界线重叠。
如果一个点在多边形外部,任意角度做射线要么与多边形有一个交点,要么有两个交点,要么没有交点,要么有与多边形边界线重叠。
利用上面的结论,我们只要判断这个点与多边形的交点个数,就可以判断出点与多边形的位置关系了。
首先罗列下注意事项:
l 射线跟多边形的边界线重叠的情况
l 区别内部点和外部点的射线在有一个交点时的情况
对于第一个注意事项,可以将射线角度设为零度,这样子只需要判断两个相邻顶点的Y值是否相等即可。然后再判断这个点的X值方位。
对于第二个注意事项,网上许多文章都说到做射线以后交点为奇数则表示在多边形内部,这是一个错误的观点,不仅对于凹多边形不成立,对于凸多边形也不成立。
例如:从外部点做射线刚好经过一顶点,这样子交点个数就为奇数,但是该点却不在多边形内部。
至于要如何区分这两种情况呢,我能想到了一个不完美的方法,外部点的射线跟多边形有一个交点的时候,该交点肯定为顶点,如果该射线上移一位或者下移一位,要么变成有两个交点要么没有交点。当然为了安全起见,这里把射线尽量往相邻点中心移动。这样子就能够判断出是外部点的射线跟多边形有一个交点。
不过这个方法并不完美,因为有了移位操作,可能会把内部点移动出外部。而且如果判断点在(60,30)位置,判断的时候先遇到(20,30),然后移位操作,就判断就出错了。
为了解决这些问题,我在起初先扫描一次判断点是否在顶点上虽然影响了一点效率,而且当判定点距离多边形一个单位时,判断可能会有误。不过只要不是需要高精度的话,这个方法还是很有效的。
代码如下:
bool InsidePolygon1( POINTD *polygon,int N,POINTD pt )
{
int i,j;
bool inside,redo;
redo = true;
for (i = 0;i < N;++i)
{
if (polygon[i].x == pt.x && // 是否在顶点上
polygon[i].y == pt.y )
{
redo = false;
inside = true;
break;
}
}
while (redo)
{
redo = false;
inside = false;
for (i = 0,j = N - 1;i < N;j = i++)
{
if ( (polygon[i].y < pt.y && pt.y < polygon[j].y) ||
(polygon[j].y < pt.y && pt.y < polygon[i].y) )
{
if (pt.x <= polygon[i].x || pt.x <= polygon[j].x)
{
double _x = (pt.y-polygon[i].y)*(polygon[j].x-polygon[i].x)/(polygon[j].y-polygon[i].y)+polygon[i].x;
if (pt.x < _x) // 在线的左侧
inside = !inside;
else if (pt.x == _x) // 在线上
{
inside = true;
break;
}
}
}
else if ( pt.y == polygon[i].y)
{
if (pt.x < polygon[i].x) // 交点在顶点上
{
polygon[i].y > polygon[j].y ? --pt.y : ++pt.y;
redo = true;
break;
}
}
else if ( polygon[i].y == polygon[j].y && // 在水平的边界线上
pt.y == polygon[i].y &&
( (polygon[i].x < pt.x && pt.x < polygon[j].x) ||
(polygon[j].x < pt.x && pt.x < polygon[i].x) ) )
{
inside = true;
break;
}
}
}
return inside;
}
3.角度和判别法
角度和判别法就是判定点与多边形所有相邻顶点组成的角的角度和来判断点与多边形的位置关系。
如果点在多边形内部,只要该点不在边界线或者顶点上,则角度和为三百六十度。
如果在边界线上,则角度和为一百八十度。
如果点在多边形外部,则角度和无法达到三百六十度。但是角度和可能会达到一百八十度,比如判断点在(60,10)。
代码如下:
// 根据需要不判断顶点
bool IsPointInLine( const POINTD &pt,const POINTD &pt1,const POINTD &pt2 )
{
bool inside = false;
if (pt.y == pt1.y &&
pt1.y == pt2.y &&
((pt1.x < pt.x && pt.x < pt2.x) ||
(pt2.x < pt.x && pt.x < pt1.x)) )
{
inside = true;
}
else if (pt.x == pt1.x &&
pt1.x == pt2.x &&
((pt1.y < pt.y && pt.y < pt2.y) ||
(pt2.y < pt.y && pt.y < pt1.y)) )
{
inside = true;
}
else if ( ((pt1.y < pt.y && pt.y < pt2.y) ||
(pt2.y < pt.y && pt.y < pt1.y)) &&
((pt1.x < pt.x && pt.x < pt2.x) ||
(pt2.x < pt.x && pt.x < pt1.x)) )
{
if (0 == (pt.y-pt1.y)/(pt2.y-pt1.y)-(pt.x - pt1.x) / (pt2.x-pt1.x))
{
inside = true;
}
}
return inside;
}
bool InsidePolygon2( POINTD *polygon,int N,POINTD p )
{
int i,j;
double angle = 0;
bool inside = false;
for (i = 0,j = N - 1;i < N;j = i++)
{
if (polygon[i].x == p.x && // 是否在顶点上
polygon[i].y == p.y)
{
inside = true;
break;
}
else if (IsPointInLine(p,polygon[i],polygon[j])) // 是否在边界线上
{
inside = true;
break;
}
double x1,y1,x2,y2;
x1 = polygon[i].x - p.x;
y1 = polygon[i].y - p.y;
x2 = polygon[j].x - p.x;
y2 = polygon[j].y - p.y;
double radian = atan2(y1,x1) - atan2(y2,x2);
radian = abs(radian);
if (radian > M_PI) radian = 2* M_PI - radian;
angle += radian; // 计算角度和
}
if ( fabs(6.28318530717958647692 - angle) < 1e-7 )
inside = true;
return inside;
}
有的人管角度和判别法叫做转角法,还有一个类似的方法叫弧长法。
http://haibarali.blog.163.com/blog/static/23474013200722813356858/
4.面积和判别法
根据角度和判别法,我们可以继续演化,可以得出如下结论:
如果一个点在多边形内部,该点与多边形所有相邻顶点组成的三角形面积和为多边形面积。反之不成立。
求三角形面积:
已知三角形A(x1,y1)B(x2,y2)C(x3,y3),则面积公式为:
|x1 x2 x3|
S(A,B,C) = |y1 y2 y3| * 0.5 (当三点为逆时针时为正,顺时针则为负的)
|1 1 1 |
= ( x1*y2 - x1*y3 - x2*y1 + x2*y3 + x3*y1 - x3*y2 ) * 0.5
求多边形面积:
根据上面求三角形的方法,
已知多边形A1A2A3...An(顺或逆时针都可以),设平面上有任意的一点P,则面积公式为:
S(A1,A2,A3...An)
= S(P,A1,A2)+ S(P,A2,A3)+...+S(P,An,A1)
既然P是可以任取,为了方便用(0,0)好了。
代码如下:
bool InsidePolygon3( POINTD *polygon,int N,POINTD pt )
{
int i,j;
bool inside = false;
double polygon_area = 0;
double trigon_area = 0;
for (i = 0,j = N - 1;i < N;j = i++)
{
polygon_area += polygon[i].x * polygon[j].y - polygon[j].x * polygon[i].y;
trigon_area += abs(
pt.x * polygon[i].y -
pt.x * polygon[j].y -
polygon[i].x * pt.y +
polygon[i].x * polygon[j].y +
polygon[j].x * pt.y -
polygon[j].x * polygon[i].y
);
}
trigon_area *= 0.5;
polygon_area = abs(polygon_area * 0.5);
if ( fabs(trigon_area - polygon_area) < 1e-7 )
inside = true;
return inside;
}
多边形面积计算
http://www.cppblog.com/zyzx/archive/2009/04/27/81241.html
行列式如何展开
http://www.tongji.edu.cn/~math/xxds/kcja/kcja_b/1-6.htm
5.点线判别法
经网友glshader提示,添加了这个方法。
如果判断点在所有边界线的同侧,就能判定该点在多边形内部。
判断方法就是判断两条同起点射线斜率差。
代码如下:
bool InsidePolygon4( POINTD *polygon,int N,POINTD p )
{
int i,j;
bool inside = false;
int count1 = 0;
int count2 = 0;
for (i = 0,j = N - 1;i < N;j = i++)
{
double value = (p.x - polygon[j].x) * (polygon[i].y - polygon[j].y) - (p.y - polygon[j].y) * (polygon[i].x - polygon[j].x);
if (value > 0)
++count1;
else if (value < 0)
++count2;
}
if (0 == count1 ||
0 == count2)
{
inside = true;
}
return inside;
}
1. 叉乘判别法(只适用于凸多边形)
想象一个凸多边形,其每一个边都将整个2D屏幕划分成为左右两边,连接每一边的第一个端点和要测试的点得到一个矢量v,将两个2维矢量扩展成3维的,然后将该边与v叉乘,判断结果3维矢量中Z分量的符号是否发生变化,进而推导出点是否处于凸多边形内外。这里要注意的是,多边形顶点究竟是左手序还是右手序,这对具体判断方式有影响。
2. 面积判别法(只适用于凸多边形)
第四点分别与三角形的两个点组成的面积分别设为S1,S2,S3,只要S1+S2+S3>原来的三角形面积就不在三角形范围中.可以使用海伦公式 。推广一下是否可以得到面向凸多边形的算法?(不确定)
3. 角度和判别法(适用于任意多边形)
double angle = 0;
realPointList::iterator iter1 = points.begin();
for (realPointList::iterator iter2 = (iter1 + 1); iter2 < points.end(); ++iter1, ++iter2)
{
double x1 = (*iter1).x - p.x;
double y1 = (*iter1).y - p.y;
double x2 = (*iter2).x - p.x;
double y2 = (*iter2).y - p.y;
angle += angle2D(x1, y1, x2, y2);
}
if (fabs(angle - span::PI2) < 0.01) return true;
else return false;
另外,可以使用bounding box来加速。
if (p.x < (*iter)->boundingBox.left ||
p.x > (*iter)->boundingBox.right ||
p.y < (*iter)->boundingBox.bottom ||
p.y > (*iter)->boundingBox.top) 。。。。。。
对于多边形来说,计算bounding box非常的简单。只需要把水平和垂直方向上的最大最小值找出来就可以了。
对于三角形:第四点分别与三角形的两个点的交线组成的角度分别设为j1,j2,j3,只要j1+j2+j3>360就不在三角形范围中。
4. 水平/垂直交叉点数判别法(适用于任意多边形)
注意到如果从P作水平向左的射线的话,如果P在多边形内部,那么这条射线与多边形的交点必为奇数,如果P在多边形外部,则交点个数必为偶数(0也在内)。所以,我们可以顺序考虑多边形的每条边,求出交点的总个数。还有一些特殊情况要考虑。假如考虑边(P1,P2),
1)如果射线正好穿过P1或者P2,那么这个交点会被算作2次,处理办法是如果P的从坐标与P1,P2中较小的纵坐标相同,则直接忽略这种情况
2)如果射线水平,则射线要么与其无交点,要么有无数个,这种情况也直接忽略。
3)如果射线竖直,而P0的横坐标小于P1,P2的横坐标,则必然相交。
4)再判断相交之前,先判断P是否在边(P1,P2)的上面,如果在,则直接得出结论:P再多边形内部。