点与多边形的包含判定，代码

[homerlu]

点与多边形的包含判定，代码
--------------------------------------------------------------------
有朋友问起这个问题的处理，就总是回复说，射线求交，可自己从来没真正下手写过
到自己需要了，竟然到处都找不到，网上也搜不见，真是奇怪，讨论的是挺多的。难道大家都停留在讨论阶段，还是windows api中的PtInRegion就够用了呢。我也挺懒，有的东西向来喜欢拿来主义，这次写写看吧，不过感觉例外情况很多，处理繁琐。比如这个代码没有处理相邻两边在同一条直线且跟射线也共线的情况。---好像说起来满别扭的.
发上来大家或许可以借鉴一下，并恳请各位帮忙检查一下还有哪些错误的地方。有劳，谢了！

一点说明，多边形定点值数组为pllPoint，nPoiCount是其定点数或者说边数

/*
* function: PtInRegion
* ascertain that a point in region or not.
* parameter: coordinate of point.
* return value:
* 0 --- out of the region;
* 1 --- in;
* 2 --- is climax;
* 3 --- in edge.
* code by Homer Lu, 2003-11-13
*/
int CRegion::PtInRegion(double fX, double fY)
{
int nStartP, nEndP;
double fX0, fY0 = fY;

int nCrossCount = 0;

for (int i = 0; i < nPoiCount; i++)
{
nStartP = i;
if(i + 1 == nPoiCount)
nEndP = 0;
else
nEndP = i + 1;

if(fY > pllPoint[nStartP].lat && fY > pllPoint[nEndP].lat) //on top of
continue ;
else if(fY < pllPoint[nStartP].lat && fY < pllPoint[nEndP].lat) //foot
continue ;
else if(fY == pllPoint[nStartP].lat && fX == pllPoint[nStartP].lon)
return 2;
else if(fY == pllPoint[nStartP].lat && fY == pllPoint[nEndP].lat)
{
if(fX > max(pllPoint[nStartP].lon, pllPoint[nEndP].lon))
continue ;
else if(fX < max(pllPoint[nStartP].lon, pllPoint[nEndP].lon) && fX > min(pllPoint[nStartP].lon, pllPoint[nEndP].lon))
return 3;
else //fX <= min(pllPoint[i].lon, pllPoint[i + 1].lon)
{
if(i - 1 < 0)
nStartP = nPoiCount - 1;
else
nStartP = i - 1;
if(i + 2 >= nPoiCount)
nEndP = i + 2 - nPoiCount;
else
nEndP = i + 2;
i++;
if(fY > max(pllPoint[nStartP].lat, pllPoint[nEndP].lat) ||
fY < min(pllPoint[nStartP].lat, pllPoint[nEndP].lat))
{
nCrossCount ++ ;
continue ;
}
else
continue ;
}
}
else if(fY == pllPoint[nStartP].lat)
{
if(fX < pllPoint[nStartP].lon)
nCrossCount ++ ;
else
continue ;
}
else if(fY == pllPoint[nEndP].lat)
{
nStartP = i;
if(i + 2 >= nPoiCount)
nEndP = i + 2 - nPoiCount;
else
nEndP = i + 2;

if(fY > max(pllPoint[nStartP].lat, pllPoint[nEndP].lat) ||
fY < min(pllPoint[nStartP].lat, pllPoint[nEndP].lat))
{
nCrossCount ++ ;
continue ;
}
else
continue ;
}
else //if (fY < max(pllPoint[nStartP].lat, pllPoint[nEndP].lat) &&
// (fY > min(pllPoint[nStartP].lat, pllPoint[nEndP].lat)))
{
// get cross-point. formula: x0 = x1 + (x2 - x1) * (y0 - y1) / (y2 - y1)
fX0 = pllPoint[nStartP].lon + (pllPoint[nEndP].lon - pllPoint[nStartP].lon) * (fY - pllPoint[nStartP].lat) / (pllPoint[nEndP].lat - pllPoint[nStartP].lat);
if(fX0 == fX)
return 3;
else if(fX0 > fX)
{
nCrossCount ++ ;
continue ;
}
else // fX0 < fX
continue ;
}
} // end for (int i = 0; i < nPoiCount - 1; i++)

return nCrossCount % 2;
}

