Hi, does anybody know where can I find this function in FORTRAN?

Thanks Douglas!

I was dealing with interpolation problem when have replicate points, this function was a great alternative! Great work!

Thank you so much! This is perfect.

Perfect for my application. Thanks a lot!!!

Thanks Doug, your program is impressive and has saved me a lot of time. On the other hand, I encounter a difficulty in applying it to the following case:
% Segment [AB]
x1=[0.119316419331154;0.167903096090215];
y1=[0.156890140334954;0.181466113513819];
%Segment [CD]
x2=[0;0.133106852905348];
y2=[0.220693954988999;0.163865577805703];
intersections(x1,y1,x2,y2,true) %Empty matrix: 0-by-1
The point D(x2(2),y2(2)) has been determined by the same function and should be the intersection point of the two segments.
when I calculate the gap between the point D and the segment [AB] I find 2.2999e-16
Could we consider a certain tolerance in the determination of the point of intersection?
Thanks!

An absolute lifesaver, thank you a bunch.

Thanks, Douglas. I ended up using the contourcs function (https://www.mathworks.com/matlabcentral/fileexchange/28447-contourcs-to-obtain-contourc-output-as-a-struct-array), to do what you achieved with the splitcontours function. But, overall, this is as awesome utility. Thanks!

Marta, here is the function I used to make the demo plot. Perhaps it will help. -Doug

function intersections_demo
np = 201;
p = peaks(np);
xy = linspace(-3,3,np);
[dx,dy] = gradient(p,xy,xy);
cx = contourc(xy,xy,dx,[0 0]);
cy = contourc(xy,xy,dy,[0 0]);
[~,vx] = splitcontours(cx);
[~,vy] = splitcontours(cy);
[xx,yy] = intersections(vx(:,1),vx(:,2),vy(:,1),vy(:,2));

contour(xy,xy,p,20,'Color',[0.85 0.85 0.85])
axis square
orient portrait
set(gcf,'paperpos',[0 0 5 5])
set(gca,'ytick',-3:3)
set(gca,'xtick',-3:3)
hold on
hh = plot(vx(:,1),vx(:,2),vy(:,1),vy(:,2),xx,yy,'ok');
legend(hh,{'curve 1','curve 2','intersections'})
title('Local Extrema and Saddle Points of peaks Function')
hold off

function [s,v] = splitcontours(con)
%splitcontours: Convert contour matrix to structure array.

% written by Douglas M. Schwarz.

len_con = size(con,2);
val_flag = false(1,len_con);
ptr = 1;
while ptr < len_con
val_flag(ptr) = true;
ptr = ptr + con(2,ptr) + 1;
end
num_cons = sum(val_flag);
val_locs = [find(val_flag),len_con+1];
nvals = diff(val_locs) - 1;
s = repmat(struct('value',[],'x',[],'y',[],'n',[]),1,num_cons);
for i = 1:num_cons
s(i).value = con(1,val_locs(i));
s(i).n = nvals(i);
i1 = val_locs(i) + 1;
i2 = val_locs(i+1) - 1;
s(i).x = con(1,i1:i2).';
s(i).y = con(2,i1:i2).';
end

if nargout > 1
v = con.';
v(val_flag,:) = NaN;
v(1,:) = [];
end

Armindo, n = 90000 points is a lot. Intermediate matrices of size n^2 have to be formed. As is, I'm afraid my function is not going to help you. It is possible to partition your points. For example, if you make four quadrants with 45000 points per curve in each, now you only have to form four 45000^2 matrices. Continuing like that might lead to something that can be solved in a reasonable amount of time. -Doug

Hello,
i need your Help,i am still a beginner in Matlab i had a task to determine intersections of polygons and curves, i tried to write a code with the Polyxpoly function, but i have a problem, i don't know how to define and minimize the total time and the time increase by calling the function

i need your help

Hi Douglas,

I'm trying to calculate stationary points in a numerically calculated array; this is done using the 2D gradient function: [Fx,Fy]=gradient(fxy,dx,dy), following which I plot the zeros of Fx and Fy using [Cx,~] = contour(x,y,Fx,[0 0]), and corresponding Cy definition, and was hoping to use the intersections function to calculate the intersection points of these zeros. Is there an easy way to do this using the intersections function in conjunction with the contour definition. The graph you have of the saddle points of the peaks function seems to do this, but are there intermediate steps I'm missing? This may also be complicated because my function is numerically calculated and not derived from an algebraic expression.

Thanks!

I can't remember how many times I've recommended people to use this code.

How to get intersections from two curves with 90000 points? Is there any any to use this function with so many points?

Used this for a complicated aerodynamics problem. Worked perfectly. Thank you very much!

Nadish,
You certainly have the right idea, but without your code I can't tell if you're doing it right. Your code should look something like this:
[x0,y0] = intersections(x,f(x),x([1 end]),[0 0]);
In other words, make sure you express the zero line with two points.
Doug

Hello Doug,

I am trying to find the roots of a function with your code. I was trying to use fsolve etc, but I did not get proper convergence and fsolve did not capture multiple roots. I was wondering if I can make the second curve as y = 0 and the first one as y = f(x), the roots must be the intercepts of the equation. I tried it, but I am getting an empty set. I think I might be doing something wrong with either in using intercepts function or in passing the values to it. Can you please help me with this?

Sincerely,

Nadish

has anyone mexed this function?

Trisha, it does work with R2018a. What error are you getting? -Doug

Hi the intersections function doesn't work on Matlab2018a is there another alternative?

Hi Ying Tang,
There are a couple of issues. 1) Because of floating point arithmetic, an intersection exactly at one of the segment end points may or may not be included. 2) An "intersection" at a tangent point is not really an intersection in that the curves don't cross. It depends on your interpretation. There is no bug in the code; some things are just unavoidable. Perhaps some other tool would be more suitable for your purpose.
Doug

Hi Doug,
Thank you for the code. but I found out that there is a issue with the tangent line of a curve. In this case it seems that the code is not able to find the point.
for instance I want to calculate the intersection of unit circle and line x = 1:
x = linspace(-10,10,5000);
y1 = sqrt(1 - x.^2);

x_filted = x(imag(y1)==0);
y_filted = y1(imag(y1)==0);

x_all = [x_filted,fliplr(x_filted)];
y_all = [y_filted,-fliplr(y_filted)];

figure,hold on,axis equal;
plot(x_all,y_all);
plot([0.5,0.5],[-5,5]);

[x0,y0] = intersections([0.5,0.5],[-5,5],x_all,y_all,0)

the output does not give a solution.

is there anyway to solve this issue?

