-2

I have two ordered arrays of x(xcor) and y(ycor) values. Joining the first and last points gives a line. I want to compute the perpendicular distances of all points from this line. This is similar to least square distance. Is there a direct way to do this in matlab? Please also note that sign of the distance should represent the side on the points lie of the line. enter image description here

 xy =
     -121.9067  -53.5483
     -122.0750  -53.5475
     -122.4750  -53.5243
     -123.0975  -53.4835
     -123.9050  -53.4168
     -124.8050  -53.3235
     -125.7025  -53.2467
     -126.5675  -53.1800
     -127.3825  -53.1215
     -128.1500  -53.0798
     -128.8825  -53.0468
     -129.6000  -53.0452
     -130.3150  -53.1133
     -131.0400  -53.2532
     -131.7850  -53.4513
     -132.5525  -53.6877
     -133.3425  -53.9345
     -134.1600  -54.1758
     -135.0075  -54.4115
     -135.8675  -54.6480
     -136.7375  -54.9040
     -137.5075  -55.2635
     -138.1875  -55.7435
     -138.7775  -56.3333
     -139.2850  -57.0665
     -139.8450  -57.9285
     -140.4550  -58.9492
     -141.1575  -60.0988
     -141.9825  -61.3415
     -142.9275  -62.6172
     -144.0050  -63.8517
     -145.2125  -65.0523
     -146.5450  -66.1715
     -147.9950  -67.1727
     -149.5575  -68.0570
     -151.2225  -68.8152
     -152.9925  -69.4493
     -154.8625  -69.9500
     -156.8300  -70.3063
     -158.8700  -70.5280
     -160.9050  -70.6017
     -162.8550  -70.6287
     -164.6525  -70.7372
     -165.5367  -70.7550
     -166.3450  -70.8620
Abhishek Bhatia
  • 9,404
  • 26
  • 87
  • 142
  • This is a math problem. As I understand this problem, you have a set of points and a line. For each point, you want to find the shortest distance to the line. What you want to do is, for each point A, find a point on the line B such that segment AB is perpendicular to the line. Then everything else is trivial. – Matthew Gunn Nov 27 '15 at 06:06
  • Note that you're calculating total least square distance here. – Jonas Nov 27 '15 at 07:02
  • @MatthewGunn: many questions asked here are math problems. Often, Matlab has a built-in function to do that for you. Unfortunately not in this case. – Jonas Nov 27 '15 at 07:04
  • @Jonas yeah you are right. – Abhishek Bhatia Nov 27 '15 at 07:07
  • @Shai: How is this a duplicate of a Mathematica question? The geometrical problem is similar, but the language not quite. – Jonas Dec 02 '15 at 07:42
  • What is a positive, what is a negative sign? – Jonas Dec 05 '15 at 20:53
  • It doesn't matter but the same convention should be followed on all points. – Abhishek Bhatia Dec 05 '15 at 21:09

2 Answers2

8

If you have a vector AB, the distance from a point C to that vector can be calculated as follows:

  • Normalize the vector AB
  • Calculate the vector AC
  • Project the vector AC onto AB
  • Subtract the projection from AC
  • Calculate the length of the result

In other words, you split AC into a component that is parallel to AB and a component that is perpendicular, and you calculate the length of the latter.

If you have arrays x and y, you can do the following

xy = [x(:),y(:)];
abVector = xy(end,:) - xy(1,:); %# a is the first, b the last point
abVectorNormed = abVector./norm(abVector);

acVector = bsxfun(@minus, xy, xy(1,:));

acParallelLength = sum(bsxfun(@times, acVector , abVectorNormed ),2);
acParallelVector = bsxfun(@times, acParallelLength, abVectorNormed );

perpendicularVector = acVector - acParallelVector;

perpendicularDistance = sqrt(sum(perpendicularVector.^2,2));

EDIT You asked for figures because the code "does not work" in your hands. See below the figures (top: raw data; bottom: perpendicular distance) and the command to plot them; the data looks fairly reasonable in my eyes.

enter image description here