[zhouqingyuan]

up

[homerlu]

改两个bug
/*
* function: PtInRegion
* ascertain that a point in region or not.
* parameter: coordinate of point.
* return value:
* 0 --- out of the region;
* 1 --- in;
* 2 --- is climax;
* 3 --- in edge.
* code by Homer Lu, 2003-11-13
*/
int CRegion::PtInRegion(double fX, double fY)
{
int nStartP, nEndP;
double fX0, fY0 = fY;

int nCrossCount = 0;

for (int i = 0; i < nPoiCount; i++)
{
nStartP = i;
if(i + 1 == nPoiCount)
nEndP = 0;
else
nEndP = i + 1;

if(fY > pllPoint[nStartP].lat && fY > pllPoint[nEndP].lat) //on top of
continue ;
else if(fY < pllPoint[nStartP].lat && fY < pllPoint[nEndP].lat) //foot
continue ;
else if((fY == pllPoint[nStartP].lat && fX == pllPoint[nStartP].lon) ||
(fY == pllPoint[nEndP].lat && fX == pllPoint[nEndP].lon))
return 2;
else if(fY == pllPoint[nStartP].lat && fY == pllPoint[nEndP].lat)
{
if(fX > max(pllPoint[nStartP].lon, pllPoint[nEndP].lon))
continue ;
else if(fX < max(pllPoint[nStartP].lon, pllPoint[nEndP].lon) &&
fX > min(pllPoint[nStartP].lon, pllPoint[nEndP].lon))
return 3;
else //fX <= min(pllPoint[i].lon, pllPoint[i + 1].lon)
{
if(i - 1 < 0)
nStartP = nPoiCount - 1;
else
nStartP = i - 1;
if(i + 2 >= nPoiCount)
nEndP = i + 2 - nPoiCount;
else
nEndP = i + 2;
i++;
if(fY < max(pllPoint[nStartP].lat, pllPoint[nEndP].lat) ||
fY > min(pllPoint[nStartP].lat, pllPoint[nEndP].lat))
{
nCrossCount ++ ;
continue ;
}
else
continue ;
}
}
else if(fY == pllPoint[nStartP].lat)
{
if(fX < pllPoint[nStartP].lon)
nCrossCount ++ ;
else
continue ;
}
else if(fY == pllPoint[nEndP].lat)
{
nStartP = i;
if(i + 2 >= nPoiCount)
nEndP = i + 2 - nPoiCount;
else
nEndP = i + 2;

if(fY > max(pllPoint[nStartP].lat, pllPoint[nEndP].lat) ||
fY < min(pllPoint[nStartP].lat, pllPoint[nEndP].lat))
{
nCrossCount ++ ;
continue ;
}
else
continue ;
}
else //if (fY < max(pllPoint[nStartP].lat, pllPoint[nEndP].lat) &&
// (fY > min(pllPoint[nStartP].lat, pllPoint[nEndP].lat)))
{
// get cross-point. formula: x0 = x1 + (x2 - x1) * (y0 - y1) / (y2 - y1)
fX0 = pllPoint[nStartP].lon + (pllPoint[nEndP].lon - pllPoint[nStartP].lon) * (fY - pllPoint[nStartP].lat) / (pllPoint[nEndP].lat - pllPoint[nStartP].lat);
if(fX0 == fX)
return 3;
else if(fX0 > fX)
{
nCrossCount ++ ;
continue ;
}
else // fX0 < fX
continue ;
}
} // end for (int i = 0; i < nPoiCount - 1; i++)

return nCrossCount % 2;
}

[homerlu]

up

[ameba]

扫描线法。
可以参考计算机图形

[BlueSky2008]

