一个判断点是否在浮点坐标边界的多边形之内的新算法 (50分)

前两天接到一个网友的来信，提出要解决判断一个点是否在单连通多边形区域之内的问题，
我搜索了一下以前的帖子，看见论坛上有很多这方面的问题，答案也有很多。
（大家可以看一看下面的帖子）
http://www.delphibbs.com/delphibbs/dispq.asp?lid=0997458
http://www.delphibbs.com/delphibbs/dispq.asp?lid=0093903
另：关于多边形的面积，大家可以看
http://www.delphibbs.com/delphibbs/dispq.asp?lid=0170597 其中yysun老师的回答
使用API函数PtInRegion无疑是非常简单的方法，但是它仅能用于座标为整数的情况，有
局限性。网友们还提出了很多算法，其中最普遍的是“跳栏”法——即根据从待考察点发出
的一条射线与多边形的边相交的次数进行判断，我研究了一下，似乎还有商榷的余地。
我在此自创了一种新的算法。源代码及说明如下：
原来写的算法在判断交点位于多边形顶点的情况时会出现错误，究其原因，还是没有考虑
周全造成的。我们现在再次进行详细的分析：
首先，判断点是否在多边形的边界线上（包括顶点），如果在，则认为点在多边形内部。
如果点不在多边形的边界上，我们就过待考查点做一条正向水平射线，并对它与多边形的
边界相交的情况进行分析：
（需要注意的是：由于做的是水平线，所以所有水平方向的边界都不用参加计算——在上面
已经将位于边界上的情况过滤掉了，对于不在水平边界上、但是y坐标又和水平边界一致的
情况，在下面的计算中不用考虑）
1.没有交点
说明点在多边形之外
2.有交点
找到距离待考查点最近的那个交点；
2.1.该交点位于直线上
可以直接根据直线的穿越方向进行判断：由于多边形区域始终处于边界前进方向的右侧，
如果边界自下而上穿越水平线，那么说明点在多边形之外；反之则说明在内部。
2.2.该交点位于多边形的顶点上
假设这个顶点和位于它两边的顶点按照多边形顶点的顺时针的排列顺序分别为A,B,C；
如果线段AB和BC都是向上穿越水平线，那么说明点在多边形之外；反之，如果它们都向
下穿越水平线，那么说明点在多边形之内；
如果在上述这两种情况之外，那么就要考察下面这四种情况：（假设待考查点为O）
A C
C A
O - - B - - > O - - B - - >
O - - B - - > O - - B - - >
C
A C A
我们可以发现，可以根据直线AB以及BC的斜率关系判断点B是多边形的局部“凸”点还
是局部“凹”点。如果是凸点，那么说明O点必在多边形之外；反之则说明在多边形之内。
type
TBorderRec=Record //点的结构
X,Youble;
end;
TPolygon=class
private
FData:array of TBorderRec;
FPointCount: Integer;
FMinX:do
uble;
FMinY:do
uble;
FMaxX:do
uble;
FMaxY:do
uble;
function Get(Index: Integer): TBorderRec;
procedure Put(Index: Integer;
const Value: TBorderRec);
public
property MaxXouble read FMaxX;
property MaxYouble read FMaxY;
property MinXouble read FMinX;
property MinYouble read FMinY;
property PointCount:Integer read FPointCount;
property Items[Index:Integer]:TBorderRec read Get write Put;
default;
procedure Clear;
function PointInside(x,youble):Boolean;
procedure LoadFromFile(const FileName:String);
constructor Create;
destructor Destroy;
override;
end;

{ TPolygon }
procedure TPolygon.Clear;
begin
SetLength(FData,0);
FPointCount:=0;
FMinX:=0;
FMinY:=0;
FMaxX:=0;
FMaxY:=0;
end;

constructor TPolygon.Create;
begin
Clear;
end;

destructor TPolygon.Destroy;
begin
Clear;
inherited;
end;

function TPolygon.Get(Index: Integer): TBorderRec;
begin
if (Index>=0) and (Index<FPointCount) then
Result:=FData[Index]
else
raise Exception.Create('数组访问越界!');
end;

