V_POLYGONXLINE Find where a line crosses a polygon [xc,ec,tc,xy0]=(p,l) Inputs P(n,2) gives the polygon vertices L(1,3) gives the line in the form L*[X; Y; 1]=0 Outputs xc(k,2) gives (x,y) coordinates of crossing points the crossings are between vertex pairs ec(k,1) and ec(k,1)+1 tc(k,1) gives parametric position of crossing points (x,y)=(x0,y0)+tc(-l(2),l(1)) xy0(1,2) gives the starting point, (x0,y0) of the parametric line
0001 function [xc,ec,tc,xy0]=v_polygonxline(p,l) 0002 %V_POLYGONXLINE Find where a line crosses a polygon [xc,ec,tc,xy0]=(p,l) 0003 % Inputs 0004 % P(n,2) gives the polygon vertices 0005 % L(1,3) gives the line in the form L*[X; Y; 1]=0 0006 % 0007 % Outputs 0008 % xc(k,2) gives (x,y) coordinates of crossing points 0009 % the crossings are between vertex pairs ec(k,1) and ec(k,1)+1 0010 % tc(k,1) gives parametric position of crossing points (x,y)=(x0,y0)+tc(-l(2),l(1)) 0011 % xy0(1,2) gives the starting point, (x0,y0) of the parametric line 0012 0013 % Copyright (C) Mike Brookes 2009 0014 % Version: $Id: v_polygonxline.m 10865 2018-09-21 17:22:45Z dmb $ 0015 % 0016 % VOICEBOX is a MATLAB toolbox for speech processing. 0017 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0018 % 0019 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0020 % This program is free software; you can redistribute it and/or modify 0021 % it under the terms of the GNU General Public License as published by 0022 % the Free Software Foundation; either version 2 of the License, or 0023 % (at your option) any later version. 0024 % 0025 % This program is distributed in the hope that it will be useful, 0026 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0027 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0028 % GNU General Public License for more details. 0029 % 0030 % You can obtain a copy of the GNU General Public License from 0031 % http://www.gnu.org/copyleft/gpl.html or by writing to 0032 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0033 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0034 0035 n=size(p,1); 0036 q=ones(n+1,3); 0037 q(1:n,1:2)=p; 0038 q(end,:)=q(1,:); % duplicate the final point 0039 cdist=q*l'; 0040 cside=cdist>0; % find which side each vertex lies 0041 cside(end)=cside(1); % just in case 0042 ec=find(cside(2:end)~=cside(1:n)); % sides of boundary that are crossed 0043 nc=numel(ec); 0044 if ~nc 0045 xc=[]; 0046 tc=[]; 0047 ec=[]; 0048 xy0=[]; 0049 else 0050 l2=l(1:2); 0051 l2m=l2*l2'; 0052 l3=l(3); 0053 a=[-l(2) l(1)]; 0054 xy0=-l3/(l2*l2')*l2; % point on line closest to origin 0055 tn=(q(:,1:2)-repmat(xy0,n+1,1))*a'/l2m; 0056 tc=(cdist(ec+1).*tn(ec)-cdist(ec).*tn(ec+1))./(cdist(ec+1)-cdist(ec)); 0057 [tc,ix]=sort(tc); 0058 ec=ec(ix); 0059 tm=tc(2:end)==tc(1:end-1); % check for duplicated points 0060 tm=[0;tm]|[tm;0]; 0061 tc(tm)=[]; % remove duplicate pairs 0062 ec(tm)=[]; 0063 nc=numel(ec); 0064 xc=repmat(xy0,nc,1)+tc*a; 0065 end 0066 if ~nargout && nc>0 0067 plot(q(:,1),q(:,2),'-r',xc(:,1),xc(:,2),'-xb',xc(1,1),xc(1,2),'ob') 0068 end