V_POLYGONWIND Test if points are inside a polygon Inputs: P(n,2) polygon vertices X(m,2) points to test Outputs: W(m) winding number for each point For a normal polygon whose vertices are listed anti-clockwise, the winding number is 0 or 1 according to whether the point is outside or inside the polygon. The winding number will be -1 if the polygon vertices go colckwise around the point.
0001 function w=v_polygonwind(p,x) 0002 %V_POLYGONWIND Test if points are inside a polygon 0003 % Inputs: 0004 % P(n,2) polygon vertices 0005 % X(m,2) points to test 0006 % 0007 % Outputs: 0008 % W(m) winding number for each point 0009 % 0010 % For a normal polygon whose vertices are listed 0011 % anti-clockwise, the winding number is 0 or 1 according to whether the 0012 % point is outside or inside the polygon. The winding number will 0013 % be -1 if the polygon vertices go colckwise around the point. 0014 0015 0016 % Copyright (C) Mike Brookes 2009 0017 % Version: $Id: v_polygonwind.m 10865 2018-09-21 17:22:45Z dmb $ 0018 % 0019 % VOICEBOX is a MATLAB toolbox for speech processing. 0020 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0021 % 0022 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0023 % This program is free software; you can redistribute it and/or modify 0024 % it under the terms of the GNU General Public License as published by 0025 % the Free Software Foundation; either version 2 of the License, or 0026 % (at your option) any later version. 0027 % 0028 % This program is distributed in the hope that it will be useful, 0029 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0030 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0031 % GNU General Public License for more details. 0032 % 0033 % You can obtain a copy of the GNU General Public License from 0034 % http://www.gnu.org/copyleft/gpl.html or by writing to 0035 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0036 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0037 n=size(p,1); 0038 m=size(x,1); 0039 q=zeros(2,n+1); 0040 q(:,1:n)=p'; 0041 q(:,n+1)=q(:,1); % append an extra point 0042 i=1:n; 0043 j=2:n+1; 0044 ym=repmat(2,m,1); 0045 yn=repmat(2,1,n); 0046 w=sum((2*((repmat(q(1,i).*q(2,j)-q(2,i).*q(1,j),m,1)+x(:,1)*(q(2,i)-q(2,j))+x(:,2)*(q(1,j)-q(1,i)))>0)-1).*abs((q(ym,j)>x(:,yn))-(q(ym,i)>x(:,yn))),2)/2; 0047 if ~nargout 0048 w0=w==0; 0049 wp=w>0; 0050 wn=w<0; 0051 plot(q(1,:),q(2,:),'k-',x(w0,1),x(w0,2),'go',x(wp,1),x(wp,2),'r+',x(wn,1),x(wn,2),'bx'); 0052 mnx=[1.05 -0.05;-0.05 1.05]*[min([p; x]);max([p; x])]; 0053 set(gca,'xlim',mnx(:,1)','ylim',mnx(:,2)'); 0054 title('Winding numbers: o=0, +=pos, x=neg'); 0055 end