利用Matlab判断某些点是否在多边形区域内

目录

 1. 代码作用及实现效果

 2. 技术分析：        

  3. 程序

1. 代码作用及实现效果

        先给各位看看具体效果，如下所示，其中红色的点表示需要判断的点，是否在蓝色区域内，从图中可知，有两个点在蓝色区域内，一个点在边界线上，一个在蓝色区域外。

 2. 技术分析：        

       

原理：射线法；从该点假想一条通向无穷远的射线，通过判断射线与多边形交点的个数奇偶来确定点与多边形的位置关系。具体原理详解可自行百度；

 可以利用Matlab内置的函数inpolygon，也可以写一个函数，加强自己的理解；

inpolygon函数用法：

   in = inpolygon(xq,yq,xv,yv)
   [in,on] = inpolygon(xq,yq,xv,yv)

  3. 程序

        下面是我根据射线法写的代码：该代码可以判断点是否在这个区域（包含边界线）上，是则输出1；反之则输出0。然后提取输出为1的索引，并找出对应的点坐标；

function Flags = inPoly(Site_Ponit,Area)
%==================================说明====================================
%   判定点目标是否在区域内函数
%   函数作用：判断点坐标是否在所需的区域内
%
%   作者：胡礼珍
%   单位：厦门大学联合遥感接收站
%   邮件：hulizhen@xmu.edu.cn
%
%   输入：
%       Site_Ponit      需要判断的点的位置坐标，n×2矩阵；
%       Area            面状区域，多边形各个顶点位置坐标(可以是经纬度)，n×2矩阵；
%   输出：
%       OutImage    输出一个二维矩阵
%       ImageName   将输出的矩阵保存为文件
%
%   语法：
%   flags = inPoly(Site_Ponit,Area)：
%           Site_Ponit表示需要判断的点的位置坐标，n×2矩阵；Area表示面状区域，
%           多边形各个顶点位置坐标(可以是经纬度)，m×2矩阵；
%--------------------------------------------------------------------------
%   参考博客：https://blog.csdn.net/ye_xiao_yu/article/details/72765726
%==========================================================================
if ~(Area(1,:) == Area(end,:))
    errordlg('首尾点不相等，请确认区域是否为面状区域。如果是，请忽略；否则，请核对重新运行；','警告：输入可能有误')
    Area = [Area;Area(1,:)];
end

Num_Piont = size(Site_Ponit,1);     % Site_Ponit的行数，也表示检测点的个数
Num_Area = size(Area,1);            % Area的行数，也表示多边形顶点数+1
Flags = zeros(1,Num_Piont);         % 新建存储空间，存放判定结果0和1，初始全为0
% h = waitbar(0,'运行……');          % 新建进度条
for i=1:Num_Piont
    % 添加一个进度条，显示进展
%    per = fix(i/Num_Piont*10000)/100;                  % 完成度
%    str=['点目标判断函数-完成： ',num2str(per),' % '];  % 进度条上显示的内容
%    waitbar(i/Num_Piont,h,str);
    if ~isempty(find(Area(:,1)==Site_Ponit(i,1)& ...
            Area(:,2)==Site_Ponit(i,2), 1))             % 判定点是否在顶点上
        Flags(i) = 1;               % 符合要求，判定为1
        continue;                   % 结束本次循环，进行下一次循环
    end
    % 射线法：从该点假想一条通向无穷远的射线，通过判断射线与多边形交点的个数奇偶来确定点与多边形的位置关系。
    % 参考：https://blog.csdn.net/bitcarmanlee/article/details/60591192
    for j=2:Num_Area
        % 条件A= 条件01 或 条件02
        %      = (条件11 与 条件12)
                % 条件01 = (条件11 与 条件12) 或
                %        = (Area(j,2)<=Site_Ponit(i,2)) && (Site_Ponit(i,2) < Area(j-1,2))
                    % 条件11：(Area(j,2)<=Site_Ponit(i,2))
                    % 条件12：(Site_Ponit(i,2) < Area(j-1,2))
                % 条件02 = (条件21 与 条件22)
                %        = ((Area(j-1,2) <= Site_Ponit(i,2)) && (Site_Ponit(i,2) < Area(j,2)))
                    % 条件21：(Area(j-1,2) <= Site_Ponit(i,2))
                    % 条件22：(Site_Ponit(i,2) < Area(j,2))
        % 条件B：Site_Ponit(i,1) < (Area(j-1,1) - Area(j,1)) * (Site_Ponit(i,2) - Area(j,2))/(Area(j-1,2)-Area(j,2)) + Area(j,1);
        if (...
                (...        % 判断：条件A = 条件01 或 条件02 
                    (...        % 判断：条件01 =(条件11 与 条件12)
                        (Area(j,2)<=Site_Ponit(i,2)) && (Site_Ponit(i,2) < Area(j-1,2))...
                     )...
                    ||...   % ********************“或”*******************
                    (...        % 判断：条件02 = (条件21 与 条件22)
                        (Area(j-1,2) <= Site_Ponit(i,2)) && (Site_Ponit(i,2) < Area(j,2))...
                     )...
                 )...   % ----------------------条件A----------------------
                && ...  % **********************“与”*********************
                (...
                    Site_Ponit(i,1) < ...
                        (Area(j-1,1) - Area(j,1)) * (Site_Ponit(i,2) - Area(j,2))/...
                        (Area(j-1,2)-Area(j,2)) +...
                        Area(j,1)...
                 )...   % ----------------------条件B----------------------
             )
            Flags(i) = Flags(i) + 1;    % 交点数累加，最后根据奇偶判定是否在其中
        end
    end
end
% close(h);               % 关闭进度条
Flags = mod(Flags,2);   % 分裂后剩余，判断奇数偶数，奇数则表示在，偶数则表示不在
end

        调用验证上面函数的代码如下：

clear;close all;clc;
% 设置3个点，第1个点和第二个点在区域内，第三个点不在区域内，第四个点在边界线上
P = [116.5412066049,25.2741767057;
     116.2906735589,25.4947457375;
     115.5006735589,24.4947457375;
     117.9292952035,26.1752619468];      % 需要检测的点
Poly = [118.0562474131,25.521042552;
        117.9292952035,26.1752619468;
        115.4280709351,25.7592188999;
        115.7898649526,24.0756121845;
        117.7748394579,24.4133199012;
        117.5551745295,25.4373847646;
        118.0562474131,25.521042552];   % 判定区域的范围  
% 第一种方法，自己编写函数
tic
flags = inPoly(P,Poly);                 % 判定点目标是否在区域内函数
In_P = P(flags== 1,:)            % 根据判定，找出在判定区域内的点坐标
toc

figure('Name','效果图1')

fill(Poly(:,1),Poly(:,2),'b');
title('判断点目标是否在该区域范围内')
hold on
scatter(P(:,1),P(:,2),'r.')
hold off

tic
% 第二种方法，用内置函数 inpolygon
flags = inpolygon(P(:,1),P(:,2),Poly(:,1),Poly(:,2));
In_P = P(flags== 1,:)            % 根据判定，找出在判定区域内的点坐标
toc                              % 
% 自己写的，速度差不多,说明还是成功的

不足之处，敬请斧正！

转载请说明出处！

路漫漫其修远兮，吾将上下而求索