//从文件获得数据，每行存放一个座标值，x和y之间用逗号或空格分隔
//点应该呈顺时针方向，且最后一个点应该和第一个点一致，以封闭多边形
procedure TPolygon.LoadFromFile(const FileName: String);
var
i,n:Integer;
x,y,MaxX,MaxY,MinX,MinYouble;
SL,msl:TStringList;
begin
Clear;
SL:=TStringList.Create;
msl:=TStringList.Create;
MaxX:=0;
MinX:=0;
MaxY:=0;
MinY:=0;
try
SL.LoadFromFile(FileName);
n:=SL.Count;
SetLength(FData,n);
for i:=0 to n-1do
begin
msl.CommaText:=SL;
if msl.Count<>2 then
raise Exception.Create('无效的数据格式!');
x:=StrToFloat(msl[0]);
y:=StrToFloat(msl[1]);
if i=0 then
begin
MaxX:=x;
MinX:=x;
MaxY:=y;
MinY:=y;
end
else
begin
if x>MaxX then
MaxX:=x
else
if x<MinX then
MinX:=x;
if y>MaxY then
MaxY:=y
else
if y<MinY then
MinY:=y;
end;
FData.X:=x;
FData.Y:=y;
end;
FPointCount:=n;
FMaxX:=MaxX;
FMaxY:=MaxY;
FMinX:=MinX;
FMinY:=MinY;
finally
SL.Free;
msl.Free;
end;
end;

function TPolygon.PointInside(x, y:do
uble): Boolean;
function PointOnLine(x1,y1,x2,y2ouble):Boolean;
var
xx,yyouble;
begin
if ((x<x1) and (x<x2)) or ((y<y1) and (y<y2))
or ((x>x1) and (x>x2)) or ((y>y1) and (y>y2)) then
begin
Result:=false;
exit;
end;
xx:=x-x1;
yy:=y-y1;
x2:=x2-x1;
y2:=y2-y1;
if x2=0 then
Result:=xx=0 //在垂线上
else
if y2=0 then
Result:=yy=0 //在水平线上
else
Result:=xx*y2=yy*x2;
end;
function Between(f,f1,f2ouble):Boolean;
begin
if f1<f2 then
Result:=(f>=f1) and (f<=f2)
else
Result:=(f<=f1) and (f>=f2);
end;
function IntersectionDistance(x1,y1,x2,y2ouble)ouble;
begin
//获得直线与水平线的交点到待考查点的距离——带方向（右侧为正）
x1:=x1-x;
x2:=x2-x;
if x1=x2 then
begin
//垂直穿越
Result:=x1;
exit;
end;
y1:=y1-y;
y2:=y2-y;
Result:=x2+y2*(x2-x1)/(y1-y2);
end;
var
i,MinPoint,SameMinPoint:Integer;
MinDistance,tmpfouble;
x1,y1,x2,y2ouble;
begin
Result:=false;
if ((x<FMinX) or (x>FMaxX)) and ((y<FMinY) or (y>FMaxY)) then
exit;
//判断是否落在多边形的边上（包含顶点）
for i:=0 to FPointCount-2do
if PointOnLine(FData.X,FData.Y,FData[i+1].X,FData[i+1].Y) then
begin
Result:=true;
exit;
end;
MinDistance:=0;
MinPoint:=-1;
SameMinPoint:=-1;
for i:=0 to FPointCount-2do
if ((FData.X>=x) or (FData[i+1].X>=x)) and
(FData.Y<>FData[i+1].Y) then
//过滤整体位于左侧的线段以及水平线
if Between(y,FData.Y,FData[i+1].Y) then
begin
tmpf:=IntersectionDistance(FData.X,FData.Y,FData[i+1].X,FData[i+1].Y);
if tmpf>0 then
//因为点不在直线上，所以交点的距离不可能等于0
if MinPoint=-1 then
begin
MinPoint:=i;
MinDistance:=tmpf;
end
else
if tmpf<MinDistance then
//找到一个距离更小的
begin
MinPoint:=i;
SameMinPoint:=-1;
MinDistance:=tmpf;
end
else
if tmpf=MinDistance then
//遇到距离同样近的情况，记下前一个的位置
begin
SameMinPoint:=MinPoint;
MinPoint:=i;
end;
end;
if (MinPoint>=0) and (SameMinPoint<0) then
//有交点，且最近交点不是多边形的顶点
Result:=FData[MinPoint].Y>FData[MinPoint+1].Y //根据穿越方向进行判别
else
if SameMinPoint>=0 then
//有交点，但是是多边形的顶点
begin
if SameMinPoint>MinPoint then
begin
i:=SameMinPoint;
SameMinPoint:=MinPoint;
MinPoint:=i;
end;
y1:=FData[SameMinPoint].Y-FData[MinPoint].Y;
y2:=FData[MinPoint].Y-FData[MinPoint+1].Y;
//由于已经过滤了水平重合的线段...
if (y1>0) and (y2>0) then
//两条线都是自上而下穿越水平线
Result:=true
else
if (y1<0) and (y2<0) then
//两条线都是自下而上穿越水平线
Result:=false
else
begin
//两条直线有升有降——位于水平线的同一侧
x1:=FData[SameMinPoint].X-FData[MinPoint].X;
x2:=FData[MinPoint].X-FData[MinPoint+1].X;
Result:=x1*y2>x2*y1;
//根据两条直线的斜率大小进行判断
end;
end;
end;