up.
点与多边形的包含判定确实问的很多.但真正给出代码的很少。非常感谢楼主！
另外还有转角法和面积法，欢迎大家都来参与，把这个问题的代码搞全，做成一个faq.尤其欢迎HUTTON和LeeMaRs!

[HUNTON]

to BlueSky2008:写错了，是HUNTON。对了上次我的那题你说用什么IEEE定义来做的，我怎么也不知道你是怎么拟合的，希望您再看一下。

面积法只适合凸多边形，我以前的是用浮点型做的，下面的代码是我用整型现写的，还没试过，不知有没有什么问题。

//功能：根据三顶点坐标求三角形面积的2倍
//参数：三点的坐标
//返回：该三角形的面积的2倍
int CalculateS(POINT P1, POINT P2, POINT P3)
{
int S = 0;

S = (P1.x -P3.x) * (P2.y - P3.y) - (P2.x - P3.x) * (P1.y - P3.y);
return abs(S);
}


//功能：判断一点与凸n多边形的位置关系
//参数：判别点、多边形顶点数组、多边形顶点个数
//返回：
// 1 在内部
// 0 在边上
// -1 在外部
int IsPointInTu(POINT P0, POINT *P, int n)
{
int S = 0; //存放多边形的面积的2倍
int S0 = 0; //存放该点与该凸多边形各边构成的三角形的面积之和的2倍
int nS0 = 0, temp;
int i = 0;

for(i = 0; i <= n - 3; i++){
S += CalculateS(P[0], P[i + 1], P[i + 2]);
}
for(i = 0 ; i <= n - 2 ; i ++){
temp = CalculateS(P0, P[i], P[i + 1]);
S0 += temp;
if(temp == 0) nS0++;
}
temp = CalculateS(P0, P[0], P[n - 1]);
if(temp == 0) nS0++;
S0 += temp;
if(S0 == S){
if(nS0 > 0) return 0;
else return 1;
}
else return -1;
}

[LeeMaRS]

这个是角度和法的，我写在ZJU 1081的程序里面了, 就一起帖上来吧..

{$MODE FPC}
Program T_1081;
Const
EPS = 1e-6;

Type
TPoint = Record
x, y : Double;
End;

TLine = Record
P1, P2 : TPoint;
End;

TPolygon = ^TPoint;

Var
Polygon : TPolygon;
P : TPoint;
N, M, i, t : LongInt;

Function Zero(n : Double) : Boolean;
Begin
Zero := Abs(n) <= EPS;
End;

Function Norm(P : TPoint) : Double;
Begin
Norm := Sqrt(Sqr(P.x) + Sqr(P.y));
End;

Function Sub(P1, P2 : TPoint) : TPoint;
Begin
Sub.x := P1.x - P2.x; Sub.y := P1.y - P2.y;
End;

Function Multiply(P1, P2 : TPoint) : Double;
Begin
Multiply := (P1.x * P2.y) - (P2.x * P1.y);
End;

Function DotMultiply(P1, P2 : TPoint) : Double;
Begin
DotMultiply := (P1.x * P2.x) + (P1.y * P2.y);
End;

Function OnLineSeg(P : TPoint; L : TLine) : Boolean;
Begin
OnLineSeg := Zero(Multiply(Sub(L.P1, P), Sub(L.P2, P))) AND
((P.x - L.P1.x) * (P.x - L.P2.x) <= 0) AND
((P.y - L.P1.y) * (P.y - L.P2.y) <= 0);
End;

Function ArcCos(n : Double) : Double;
Begin
If n > 1 Then n := 1 Else If n < -1 Then n := -1;

If Zero(n) Then
ArcCos := PI / 2
Else
Begin
ArcCos := ArcTan(Sqrt(1 - Sqr(n)) / n);
If n < 0 Then
ArcCos := ArcCos + PI;
End;
End;