subplot(2,1,1),plot(xy(:,1),xy(:,2),'or')
hold on, plot([xy(1,1),xy(1,1)+abVector(1)],[xy(1,2),xy(1,2)+abVector(2)],'b')
hold on, plot([xy(1,1)+acParallelVector(:,1),xy(:,1)]',[xy(1,2)+acParallelVector(:,2),xy(:,2)]','r')
axis equal %# important to see right angles as such
subplot(2,1,2),stem(xy(:,1),perpendicularDistance,'r')
ylabel('perpendicular distance')
Jonas
  • 74,690
  • 10
  • 137
  • 177
  • It shows me this error: `Error using sqrt Too many input arguments.` perpendicularVectoris 45x 2. – Abhishek Bhatia Nov 27 '15 at 07:11
  • @AbhishekBhatia: Sorry, typo. Fixed now. – Jonas Nov 27 '15 at 07:17
  • Thanks for edit! What does AB represent here? Does it represent the line? – Abhishek Bhatia Nov 27 '15 at 08:00
  • @AbhishekBhatia the line goes from point A to point B, thus I call it "abVector" – Jonas Nov 27 '15 at 08:38
  • Can you plot the line computed and points. I am getting an incorrect answer, I just wanted to be sure. I have added points in the question for reproducing it. – Abhishek Bhatia Nov 27 '15 at 08:50
  • I am unsure how you are representing the line with two parameters a and b instead of 3. – Abhishek Bhatia Nov 27 '15 at 08:56
  • 2
    @AbhishekBhatia: I've plotted the figures with all the lines. Result looks reasonable to me. Why do you think it's wrong? – Jonas Nov 27 '15 at 11:40
  • Thanks so much! Sorry there some problem with my dataset. Can you please possible. Add sign to distance, that is, a point on right has +ve distance from the line and point on left has a negative distance from the line. – Abhishek Bhatia Dec 01 '15 at 15:05
  • @AbhishekBhatia: augment `perpendicularVector` and `acParallelVector` with a column of zeros, take the cross product (`cross`) and look at the sign of the third column of the result. – Jonas Dec 02 '15 at 07:52
  • Please make changes in the answer. – Abhishek Bhatia Dec 02 '15 at 22:11
  • @AbhishekBhatia: Why? The question has not changed. – Jonas Dec 03 '15 at 09:15
  • Okay I added it now. Please do so in the answer if possible. – Abhishek Bhatia Dec 04 '15 at 12:26
  • Please update the answer that considers the sign soon. – Abhishek Bhatia Dec 23 '15 at 13:18
  • @AbhishekBhatia are you having difficulties implementing [Jonas's answer](http://stackoverflow.com/questions/33950784/least-square-distances#comment55829872_33951857), or are you just being needy? – Andras Deak -- Слава Україні Dec 23 '15 at 16:46
  • I tried this.`sign_vector=cross([perpendicularVector,zeros(size(perpendicularVector,1),1)],[acParallelVector,zeros(size(acParallelVector,1),1)]); signed_log=sign_vector(:,3)./abs(sign_vector(:,3)); perpendicularDistance=perpendicularDistance.*signed_log; perpendicularDistance=perpendicularDistance(~isnan(perpendicularDistance));` It would be good to verify and accept the answer I guess. – Abhishek Bhatia Dec 23 '15 at 17:57
  • 1
    @AbhishekBhatia it very well should be. Note that you have to assign the bounty yourself to Jonas. Otherwise he won't get the bounty, since his answer is older than your bounty. (And even if the answer was newer, he'd only get *half* the assigned bounty automatically, which wouldn't be fair: he fully solved your problem, including additional constraints.). Also, note that you can just use `sign(sign_vector)` to avoid `NaN`s. – Andras Deak -- Слава Україні Dec 23 '15 at 19:34
0
function tot_distance=compute_distance2(x,y)
xA=x(1);xB=x(end);yA=y(1);yB=y(end);
a=(yA-yB); 
b=(xB-xA);
c=xA*yB-xB*yA;
d=0;
for p=2:numel(x)-1,
    d=d+(a*x(p)+b*y(p)+c)/sqrt(a^2+b^2);
end
tot_distance=abs(d);
end

Easier way.

Abhishek Bhatia
  • 9,404
  • 26
  • 87
  • 142