intsecpl

PURPOSE ^

INTSECPL - Intersection of a polygon and a line.

SYNOPSIS ^

function [xo,yo] = intsecpl(xv,yv,xl,yl,trace)

DESCRIPTION ^

INTSECPL - Intersection of a polygon and a line.

    [XI,YI] = INTSECPL(XV,YV,XL,YL) calculates
    intersections XI, YI of a polygon with vertices XV,
    YV and a line specified by pairs of end coordinates
    XL = [XL0 XL1], YL = [YL0 YL1]. Line is assumed to
    continue beyond the range of end points.
    INTSECPL(XV,YV,[A B]) uses another specification for
    a line: Y = A*X+B.

    If a line does not intersect polygon, returns empty
    XI, YI.
    For convex polygon maximum number of intersections is
    2, for non-convex polygons multiple intersections are 
    possible.

    INTSECPL(XV,YV,XL,YL)  by itself or
    [XI,YI] = INTSECPL(XV,YV,XL,YL,1) plots polygon,
    a line segment and intersection segment(s) 
    (part(s) of the same line inside the polygon).

  Copyright (c) 1995 by Kirill K. Pankratov,
       kirill@plume.mit.edu.
       06/25/95, 08/27/95, 09/27/95

  Calls ISCROSS, INTSECL programs.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function  [xo,yo] = intsecpl(xv,yv,xl,yl,trace)
