octave_packages/m/geometry/inpolygon.m

   1 ## Copyright (C) 2006-2012 Frederick (Rick) A Niles
   2 ##               and S�ren Hauberg
   3 ##
   4 ## This file is part of Octave.
   5 ##
   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.
  10 ##
  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.
  15 ##
  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/>.
  19
  20 ## -*- texinfo -*-
  21 ## @deftypefn {Function File} {[@var{in}, @var{on}] =} inpolygon (@var{x}, @var{y}, @var{xv}, @var{yv})
  22 ##
  23 ## For a polygon defined by vertex points @code{(@var{xv}, @var{yv})}, 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.
  27 ##
  28 ## @end deftypefn
  29
  30 ## Author: Frederick (Rick) A Niles <niles@rickniles.com>
  31 ## Created: 14 November 2006
  32
  33 ## Vectorized by S�ren Hauberg <soren@hauberg.org>
  34
  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.
  39
  40 function [in, on] = inpolygon (x, y, xv, yv)
  41
  42   if (nargin != 4)
  43     print_usage ();
  44   endif
  45
  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");
  52   endif
  53
  54   npol = length (xv);
  55   do_boundary = (nargout >= 2);
  56
  57   in = zeros (size(x), "logical");
  58   if (do_boundary)
  59     on = zeros (size(x), "logical");
  60   endif
  61
  62   j = npol;
  63   for i = 1 : npol
  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;
  68     ##
  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);
  74
  75     ## Check if (x,y) are actually on the boundary of the polygon.
  76     if (do_boundary)
  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));
  80        on (idx2) = true;
  81     endif
  82     j = i;
  83   endfor
  84
  85 endfunction
  86
  87 %!demo
  88 %!  xv=[ 0.05840, 0.48375, 0.69356, 1.47478, 1.32158, \
  89 %!       1.94545, 2.16477, 1.87639, 1.18218, 0.27615, \
  90 %!       0.05840 ];
  91 %!  yv=[ 0.60628, 0.04728, 0.50000, 0.50000, 0.02015, \
  92 %!       0.18161, 0.78850, 1.13589, 1.33781, 1.04650, \
  93 %!       0.60628 ];
  94 %! xa=[0:0.1:2.3];
  95 %! ya=[0:0.1:1.4];
  96 %! [x,y]=meshgrid(xa,ya);
  97 %! [in,on]=inpolygon(x,y,xv,yv);
  98 %!
  99 %! inside=in & !on;
 100 %! plot(xv,yv)
 101 %! hold on
 102 %! plot(x(inside),y(inside),"@g")
 103 %! plot(x(~in),y(~in),"@m")
 104 %! plot(x(on),y(on),"@b")
 105 %! hold off
 106 %! disp("Green points are inside polygon, magenta are outside,");
 107 %! disp("and blue are on boundary.");
 108
 109 %!demo
 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, \
 117 %!       0.82096, 0.60628];
 118 %! xa=[0:0.1:2.3];
 119 %! ya=[0:0.1:1.4];
 120 %! [x,y]=meshgrid(xa,ya);
 121 %! [in,on]=inpolygon(x,y,xv,yv);
 122 %!
 123 %! inside=in & ~ on;
 124 %! plot(xv,yv)
 125 %! hold on
 126 %! plot(x(inside),y(inside),"@g")
 127 %! plot(x(~in),y(~in),"@m")
 128 %! plot(x(on),y(on),"@b")
 129 %! hold off
 130 %! disp("Green points are inside polygon, magenta are outside,");
 131 %! disp("and blue are on boundary.");
 132
 133 %!error inpolygon ();
 134 %!error inpolygon (1, 2);
 135 %!error inpolygon (1, 2, 3);
 136
 137 %!error inpolygon (1, [1,2], [3, 4], [5, 6]);
 138 %!error inpolygon ([1,2], [3, 4], [5, 6], 1);
 139
 140 %!test
 141 %! [in, on] = inpolygon ([1, 0], [1, 0], [-1, -1, 1, 1], [-1, 1, 1, -1]);
 142 %! assert (in, [false, true]);
 143 %! assert (on, [true, false]);