procedure TPolygon.Put(Index: Integer;
const Value: TBorderRec);
begin
if (Index>=0) and (Index<FPointCount) then
FData[Index]:=Value
else
raise Exception.Create('数组访问越界!');
end;

给大家一组测试用数据：
0,-0.5
0,0
0.5,-1
1,0
1,-1
0,-2
1,-1.5
2,-3
1,-4
1.25,-3.5
1,-2.5
0,-2.5
-0.5,-3.5
-1,-2.5
-1.5,-3.1
-2,-2.5
-1.5,-2.8
-1,-1
0,-0.5
欢迎查错、提出改进意见！（分数可以另外给）

呵学习有机会拜你为师！·

困扰我3年的问题被你解决了，再次非常感谢。例外我向你求教两个问题，是该问题的延伸：
1、由直线、（圆、椭圆）弧线任意组合（边与便不相交）成的封闭多边形，应该如何判断
2、由第一个问题组成的多边形，如何计算面积。
非常希望你能解决，我的E-mail：llh0223@xinhuanet.com

多谢，学习了。

碰到终结版都很“谦虚”！！！！！！！！
我非常佩服作者缜密而复杂的思维。但是还是忍不住
废话几句。
首先以上的算法仅适用于“凸”多边组。
其次，既然是对凸多边形所做的算法，那么有更简单的
思路。
大家可以考虑一下，对于任意的凸多边形，取其内部的
任意一点。通过该点画一条任何方向的直线。那么至少和多
边形的两条边相交。（如果将多边形的顶点只作为边的终结
点或起始点，那么所画的直线一定和多边形的两条边相交，
如果我们将所取的内部点作为原点，在所画的直线上标明方向，
那么这两点分别位于下负轴。）

综合所述，算法也不难实现，不考虑多边形边的情况下，
就按作者的“水平扫描法”可如下进行。当然也可以考虑将
边的情况揉进去。
var
i:integer;
HaveRightBorder:Boolean=false;
HaveLeftBorder:Boolean=false;
begin
for i:=0 to FPointCount-2do
begin
if not Between(y,FData.Y,FData[i+1].Y) then
continue;
//这种边不可能相交
if (IntersectionDistance(FData.X,FData.Y,FData[i+1].X,FData[i+1].Y)>0)
then
HaveRightborder:=true //正方向有边相交,只要是两个方向，不管熟正熟负。
else
HaveLeftborder:=true;
//负方向有边相交
if (HaveRightborder and HaveLeftborder) then
Break;
//肯定在内部；
end;
Result:=(HaveRightborder and HaveLeftborder);
end;

to jsxjd:
>>首先以上的算法仅适用于“凸”多边组。
您没有完全理解算法的用途，我在上面已经评论过“跳栏法”的不足——如果是凸多边形，
那么穿越的次数只可能是0,1,2，不用扩大到奇数、偶数。同样，我的程序弥补了它的不足，
当然可以用于任意单连通多边形区域。
估计使您产生误解的是“多边形的顶点已经事先按顺时针方向排列”一句。请看：
J A D
B
I H C
E
G
F
上面这个 A..B..C.....I..J..A 组成的封闭多边形，我们同样认为它的边界点的排列顺序
是“顺时针”的——可它不是一个凸多边形。
——我的算法已经针对湖南省的边界数据（892条边）进行了实验，完全达到了要求。

感谢提出不同意见！请大家继续。

to creation-zy:
我没有仔细看你的代码,可能有些不全.但是感到似乎有点问题：
您的凹多变形好像最多只凹一次，但是实际上可能是这样的：
A B
F E
I J
D C
G H

L K
象这个多变形ABCDEFGHIJKLA你的凸凹判定好像有问题。

