Projecting a point onto a line

Hi,

I hope you can help.

I have some x,y data that has a shoulder I would like to identify with a signature of a transition as a function of some variable. I have many data files and would like do this in a systematic way. I figured the best way would be to draw a chord across the shoulder and find the maximum distance between the chord and the data. The line between the chord and the data point should be orthogonal to the chord. The chord does not go through the origin.

I have tried two ways and I can't understand why they don't work, I suspect it is ignorance of the maths. Which I stole from here: http://cs.nyu.edu/~yap/classes/visual/03s/hw/h2/math.pdf

% write function that projects the  point (q = X,Y) on a vector
% which is composed of two points - vector = [p0x p0y; p1x p1y]. 
% i.e. vector is the line between point p0 and p1. 
%
% The result is a point qp = [x y] and the length [length_q] of the vector drawn 
% between the point q and qp . This resulting vector between q and qp 
% will be orthogonal to the original vector between p0 and p1. 
% 
% This uses the maths found in the webpage:
% http://cs.nyu.edu/~yap/classes/visual/03s/hw/h2/math.pdf
%

function [ProjPoint, length_q] = ProjectPoint(vector, q)

      p0 = vector(1,:);
      p1 = vector(2,:);
      length_q = 1; %ignore for now

      a = [p1(1) - p0(1), p1(2) - p0(2); p0(2) - p1(2), p1(1) - p0(1)];
      b = [q(1)*(p1(1) - p0(1)) + q(2)*(p1(2) - p0(2)); ...
          p0(2)*(p1(1) - p0(1)) - p0(1)*(p1(2) - p0(2))] ;

      ProjPoint = a\b;
  end

This 'works' but the point on the line that it returns is not at the point on the chord that would be form a line with the data point orthogonal to the chord. It is orthogonal to the x-axis.

I won't post the second way which I tried because it gives the same answer, I suspect that once I understand what is wrong with this the other one can be fixed too.

Thanks!

[Ajithabh K S]

Hi Ali,

Please try the below code

% vector = [1,1; 4,4];
% q = [3,2];
% I used above values for calling function

function [ProjPoint] = proj(vector, q)
p0 = vector(1,:);
p1 = vector(2,:);
a = [-q(1)*(p1(1)-p0(1)) - q(2)*(p1(2)-p0(2)); ...
    -p0(2)*(p1(1)-p0(1)) + p0(1)*(p1(2)-p0(2))]; 
b = [p1(1) - p0(1), p1(2) - p0(2);...
    p0(2) - p1(2), p1(1) - p0(1)];
ProjPoint = -(b\a);
end

See the plot

%

[Saneesh Ali]

%%% Try using this code
%%% vector = [1,1; 4,4];
%%% p = [3,2];
%%% a = vector(1,:)
%%% b = vector(2,:)

function result = point_on_line(a,b,p)
ap = p-a;
ab = b-a;
result = a+(dot(ap,ab)/dot(ab,ab))*ab;
end

Answer 1

[Ali]

Note that this was working fine. The problem was that the scales of the two axes were different and so the line connecting the projected point onto the line look didn't look like it was at a right angle.

Answer 2

[Pedro Souza]

You can set the scales of the two axis to equal with

 axis equal