From octave-sources-request at bevo dot che dot wisc dot edu  Fri Mar 26 13:02:02 1999
Subject: inpoly.m clean up
From: David Doolin <doolin at cs dot utk dot edu>
To: octave-sources at bevo dot che dot wisc dot edu
cc: doolin at cs dot utk dot edu
Date: Fri, 26 Mar 1999 14:01:17 -0500


Code clean, documentation clean, minor 
corrections.


----------------snip snap------------

% [ 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