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