点在多边形内判断方法matlab实现（点在凸包内、任意形状内判定）

惯例声明：本人没有相关的工程应用经验，只是纯粹对相关算法感兴趣才写此博客。所以如果有错误，欢迎在评论区指正，不胜感激。本文主要关注于算法的实现，对于实际应用等问题本人没有任何经验，所以也不再涉及。

0 前言

本文大概说一下目前常用计算机判断点在多边形内的方法。这个实际应用还是比较普遍的，比如单纯的几何意义上的判断，或者选定一个范围，看看数据是否超出这个范围之类的。

咱们人在判断点是否在多边形内，通常是先识别多边形图形，然后看看点是否在这个识别的图形范围内。但是对于很复杂的图形，就不容易一下子得到结论，而且这种基于图形的算法实际上效率并不高，不适用于计算机计算。

常用的判断点是否在多边形内的算法，与多边形的变数n相关，通常都是O(n)量级的，所以其实大家都相差不大，就看谁更能优化。但是对于凸多边形来说，可以用二分法算法，将计算复杂度降到O(logn)量级，这个在后面会涉及到。

之后如果有空，准备试一试写个判断点是否在三维以及高维多面体中的文章（又挖坑不填 ）。

本文参考的博客和论文如下：

[1] 凸多边形（计算几何，判断点是否在多边形内，二分法）
https://www.cnblogs.com/yym2013/p/3673616.html
[2] 判断点在多边形内的算法（Winding Number详解）
https://www.codenong.com/cs106265948/
[3]Hormann K, Agathos A. The point in polygon problem for arbitrary polygons[J]. Computational geometry, 2001, 20(3): 131-144.

1 射线交叉点法 （Crossing Number）

该方法称为射线法或者奇偶规则法。原理就是从该点P出发向外射出一条射线，如果射线与多边形的交点数量为奇数，则点P在多边形内。

可以想象一个肥皂泡，如果想从肥皂泡里冲出去，则一定要穿过一次肥皂泡。如果穿过了两次肥皂泡，则一进一出等于没有穿过。

通常为了简单起见，射线直接选择为+x方向，也就是向右的方向。

一般优化的思路有：
1大概定义一个多边形的范围，如果点超过这个范围，则无需再判断，直接判断为不在图形内。
2与+x方向射线相交的线段，其y的范围一定要在射线两端。也就是说如果线段在射线的正上方，或者线段在射线的正下方，都不需要再计算交点，因为不可能相交。
3与+x方向射线相交的线段，不可能在射线的左侧。

这个方法看着很简单，但是还有一些细节需要扣。比如当边线正好也是水平线，怎么计算交点？当点落在边线上，怎么计算交点？当边线正好是水平线，而且点也落在了这个边线上，交叉点算几个？当点的射线正好穿过多边形顶点，会不会得到与两个边相交的bug？

关于这几个问题，我大概解决方法为：设置一个边线点的概念，专门用来判断点是否在边线上。如果点刚好落在水平线的边线上，则不再计算是否相交，直接当做点在边线上处理。交叉点的判断，用一头小于等于，一头小于的方式，避免重复计算。而且因为计算中存在误差，所以对于计算得到的交点xc，需要加一个误差带来防止判断失误。

具体不多说，matlab程序如下：

clear
clc
close all
%多边形定义（连线必须按照首尾相接的顺序）
BD=[0,0;0,-1;1,0;2,-1;2,0;2,1;2,2;1,2;1,1;0,1];%凹多边形
%BD=[0,-1;0,-2;1,-1;2,-2;2,-1;2,0;2,3;-1,3;-1,1;...
%    3,1;3,1.5;0,1.5;0,2.5;1,2.5;1,0;0,0];%凹多边形，而且自相交
%BD=[0,1;-1,0;-1,-1;0,-1;1,-1;1,0];%凸多边形
%要判断的点
[X,Y]=meshgrid(-5:0.1:5,-5:0.1:5);
xy2=[X(:),Y(:)];

[IsInPoly,IsOnBD]=IfInPoly1(BD,xy2);

figure()
hold on
plot(BD(:,1),BD(:,2))
scatter(xy2(:,1),xy2(:,2),24,IsInPoly+2*IsOnBD,'Marker','.')
%scatter(xy2(:,1),xy2(:,2),24,or(IsInPoly,IsOnBD),'Marker','.');%不显示边线

%后置函数
function [IsInPoly,IsOnBD]=IfInPoly1(BD,xy2)
%输出逻辑索引（IsInPoly表示在内部，IsOnBD表示在多边形边界上）
%BD是多边形边界，存在顺序，两列。xy2是点的坐标，两列。
%方法1，射线法
if (BD(1,1)~=BD(end,1)) || (BD(1,2)~=BD(end,2))
    BD2=[BD;BD(1,:)];
else %给出的边界已经收尾相接
    BD2=BD;
end
%删除边界中相邻重复的点
IsSame=and(BD2(1:end-1,1)==BD2(2:end,1),BD2(1:end-1,2)==BD2(2:end,2));
BD2([IsSame;false],:)=[];
NB=size(BD2,1)-1;
NP=size(xy2,1);

