Mercurial > octave-antonio

comparison scripts/geometry/inpolygon.m @ 20133:bef9086a6725

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
inpolygon.m: Overhaul function and fix bug #40970. * NEWS: Announce changes in inpolygon definition. * inpolygon.m: Use 'in |= on' to include points on the boundary since this is what Matlab does. Rewrite docstring with new definition. Change error() messages to be more explicit. Add more input validation tests.
author Rik <rik@octave.org>
date Sat, 25 Apr 2015 20:50:36 -0700
parents 9fc020886ae9
children 7503499a252b
comparison
equal deleted inserted replaced
20132:d7bea5b11fc3 20133:bef9086a6725
16 ## You should have received a copy of the GNU General Public License 16 ## You should have received a copy of the GNU General Public License
17 ## along with Octave; see the file COPYING. If not, see 17 ## along with Octave; see the file COPYING. If not, see
18 ## <http://www.gnu.org/licenses/>. 18 ## <http://www.gnu.org/licenses/>.
19 19
20 ## -*- texinfo -*- 20 ## -*- texinfo -*-
21 ## @deftypefn {Function File} {[@var{in}, @var{on}] =} inpolygon (@var{x}, @var{y}, @var{xv}, @var{yv}) 21 ## @deftypefn {Function File} {@var{in} =} inpolygon (@var{x}, @var{y}, @var{xv}, @var{yv})
22 ## @deftypefnx {Function File} {[@var{in}, @var{on}] =} inpolygon (@var{x}, @var{y}, @var{xv}, @var{yv})
22 ## 23 ##
23 ## For a polygon defined by vertex points @code{(@var{xv}, @var{yv})}, determine 24 ## For a polygon defined by vertex points @code{(@var{xv}, @var{yv})}, return
24 ## if the points @code{(@var{x}, @var{y})} are inside or outside the polygon. 25 ## true if the points @code{(@var{x}, @var{y})} are inside (or on the boundary)
26 ## of the polygon; Otherwise, return false.
27 ##
25 ## The variables @var{x}, @var{y}, must have the same dimension. The optional 28 ## 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. 29 ## output @var{on} returns true if the points are exactly on the polygon
27 ## 30 ## edge, and false otherwise.
31 ## @seealso{delaunay}
28 ## @end deftypefn 32 ## @end deftypefn
29 33
30 ## Author: Frederick (Rick) A Niles <niles@rickniles.com> 34 ## Author: Frederick (Rick) A Niles <niles@rickniles.com>
31 ## Created: 14 November 2006 35 ## Created: 14 November 2006
32 36
33 ## Vectorized by Søren Hauberg <soren@hauberg.org> 37 ## Vectorized by Søren Hauberg <soren@hauberg.org>
34 38
35 ## The method for determining if a point is in in a polygon is based on 39 ## The method for determining if a point is in in a polygon is based on
36 ## the algorithm shown on 40 ## the algorithm shown on
37 ## http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/ and is 41 ## http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/
38 ## credited to Randolph Franklin. 42 ## and is credited to Randolph Franklin.
39 43
40 function [in, on] = inpolygon (x, y, xv, yv) 44 function [in, on] = inpolygon (x, y, xv, yv)
41 45
42 if (nargin != 4) 46 if (nargin != 4)
43 print_usage (); 47 print_usage ();
44 endif 48 endif
45 49
46 if (! (isreal (x) && isreal (y) && ismatrix (y) && ismatrix (y) 50 if (! (isreal (x) && isreal (y) && ismatrix (y) && ismatrix (y)
47 && size_equal (x, y))) 51 && size_equal (x, y)))
48 error ("inpolygon: first two arguments must be real matrices of same size"); 52 error ("inpolygon: X and Y must be real matrices of the same size");
49 elseif (! (isreal (xv) && isreal (yv) && isvector (xv) && isvector (yv) 53 elseif (! (isreal (xv) && isreal (yv) && isvector (xv) && isvector (yv)
50 && size_equal (xv, yv))) 54 && size_equal (xv, yv)))
51 error ("inpolygon: last two arguments must be real vectors of same size"); 55 error ("inpolygon: XV and YV must be real vectors of the same size");
52 endif 56 endif
53 57
54 npol = length (xv); 58 npol = length (xv);
55 do_boundary = (nargout >= 2);
56 59
57 in = zeros (size (x), "logical"); 60 in = on = false (size (x));
58 if (do_boundary)
59 on = zeros (size (x), "logical");
60 endif
61 61
62 j = npol; 62 j = npol;
63 for i = 1 : npol 63 for i = 1 : npol
64 delta_xv = xv(j) - xv(i); 64 delta_xv = xv(j) - xv(i);
65 delta_yv = yv(j) - yv(i); 65 delta_yv = yv(j) - yv(i);
66 ## distance = [distance from (x,y) to edge] * length(edge) 66 ## distance = [distance from (x,y) to edge] * length(edge)
67 distance = delta_xv .* (y - yv(i)) - (x - xv(i)) .* delta_yv; 67 distance = delta_xv .* (y - yv(i)) - (x - xv(i)) .* delta_yv;
68 ## 68
69 ## is y between the y-values of edge i,j 69 ## is y between the y-values of edge i,j AND (x,y) on the left of the edge?
70 ## AND (x,y) on the left of the edge ?
71 idx1 = (((yv(i) <= y & y < yv(j)) | (yv(j) <= y & y < yv(i))) 70 idx1 = (((yv(i) <= y & y < yv(j)) | (yv(j) <= y & y < yv(i)))
72 & 0 < distance.*delta_yv); 71 & 0 < distance.*delta_yv);
73 in(idx1) = ! in(idx1); 72 in(idx1) = ! in(idx1);
74 73
75 ## Check if (x,y) are actually on the boundary of the polygon. 74 ## Check if (x,y) are actually on the boundary of the polygon.
76 if (do_boundary) 75 idx2 = (((yv(i) <= y & y <= yv(j)) | (yv(j) <= y & y <= yv(i)))
77 idx2 = (((yv(i) <= y & y <= yv(j)) | (yv(j) <= y & y <= yv(i))) 76 & ((xv(i) <= x & x <= xv(j)) | (xv(j) <= x & x <= xv(i)))
78 & ((xv(i) <= x & x <= xv(j)) | (xv(j) <= x & x <= xv(i))) 77 & (0 == distance | !delta_xv));
79 & (0 == distance | !delta_xv)); 78 on(idx2) = true;
80 on(idx2) = true; 79
81 endif
82 j = i; 80 j = i;
83 endfor 81 endfor
82
83 ## Matlab definition include both in polygon and on polygon points.
84 in |= on;
84 85
85 endfunction 86 endfunction
86 87
87 88
88 %!demo 89 %!demo
117 %! 0.18161, 0.78850, 1.13589, 1.33781, 1.04650, ... 118 %! 0.18161, 0.78850, 1.13589, 1.33781, 1.04650, ...
118 %! 0.60628, 0.82096, 0.67155, 0.96114, 1.14833, ... 119 %! 0.60628, 0.82096, 0.67155, 0.96114, 1.14833, ...
119 %! 0.82096, 0.60628]; 120 %! 0.82096, 0.60628];
120 %! xa = [0:0.1:2.3]; 121 %! xa = [0:0.1:2.3];
121 %! ya = [0:0.1:1.4]; 122 %! ya = [0:0.1:1.4];
122 %! [x,y] = meshgrid (xa, ya); 123 %! [x, y] = meshgrid (xa, ya);
123 %! [in,on] = inpolygon (x, y, xv, yv); 124 %! [in, on] = inpolygon (x, y, xv, yv);
124 %! inside = in & !on; 125 %! inside = in & !on;
125 %! 126 %!
126 %! clf; 127 %! clf;
127 %! plot (xv, yv); 128 %! plot (xv, yv);
128 %! hold on; 129 %! hold on;
132 %! hold off; 133 %! hold off;
133 %! disp ("Green points are inside polygon, magenta are outside,"); 134 %! disp ("Green points are inside polygon, magenta are outside,");
134 %! disp ("and blue are on boundary."); 135 %! disp ("and blue are on boundary.");
135 136
136 %!test 137 %!test
137 %! [in, on] = inpolygon ([1, 0], [1, 0], [-1, -1, 1, 1], [-1, 1, 1, -1]); 138 %! [in, on] = inpolygon ([1, 0, 2], [1, 0, 0], [-1, -1, 1, 1], [-1, 1, 1, -1]);
138 %! assert (in, [false, true]); 139 %! assert (in, [true, true, false]);
139 %! assert (on, [true, false]); 140 %! assert (on, [true, false, false]);
140 141
141 ## Test input validation 142 ## Test input validation
142 %!error inpolygon () 143 %!error inpolygon ()
143 %!error inpolygon (1, 2) 144 %!error inpolygon (1, 2)
144 %!error inpolygon (1, 2, 3) 145 %!error inpolygon (1, 2, 3)
145 %!error inpolygon (1, [1,2], [3, 4], [5, 6]) 146 %!error inpolygon (1, 2, 3, 4, 5)
146 %!error inpolygon ([1,2], [3, 4], [5, 6], 1) 147 %!error <X and Y must be real matrices> inpolygon (1i, 1, [3, 4], [5, 6])
148 %!error <X and Y must be real matrices> inpolygon (1, {1}, [3, 4], [5, 6])
149 %!error <X and Y must be .* the same size> inpolygon (1, [1,2], [3, 4], [5, 6])
150 %!error <XV and YV must be real vectors> inpolygon (1, 1, [3i, 4], [5, 6])
151 %!error <XV and YV must be real vectors> inpolygon (1, 1, [3, 4], {5, 6})
152 %!error <XV and YV must .* the same size> inpolygon ([1,2], [3, 4], [5, 6], 1)
147 153