Finding the point of intersection on a 3D line per

Question

I have a line and a point and I wanted to find a point (x,y,z) on the line that is 90 degrees or perpendicular if I were to draw a line from this point of intersection with another point.

So far I can create a line with this code and i have another code that calculates angle between three points, but that doesn't really apply here:

a = [1 1 2]; %line
b = [20 28 90]; % line

c = [50 30 67]; %point 

ab = b - a;


n = max(abs(ab)) + 1;

s = repmat(linspace(0, 1, n)', 1, 3);
for d = 1:3
    s(:, d) = s(:, d) * ab(d) + a(d);
end


s = round(s);


Z = 100; 
N = 100;
X = zeros(N, N, Z);

X(sub2ind(size(X), s(:, 1), s(:, 2), s(:, 3))) = 1;


x = c(:,1);


clf

plot3(s(:, 1), s(:, 2), s(:, 3), 'r.-')


axis(N * [0 1 0 1 0 1])
grid on

Answer 1

This will require a bit of mathematics to determine that analytically. By "90 degrees", I'm assuming you want to find the point on this 3D line that would be perpendicular to this line if you extended a line from this point of intersection to the desired point.

I'm assuming the two points a and b denote the coordinates in 3D space where a line can join them and that c is the point of interest. Here's a better diagram of what I'm talking about:

^{Source: MathWorld}

In your case, x1 and x2 denote a and b in your code and x0 denotes c. The distance d would be the distance on the line that would allow the point to be perpendicular to this line if you extended a line from the point of intersection to the point c.

You can define a parametric equation that describes the line between x1 and x2 like so:

A point on this line between x1 and x2 can be described by taking each of the (x,y,z) values for x1 and x2 and writing it in the above parametric form and varying the parameter t, which goes from [0,1]. Therefore t=0 would give you the first point x1 or a and t=1 would give you the second point x2 or b. Any value of t between [0,1] would give you a point along the line. The objective is to find the value t that would minimize the distance from x0 or c to this line. Like I said before, we all know from geometry that the smallest distance from a point to a line would make the angle of intersection perpendicular / 90 degrees if you extended a line from this point of intersection to the point x0 or c.

Therefore, all you have to do is find this value of t, then substitute this into the above parametric equation to find your desired point. In other words, we want to minimize the distance between the point and the line, and the distance can be described like so:

To find the minimum distance, you'd find the parameter t that minimized the above equation by finding the derivative with respect to t and setting it equal to 0. Logistically, you would take the square root of the equation so that you would be minimizing the distance, and not the distance squared. However, it's actually easier to minimize the distance squared, which is why the equation above is represented like so. It makes sense because if you minimize the distance squared.... the distance would also be minimized as you'd simply just place a square root on the answer to get what you asked for. Eliminating the square root from the equation would make calculating the derivative easier.

If you do that, and solve for t, we get this equation:

Therefore, find the difference between a and c, take this with the dot product with the difference between b and a, then divide this by the magnitude squared of the difference between b and a. This solves for t, and then you'd substitute this into the above parametric equation to find your point.

In MATLAB code, it could look something like this:

a = [1 1 2]; %line - x1
b = [20 28 90]; % line - x2

c = [50 30 67]; %point - x0

ab = b - a; %// Find x2 - x1

%// -(x1 - x0).(x2 - x1) / (|x2 - x1|^2)
t = -(a - c)*(ab.') / (ab*ab.'); %// Calculate t

%// Find point of intersection
Xinter = a + (b - a)*t;

The code for t I took advantage of using matrix multiplication. The dot product can be found by multiplying a row array by a column array and in a similar fashion, if the row array and column array have the same coefficients, this results in the magnitude squared of a vector.

For your example, we get:

Xinter =

   16.9889   23.7211   76.0539

To show that this is right, let's plot the line, and the point and the point of intersection:

The code to produce the above figure is:

figure;

%// Plot line
plot3([a(1) b(1)], [a(2) b(2)], [a(3) b(3)]);
hold on;

%// Plot point of interest in red
plot3(c(1), c(2), c(3), 'r.');

%// Plot intersection point in green
plot3(Xinter(1), Xinter(2), Xinter(3), 'g.');

%// Plot line from intersection point to point of interest in black
plot3([c(1) Xinter(1)], [c(2) Xinter(2)],  [c(3) Xinter(3)], 'k');

%// Turn on a grid
grid;