--- /dev/null
+## Copyright (C) 2006-2012 Frederick (Rick) A Niles
+## and S�ren Hauberg
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{in}, @var{on}] =} inpolygon (@var{x}, @var{y}, @var{xv}, @var{yv})
+##
+## For a polygon defined by vertex points @code{(@var{xv}, @var{yv})}, determine
+## if the points @code{(@var{x}, @var{y})} are inside or outside the polygon.
+## The variables @var{x}, @var{y}, must have the same dimension. The optional
+## output @var{on} gives the points that are on the polygon.
+##
+## @end deftypefn
+
+## Author: Frederick (Rick) A Niles <niles@rickniles.com>
+## Created: 14 November 2006
+
+## Vectorized by S�ren Hauberg <soren@hauberg.org>
+
+## The method for determining if a point is in in a polygon is based on
+## the algorithm shown on
+## http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/ and is
+## credited to Randolph Franklin.
+
+function [in, on] = inpolygon (x, y, xv, yv)
+
+ if (nargin != 4)
+ print_usage ();
+ endif
+
+ if (! (isreal (x) && isreal (y) && ismatrix (y) && ismatrix (y)
+ && size_equal (x, y)))
+ error ("inpolygon: first two arguments must be real matrices of same size");
+ elseif (! (isreal (xv) && isreal (yv) && isvector (xv) && isvector (yv)
+ && size_equal (xv, yv)))
+ error ("inpolygon: last two arguments must be real vectors of same size");
+ endif
+
+ npol = length (xv);
+ do_boundary = (nargout >= 2);
+
+ in = zeros (size(x), "logical");
+ if (do_boundary)
+ on = zeros (size(x), "logical");
+ endif
+
+ j = npol;
+ for i = 1 : npol
+ delta_xv = xv(j) - xv(i);
+ delta_yv = yv(j) - yv(i);
+ ## distance = [distance from (x,y) to edge] * length(edge)
+ distance = delta_xv .* (y - yv(i)) - (x - xv(i)) .* delta_yv;
+ ##
+ ## is y between the y-values of edge i,j
+ ## AND (x,y) on the left of the edge ?
+ idx1 = (((yv(i) <= y & y < yv(j)) | (yv(j) <= y & y < yv(i)))
+ & 0 < distance.*delta_yv);
+ in (idx1) = !in (idx1);
+
+ ## Check if (x,y) are actually on the boundary of the polygon.
+ if (do_boundary)
+ idx2 = (((yv(i) <= y & y <= yv(j)) | (yv(j) <= y & y <= yv(i)))
+ & ((xv(i) <= x & x <= xv(j)) | (xv(j) <= x & x <= xv(i)))
+ & (0 == distance | !delta_xv));
+ on (idx2) = true;
+ endif
+ j = i;
+ endfor
+
+endfunction
+
+%!demo
+%! xv=[ 0.05840, 0.48375, 0.69356, 1.47478, 1.32158, \
+%! 1.94545, 2.16477, 1.87639, 1.18218, 0.27615, \
+%! 0.05840 ];
+%! yv=[ 0.60628, 0.04728, 0.50000, 0.50000, 0.02015, \
+%! 0.18161, 0.78850, 1.13589, 1.33781, 1.04650, \
+%! 0.60628 ];
+%! xa=[0:0.1:2.3];
+%! ya=[0:0.1:1.4];
+%! [x,y]=meshgrid(xa,ya);
+%! [in,on]=inpolygon(x,y,xv,yv);
+%!
+%! inside=in & !on;
+%! plot(xv,yv)
+%! hold on
+%! plot(x(inside),y(inside),"@g")
+%! plot(x(~in),y(~in),"@m")
+%! plot(x(on),y(on),"@b")
+%! hold off
+%! disp("Green points are inside polygon, magenta are outside,");
+%! disp("and blue are on boundary.");
+
+%!demo
+%! xv=[ 0.05840, 0.48375, 0.69356, 1.47478, 1.32158, \
+%! 1.94545, 2.16477, 1.87639, 1.18218, 0.27615, \
+%! 0.05840, 0.73295, 1.28913, 1.74221, 1.16023, \
+%! 0.73295, 0.05840 ];
+%! yv=[ 0.60628, 0.04728, 0.50000, 0.50000, 0.02015, \
+%! 0.18161, 0.78850, 1.13589, 1.33781, 1.04650, \
+%! 0.60628, 0.82096, 0.67155, 0.96114, 1.14833, \
+%! 0.82096, 0.60628];
+%! xa=[0:0.1:2.3];
+%! ya=[0:0.1:1.4];
+%! [x,y]=meshgrid(xa,ya);
+%! [in,on]=inpolygon(x,y,xv,yv);
+%!
+%! inside=in & ~ on;
+%! plot(xv,yv)
+%! hold on
+%! plot(x(inside),y(inside),"@g")
+%! plot(x(~in),y(~in),"@m")
+%! plot(x(on),y(on),"@b")
+%! hold off
+%! disp("Green points are inside polygon, magenta are outside,");
+%! disp("and blue are on boundary.");
+
+%!error inpolygon ();
+%!error inpolygon (1, 2);
+%!error inpolygon (1, 2, 3);
+
+%!error inpolygon (1, [1,2], [3, 4], [5, 6]);
+%!error inpolygon ([1,2], [3, 4], [5, 6], 1);
+
+%!test
+%! [in, on] = inpolygon ([1, 0], [1, 0], [-1, -1, 1, 1], [-1, 1, 1, -1]);
+%! assert (in, [false, true]);
+%! assert (on, [true, false]);
git://git.creatis.insa-lyon.fr / CreaPhase.git / blobdiff