Thread: Is a point inside or outside a polygon?

Hi All,

My first post, any help greatly appreciated.

In need to calculate if a point is inside or outside a polygon.

The application is to calculate which oceanic polygon a ship is in, or conversly, which ships are within a certain polygon.

The vertices of the polygon will be in degrees lat and long, but this should translate fine to a simple 2D polygon on a cartesian diagram. The ship's position will also be in degrees lat/long, again should corresponde fine to cartesian coordinates.

I realise that it's actually a polygon on a sphere i should be calculating, but for this i'm happy to assume the edges are straight.

I've added a diagram to show what I mean. In the first case, the point P is inside the polygon, in the second case it's outside.

How do I determine which?? Thanks everyone,

Luke

Originally Posted by LukeCoverdale

Hi All,

My first post, any help greatly appreciated.

In need to calculate if a point is inside or outside a polygon.

The application is to calculate which oceanic polygon a ship is in, or conversly, which ships are within a certain polygon.

The vertices of the polygon will be in degrees lat and long, but this should translate fine to a simple 2D polygon on a cartesian diagram. The ship's position will also be in degrees lat/long, again should corresponde fine to cartesian coordinates.

I realise that it's actually a polygon on a sphere i should be calculating, but for this i'm happy to assume the edges are straight.

I've added a diagram to show what I mean. In the first case, the point P is inside the polygon, in the second case it's outside.

How do I determine which?? Thanks everyone,

Luke

If you have access to Matlab there is a function inpolygon that will do this, and similar in Octave the inpoly function should do this.

RonL

Thanks for your help.

Do you know of any actual code that would do it, like a function in Java or something?

The reason I ask is that I would need to do this calculation within my program. I wouldn't be able to do open Matlab or Octave each time.

Thanks,
Luke

Originally Posted by LukeCoverdale

Thanks for your help.

Do you know of any actual code that would do it, like a function in Java or something?

The reason I ask is that I would need to do this calculation within my program. I wouldn't be able to do open Matlab or Octave each time.

Thanks,
Luke

The Octave code is open source and should be easy to translate into any other language

Now I have not tried this out myself, and this looks different from my Euler code to do this but:

Code:

% [ inside ] = inpoly(polygon, point)
%
% Determines whether a given point is inside a polygon.
% The polygon is determined by point pairs on its 
% boundary.  There is no error checking, and to be safe
% only use convex polygons.
%
% input:
%       polygon: column vector of (x,y) pairs representing the vertex
%                points of the polygon.
%       point: (x,y) pair representing point that is either
%              in or out of the polygon
%
% output:
%       inside: boolean true for inside, false otherwise.
%
% todo:  
%       Implement the ability to handle an array of points
%       and return an array of boolean values.
%
% bugs and limitations:  
%        
%         Need to check whether convexity is a necessary 
%         condition for this algorithm (don't think it is).
%
%         Behavior with multiple loops, i.e., if sides cross
%         over, is unknown.  
%          
%         Algorithm implemented in integer, first quadrant
%         numbers.  Need to check generality of algorithm.
%
%         Points on the boundaries of polygons may not be handled
%         correctly or consistently.

%
%         This code was derived from a public domain code
%         copyright 1995-1996 Galacticomm, Inc., modifications
%         allowed for any purpose provided redistribution is 
%         not restricted
%
% This script properly belongs in octave/compgeom/
% (computational geometry).
%
% This modified code is copyright David M. Doolin and placed in
% the public domain, with the exception that redistribution 
% is not restricted in accordance with the copyright of unmodified
% (original) C code.
%
% $Author: doolin $  doolin at ce dot berkeley dot edu
% $Date: 1999/03/26 18:42:00 $
% $Source: /shag/homes/doolin/cvsroot/octave/compgeom/inpoly.m,v $
% $Revision: 1.2 $


function [ inside ] = inpoly(polygon, point)

% Check for the correct number of arguments

% Check for argument validation
% If argument validation required, validate arguments.

npoints = rows(polygon);
inside = 0;

xt = point(1);
yt = point(2);

xold = polygon(npoints,1);
yold = polygon(npoints,2);

for i = 1:1:npoints

   xnew = polygon(i,1);
   ynew = polygon(i,2);

      if (xnew > xold) 
         x1=xold;
         x2=xnew;
         y1=yold;
         y2=ynew;
      else 
         x1=xnew;
         x2=xold;
         y1=ynew;
         y2=yold;
       endif

      if ((xnew < xt) == (xt <= xold)         %/* edge "open" at left end */
           && (yt-y1)*(x2-x1) < (y2-y1)*(xt-x1) ) 
               inside=!inside;
      endif

      xold=xnew;
      yold=ynew;

end  % for loop

endfunction

RonL

Code:

function [ inside ] = inpoly(polygon,xt,yt)

rows = size(polygon);
npoints = rows(1);
disp (npoints);
inside = 0;

xold = polygon(npoints,1);
yold = polygon(npoints,2);

for i = 1:1:npoints

   xnew = polygon(i,1);
   ynew = polygon(i,2);

      if (xnew > xold) 
         x1=xold;
         x2=xnew;
         y1=yold;
         y2=ynew;
      else 
         x1=xnew;
         x2=xold;
         y1=ynew;
         y2=yold;
      end

      if ((xnew < xt) == (xt <= xold) & (yt-y1)*(x2-x1) < (y2-y1)*(xt-x1) ) 
               inside=~inside;
      end

      xold=xnew;
      yold=ynew;

end

endfunction