1## Copyright (C) 2006, 2007, 2008, 2009 Frederick (Rick) A Niles
4## This file is part of Octave.
6## Octave is free software; you can redistribute it and/or modify it
7## under the terms of the GNU General Public License as published by
8## the Free Software Foundation; either version 3 of the License, or (at
9## your option) any later version.
11## Octave is distributed in the hope that it will be useful, but
12## WITHOUT ANY WARRANTY; without even the implied warranty of
13## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14## General Public License for more details.
16## You should have received a copy of the GNU General Public License
17## along with Octave; see the file COPYING. If not, see
18## <http://www.gnu.org/licenses/>.
21## @deftypefn {Function File} {[@var{in}, @var{on}] =} inpolygon (@var{x}, @var{y}, @var{xv}, @var{xy})
23## For a polygon defined by @code{(@var{xv}, @var{yv})} points, determine
24## if the points @code{(@var{x}, @var{y})} are inside or outside the polygon.
25## The variables @var{x}, @var{y}, must have the same dimension. The optional
26## output @var{on} gives the points that are on the polygon.
30## Author: Frederick (Rick) A Niles <niles@rickniles.com>
31## Created: 14 November 2006
33## Vectorized by S�ren Hauberg <soren@hauberg.org>
35## The method for determining if a point is in in a polygon is based on
36## the algorithm shown on
37## http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/ and is
38## credited to Randolph Franklin.
40function [IN, ON] = inpolygon (X, Y, xv, yv)
46 if (! (isreal (X) && isreal (Y) && ismatrix (Y) && ismatrix (Y)
47 && size_equal (X, Y)))
48 error ("inpolygon: first two arguments must be real matrices of same size");
49 elseif (! (isreal (xv) && isreal (yv) && isvector (xv) && isvector (yv)
50 && size_equal (xv, yv)))
51 error ("inpolygon: last two arguments must be real vectors of same size");
55 do_boundary = (nargout >= 2);
57 IN = zeros (size(X), "logical");
59 ON = zeros (size(X), "logical");
64 delta_xv = xv(j) - xv(i);
65 delta_yv = yv(j) - yv(i);
66 ## distance = [distance from (X,Y) to edge] * length(edge)
67 distance = delta_xv .* (Y - yv(i)) - (X - xv(i)) .* delta_yv;
69 ## is Y between the y-values of edge i,j
70 ## AND (X,Y) on the left of the edge ?
71 idx1 = (((yv(i) <= Y & Y < yv(j)) | (yv(j) <= Y & Y < yv(i)))
72 & 0 < distance.*delta_yv);
73 IN (idx1) = !IN (idx1);
75 ## Check if (X,Y) are actually ON the boundary of the polygon.
77 idx2 = (((yv(i) <= Y & Y <= yv(j)) | (yv(j) <= Y & Y <= yv(i)))
78 & ((xv(i) <= X & X <= xv(j)) | (xv(j) <= X & X <= xv(i)))
79 & (0 == distance | !delta_xv));
88%! xv=[ 0.05840, 0.48375, 0.69356, 1.47478, 1.32158, \
89%! 1.94545, 2.16477, 1.87639, 1.18218, 0.27615, \
91%! yv=[ 0.60628, 0.04728, 0.50000, 0.50000, 0.02015, \
92%! 0.18161, 0.78850, 1.13589, 1.33781, 1.04650, \
96%! [x,y]=meshgrid(xa,ya);
97%! [IN,ON]=inpolygon(x,y,xv,yv);
102%! plot(x(inside),y(inside),"@g")
103%! plot(x(~IN),y(~IN),"@m")
104%! plot(x(ON),y(ON),"@b")
106%! disp("Green points are inside polygon, magenta are outside,");
107%! disp("and blue are on boundary.");
110%! xv=[ 0.05840, 0.48375, 0.69356, 1.47478, 1.32158, \
111%! 1.94545, 2.16477, 1.87639, 1.18218, 0.27615, \
112%! 0.05840, 0.73295, 1.28913, 1.74221, 1.16023, \
113%! 0.73295, 0.05840 ];
114%! yv=[ 0.60628, 0.04728, 0.50000, 0.50000, 0.02015, \
115%! 0.18161, 0.78850, 1.13589, 1.33781, 1.04650, \
116%! 0.60628, 0.82096, 0.67155, 0.96114, 1.14833, \
120%! [x,y]=meshgrid(xa,ya);
121%! [IN,ON]=inpolygon(x,y,xv,yv);
126%! plot(x(inside),y(inside),"@g")
127%! plot(x(~IN),y(~IN),"@m")
128%! plot(x(ON),y(ON),"@b")
130%! disp("Green points are inside polygon, magenta are outside,");
131%! disp("and blue are on boundary.");