Doug,
Thank you! This works exactly as I'd hoped.

Doug,

Thank you very much for this function! It worked perfectly! I wish I could attached some jpegs. I spent several hours trying to locate intersection of two curves when trying to find yield strength with no success. I found your function and it worked seamlessly - intersections found! Cheers!

Amazing! Thanks, really helpful!

Hi, Doug. Greetings from a former colleague. Thanks for posting this excellent function - it has been of great use.

I found a case where the two curves intersect at a common endpoint and the function returns empty (x0, y0) because of the finite precision arithmetic when calculating T (e.g. coordinate 1 versus coordinate 1 + 1e-16). I believe that adding a decimalPrecision input argument and inserting T(:,k) = round(T(:,k)*10^decimalPrecision)/10^decimalPrecision; after line 266 (T(:,k) = AA(:,:,k)\B(:,k);) fixes the problem of the intersection point failing the in-range condition (i.e. in_range = (T(1,:) >= 0 & T(2,:) >= 0 & T(1,:) <= 1 & T(2,:) <= 1).';) due to the finite arithmetic precision. I haven't thought about this long enough to say that this is definitely the best solution, but I thought I'd let you know how I am dealing with it for my data set. - Efrain

Is there a 3d version of this?

Carl, at this point I can't change the output arguments as doing so would break backward compatibility. I generally prefer to have separate arguments for x- and y-values which is why I did it that way. I find that I usually just have to separate them at some point anyway. You can, of course, change your copy of the function any way you desire. -Doug

I strongly recommend changing the output to a N-by-2 array. There's no reason to require two separate objects (with two separate names) for such closely-linked data. (Plus it's a royal pain to have to type " [ name ,name ]" ) It is far easier to index into an array than to have to index into two separate objects to retrieve a point's coordinates.

Hi Doug:
Thanks for the response; I'll check out your new version. My change for closed polygons does not change the meaning of NaNs. It just means that if you have a sequence of points in your curve like:

NaN A B C D A NaN

that no intersection is ever reported for the A-B segment with the D-A segment. With your old code this would happen very often, reporting point A as an intersection. Your self-intersection code (as I read it) already eliminates checks for consecutive segments, i.e., it will never check A-B against B-C due to your test "j <= i + 1" forcing the second segment of a test to be at least two greater than the first. My addition adds an extra test "polyidx(j) == i" to catch these cases. Here's a minimal example you can try to illustrate the problem. Your code will report that the unit square is self-intersecting.

polygonx = [0 1 1 0 0];
polygony = [0 0 1 1 0];
[x0,y0,iout,jout] = intersections(polygonx, polygony, 1)

If you're happy for me to do so, I'll submit a renamed version, keeping your name with your code.

P.S.: I just did a quick test with my 11,000-point curve and my recursive version is still over ten times faster than your new version.

James, I just submitted an update of the function. It should work much better when the number of line segments is large. My algorithm is different from your suggestion so I would be interested in a comparison. However, I did not include your code to split closed polygons with NaNs as my original function used NaNs to split a curve into separate sections. Changing the meaning of NaNs would break existing code. You might want to change the name of the function and submit your version to the FEX. Thanks for the incentive to post this -- I had worked on it a couple of months ago to address the concerns of Far Bod, but hadn't gotten around to finishing up the details. -Doug

Hi Doug:

Two issues I found and fixed:
1) My curves are sets of polygons which are properly closed, i.e., NaNs separate the polygons and, for each, the last point is identical to the first. This is common for GIS work. If I test such a curve for self-intersection it often finds that the first segment of a polygon intersects with the last. I've come up with a few lines to eliminate testing for these cases.

2) intersections.m uses n^2 space so it won't work once your curves go far past 10,000 points, which is not uncommon for GIS work (see Far Bod's enquiry below). What you can do is recursively divide the problem up into four smaller sub-problems; it's essentially a quad-tree approach. This *massively* speeds up the code for large problems, as well as massively reducing the memory requirements. Testing my 11,000-point map for self-intersections went from 3.0 seconds to 0.02 seconds. (Poor old polyxpoly takes 22 seconds).

I'm pasting your code with my updates below. Don't know if that will work or whether the size or formatting will cause issues. Feel free to make what use you like of my additions.

function [x0,y0,iout,jout] = intersections(x1,y1,x2,y2,robust)
%INTERSECTIONS Intersections of curves.
% Computes the (x,y) locations where two curves intersect. The curves
% can be broken with NaNs or have vertical segments.
%
% Example:
% [X0,Y0] = intersections(X1,Y1,X2,Y2,ROBUST);
%
% where X1 and Y1 are equal-length vectors of at least two points and
% represent curve 1. Similarly, X2 and Y2 represent curve 2.
% X0 and Y0 are column vectors containing the points at which the two
% curves intersect.
%
% ROBUST (optional) set to 1 or true means to use a slight variation of the
% algorithm that might return duplicates of some intersection points, and
% then remove those duplicates. The default is true, but since the
% algorithm is slightly slower you can set it to false if you know that
% your curves don't intersect at any segment boundaries. Also, the robust
% version properly handles parallel and overlapping segments.
%
% The algorithm can return two additional vectors that indicate which
% segment pairs contain intersections and where they are:
%
% [X0,Y0,I,J] = intersections(X1,Y1,X2,Y2,ROBUST);
%
% For each element of the vector I, I(k) = (segment number of (X1,Y1)) +
% (how far along this segment the intersection is). For example, if I(k) =
% 45.25 then the intersection lies a quarter of the way between the line
% segment connecting (X1(45),Y1(45)) and (X1(46),Y1(46)). Similarly for
% the vector J and the segments in (X2,Y2).
%
% You can also get intersections of a curve with itself. Simply pass in
% only one curve, i.e.,
%
% [X0,Y0] = intersections(X1,Y1,ROBUST);
%
% where, as before, ROBUST is optional.

% Version: 1.12, 27 January 2010
% Author: Douglas M. Schwarz
% Email: dmschwarz=ieee*org, dmschwarz=urgrad*rochester*edu
% Real_email = regexprep(Email,{'=','*'},{'@','.'})

