function [x,y] = curveintersect(varargin)
%
% Curve Intersections.
% [X,Y]=CURVEINTERSECT(H1,H2) or [X,Y]=CURVEINTERSECT([H1 H2]) finds the
% intersection points of the two curves on the X-Y plane identified
% by the line or lineseries object handles H1 and H2.
%
% [X,Y]=CURVEINTERSECT(X1,Y1,X2,Y2) finds the intersection points of the
% two curves described by the vector data pairs (X1,Y1) and (X2,Y2).
%
% X and Y are empty if no intersection exists.
%
% Example
% -------
% x1=rand(10,1); y1=rand(10,1); x2=rand(10,1); y2=rand(10,1);
% [x,y]=curveintersect(x1,y1,x2,y2);
% plot(x1,y1,'k',x2,y2,'b',x,y,'ro')
%
% Original Version (-> curveintersect_local)
% ---------------------------------------
% D.C. Hanselman, University of Maine, Orono, ME 04469
% Mastering MATLAB 7
% 2005-01-06
%
% Improved Version (-> this function)
% -----------------------------------
% S. Hlz, TU Berlin, Germany
% v 1.0: October 2005
% v 1.1: April 2006 Fixed some minor bugs in function 'mminvinterp'
x=[]; y=[];
[x1,y1,x2,y2]=local_parseinputs(varargin{:});
ind_x1=sign(diff(x1)); ind_x2=sign(diff(x2));
ind1=1;
while ind1<length(x1)
ind_max = ind1+min(find(ind_x1(ind1:end)~=ind_x1(ind1)))-1;
if isempty(ind_max) | ind_max==ind1; ind_max=length(x1); end
ind1=ind1:ind_max;
ind2=1;
while ind2<length(x2)
ind_max = ind2+min(find(ind_x2(ind2:end)~=ind_x2(ind2)))-1;
if isempty(ind_max) | ind_max==ind2; ind_max=length(x2); end
ind2=ind2:ind_max;
% Fallunterscheidung
if ind_x1(ind1(1))==0 & ind_x2(ind2(1))~=0
x_loc=x1(ind1(1));
y_loc=interp1(x2(ind2),y2(ind2),x_loc);
if ~(y_loc>=min(y1(ind1)) && y_loc<=max(y1(ind1))); y_loc=[]; x_loc=[]; end
elseif ind_x2(ind2(1))==0 & ind_x1(ind1(1))~=0
x_loc=x2(ind2(1));
y_loc=interp1(x1(ind1),y1(ind1),x_loc);
if ~(y_loc>=min(y2(ind2)) && y_loc<=max(y2(ind2))); y_loc=[]; x_loc=[]; end
elseif ind_x2(ind2(1))~=0 & ind_x1(ind1(1))~=0
[x_loc,y_loc]=curveintersect_local(x1(ind1),y1(ind1),x2(ind2),y2(ind2));
elseif ind_x2(ind2(1))==0 & ind_x1(ind1(1))==0
[x_loc,y_loc]=deal([]);
end
x=[x; x_loc(:)];
y=[y; y_loc(:)];
ind2=ind2(end);
end
ind1=ind1(end);
end
% ----------------------------------------------
function [x,y]=curveintersect_local(x1,y1,x2,y2)
if ~isequal(x1,x2)
xx=unique([x1 x2]); % get unique data points
xx=xx(xx>=max(min(x1),min(x2)) & xx<=min(max(x1),max(x2)));
if numel(xx)<2
x=[];
y=[];
return
end
yy=interp1(x1,y1,xx)-interp1(x2,y2,xx);
else
xx=x1;
yy=y1-y2;
end
x=mminvinterp(xx,yy,0); % find zero crossings of difference
if ~isempty(x)
y=interp1(x1,y1,x);
else
x=[];
y=[];
end
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
function [xo,yo]=mminvinterp(x,y,yo)
%MMINVINTERP 1-D Inverse Interpolation. From the text "Mastering MATLAB 7"
% [Xo, Yo]=MMINVINTERP(X,Y,Yo) linearly interpolates the vector Y to find
% the scalar value Yo and returns all corresponding values Xo interpolated
% from the X vector. Xo is empty if no crossings are found. For
% convenience, the output Yo is simply the scalar input Yo replicated so
% that size(Xo)=size(Yo).
% If Y maps uniquely into X, use INTERP1(Y,X,Yo) instead.
%
% See also INTERP1.
if nargin~=3
error('Three Input Arguments Required.')
end
n = numel(y);
if ~isequal(n,numel(x))
error('X and Y Must have the Same Number of Elements.')
end
if ~isscalar(yo)
error('Yo Must be a Scalar.')
end
x=x(:); % stretch input vectors into column vectors
y=y(:);
if yo<min(y) || yo>max(y) % quick exit if no values exist
xo = [];
yo = [];
else % find the desired points
below = y<yo; % True where below yo
above = y>yo; % True where above yo
on = y==yo; % True where on yo
kth = (below(1:n-1)&above(2:n)) | (above(1:n-1)&below(2:n)); % point k
kp1 = [false; kth]; % point k+1
xo = []; % distance between x(k+1) and x(k)
if any(kth);
alpha = (yo - y(kth))./(y(kp1)-y(kth));
xo = alpha.*(x(kp1)-x(kth)) + x(kth);
end
xo = sort([xo; x(on)]); % add points, which are directly on line
yo = repmat(yo,size(xo)); % duplicate yo to match xo points found
end
%--------------------------------------------------------------------------
function [x1,y1,x2,y2]=local_parseinputs(varargin)
if nargin==1 % [X,Y]=CURVEINTERSECT([H1 H2])
arg=varargin{1};
if numel(arg)==2 && ...
all(ishandle(arg)) && all(strcmp(get(arg,'type'),'line'))
data=get(arg,{'XData','YData'});
[x1,x2,y1,y2]=deal(data{:});
else
error('Input Must Contain Two Handles to Line Objects.')
end
elseif nargin==2 % [X,Y]=CURVEINTERSECT(H1,H2)
arg1=varargin{1};
arg2=varargin{2};
if numel(arg1)==1 && ishandle(arg1) && strcmp(get(arg1,'type'),'line')...
&& numel(arg2)==1 && ishandle(arg2) && strcmp(get(arg2,'type'),'line')
data=get([arg1;arg2],{'XData','YData'});
[x1,x2,y1,y2]=deal(data{:});
else
error('Input Must Contain Two Handles to Line Objects.')
end
elseif nargin==4
[x1,y1,x2,y2]=deal(varargin{:});
if ~isequal(numel(x1),numel(y1))
error('X1 and Y1 Must Contain the Same Number of Elements.')
elseif ~isequal(numel(x2),numel(y2))
error('X2 and Y2 Must Contain the Same Number of Elements.')
end
x1=reshape(x1,1,[]); % make data into rows
x2=reshape(x2,1,[]);
y1=reshape(y1,1,[]);
y2=reshape(y2,1,[]);
else
error('Incorrect Number of Input Arguments.')
end
if min(x1)>max(x2) | min(x2)>max(x1) | min(y2)>max(y1) | min(y1)>max(y2) % Polygons can not have intersections
x1=[]; y1=[]; x2=[]; y2=[]; return
end
if numel(x1)<2 || numel(x2)<2 || numel(y1)<2 || numel(y2)<2
error('At Least Two Data Points are Required for Each Curve.')
end