Function AngleBetween(P1, P2 : TPoint) : Double;
Var
t : Double;
Begin
t := DotMultiply(P1, P2) / (Norm(P1) * Norm(P2));
If Multiply(P1, P2) > 0 Then
AngleBetween := ArcCos(t)
Else
AngleBetween := - ArcCos(t);
End;

Function InsidePolygon(P : TPoint; Polygon : TPolygon) : Boolean;
Var
L : TLine;
Angle : Double;
i : LongInt;
Begin
If N = 1 Then
Begin
InsidePolygon := Zero(P.x - Polygon[0].x) AND Zero(P.y - Polygon[0].y);
Exit;
End;

Angle := 0;

For i := 0 To N - 1 Do
Begin
L.P1 := Polygon[i]; L.P2 := Polygon[(i+1) MOD N];

If OnLineSeg(P, L) Then
Begin
InsidePolygon := True; Exit;
End;

Angle := Angle + AngleBetween(Sub(L.P1, P), Sub(L.P2, P));
End;

InsidePolygon := NOT Zero(Angle);
End;

Begin
t := 0; Read(N);
While (N <> 0) Do
Begin
Readln(M);
GetMem(Polygon, SizeOf(TPoint) * N);
For i := 1 To N Do
With Polygon[i-1] Do
Readln(x, y);

Inc(t);
Writeln('Problem ', t, ':');
For i := 1 To M Do
Begin
With P Do
Readln(x, y);
If InsidePolygon(P, Polygon) Then
Writeln('Within')
Else
Writeln('Outside');
End;

FreeMem(Polygon, SizeOf(TPoint) * N);

Read(N);
If N <> 0 Then Writeln;
End;
End.

[ZhangYv]

http://expert.csdn.net/Expert/topic/2463/2463557.xml?temp=.2284967
我整理了旧的精华贴 ：）

[homerlu]

修改：
/*
* function: PtInRegion
* ascertain that a point in region or not.
* parameter: coordinate of point.
* return value:
* 0 --- out of the region;
* 1 --- in;
* 2 --- is climax;
* 3 --- in edge.
* code by Homerlu 2003-11-13
* HomerLu 2003-11-16
*/
int CRegion::PtInRegion(double fX, double fY)
{
int nStartP, nEndP;
double fX0, fY0 = fY;

int nCrossCount = 0;

for (int i = 0; i < nPoiCount; i++)
{
nStartP = i;
if(i + 1 == nPoiCount)
nEndP = 0;
else
nEndP = i + 1;

if(fY > pllPoint[nStartP].lat && fY > pllPoint[nEndP].lat) //on top of
continue ;
else if(fY < pllPoint[nStartP].lat && fY < pllPoint[nEndP].lat) //foot
continue ;
else if((fY == pllPoint[nStartP].lat && fX == pllPoint[nStartP].lon) ||
(fY == pllPoint[nEndP].lat && fX == pllPoint[nEndP].lon))
return 2;
else if(fY == pllPoint[nStartP].lat && fY == pllPoint[nEndP].lat)
{
if(fX > max(pllPoint[nStartP].lon, pllPoint[nEndP].lon))
continue ;
else if(fX < max(pllPoint[nStartP].lon, pllPoint[nEndP].lon) &&
fX > min(pllPoint[nStartP].lon, pllPoint[nEndP].lon))
return 3;
else //fX <= min(pllPoint[i].lon, pllPoint[i + 1].lon)
{
if(i - 1 < 0)
nStartP = nPoiCount - 1;
else
nStartP = i - 1;
if(i + 2 >= nPoiCount)
nEndP = i + 2 - nPoiCount;
else
nEndP = i + 2;
i++;
if(fY < max(pllPoint[nStartP].lat, pllPoint[nEndP].lat) ||
fY > min(pllPoint[nStartP].lat, pllPoint[nEndP].lat))
{
nCrossCount ++ ;
continue ;
}
else
continue ;
}
}
else if(fY == pllPoint[nStartP].lat)
{
if(fX < pllPoint[nStartP].lon)
nCrossCount ++ ;
else
continue ;
}
else if(fY == pllPoint[nEndP].lat)
{
nStartP = i;
if(i + 2 >= nPoiCount)
nEndP = i + 2 - nPoiCount;
else
nEndP = i + 2;

if(fY > max(pllPoint[nStartP].lat, pllPoint[nEndP].lat) ||
fY < min(pllPoint[nStartP].lat, pllPoint[nEndP].lat))
{
nCrossCount ++ ;
continue ;
}
else
continue ;
}
else //if (fY < max(pllPoint[nStartP].lat, pllPoint[nEndP].lat) &&
// (fY > min(pllPoint[nStartP].lat, pllPoint[nEndP].lat)))
{
// get cross-point. formula: x0 = x1 + (x2 - x1) * (y0 - y1) / (y2 - y1)
fX0 = pllPoint[nStartP].lon + (pllPoint[nEndP].lon - pllPoint[nStartP].lon) * (fY - pllPoint[nStartP].lat) / (pllPoint[nEndP].lat - pllPoint[nStartP].lat);
if(fX0 == fX)
return 3;
else if(fX0 > fX)
{
nCrossCount ++ ;
continue ;
}
else // fX0 < fX
continue ;
}
} // end for (int i = 0; i < nPoiCount - 1; i++)

return nCrossCount % 2;
}