% Theory of operation:
%
% Given two line segments, L1 and L2,
%
% L1 endpoints: (x1(1),y1(1)) and (x1(2),y1(2))
% L2 endpoints: (x2(1),y2(1)) and (x2(2),y2(2))
%
% we can write four equations with four unknowns and then solve them. The
% four unknowns are t1, t2, x0 and y0, where (x0,y0) is the intersection of
% L1 and L2, t1 is the distance from the starting point of L1 to the
% intersection relative to the length of L1 and t2 is the distance from the
% starting point of L2 to the intersection relative to the length of L2.
%
% So, the four equations are
%
% (x1(2) - x1(1))*t1 = x0 - x1(1)
% (x2(2) - x2(1))*t2 = x0 - x2(1)
% (y1(2) - y1(1))*t1 = y0 - y1(1)
% (y2(2) - y2(1))*t2 = y0 - y2(1)
%
% Rearranging and writing in matrix form,
%
% [x1(2)-x1(1) 0 -1 0; [t1; [-x1(1);
% 0 x2(2)-x2(1) -1 0; * t2; = -x2(1);
% y1(2)-y1(1) 0 0 -1; x0; -y1(1);
% 0 y2(2)-y2(1) 0 -1] y0] -y2(1)]
%
% Let's call that A*T = B. We can solve for T with T = A\B.
%
% Once we have our solution we just have to look at t1 and t2 to determine
% whether L1 and L2 intersect. If 0 <= t1 < 1 and 0 <= t2 < 1 then the two
% line segments cross and we can include (x0,y0) in the output.
%
% In principle, we have to perform this computation on every pair of line
% segments in the input data. This can be quite a large number of pairs so
% we will reduce it by doing a simple preliminary check to eliminate line
% segment pairs that could not possibly cross. The check is to look at the
% smallest enclosing rectangles (with sides parallel to the axes) for each
% line segment pair and see if they overlap. If they do then we have to
% compute t1 and t2 (via the A\B computation) to see if the line segments
% cross, but if they don't then the line segments cannot cross. In a
% typical application, this technique will eliminate most of the potential
% line segment pairs.

% Extra function by James Ashton, Canberra, Australia
% Keeping this as a separate function means Matlab can automatically
% free up all the temporary space once we're done.
function [i2,j2] = bBoxOverlaps()

% recursive function for searching in sub-rectangles
function [i3,j3] = quadtree(bb1, in1, bb2, in2, bb, outerprobsize)
nn1 = numel(in1);
nn2 = numel(in2);
probsize = nn1 * nn2;
% If there are many intersections to consider, divide the
% rectangle under consideration into four equal smaller
% rectangles and solve for each separately. Don't divide if
% doing so isn't reducing the problem (many long overlapping lines).
if (probsize > 3000 & probsize < outerprobsize)
midx = (bb(1) + bb(3)) * 0.5;
midy = (bb(2) + bb(4)) * 0.5;
bb = [bb midx midy];

ix1 = find(bb1(:,1) <= midx & bb1(:,2) <= midy);
ix2 = find(bb2(:,1) <= midx & bb2(:,2) <= midy);
[i3, j3] = quadtree(bb1(ix1,:), in1(ix1), bb2(ix2,:), in2(ix2), bb([1 2 5 6]), probsize);

ix1 = find(bb1(:,3) >= midx & bb1(:,2) <= midy);
ix2 = find(bb2(:,3) >= midx & bb2(:,2) <= midy);
[i4, j4] = quadtree(bb1(ix1,:), in1(ix1), bb2(ix2,:), in2(ix2), bb([5 2 3 6]), probsize);
i3 = [i3;i4];
j3 = [j3;j4];

ix1 = find(bb1(:,1) <= midx & bb1(:,4) >= midy);
ix2 = find(bb2(:,1) <= midx & bb2(:,4) >= midy);
[i4, j4] = quadtree(bb1(ix1,:), in1(ix1), bb2(ix2,:), in2(ix2), bb([1 6 5 4]), probsize);
i3 = [i3;i4];
j3 = [j3;j4];

