|
6847
|
1 ## Copyright (C) 2006 Frederick (Rick) A Niles, S�ren Hauberg |
|
|
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 |
|
|
41 if ( !(isreal(X) && isreal(Y) && ismatrix(Y) && ... |
|
|
42 ismatrix(Y) && size_equal(X,Y)) ) |
|
|
43 error ("inpolygon: first two arguments must be real matrices of same size"); |
|
|
44 elseif ( !(isreal(xv) && isreal(yv) && isvector(xv) && ... |
|
|
45 isvector(yv) && size_equal(xv,yv)) ) |
|
|
46 error ("inpolygon: last two arguments must be real vectors of same size"); |
|
|
47 endif |
|
|
48 |
|
|
49 npol = length(xv); |
|
|
50 do_boundary = (nargout >= 2); |
|
|
51 |
|
|
52 IN = zeros (size(X), "logical"); |
|
|
53 if (do_boundary) |
|
|
54 ON = zeros (size(X), "logical"); |
|
|
55 endif |
|
|
56 |
|
|
57 j = npol; |
|
|
58 for i = 1 : npol |
|
|
59 delta_xv = xv(j) - xv(i); |
|
|
60 delta_yv = yv(j) - yv(i); |
|
|
61 ## distance = [distance from (X,Y) to edge] * length(edge) |
|
|
62 distance = delta_xv .* (Y - yv(i)) - (X - xv(i)) .* delta_yv; |
|
|
63 ## |
|
|
64 ## is Y between the y-values of edge i,j |
|
|
65 ## AND (X,Y) on the left of the edge ? |
|
|
66 idx1 = ((yv(i) <= Y & Y < yv(j)) | (yv(j) <= Y & Y < yv(i)) ) & ... |
|
|
67 0 < distance.*delta_yv; |
|
|
68 IN (idx1) = !IN (idx1); |
|
|
69 |
|
|
70 ## Check if (X,Y) are actually ON the boundary of the polygon. |
|
|
71 if (do_boundary) |
|
|
72 idx2 = ((yv(i) <= Y & Y <= yv(j)) | (yv(j) <= Y & Y <= yv(i))) & ... |
|
|
73 ((xv(i) <= X & X <= xv(j)) | (xv(j) <= X & X <= xv(i))) & ... |
|
|
74 (0 == distance | !delta_xv); |
|
|
75 ON (idx2) = true; |
|
|
76 endif |
|
|
77 j = i; |
|
|
78 endfor |
|
|
79 endfunction |
|
|
80 |
|
|
81 %!demo |
|
|
82 %! xv=[ 0.05840, 0.48375, 0.69356, 1.47478, 1.32158, \ |
|
|
83 %! 1.94545, 2.16477, 1.87639, 1.18218, 0.27615, \ |
|
|
84 %! 0.05840 ]; |
|
|
85 %! yv=[ 0.60628, 0.04728, 0.50000, 0.50000, 0.02015, \ |
|
|
86 %! 0.18161, 0.78850, 1.13589, 1.33781, 1.04650, \ |
|
|
87 %! 0.60628 ]; |
|
|
88 %! xa=[0:0.1:2.3]; |
|
|
89 %! ya=[0:0.1:1.4]; |
|
|
90 %! [x,y]=meshgrid(xa,ya); |
|
|
91 %! [IN,ON]=inpolygon(x,y,xv,yv); |
|
|
92 %! |
|
|
93 %! inside=IN & !ON; |
|
|
94 %! plot(xv,yv) |
|
|
95 %! hold on |
|
|
96 %! plot(x(inside),y(inside),"@g") |
|
|
97 %! plot(x(~IN),y(~IN),"@m") |
|
|
98 %! plot(x(ON),y(ON),"@b") |
|
|
99 %! hold off |
|
|
100 %! disp("Green points are inside polygon, magenta are outside,"); |
|
|
101 %! disp("and blue are on boundary."); |
|
|
102 |
|
|
103 %!demo |
|
|
104 %! xv=[ 0.05840, 0.48375, 0.69356, 1.47478, 1.32158, \ |
|
|
105 %! 1.94545, 2.16477, 1.87639, 1.18218, 0.27615, \ |
|
|
106 %! 0.05840, 0.73295, 1.28913, 1.74221, 1.16023, \ |
|
|
107 %! 0.73295, 0.05840 ]; |
|
|
108 %! yv=[ 0.60628, 0.04728, 0.50000, 0.50000, 0.02015, \ |
|
|
109 %! 0.18161, 0.78850, 1.13589, 1.33781, 1.04650, \ |
|
|
110 %! 0.60628, 0.82096, 0.67155, 0.96114, 1.14833, \ |
|
|
111 %! 0.82096, 0.60628]; |
|
|
112 %! xa=[0:0.1:2.3]; |
|
|
113 %! ya=[0:0.1:1.4]; |
|
|
114 %! [x,y]=meshgrid(xa,ya); |
|
|
115 %! [IN,ON]=inpolygon(x,y,xv,yv); |
|
|
116 %! |
|
|
117 %! inside=IN & ~ ON; |
|
|
118 %! plot(xv,yv) |
|
|
119 %! hold on |
|
|
120 %! plot(x(inside),y(inside),"@g") |
|
|
121 %! plot(x(~IN),y(~IN),"@m") |
|
|
122 %! plot(x(ON),y(ON),"@b") |
|
|
123 %! hold off |
|
|
124 %! disp("Green points are inside polygon, magenta are outside,"); |
|
|
125 %! disp("and blue are on boundary."); |
|
|
126 |