0002 %INTSECPL - Intersection of a polygon and a line.
0003 %
0004 %    [XI,YI] = INTSECPL(XV,YV,XL,YL) calculates
0005 %    intersections XI, YI of a polygon with vertices XV,
0006 %    YV and a line specified by pairs of end coordinates
0007 %    XL = [XL0 XL1], YL = [YL0 YL1]. Line is assumed to
0008 %    continue beyond the range of end points.
0009 %    INTSECPL(XV,YV,[A B]) uses another specification for
0010 %    a line: Y = A*X+B.
0011 %
0012 %    If a line does not intersect polygon, returns empty
0013 %    XI, YI.
0014 %    For convex polygon maximum number of intersections is
0015 %    2, for non-convex polygons multiple intersections are
0016 %    possible.
0017 %
0018 %    INTSECPL(XV,YV,XL,YL)  by itself or
0019 %    [XI,YI] = INTSECPL(XV,YV,XL,YL,1) plots polygon,
0020 %    a line segment and intersection segment(s)
0021 %    (part(s) of the same line inside the polygon).
0022 %
0023 %  Copyright (c) 1995 by Kirill K. Pankratov,
0024 %       kirill@plume.mit.edu.
0025 %       06/25/95, 08/27/95, 09/27/95
0026 %
0027 %  Calls ISCROSS, INTSECL programs.
0028 
0029 
0030  % Defaults and parameters .................................
0031 tol = 1e-14;  % Tolerance
0032 marg = tol;   % Margins for polygon frame
0033 is_ab = 0;    % Default A*X+B  mode
0034 
0035  % Handle input ............................................
0036 if nargin==0, help intsecpl, return, end
0037 if nargin < 3
0038   error(' Not enough input arguments')
0039 end
0040 if nargin<5, trace = 0; end
0041 if nargin==4  % Check if 4-th arg is trace
0042   if max(size(yl))==1, trace = yl;  is_ab = 1; end
0043 end
0044 if nargin==3, is_ab = 1; end
0045 trace = trace | nargin<2;
0046 if length(xv)~=length(yv)
0047   error(' Vectors X, Y must have the same size')
0048 end
0049 
0050  % Auxillary ...........
0051 xv = [xv(:); xv(1)];
0052 yv = [yv(:); yv(1)];
0053 ii = find(abs(diff(xv))<tol & abs(diff(yv))<tol);
0054 xv(ii) = []; yv(ii) = [];
0055 nv = length(xv);
0056 ov = ones(nv-1,1);
0057 
0058 
0059   % Polygon frame
0060 lim = [min(xv)-marg max(xv)+marg min(yv)-marg max(yv)+marg];
0061  % Estimate for diameter
0062 d = sqrt((lim(2)-lim(1)).^2+(lim(4)-lim(3)).^2);
0063 
0064 
0065  % Form line segment depending on how line is specified
0066 if is_ab       % A*X+B  mode ...................
0067 
0068   xl = xl(:);
0069   if length(xl)<2
0070     error(' Line is specified by at least two parameters')
0071   end
0072 
0073   a = xl(1); b = xl(2);
0074   xl = [lim(1)-1 lim(2)+1];
0075   yl = a*xl+b;
0076 
0077 else          % [X1 X2],[Y1 Y2]  mode ..........
0078 
0079   x0 = (xl(1)+xl(2))/2;
0080   y0 = (yl(1)+yl(2))/2;
0081   dx = xl(2)-x0;
0082   dy = yl(2)-y0;
0083   dl = max(abs([lim(1:2)-x0 lim(3:4)-y0]));
0084   d = max(d,dl);
0085   dl = sqrt(dx.^2+dy.^2);
0086   dl = max(d,d/dl);
0087   dx = dx*dl; dy = dy*dl;
0088   xl = [x0-dx x0+dx];
0089   yl = [y0-dy y0+dy];
0090 
0091 end
0092 
0093 
0094  % Find intersecting segments ..............................
0095 is = iscross([xv(1:nv-1) xv(2:nv)]',[yv(1:nv-1) yv(2:nv)]',...
0096              xl(ov,:)',yl(ov,:)',0);
0097 
0098 
0099  % Quick exit if no intersections .........................
0100 if ~any(is)
0101   if trace
0102     % Intersection with polygon frame
0103     [xl,yl] = intsecpl(lim([1 2 2 1]),lim([3 3 4 4]),xl,yl);
0104     plot(xv,yv,xl,yl) % Plotting itself
0105   end
0106   xo = []; yo = []; 
0107   return
0108 end
0109 
0110 
0111  % For segments touching the line (is==.5) find whether pairs of
0112  % successive segments are on the same side
0113 ii = find(is==.5)';
0114 if ~isempty(ii)
0115   xo = [ii-1 ii+1];
0116   xo = xo+(nv-1)*(xo<1);
0117   yo = iscross([xv(xo(:,1)) xv(xo(:,2))]',[yv(xo(:,1)) yv(xo(:,2))]',...
0118        xl,yl,tol);
0119   ii = ii(find(yo==1));
0120   is(ii) = zeros(size(ii));
0121 end
0122 
0123 
0124  % Calculate intersection coordinates ......................
0125 ii = find(is);
0126 oi = ones(size(ii));
0127 [xo,yo] = intsecl([xv(ii) xv(ii+1)]',[yv(ii) yv(ii+1)]',...
0128                   xl(oi,:)',yl(oi,:)');
0129 
0130 dx = find(~finite(xo));
0131 xo(dx) = []; yo(dx) = [];
0132 ii(dx) = []; oi(dx) = [];
0133 
0134 
0135  % Sort intersections along the line ..........
0136 xo = xo(:); yo = yo(:);
0137 if any(diff(xo))
0138   [xo,ii] = sort(xo);
0139   yo = yo(ii);
0140 
0141 else
0142 
0143   [yo,ii] = sort(yo);
0144   xo = xo(ii);
0145 end
0146 
0147  % Exclude repeated points (degenerate cases)
0148  % ///////// It seems that this is not needed, figure this
0149  % later ///////////////////
0150 if length(ii)>1 & 0  % Do not execute
0151   xx = [xo yo];
0152   yy = diff(xx)';
0153   ii = [1 find(any(abs(yy)>tol))+1];
0154   xo = xx(ii,1); yo = xx(ii,2);
0155   oi = ones(size(xo));
0156 end
0157 
0158 
0159  % Plotting ................................................
0160 if trace
0161   oi(3:2:length(oi)) = oi(3:2:length(oi))+1;
0162   oi = cumsum(oi);
0163   len = max(oi);
0164   xp = nan*ones(len,1); yp = xp;
0165   xp(oi) = xo;
0166   yp(oi) = yo;
0167 
0168   % Intersection with polygon frame
0169   [xl,yl] = intsecpl(lim([1 2 2 1]),lim([3 3 4 4]),xl,yl);
0170 
0171   plot(xv,yv,xl,yl,xp,yp) % Plotting itself
0172 end
0173

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003