0001 function [xo,yo] = intsecpl(xv,yv,xl,yl,trace)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031 tol = 1e-14;
0032 marg = tol;
0033 is_ab = 0;
0034
0035
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
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
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
0060 lim = [min(xv)-marg max(xv)+marg min(yv)-marg max(yv)+marg];
0061
0062 d = sqrt((lim(2)-lim(1)).^2+(lim(4)-lim(3)).^2);
0063
0064
0065
0066 if is_ab
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
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
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
0100 if ~any(is)
0101 if trace
0102
0103 [xl,yl] = intsecpl(lim([1 2 2 1]),lim([3 3 4 4]),xl,yl);
0104 plot(xv,yv,xl,yl)
0105 end
0106 xo = []; yo = [];
0107 return
0108 end
0109
0110
0111
0112
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
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
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
0148
0149
0150 if length(ii)>1 & 0
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
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
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)
0172 end
0173