Find point of intersection between two vectors in MATLAB

Solution 1:

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.

Solution 2:

Well, you really have two points on two different lines, and you want to find the intersection. The easiest way is to find the equations of the two lines and then calculate the intersection.

The equation of a line is given by y = mx + b where m is the slope and b is the y-intercept. For one line you have two points which gives two equations. So, you can solve for the constants m and b. This gives the following two equations:

 0 = 0*m + 1*b  % Using the first point x=y=0 into y=m*x+b
 6 = 6*m + 1*b  % Using the second point x=y=6

Or in matrix form:

 [ 0 ] = [ 0 1 ]* [ m ]
 [ 6 ]   [ 6 1 ]  [ b ]

For the first line the constants can be calculated in MATLAB by

 C1 = inv([0 1;6 1]*[1;0]; % m=C1(1) and b=C(2)

Now that you have the equation for the two lines you can solve for the intersection by solving the following system of equations (which are obtained by manipulating the equation for a line):

 m_1*x-y = -b_1
 m_2*x-y = -b_2

All that is left is to write the above system of equations in matrix form and solve:

 [x] = inv [m_1 -1] * [-b_1]
 [y]       [m_2 -1]   [-b_2]

Or in MATLAB syntax:

 I = inv([m_1 -1; m_2 -1])*[-b_1;-b_2]; % I is the intersection.

Notes

• As per gnovice's comment if the lines are actually line segments you need to check if the intersection is between the end points of the line segments.

• If the two slopes are equal, m_1 = m_2, then there will either be no intersection or infinitely many intersections.

Solution 3:

For a general multi-dimensional solution, what you're actually doing is solving a series of linear systems.

First you need to reduce the equations to linear form: Ax+By=C (expand dimensions as necessary)

For two-points:

y - y1 = (y2 - y1) / (x2 - x1) * (x - x1)
y - y1 = (y2 - y1) / (x2 - x1) * x - (y2 - y1) / (x2 - x1) * x1
(y1 - y2) / (x2 - x1) * x + y - y1 = (y1 - y2) / (x2 - x1) * x1
(y1 - y2) / (x2 - x1) * x + y = (y1 - y2) / (x2 - x1) * x1 + y1
(y1 - y2) * x + (x2 - x1) * y = (y1 - y2) * x1 + (x2 - x1) * y1

A = (y1 - y2)
B = (x2 - x1)
C = (y1 - y2) * x1 + (x2 - x1) * y1 = A * x1 + B * y1

For your example:

x1 = 0, x2 = 6, y1 = 0, y2 = 6
A1 = (0 - 6) = -6
B1 = (6 - 0) = 6
C1 = A1 * 0 + B1 * 0 = 0

x1 = 0, x2 = 6, y1 = 6, y2 = 0
A2 = (6 - 0) = 6
B2 = (6 - 0) = 6
C2 = A2 * 0 + B2 * 6 = 6 * 6 = 36

Then, form a matrix, with A B and C in rows:

[A1 B1 C1]
[A2 B2 C2]

[-6 6 0]
[ 6 6 36]

Now reduce to reduced echelon form using the Matlab function rref(matrix):

[ 1 0 3]
[ 0 1 3]

As you can guess, the last column is your intersection point. This is expandable to as many dimensions as necessary. If your reduced echelon form has something other than the identity matrix for the front part of it, your vectors either do not have a unique intersection point, or have no intersection point, depending on the form of the matrix.

dim = 2;

% Do other stuff, ending with rref(matrix)

if (matrix(:,1:dim) == eye(dim))
    % Matrix has unique solution.
    solution = (matrix(:,dim+1))'
else
    % No unique solution.
end

In two dimensions, the variations are:

[ 1 B C]
[ 0 0 0]

[ 1 B C1]
[ 0 0 C2]

Solution 4:

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.