v_polygonxline

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