[homerlu]

多谢各位整理的和关注的朋友

[homerlu]

ZhangYv(网络故障中)
http://expert.csdn.net/Expert/topic/2463/2463557.xml?temp=.2284967
我整理了旧的精华贴 ：）

先生整理的东东实在太好了，多谢

[BlueSky2008]

Sorry，HUNTON,把你的名字拼错了。你那个题，说用IEEE定义来做，是Happy_888说的啊，不是我说的啊。你的题我还是不太理解，和一般的一维插值，二维插值有什么区别吗？

多谢HONTON and LeeMaRS 支持，楼主要是没有什么问题，就赶快结贴吧，你的程序最好再加点注释，前面的说明可以就用ZhangYv的了：）

[homerlu]

不知道为什么，我2m的带宽访问这里很慢很慢，向来很难结贴，至今只结过两三个，其他的几次不成功时候我就放弃了。所以我的信誉分都被扣了
这个也不定能成功呢

哎，以后啊，
结就一个字，我只结一次

[homerlu]

最近在作数据处理时，数据量大了，很慢。而且有判断错误（真是失败，改了三次还是有错：（）
重新检查修改。
看到斑竹对上面的代码作了总结处理。
http://expert.csdn.net/expert/FAQ/FAQ_Index.asp?id=184334
赶快把更新的代码贴上来，还烦请斑竹处理一下
谢谢