IsInPoly=false(NP,1);%true(NP,1);%false(NP,1);
IsOnBD=false(NP,1);%是否在边线上
%判断整个边线的大概范围
min_X_BD=min(BD(:,1));
max_X_BD=max(BD(:,1));
min_Y_BD=min(BD(:,2));
max_Y_BD=max(BD(:,2));
%做+x方向的射线
for kp=1:NP
    %如果这个点的xy超过整个边线的xy，则肯定不在边线内
    xy_k=xy2(kp,:);
    xp=xy_k(1);
    yp=xy_k(2);
    if xp<min_X_BD || xp>max_X_BD || yp<min_Y_BD || yp>max_Y_BD
        %IsInPoly(kp)=false;
        continue
    end
    %循环每一条边
    NCross=0;%初始化交点数量
    for kB=1:NB
        xB1=BD2(kB,1);xB2=BD2(kB+1,1);
        yB1=BD2(kB,2);yB2=BD2(kB+1,2); 
        %如果在+x方向上相交，在点一定在线的左边
        if max([xB1,xB2])<xp
            continue %所以点在线右侧的情况无需计算
        end
        %判断是否在水平的边线上
        if yB1==yp && yB2==yp && ( min([xB1,xB2])<=xp && max([xB1,xB2])>=xp)
            IsOnBD(kp)=true;%如果点在水平线段上，则证明点在边缘
            break %停止循环
        end
        %判断是否是某个顶点
        if (xp==xB1 && yp==yB1) || (xp==xB2 && yp==yB2) 
            IsOnBD(kp)=true;%如果点是某个边缘顶点，则证明点在边缘
            break
        end
        %如果射线穿过这个边，则y值一定介于这个边的两个y值中间
        if (yB1<=yp && yp<yB2) || (yB2<=yp && yp<yB1) 
            %如果在+x方向上相交，则交点一定在xp的右边
            yc=yp;
            xc=xB1+(xB2-xB1)/(yB2-yB1)*(yc-yB1);%计算射线与边的交点(xc,yc)
            if (xc-xp)>4*eps %由于计算会导致xc有一定的误差，所以保守估算为4*eps
                NCross=NCross+1;%证明的确相交，交点+1
            end
            %再次判断是否在边线上
            if abs(xc-xp)<=4*eps && abs(yc-yp)<=4*eps %这里因为也涉及到xc，所以也设置了一个误差带
                IsOnBD(kp)=true;%如果交点xc就是点xp，则证明点在线上
                break
            end
        end
        
    end
    %根据相交点数的奇偶性判断是否在多边形内
    if ~IsOnBD(kp) && mod(NCross,2)
        IsInPoly(kp)=true;%如果是奇数，则证明在多边形内
    end
end
end

注：这里给出了点在边界内IsInPoly，点在多边形边界上IsOnBD两个结果。如果需要合并，可以用or(IsInPoly,IsOnBD)进行合并结果。

计算结果如下：

2 环绕数法（Winding Number）

环绕数法和相交数法类似，都是引出一条射线。但是环绕数不再计算具体交点，而是判断相交的线段，是向上穿过这个射线还是向下穿过这个射线。

举个例子，还是用上一章节的图进行举例子。但是这里的线段都被赋予了方向。我们定义向上穿过射线（点在线段左侧），环绕数+1，向下穿过射线（点在线段右侧），环绕数-1。最终对环绕数相加，等于0说明点在多边形外，不等于0说明点在多边形内。

这个方法的matlab程序如下：

clear
clc
close all
%多边形定义（连线必须按照首尾相接的顺序）
BD=[0,0;0,-1;1,0;2,-1;2,0;2,1;2,2;1,2;1,1;0,1];%凹多边形
%BD=[0,-1;0,-2;1,-1;2,-2;2,-1;2,0;2,3;-1,3;-1,1;...
%    3,1;3,1.5;0,1.5;0,2.5;1,2.5;1,0;0,0];%凹多边形，而且自相交
%BD=[0,1;-1,0;-1,-1;0,-1;1,-1;1,0];%凸多边形
%要判断的点
[X,Y]=meshgrid(-5:0.1:5,-5:0.1:5);
xy2=[X(:),Y(:)];

[IsInPoly,IsOnBD]=IfInPoly2(BD,xy2);

figure()
hold on
plot(BD(:,1),BD(:,2))
scatter(xy2(:,1),xy2(:,2),24,IsInPoly+2*IsOnBD,'Marker','.')
%如果考虑边界情况，则把两个结果合并
%IsInPoly=or(IsInPoly,IsOnBD);

function [IsInPoly,IsOnBD]=IfInPoly2(BD,xy2)
%输出逻辑索引（IsInPoly表示在内部，IsOnBD表示在多边形边界上）
%BD是多边形边界，存在顺序，两列。xy2是点的坐标，两列。
%方法2 winding number 环绕数法
%将边界收尾相接
if (BD(1,1)~=BD(end,1)) || (BD(1,2)~=BD(end,2))
    NB=size(BD,1);
    BD2=[BD;BD(1,:)];
else %给出的边界已经收尾相接
    NB=size(BD,1)-1;
    BD2=BD;
end
%删除边界中相邻重复的点
IsSame=and(BD2(1:end-1,1)==BD2(2:end,1),BD2(1:end-1,2)==BD2(2:end,2));
BD2([IsSame;false],:)=[];
NB=size(BD2,1)-1;
NP=size(xy2,1);

IsInPoly=false(NP,1);%true(NP,1);%false(NP,1);
IsOnBD=false(NP,1);%是否在边线上
%判断整个边线的大概范围
min_X_BD