Find point of intersection between two vectors in MATLAB

Question

I have a very simple MATLAB question. What is the easiest way to find the point of intersection between two vectors. I am not familiar with the various MATLAB functions -- it seems like there should be one for this.

For example, if I have one vector from (0,0) to (6,6) and another vector from (0,6) to (6,0), I need to determine that they intersect at (3,3).

Answer 1

The other results are confusing, verbose and incomplete, IMO. So here's my two cents - also potentially confusing and verbose.

If you are sure that your lines are not skew-parallel or parallel, the following is all you need:

% Let each point be def as a 3x1 array
% Let points defining first line be  : p1, q1
% Let points defining second line be : p2, q2

L = p1-p2;
M = p1-q1;
N = p2-q2;
A = [M N];
T = pinv(A)*L;
h = p1-T(1)*(p1-q1); % h is a 3x1 array representing the actual pt of intersection

Yeah, the Moore-Penrose pseudoinverse is a powerful thing. The explanation for the approach is: You want to find the weights or the scaling factors of the 'direction vectors' (M and N are direction vectors), that linearly combine M and N to give L.

A full description is presented below. It presents a simple exception detection scheme, and their handling is left to the user. (The minimum distance between two line algorithms is from Wikipedia; the comparison of direction cosines (DCS) to check vector attitudes is common knowledge.)

% Let each point be def as a 3x1 array
% Let points defining first line be : p1, q1
% Let points defining second line be: p2, q2

% There are two conditions that prevent intersection of line segments/lines
% in L3 space. 1. parallel 2. skew-parallel (two lines on parallel planes do not intersect)
% Both conditions need to be identified and handled in a general algorithm.

% First check that lines are not parallel, this is done by comparing DCS of
% the line vectors
% L, M, N ARE DIRECTION VECTORS.
L = p1-p2;
M = p1-q1;
N = p2-q2;

% Calculate a normalized DCS for comparison. If equal, it means lines are parallel.
MVectorMagnitude = sqrt(sum(M.*M,2)); % The rowsum is just a generalization for N-D vectors.
NVectorMagnitude=sqrt(sum(N.*N,2)); % The rowsum is just a generalization for N-D vectors.

if isequal(M/MVectorMagnitude,N/NVectorMagnitude) % Compare the DCS for equality
     fprintf('%s\n', 'lines are parallel. End routine')
end;

% Now check that lines do not exist on parallel planes
% This is done by checking the minimum distance between the two lines. If there's a minimum distance, then the lines are skew.

a1 = dot(M,L); b1 = dot(M,M); c1 = dot(M,N);
a2 = dot(N,L); b2 = dot(N,M); c2 = dot(N,N);

s1 = -(a1*c2 - a2*c1)/(b1*c2-b2*c1);
s2 = -(a1*b2 - a2*b1)/(b1*c2-b2*c1);

Sm = (L + s1*M - s2*N);
s = sqrt(sum(Sm.*Sm,2));

if ~isequal(s,0) % If the minimum distance between two lines is not zero, then the lines do not intersect
    fprintf('%s\n','lines are skew. End routine')
end;

% Here's the actual calculation of the point of intersection of two lines.
A = [M N];
T = pinv(A)*L;
h = p1-T(1)*(p1-q1); % h is a 3x1 array representing the actual pt of intersection.

So the pinv approach will give you results even when your M and N vectors are skew (but not parallel, because inv(A'.A) is required to exist). You can use this to determine the minimum distance between two parallel lines or between two parallel planes - to do this, define k = p2+T(2)*(p2-q2), and then the required distance is h-k. Also note that h and k are the points on the lines that are closest to each other IFF lines are skew.

So the use of the pseudoinverse and projection spaces gives us a concise algorithm for:

1. Determining the point of intersection of two lines (not parallel, and not skew)
2. Determining the minimum distance between two lines (not parallel)
3. Determining the points closest to each other on two skew lines.

Concise is not the same as time-efficient. A lot depends on your exact pinv function implementation - MATLAB uses svd which solves to a tolerance. Also, some results will only be approximately accurate in higher dimensions and higher order definitions of the measurement metric (or vector norms). Besides the obvious dimension independent implementation, this can be used in statistical regression analysis and algebraically maximizing likelihood of point estimates.

Answer 2

One solution is to use the equations derived in this tutorial for finding the intersection point of two lines in 2-D (update: this is an internet archive link since the site no longer exists). You can first create two matrices: one to hold the x coordinates of the line endpoints and one to hold the y coordinates.

x = [0 0; 6 6];  %# Starting points in first row, ending points in second row
y = [0 6; 6 0];

The equations from the above source can then be coded up as follows:

dx = diff(x);  %# Take the differences down each column
dy = diff(y);
den = dx(1)*dy(2)-dy(1)*dx(2);  %# Precompute the denominator
ua = (dx(2)*(y(1)-y(3))-dy(2)*(x(1)-x(3)))/den;
ub = (dx(1)*(y(1)-y(3))-dy(1)*(x(1)-x(3)))/den;

And you can now compute the intersection point of the two lines:

xi = x(1)+ua*dx(1);
yi = y(1)+ua*dy(1);

For the example in the question, the above code gives xi = 3 and yi = 3, as expected. If you want to check that the intersection point lies between the endpoints of the lines (i.e. they are finite line segments), you just have to check that the values ua and ub both lie between 0 and 1:

isInSegment = all(([ua ub] >= 0) & ([ua ub] <= 1));

A couple more points from the tutorial I linked to above:

• If the denominator den is 0 then the two lines are parallel.
• If the denominator and numerator for the equations for ua and ub are 0 then the two lines are coincident.