/*
* function: PtInRegion()
* ascertain that a point in region or not.
* parameter: coordinate of point.
* return value:
* 0 --- out of the region;
* 1 --- in;
* 2 --- is climax;
* 3 --- on edge.
* code by Homerlu 2003-11-13
* modify 2003-12-23
*/
int CRegion::PtInRegion(double fX, double fY)
{
int nStartP, nEndP;
double fX0, fY0 = fY;

int nCrossCount = 0;

if( fX > m_mbr.llRightBottom.lon || fX < m_mbr.llLeftTop.lon ||
fY > m_mbr.llLeftTop.lat || fY < m_mbr.llRightBottom.lat)
return 0;

for (int i = 0; i < m_nPoiCount; i++)
{
nStartP = i;
if(i + 1 == m_nPoiCount)
nEndP = 0;
else
nEndP = i + 1;

if(fX > max(m_pllPoint[nStartP].lon, m_pllPoint[nEndP].lon)) //right
continue ;
else if(fY > max(m_pllPoint[nStartP].lat, m_pllPoint[nEndP].lat)) //on top of
continue ;
else if(fY < min(m_pllPoint[nStartP].lat, m_pllPoint[nEndP].lat)) //under
continue ;
else if((fY == m_pllPoint[nStartP].lat && fX == m_pllPoint[nStartP].lon) ||
(fY == m_pllPoint[nEndP].lat && fX == m_pllPoint[nEndP].lon))
return 2; //climax
else if(fY == m_pllPoint[nStartP].lat && fY == m_pllPoint[nEndP].lat)//parallel
{
// if(fX > max(m_pllPoint[nStartP].lon, m_pllPoint[nEndP].lon))
// continue ; //right
//else
if(fX < max(m_pllPoint[nStartP].lon, m_pllPoint[nEndP].lon) &&
fX > min(m_pllPoint[nStartP].lon, m_pllPoint[nEndP].lon))
return 3; //on edge
else //fX < min(m_pllPoint[i].lon, m_pllPoint[i + 1].lon)
{
if(i - 1 < 0)
nStartP = m_nPoiCount - 1;
else
nStartP = i - 1;
if(i + 2 >= m_nPoiCount)
nEndP = i + 2 - m_nPoiCount;
else
nEndP = i + 2;
i++;
if(fY > max(m_pllPoint[nStartP].lat, m_pllPoint[nEndP].lat) ||
fY < min(m_pllPoint[nStartP].lat, m_pllPoint[nEndP].lat))
{
nCrossCount ++ ;
continue ;
}
else
continue ;
}
}
else if(fY == m_pllPoint[nStartP].lat)//but fY != m_pllPoint[nEndP].lat
{
if(fX < m_pllPoint[nStartP].lon)
nCrossCount ++ ;
else
continue ;
}
else if(fY == m_pllPoint[nEndP].lat)//but fY != m_pllPoint[nStartP].lat
{
if(fX > m_pllPoint[nEndP].lon)
continue ;
nStartP = i;
if(i + 2 >= m_nPoiCount)
nEndP = i + 2 - m_nPoiCount;
else
nEndP = i + 2;

if(fY > max(m_pllPoint[nStartP].lat, m_pllPoint[nEndP].lat) ||
fY < min(m_pllPoint[nStartP].lat, m_pllPoint[nEndP].lat))
{
nCrossCount ++ ;/////////////////////////
continue ;
}
else
continue ;
}
else //if (fY < max(m_pllPoint[nStartP].lat, m_pllPoint[nEndP].lat) &&
// (fY > min(m_pllPoint[nStartP].lat, m_pllPoint[nEndP].lat)) &&
// (fX <= max(m_pllPoint[nStartP].lon, m_pllPoint[nEndP].lon))
{
if(fX <= min(m_pllPoint[nStartP].lon, m_pllPoint[nEndP].lon))
{
nCrossCount ++ ;
continue ;
}

// get cross-point. formula: x0 = x1 + (x2 - x1) * (y0 - y1) / (y2 - y1)
fX0 = m_pllPoint[nStartP].lon + (m_pllPoint[nEndP].lon - m_pllPoint[nStartP].lon) * (fY - m_pllPoint[nStartP].lat) / (m_pllPoint[nEndP].lat - m_pllPoint[nStartP].lat);
if(fX0 == fX) //edge
return 3;
else if(fX0 > fX) //right
{
nCrossCount ++ ;
continue ;
}
else // fX0 < fX //left
continue ;
}
} // end for (int i = 0; i < m_nPoiCount - 1; i++)

return nCrossCount % 2; //0 or 1
}

[homerlu]

在vc里面的排版还可以，贴到这里来显得比较乱

另外：
1、因为我做的是地理数据，数据量比较大（常常一个多边形几千个点）。
所以为每个多边形生成一个mbr,判断的时候速度快了许多
这次的代码也加上了
if( fX > m_mbr.llRightBottom.lon || fX < m_mbr.llLeftTop.lon ||
fY > m_mbr.llLeftTop.lat || fY < m_mbr.llRightBottom.lat)
return 0;
2、针对不同的数据，可能我的代码前半部分不需要（或者修改），可以归到求交点上面去。但要看你的数据什么情况的居多，考虑计算效率的问题