【原】判断点是否在多边形内(matlab版)

在网上，类似这种的理论很多，但是代码大部分是用C/C++写的。

在cnblog里面，有一个大牛用matlab实现了这个算法传送门
但是，身为程序猿的我并不满足他用脚本的方式运行。所以我改写成了函数，并且同时加入了快速判断。

function flags = inPoly(p,poly)
% 判断点是否在多边形内
% flag(i)为奇数，那么在，偶数为不在
% p : pn*2 矩阵
% poly ：polyn*2 矩阵，注意，起点和终点需要相同
% 致谢：代码部分引用自Dsp Tian的博客，原文链接
% http://www.cnblogs.com/tiandsp/p/4019880.html

if poly(1,:) ~= poly(end,:)
    poly = [poly;poly(1,:)];
end
pn = size(p,1);
polyn = size(poly,1);
flags = zeros(1,pn);
for i=1:pn
    
    t = find(poly(:,1)==p(i,1)& poly(:,2)==p(i,2), 1);
    if ~isempty(t)
        flags(i) = 1;
        continue;
    end
    for j=2:polyn
        x1=poly(j-1,1);         %多边形前后两个点
        y1=poly(j-1,2);
        x2=poly(j,1);
        y2=poly(j,2);
        
        k=(y1-y2)/(x1-x2);      %多边形一条边直线
        b=y1-k*x1;
        x=p(i,1);               %过当前点直线和多边形交点
        y=k*x+b;
        
        if min([x1 x2])<=x && x<=max([x1 x2]) && ...        %点同时在射线和多边形边上
                min([y1 y2])<=y && y<=max([y1 y2]) &&  y>=p(i,2)
            flags(i)=flags(i)+1;
        end
    end
end
flags = mod(flags,2);
end

然而，在使用中发现，这段代码有一个小小的bug。当
startPoint = [1.7,9.5]，
obj7 = [1.2,9.8;1.7,9.8;1.7,9.1;1.2,9.1;1.2,9.4;1.5,9.4;1.5,9.6;1.2,9.6;]
的时候，代码不能正常判断出这个点在多边形中(此Bug由朋友赵逸轩提出)。分析代码来看，事实上，这个点在多边形上,且这个边是垂直的，属于无斜率的特殊情况。

针对这一点，我参考了一篇博客传送门,这是用C写的一篇代码。
由此，修改后的代码为

function flags = inPoly(p,poly)
% 判断点是否在多边形内
% flag(i)为奇数，那function flags = inPoly(p,poly)
% 判断点是否在多边形内
% flag(i)为奇数，那么在，偶数为不在
% p : pn*2 矩阵
% poly ：polyn*2 矩阵，注意，起点和终点需要相同
% 致谢博客 http://www.cnblogs.com/dwdxdy/p/3230647.html


if ~(poly(1,1) == poly(end,1)&&poly(1,2) == poly(end,2))
    poly = [poly;poly(1,:)];
end

pn = size(p,1);
polyn = size(poly,1);
flags = zeros(1,pn);
for i=1:pn

    if ~isempty(find(poly(:,1)==p(i,1)& poly(:,2)==p(i,2)))%找到一个相同的点即可
        flags(i) = 1;
        continue;%%结束pn=1，进入pn=2
    end
    for j=2:polyn
        if ((((poly(j,2)<=p(i,2)) && (p(i,2) < poly(j-1,2) )) ||...
                ((poly(j-1,2) <= p(i,2)) && (p(i,2) < poly(j,2)))) && ...
                (p(i,1) < (poly(j-1,1) - poly(j,1)) * (p(i,2) - poly(j,2))/(poly(j-1,2)-poly(j,2)) + poly(j,1)))
            flags(i) = flags(i) + 1;
        end
    end
end
flags = mod(flags,2);
end

至此，代码bug完美解决。
PS：如果还有bug请各位大牛指出。欢迎讨论。