ix1 = find(bb1(:,3) >= midx & bb1(:,4) >= midy);
ix2 = find(bb2(:,3) >= midx & bb2(:,4) >= midy);
[i4, j4] = quadtree(bb1(ix1,:), in1(ix1), bb2(ix2,:), in2(ix2), bb([5 6 3 4]), probsize);
i3 = [i3;i4];
j3 = [j3;j4];
else
[i3,j3] = find(...
repmat(bb1(:,1),1,nn2) <= repmat(bb2(:,3).',nn1,1) & ...
repmat(bb1(:,3),1,nn2) >= repmat(bb2(:,1).',nn1,1) & ...
repmat(bb1(:,2),1,nn2) <= repmat(bb2(:,4).',nn1,1) & ...
repmat(bb1(:,4),1,nn2) >= repmat(bb2(:,2).',nn1,1) ...
);
i3 = reshape(in1(i3),[],1);
j3 = reshape(in2(j3),[],1);
end
end

% find the bounding boxes of each segment
% each bbox row will have: [minx, miny, maxx, maxy]
bbox1 = horzcat(...
min(x1(1:end-1), x1(2:end)), ...
min(y1(1:end-1), y1(2:end)), ...
max(x1(1:end-1), x1(2:end)), ...
max(y1(1:end-1), y1(2:end)) ...
);
indx1 = 1:n1; % keep a list of the segment indexes
del = find(isnan(dxy1(:,1)));
bbox1(del,:) = [];
indx1(del) = [];

bbox2 = horzcat(...
min(x2(1:end-1), x2(2:end)), ...
min(y2(1:end-1), y2(2:end)), ...
max(x2(1:end-1), x2(2:end)), ...
max(y2(1:end-1), y2(2:end)) ...
);
indx2 = 1:n2;
del = find(isnan(dxy2(:,1)));
bbox2(del,:) = [];
indx2(del) = [];

% find the overall bounding box where potential overlaps might
% occur.
minx = max(min(bbox1(:,1)),min(bbox2(:,1)));
miny = max(min(bbox1(:,2)),min(bbox2(:,2)));
maxx = min(max(bbox1(:,3)),max(bbox2(:,3)));
maxy = min(max(bbox1(:,4)),max(bbox2(:,4)));

[i2,j2] = quadtree(bbox1, indx1, bbox2, indx2, [minx miny maxx maxy], inf);
end

% Input checks.
narginchk(2,5);
% Adjustments when fewer than five arguments are supplied.
switch nargin
case 2
robust = true;
x2 = x1;
y2 = y1;
self_intersect = true;
case 3
robust = x2;
x2 = x1;
y2 = y1;
self_intersect = true;
case 4
robust = true;
self_intersect = false;
case 5
self_intersect = false;
end

% x1 and y1 must be vectors with same number of points (at least 2).
if sum(size(x1) > 1) ~= 1 || sum(size(y1) > 1) ~= 1 || ...
length(x1) ~= length(y1)
error('X1 and Y1 must be equal-length vectors of at least 2 points.')
end
% x2 and y2 must be vectors with same number of points (at least 2).
if sum(size(x2) > 1) ~= 1 || sum(size(y2) > 1) ~= 1 || ...
length(x2) ~= length(y2)
error('X2 and Y2 must be equal-length vectors of at least 2 points.')
end

% Force all inputs to be column vectors.
x1 = x1(:);
y1 = y1(:);
x2 = x2(:);
y2 = y2(:);

% Compute number of line segments in each curve and some differences we'll
% need later.
n1 = length(x1) - 1;
n2 = length(x2) - 1;
xy1 = [x1 y1];
xy2 = [x2 y2];
dxy1 = diff(xy1);
dxy2 = diff(xy2);

% Determine the combinations of i and j where the rectangle enclosing the
% i'th line segment of curve 1 overlaps with the rectangle enclosing the
% j'th line segment of curve 2.

% [i,j] = find(...
% repmat(min(x1(1:end-1),x1(2:end)),1,n2) <= ...
% repmat(max(x2(1:end-1),x2(2:end)).',n1,1) & ...
% repmat(max(x1(1:end-1),x1(2:end)),1,n2) >= ...
% repmat(min(x2(1:end-1),x2(2:end)).',n1,1) & ...
% repmat(min(y1(1:end-1),y1(2:end)),1,n2) <= ...
% repmat(max(y2(1:end-1),y2(2:end)).',n1,1) & ...
% repmat(max(y1(1:end-1),y1(2:end)),1,n2) >= ...
% repmat(min(y2(1:end-1),y2(2:end)).',n1,1));
% % Force i and j to be column vectors, even when their length is zero, i.e.,
% % we want them to be 0-by-1 instead of 0-by-0.
% i = reshape(i,[],1);
% j = reshape(j,[],1);

[i, j] = bBoxOverlaps();

% Find segments pairs which have at least one vertex = NaN and remove them.
% This line is a fast way of finding such segment pairs. We take
% advantage of the fact that NaNs propagate through calculations, in
% particular subtraction (in the calculation of dxy1 and dxy2, which we
% need anyway) and addition.
% At the same time we can remove redundant combinations of i and j in the
% case of finding intersections of a line with itself.
if self_intersect
% lastpoly code by James Ashton, Canberra, Australia
% for closed polygons, don't test the first and last segments
lastpoly = isnan(dxy1(:,1)); % find polygon separators
firstpoly = find(~lastpoly & [true;lastpoly(1:end-1)]); % first segment of each poly
lastpoly = find(~lastpoly & [lastpoly(2:end);true]); % last segment of each poly
% not a polygon (not closed) if first and last points aren't identical
remove = x1(firstpoly) ~= x1(lastpoly+1) | y1(firstpoly) ~= y1(lastpoly+1);
firstpoly(remove) = [];
lastpoly(remove) = [];
polyidx = (1:n1)';
polyidx(lastpoly) = firstpoly; % create a lookup array to find the first segment given the last segment
remove = isnan(sum(dxy1(i,:) + dxy2(j,:),2)) | j <= i + 1 | polyidx(j) == i;
else
remove = isnan(sum(dxy1(i,:) + dxy2(j,:),2));
end
i(remove) = [];
j(remove) = [];

% Initialize matrices. We'll put the T's and B's in matrices and use them
% one column at a time. AA is a 3-D extension of A where we'll use one
% plane at a time.
n = length(i);
T = zeros(4,n);
AA = zeros(4,4,n);
AA([1 2],3,:) = -1;
AA([3 4],4,:) = -1;
AA([1 3],1,:) = dxy1(i,:).';
AA([2 4],2,:) = dxy2(j,:).';
B = -[x1(i) x2(j) y1(i) y2(j)].';

% Loop through possibilities. Trap singularity warning and then use
% lastwarn to see if that plane of AA is near singular. Process any such
% segment pairs to determine if they are colinear (overlap) or merely
% parallel. That test consists of checking to see if one of the endpoints
% of the curve 2 segment lies on the curve 1 segment. This is done by
% checking the cross product
%
% (x1(2),y1(2)) - (x1(1),y1(1)) x (x2(2),y2(2)) - (x1(1),y1(1)).
%
% If this is close to zero then the segments overlap.

% If the robust option is false then we assume no two segment pairs are
% parallel and just go ahead and do the computation. If A is ever singular
% a warning will appear. This is faster and obviously you should use it
% only when you know you will never have overlapping or parallel segment
% pairs.

if robust
overlap = false(n,1);
warning_state = warning('off','MATLAB:singularMatrix');
% Use try-catch to guarantee original warning state is restored.
try
lastwarn('')
for k = 1:n
T(:,k) = AA(:,:,k)\B(:,k);
[~,last_warn] = lastwarn;
lastwarn('')
if strcmp(last_warn,'MATLAB:singularMatrix')
% Force in_range(k) to be false.
T(1,k) = NaN;
% Determine if these segments overlap or are just parallel.
overlap(k) = rcond([dxy1(i(k),:);xy2(j(k),:) - xy1(i(k),:)]) < eps;
end
end
warning(warning_state)
catch err
warning(warning_state)
rethrow(err)
end
% Find where t1 and t2 are between 0 and 1 and return the corresponding
% x0 and y0 values.
in_range = (T(1,:) >= 0 & T(2,:) >= 0 & T(1,:) <= 1 & T(2,:) <= 1).';
% For overlapping segment pairs the algorithm will return an
% intersection point that is at the center of the overlapping region.
if any(overlap)
ia = i(overlap);
ja = j(overlap);
% set x0 and y0 to middle of overlapping region.
T(3,overlap) = (max(min(x1(ia),x1(ia+1)),min(x2(ja),x2(ja+1))) + ...
min(max(x1(ia),x1(ia+1)),max(x2(ja),x2(ja+1)))).'/2;
T(4,overlap) = (max(min(y1(ia),y1(ia+1)),min(y2(ja),y2(ja+1))) + ...
min(max(y1(ia),y1(ia+1)),max(y2(ja),y2(ja+1)))).'/2;
selected = in_range | overlap;
else
selected = in_range;
end
xy0 = T(3:4,selected).';

% Remove duplicate intersection points.
[xy0,index] = unique(xy0,'rows');
x0 = xy0(:,1);
y0 = xy0(:,2);

% Compute how far along each line segment the intersections are.
if nargout > 2
sel_index = find(selected);
sel = sel_index(index);
iout = i(sel) + T(1,sel).';
jout = j(sel) + T(2,sel).';
end
else % non-robust option
for k = 1:n
[L,U] = lu(AA(:,:,k));
T(:,k) = U\(L\B(:,k));
end

% Find where t1 and t2 are between 0 and 1 and return the corresponding
% x0 and y0 values.
in_range = (T(1,:) >= 0 & T(2,:) >= 0 & T(1,:) < 1 & T(2,:) < 1).';
x0 = T(3,in_range).';
y0 = T(4,in_range).';

% Compute how far along each line segment the intersections are.
if nargout > 2
iout = i(in_range) + T(1,in_range).';
jout = j(in_range) + T(2,in_range).';
end
end

% Plot the results (useful for debugging).
% plot(x1,y1,x2,y2,x0,y0,'ok');
end

Hello Doug,

Thanks for amazing code, do you have any solution for finding the intersection for very large size of arrays? my 2 curves consist of 90,000 elements and repmat gives a memory error ... I cannot reduce my time increment since it affects the accuracy of method

Sara,
My function does not assume that the curves you specify are closed. To see what I mean, plot your two curves with

plot(A(:,1),A(:,2),B,C)

You'll see that the line only intersects the curve at one point. To get an intersection at the second point you will have to close the curve A by replicating the first point as the last point.

>> A2 = [A;A(1,:)];
>> [x0,y0] = intersections(A2(:,1),A2(:,2),B,C,1)
x0 =
5.9817
5.9817
y0 =
2.9941
5.9030

Now you have your two desired intersections.
-Doug

Hi Douglas,
When looking for the intersection between shape A, with vertical line that passes through BC, I get only one intersection point, whereas if you plot the data, it should be two. Can you please explain this?

A =
6.1983 2.9970
6.5888 3.1271
7.0661 3.4502
7.4663 3.7828
7.7652 4.2264
8.0063 4.6748
8.1076 5.0701
8.1172 5.5812
8.0786 5.9332
7.6736 6.1501
7.1626 6.0585
6.5310 5.9235
6.0392 5.8946
5.4124 5.9862
4.8628 6.1501
4.3903 6.1694
4.0046 5.9910
3.7876 5.5282
3.9757 5.0653
4.2987 4.3180
4.6892 3.6430
5.2051 3.2043
5.8367 2.9921
B =

5.9817
5.9817

C =

6.6694
2.4921

[x0,y0]=intersections(A(:,1),A(:,2),B,C,1)

Hi Douglas,

when sections of both curves overlap, the output is the midpoint of the overlapping region. What's missing is the fractional index of that point (third and fourth function output). Somewhere between lines 227 and 236 T(1,overlap) and T(2,overlap) have to be set, otherwise they stay NaN.
Apart from that rather special case great function!
Also please check the tags on this file, there are some commas missing.

Looks like MATLAB is moving underneath, so time for an update.

[xnew,ynew] = intersections(x,y1,x,y2)
Warning: NARGCHK will be removed in a future release. Use NARGINCHK or NARGOUTCHK instead.
> In intersections (line 91)

straight forward and fast, thanks!

Works perfectly

Working great, finds the intersection correctly and quickly even for a limited number of spaced out points; just what I needed.

Thanks so much!

Great function! It's easy to work with and very well documented.
I've used it in Matlab R2016a and I got the message bellow:
" Warning: NARGCHK will be removed in a future release. Use NARGINCHK or NARGOUTCHK instead. "

It works when you just apply the suggested fix, but maybe it would be relevant to already update the script as well

oh man this is crazy good

thanks a lot

Excellent! Thank you very much!

Very helpful in my coursework! Thank you!!

Great work. Just what I needed

Hello Douglas, great job! I was desperately looking for a function like this. However, I have to use this function in another language.

Can you please provide supporting information or documents that explain the concept of the function? This would help me in writing the code in other language.

Thanks!

works like a charm!

Amazing work, 5*

Perfect function. Thanks for sharing.

Pedro, thanks for the suggestion. I tried bsxfun when I first wrote this and it was not faster at that time. Perhaps things have changed so i will revisit it.

Thanks for this functions. It was very helpfull.
I would just remocomend the use of bsxfun instead of repmat, simpy because it's faster. The equivalent would be I believe:
[i, j] = find(bsxfun(@le, min(x1(1:end-1),x1(2:end)), max(x2(1:end-1),x2(2:end)).') & ...
bsxfun(@ge, max(x1(1:end-1),x1(2:end)), min(x2(1:end-1),x2(2:end)).') & ...
bsxfun(@le, min(y1(1:end-1),y1(2:end)), max(y2(1:end-1),y2(2:end)).') & ...
bsxfun(@ge, max(y1(1:end-1),y1(2:end)), min(y2(1:end-1),y2(2:end)).'));

Thanks again, Pedro

Ilya and Jan, because of floating point arithmetic, it's impossible to find intersections perfectly in all cases. Jan, your example has two curves that touch at a single point; some people might define this as an intersection. Your assertion that (0,0) isn't an intersection is debatable.

It also erroneously finds contact points (no real intersections).
Example: [x0, y0] = intersections([-1,0,-1], [-1,0,1], [1, 0, 1], [-1,0,1], 1);
returns point (0, 0) as intersection point although it isn't.

An excellent function! However, if both intersecting curves already include the intersection point, weird results are possible (see the provided example). This problem was already touched by John Mahoney (the results were different).

My test case: several lines pass through one point and every line must eventually include this point. However, it's not known in advance that we have this situation, it simply may appear.

x1=[-0.49313932739246, -0.02127781500161, 0.450583697389237];
y1=[-10.01, 0, 10.01];
x2=[-1.22073877122679, -0.02127781500161, 1.17818314122357];
y2=[10.01, 0, -10.01];
[x,y,seg1,seg2] = intersections(x1, y1, x2, y2)

The case is published because it may be useful to know about. Otherwise I agree that the function works fine in 99,99999999% of cases.

exactly what i'm looking for.
thanks

where to find intersection.m file

thanks to douglas sir. it is running .

hello sir i asked a question on mathwork ask the question please help me out this question.
i want to intersection point of curve by the .fig file
</matlabcentral/answers/uploaded_files/12588/Unti.jpg> this is my graph fig. i want to intersection point of efficiency & unit discharge at 100 unit speed & other unit speed. so how we use matlab code finding intersection of this.

Deepak, you must install intersections.m into a folder that is on your MATLAB path. The error message is simply informing you that MATLAB cannot find the function.

when i run the code
t= 0:pi/64:3*pi;
y = sin(t);
>> y2=cos(t);
>> [xout,yout] = intersections(t,y,t,y2,1);
plot(t,y,'linewidth',2)
set(gca,'xlim',[min(t) max(t)],'ylim',[-1.1 1.1])
hold on
plot(t,y2,'g','linewidth',2)
plot(xout,yout,'r.','markersize',18)

that problem is occured:-

"Undefined function 'intersections' for input arguments of type 'double'."

Thank you!

Just what i was looking for! thanks

Thanks.

Does what it should, no problems so far.

John,

Welcome to the fun of floating-point arithmetic. First, the duplicate points *are* removed, but this is not obvious unless you display all the digits of the result. For your third example, the four "duplicates" are actually all different, but only in the final bit (you can see it with "format hex"). I don't think there's anything I can do to make this work correctly in all cases.

Second, why should something fail if an intersection occurs at a segment boundary? Some people would like to consider this a proper intersection and others not (I know this from several discussions on the topic). As a user, you will have to evaluate the results you get for your particular application.

In many real world applications, where the two curves are formed from real data, this issue never comes up at all. It's (mostly) only from textbook type problems that you'll see this behavior.

Doug

Douglas,

Thanks for the nice function.
Could you please clarify the behavior depending on the ROBUST flag?

From the function doc:

% ROBUST (optional) set to 1 or true means to use a slight variation of the
% algorithm that might return duplicates of some intersection points, and
% then remove those duplicates. The default is true, but since the
% algorithm is slightly slower you can set it to false if you know that
% your curves don't intersect at any segment boundaries. Also, the robust
% version properly handles parallel and overlapping segments.

First, based on the first sentence, it would seem that the removal of duplicates is handled by intersections.m. This is not the case.

Second, it would seem that something should fail if the curves *do* intersect at segment boundaries (points specifying the line) and the ROBUST = 0/false is used.
At least for the example below, this is not the case.

Thanks for your help!

% Test intersections function
% Depending on whether ROBUST is true or false, and depending on how points
% defining each line do or don't lay (exactly = machine prec?) on the other
% line, the algorithm returns a different number of intersections.

clear all

% % two lines where no point is exactly on other line
% x1 = [0.1, 0.4, 0.5, 0.8];
% y1 = [0.2, 0.5, 0.7, 0.9];
% x2 = [0.1, 0.4, 0.8];
% y2 = [0.9, 0.3, 0.1];

% % first line has point on the other, but not vice versa
% x1 = [0.1, 0.25, 0.5, 0.8];
% y1 = [0.2, 0.6, 0.7, 0.9];
% x2 = [0.1, 0.4, 0.8];
% y2 = [0.9, 0.3, 0.1];

% % two lines that share a common point
% x1 = [0.1, 0.4, 0.5, 0.8];
% y1 = [0.2, 0.3, 0.7, 0.9];
% x2 = [0.1, 0.4, 0.8];
% y2 = [0.9, 0.3, 0.1];

robust = 0;
[xinter, yinter, iinter, jinter] = intersections(x1, y1, x2, y2, robust);

figure
hold on
scatter(xinter, yinter, 200, 'r','filled')
plot(x1, y1, '.b-','markersize',30)
plot(x2, y2, '.g-','markersize',20)
text(0.5, 0.5, ['found ', num2str(numel(xinter)), ' intersection/s'], 'fontsize',14)

did anyone used this code just for one curve? could you help me with the transformation? i am new in matlab, and i can't achieve a working code for just one curve. thank you so so much.

Paul,

Thank you for your comments. First of all, you should know that the definitions for I and J were specifically requested by users -- they were not included at all in early versions of the program. I have no opinion about them since I don't use them myself. Certainly, backwards compatibility is important so they are not going to change.

I would encourage you to write a simple wrapper function that runs intersections and then converts its outputs into whatever you require. I hesitate to do that for you since I'm not completely sure what it is you need, but feel free to contact me via email if you need help with it.

Doug

This is an absolutely fantastic piece of code. 5* well deserved. Thank you!

That said I really regret writing this but...

To me the returned values for I and J feel awkward. Because X0 and Y0 return the intersection points in terms of the original space in which the curves lie (X,Y), I would expect the second set of returned values to give coordinates in terms of the curves themselves. So if a curve has (Euclidean) length 4.8, then I would expect 2.4 to be halfway the curve, not at 4 tenth of the second line segment.

This way we could also compute derivatives and integrals along a curve and interpolate these values simply by looking op the value at I (or J for that matter). Furthermore, the distance along the curve is a nice linear function, whereas the currently returned values for I and J are piecewise linear and non-differentiable it the points.

I understand that such a change of I and J would not be backwards compatible, and also it might be slightly (but hardly) more involved to compute, so I see the big arguments against. Nevertheless, I would be very happy if you would at least give it a thought :-)

Excellent work

Really good, faster and as useful as polyxpoly.

Dear Sir/Madam,

I am looking for the intersection of dynamic waves and base line. SO far, I could manage to get few of the intersection points using ratio (Y data of dynamic wave/ Y data of base line) for whole of samples. Pls advise me how should I continue since I was stuck in here so long.

I used your function but mine is dynamic wave form and base straing line so the X -axis is not constant line. So the array is only (1900x1) for wave form and base line. The Y values of wave form is changing but the base line is seems like constant. Let say, there is 64 intersection points visuality but what i can get is only 47 but some are not accutare positions.

Thanks and best regards
Kyaw

Very nice work. Thanks for sharing.

I agree that the index identification is a very useful extra feature.

Very nice work. The option to identify index pairs between which the intersection points lie is a nice touch. This saves me several hours of work, and the code has worked great on my test cases so far. Thanks so much!

Jack,

It just occurred to me that you can use the 3rd and 4th output arguments to tell you in which line segments the intersections lie. See the help for more info.

Doug

Jack,

No, it's not random.

For the NON-ROBUST option:
Number the line segments of your polygon from 1 to N, as defined by the order of the x- and y-coordinates. For each P1P2 line segment, the intersections are found in order of the polygon segments. In other words, if P1P2 crosses polygon segments 12 and 25, the intersection with segment 12 will be first followed by the intersection with segment 25, because 12 < 25.

For the ROBUST option:
Because unique() is run on the intersections, they end up being sorted in increasing order of the x-coordinate (with any tie broken by sorting y).

Doug

Hi,

I've got a pollygon that encloses the origin and I have 2 points P1 and P2 outside the pollygon such that the line P1P2 intersects the pollygon in 2 points, I1 and I2, that I determine using the 'intersections' function. Then I rotate my line P1P2 by 10 degrees 36 times and I determine the intersection points with the pollygon for each angle. The results are stored in 2 vectors, x0 and y0, each with 2 components. Now,what I am trying to figure out here is how does the function choose to allocate the numerical values of the coordinates of the 2 points to x0(1) and x0(2) respectively?More explicitly, what does x0(1) represent, the x coordinate of which one of the 2 intersection points? Is it random? Cheers!

Using in combination with polymorph(), thanks!

Thank you so much for sharing this wonderful code with us.

One small issue though, I found a case where instead of detecting 3 intersections it finds 5. Perhaps it is pathological but let me know if you would like more info.

Thanks a million!

Thanks a lot bruder !!!

NURUL, please send me an example.

Hi, I have used this code to find the coordinates intersection between 2 vectors, where X1, Y1 is my latitude and longitude for vector 1, and X2,Y2 is my latitude and longitude for vector 2.
Previously this function works very well in my data, but now, I am not so sure why it is sometimes cannot find the intersections, even it is implemented on the same datasets. I am really wondering why?

Laura,

You have defined two curves consisting of just one line segment each. These two segments do not intersect so that is why my function returns empty vectors.

In order to find the intersection of the two infinite straight lines defined by the segments you can set

x1=[627 619];
y1=[338 288];
x2=[621 627];
y2=[351 339];

and then follow the calculations in the comments in intersections.m. Set up the A and B matrices and then solve with T = A\B. Since you are looking for an intersection outside the endpoints you can ignore t1 and t2 and just look at x0 and y0.

Doug

Dear Douglas,

My question: is it possible to obtain with you code the intersection of 2 curves if their intersection is outside the endpoints given?
thank you!
Laura

I have asked my question wrongly. The intersection is outside the interval of the selected points:
e.g;

C1=(627,338),(619,288 )
C2= (621,351), (627,339)

because this gives empty for me and it should intersect

Excellent code ... well written, fast and above all it works!!! Has saved me a huge amount of coding time!

Thanks for the file. Saved me lines of code and hours of computing time.

The function seems great but I get an error about the try/catch loop syntax for my version of matlab.
changing line 218 to "catch, err" avoids the error.

May I suggest an improvement to your code.

Currently the code checks for intersecting domains and then for intersecting domain checks if the lines actually cross. This op takes 16 flops (that I count) per pair of lines, with then the LU decomp and back substitution for a 4 by 4 system of equations where the domains do cross.

However if the intersection condition is stated in terms of L1 = O1 + k1*D1 and L2 = O2 + k2*D2 where L1 == L2, you are left with an equation in 2 unknowns (k1 and k2). This involves ~7 operations to find k1 and k2, with four more to find x0 and y0.

N.B. Your implementation of this was far faster than my one using matlab classes (of the order of 3000-4000 times faster) :).

Cheers

David

2Douglas Schwarz
Thank you! Thank you! Thank you!

2Truc Phan
1. Not all have the Mapping Toolbox.
2. Polyxpoly does not work when you want to find the intersection curve itself (loops in curve).

Reliable and robust, but it takes 6.6 seconds for two 10,000 point-functions with ~20 intersections, with the robust flag off. May be I can figure out how to make it less robust, but a bit faster.

Thanks for sharing.

However, can you please tell us the major difference and/or advantage between this function and the polyxpoly function in Mapping Toolbox. Thanks

Excellent work!

I have submitted a code doing similar task but faster. You can find it here:

http://www.mathworks.com/matlabcentral/fileexchange/27205

I would appreciate your feedback.

Ya thanks Doug. Didnt have your email address so contacted you with a private message on matlab central itself. If you can give me you email i can mail you our code so far & ear images as an attachment so you can get a better idea. Thanks

Rui, please email me directly with details of your problem and I will see if I can help.
Doug

Hi douglas!! i am currently doing a final year project on ear recognition.the idea of finding the feature vector to be trained involves finding the distances of a centroid from the intersection points(lines have been drawn at different angles to intersect a edge detected ear image)..but we are unable to find the coordinates of the intersections.hoping u can help us out asap

thanking u in anticipation
Rui

Excellent! It does exactly what is needed and saved so much of time. I am using this function heavily, it saved 300 lines of my code which was inefficient in finding the intersection points.

A suggestion to all: Ignore the post by M MA, this function works perfect, and is fast as stated by Douglas.

Thanks! And i hope you will help us all with such great functions in the future.

Shadab

Arrgh! I actually did exactly that and then removed it thinking (mistakenly) that it wasn't needed. I ended up using reshape(i,[],1) rather than i(:) because it is documented to force the second dimension to be 1 even if i starts out as 0-by-0. i(:) does the right thing today, but maybe not in the future. Thanks again, Rodrigo! (New version just uploaded.)

Douglas,
there is still a little bug related to mismatching dimensions of empty matrices. When we try, for instance,
[xi,yi,i,j] = intersections([0,1],[0,0,],[1,0],[1,1])
that's OK. However, if we try
[xi,yi,i,j] = intersections([0,1],[0,0,],[1,0.5,0],[1,0.5,1])
there is the error "Matrix dimensions should agree", because empty matrices should have compatible dimensions, even though they are empty entities.

So I have look at the code, and I have inserted the following code
i = i(:); j = j(:);
just after "i" and "j" definitions at line 165. In this way it seems that this little annoying bug was removed. Could you check it please?

Personally (and strangely) I am very pleased about this bug, because reading the help on empty matrices I've learned that they have a strong operational structure underneath. As a matter of fact, they behave like non-empty matrices, when we try to make math operations on them!

Congratulations and thanks for the prompt response.

Rodrigo

Rodrigo, I agree. I just fixed it and uploaded a new version -- it should appear soon. Thanks for the report.

A really really excellent function. I'm using all the time. However I have a suggestion. When the curves do not interecept, like in this example
[xi,yi,i,j] = intersections([0,1],[0,0,],[1,0],[1,1])
instead to produce an error, it would be better, in my opinion, if it produces empty matrices as an output. This is partially done when the function is called with two outputs, such as,
[xi,yi] = intersections([0,1],[0,0,],[1,0],[1,1])
but not with four. I give you five points anyway!
Rodrigo

An excellent function. Met my needs 100%. Thanks for making life easier for me.

Perfect

Excellent!
Found the link for this on the following thread <http://www.mathworks.de/matlabcentral/newsreader/view_thread/258260#672305>
It was just what I needed :)

Excellent file. It works perfectly, and extremely quickly. Just what I was looking for!

Yuri, no problem! I'm glad it's resolved and that you're still following this thread. -Doug

Doug, thanks and sorry for my stupid mistake and your time. I didn't test it yet, but believe that was the problem.

It bothered me that I never resolved Yuri's issue so I just took another look at it and realized that his code has a bug. When he gets X2 and Y2 he is using the handle H1 instead of H2. After fixing that there is no problem.

Yuri,

Please email me a mat-file and m-file that demonstrate what you are saying and I will investigate further.

Doug

Douglas, thank you for comment.
That was exactly my point. At smaller N the polygons don't touch, but the function does find an intersection. If I use robust as true, much more intersection points on the 1st polygon are found.
At high N many points from two circles are actually very close, almost ideal touch. (Although I confirm there is a small gap at very close zoom in.) But at let's say N=1000, nothing found at robust=0, and 923 points all around the 1st circle at robust=1. Why is it?
Yuri

Yuri,

If you zoom in on the supposed intersection point I think you'll find that there is no intersection at all. Remember, those circles are made up of line segments and so are actually polygons with N sides. For most values of N the two polygons won't even touch. In fact, the closest they can ever come is to touch at just one point and because of floating point round-off effects, my routine may or may not find that point.

Doug

Great function, works perfectly in almost all cases. However I found situation when it gives unpredictable result.

I draw 2 circles (using famous CIRCLE function from FEX), which touch each other.
clf
N=30;
hold on
H1=circle([0 0],5,N);
X1=get(H1,'XData');
Y1=get(H1,'YData');
H2=circle([0 10],5,N);
X2=get(H1,'XData');
Y2=get(H1,'YData');
[x,y]=intersections(X1,Y1,X2,Y2,0);
plot(x,y,'ro')
hold off
axis equal

It doesn't find intersection even with larger N and returns many warnings like
Warning: Matrix is singular to working precision.
or
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.559090e-016.

If I use N=10 or N=20 it actually finds intersection but in a weird place, specially at 20.

Thanks a lot anyway. 4 points.

Doug,

Fanstastic job on this. The inclusion of indices is a great convenience as well. Thank you for sharing this code and continuing to improve on your work.

Jeremy

Great code, works perfectly, and fast.

Excellent, thank you.

Fantastic!

Thanks so much Doug.

FYI: the reason for wanting a index is for the case where you have a set of spatial vectors (x1,y1) and (x2,y2) but you also have a third set of data (i.e. t1) which is in one-to-one correspondence with (x1,y1). The indices which you have now included make it trivial to obtain the value of t1 where (x1,y1) and (x2,y2) intersect.

Hopefully this make sense.

Thanks again,
Levi

Levi,

I'm not sure what use that would be, but it was trivial to implement and doesn't hurt the execution time so I have no objection to adding it. I have just uploaded a new version.

Doug

This is a great function. Best intersection routine out there so far. No doubt.

One thing is missing though:
Returning the indices of the intersection. I think that Tim Meyers got this piece right with "intxy" ID:15239. He used floating point indices to indicate the index (read his description). However, his routine does not return all intersections (only the first), and I would guess it is not as fast/robust as this one.

I haven't looked through this function long enough to figure out how difficult such an addition would be. However, I will say that this routine would (IMHO) become MUCH more powerful/useful if some sort of index returning were included.

I'd be happy to work with you on implementing this if you gave me an idea where to start. I've done a bit of thinking about the "intersection problem", but clearly not as much as you.

Thanks,
Levi

Works perfectly.

This routine works well and quickly.

One general comment for users: "intersections" is a general curve intersection routine, not a polygon itersection routine, although it can be used for this as well. When testing polygons, add the first point of the boundary to the end of the boundary list - "intersections" does not add the last point automatically like "fill", nor is it intended to.

Scott

Gabriel, thanks for the good example. The definition of intersection is debatable and your example only touches at the end points. Nevertheless, I have reverted to using \ in my algorithm instead of linsolve as it seems to be more accurate (at least in this case). This particular failure didn't occur in any of my test cases. The new version should be available shortly.

The function was working well, but i found no intersection for [X,Y]=intersections([1 6],[1 2],[6 8],[2 2]) when [6 2] actually is an intersection point. Other intersection tools do find this point....

The best intersection tool on FC, thanks.

A very nicely written and well documented function that is also efficient. I appreciate the error checking and handling of special cases. The comments make it straightforward to understand how the function works. I need to know which line segments contain intersection points and it was easy to add a couple of statements to return this information ( i(in_range) and j(in_range) ). I can confirm that the function is slower if A is changed from sparse to full. Does anyone know why?

This is exactly the sort of function that I hope to find whenever I browse the FEX.

Does exactly what it says on the tin. Excellent!

you're my hero.works like a dream; thanks very much

Excellent! Just what I am looking for. Thanks a lot!

Works very fine.

Excellent!

Good job. Well coded and well documented.

Doug has done a very nice job here, writing a code that is both fast and carefully constructed, watching out for degeneracies. Excellent explanations internally, in case you wish to know how he did it. Ignore the sour grapes review by M Ma. Or, if you prefer, consider the following comparison:

n = 100;
x1 = rand(1,n);
x2 = rand(1,n);
y1 = rand(1,n);
y2 = rand(1,n);

tic,[x0,y0]=intersections(x1,y1,x2,y2);toc
Elapsed time is 0.282314 seconds.

tic,[x0m,y0m]=meetpoint(x1,y1,x2,y2);toc
Elapsed time is 2.219815 seconds.

Both codes found the same points of intersections (2183 of them in total.) But intersections was considerably faster.

M MA,
I'm sorry you felt you needed to give my function such a low rating, but it is not doing quite the same thing as yours (meetpoint). Comparing test suites is a little unfair under these circumstances. Nevertheless, I have attempted to make my function behave better when presented with pathological data. Speed and accuracy with real data are my priorities. I have just uploaded a new version -- it should be available soon.

Well, the function meetpoint (http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=8997&objectType=FILE) already works with vertical lines... but the rate 2 is because I found one error at the first test!! using the example in meetpoint:
x1=[0 1.5 3 5 5];
y1=[0 1.5 1.5 2 1];
x2=[.2 .5 .5 1 1 2 2 5 5];
y2=[.2 .5 .8 1 1.3 1.3 1.8 1.8 .5];
[x,y,m1,m2]=meetpoint(x1,y1,x2,y2);
[x0,y0]=intersections(x1,y1,x2,y2);

The output of meetpoint is:
>> [x' y']
ans =
0.5000 0.5000
1.0000 1.0000
1.3000 1.3000
2.0000 1.5000
4.2000 1.8000
5.0000 1.8000

While the output of intersections is:
>> [x0 y0]
ans =
0.5000 0.5000
1.0000 1.0000
1.3000 1.3000
2.0000 1.5000
4.2000 1.8000

So, it looks like its not handling vertical lines (both) correctly as meetpoint does!
Moreover, meetpoint also returns the slopes at the intersections!

Finally! A code that will work with vertical lines. Perfect timing on posting this too! Thanks!!!!