谢谢zjczxd兄！
您认为我的多边形“最多只凹一次”，这很可能是因为误解了我在上面的 2.2 中的分析
过程。的确，在我上面的图像中，看起来水平扫描线似乎只能合一个顶点重合，请注意，我
上面的图像是基于“距离待考查点最近的那个交点”的，因为我的算法是基于最近交点的穿
越方向的，除了最近的一个交点，根本不用考虑其余的相交情况，比如：
A C A C E F
O- - -B- - - - E - -> O- - -B- - - - - - - ->
D F D
在这两种情况中，我只要考虑最近的一个交点——B的相交情况即可，更远的相交情况不会
影响最终的判定结果。
我在此给出我所谓的“顺时针方向”多边形顶点序列的定义：
1.该序列的首元素应该和尾元素一致，以封闭多边形；
2.该序列中，任意两个相邻元素不能相同，它们代表一条线段的两个端点；
3.由该序列生成的一系列线段，不能相交（当然可以“相接”——本来就是折线嘛）；
4.当我们从一个顶点走向下一个顶点时，多边形的“内部”始终处于右手方向。
请注意，最关键的是第四条定义——我的算法就是基于它的。无论从待考查点发射出去的
射线与多边形有多少个交点，我只要考察它在最近的交点所在线段的左侧还是右侧即可。

to creation-zy: 好，如果有你的第四个假设就好成立了，不过适用性就可能会比较差了。
下面我来讨论一个一般情况的解法：
实际上跳栅栏法的缺点就是：
这条线可能会进进出出好几次，而你无法知道哪一次是出哪一次是进。
当碰到一个顶点时它可能是只出（或进）了一次，还是出（或进）了两次。
所以只要讨论顶点的相邻两条直线和水平线的关系就可以了，而不应该讨论是在多边形外还是内。
A C
---B------------
A,D在水平线的同一侧就说明是在两次穿越
A
-----B------
C
A,D在水平线的异侧就是一次穿越
至于在水平线上可以认为0次穿越,即将B,C看成同一个点。
----B---C----
这样如下的情况就好办了:
A D
---B---C-----
而判断是否同侧是极为简单的,只要判断他们的纵坐标就可以了。
麻烦creation-zy帮我看看是否有不对的地方？

DarwinZhang兄，您的分析是完全正确的。Good ideal!
在加入了您的改进之后，跳栏法就可以进行准确的判断了。现在，两种算法都能完成任务
（当然，我的算法应用面没有跳栏法宽），并且两个算法的时间复杂度都是相当的（都是O(
顶点个数)）。
看来讨论应该结束了...

有没有考虑过复连通区,就是一个多边行里面又有一个洞(多边形套在里面),
或洞里面还有洞...
在数学(复变函数)中内部好像是这样定义的:沿边界逆时针旋转的左侧
而判断时针方向是用右手螺旋的法方向
这个问题我也想过很九,没写出安这个定义的算法

对于多连通区域，可以有多个边界列表，只要这些列表均满足我在上面提到的4点定义，
我的算法是完全可以进行正确判断的（跳栏法也可以）。

我以为跳栅栏法还有一个奇点的问题需要讨论:就是区域内的某点不包含在该区域内.
这时需要考虑恰好穿过该点时认为是穿过两次.

DarwinZhang兄提出的“奇点”问题的确存在。我们是不是可以这么分析它：根据我在上面
提到的几点条件中第一点的的要求，要定义一个点，就必须要求序列中有两个元素——起始
点和终止点（在这种特殊情况下，第二个条件可以忽略）——这样一来，“奇点”可以被看
成是由两个长度为0的水平线段所围成的区域——而水平线在两种算法中都是被忽略的。但是
我们还会发现，“奇点”本身因为位于线段的边界，也可以被认为处于多边形之内。
我们发现了一个问题，为了精确的判断，我们必须确定刚好位于边界上的点到底属于哪一
方。为了实现这个功能，我们可以在定义各个线段序列的同时规定一个布尔值——如果为真，
那么就意味着线段上的点位于多边形之内；反之，则在多边形之外。

creation-zy,
你好，由于不知道你的邮箱，我不知道你能否收到我的信件，我要问一下关于你在
http://www.delphibbs.com/delphibbs/dispq.asp?lid=990544中提到的穷举算法。能否
给一个实例给我。行吗？guoyunzhi@yeah.net 我写的得程序运行时间太长了。是从M个数
中选N个数，列举所有的组合。能否赐教。我的QQ是8185606.请联系我，好吗？
谢谢

谢谢大家了！

这个问题以前也想过，就是想不出，今天竟然在这里找到了，高手！！！