1

I am trying to implement a script to analyze some of my data. I have position information for three points (p1,p2,p3). I would like to find the angular displacement of point p3 from the vector p1p2, as in this picture:

angles

p3a, p3b, p3c, p3d show possible relative positions for p3. As shown, I would like the sign of the output angle to describe its relative position to vector p1p2.

The code I am using is as follows (adapted to the diagram):

v1 = p2 - p1;
x1 = v1(1);
y1 = v1(2);
v2 =  p1 - p3;
x2 = v2(1);
y2 = v2(2);
angle = atan2d(x1*y2-y1*2,x1*x2+y1*y2);

This works as desired when p3 is at p3a, giving a negative angle of the right size (-77 degrees). However, when p3 is at p3d, it outputs a large positive angle (+150 degrees), as opposed to the desired large negative angle.

  • 3
    should you have `v1=p2-p1` and `v2=p3-p1`? – David Feb 25 '15 at 22:39
  • I have tried this approach, but it in fact makes all of the outputs. incorrect. When p3 is at position a, this version returns a large positive value (103degrees), instead of a smaller negative value (-77) and when p3 is at position d, it returns a small negative value (-16) as opposed to a large negative value. – Sebastian Echeverri Feb 25 '15 at 23:29
  • Are you trying to find the angle between the lines P1->P2 and P1->P3, or between P2->P1 and P1->P3? – David Feb 25 '15 at 23:49
  • I am trying to find the angle between P2->P1 and P1->P3. I apologize for the confusion in naming them. I also believe that this is causing a problem for the calculation, as the angle between P2->P1 and P1->P3 is indeed my desired output when P3 is "in front" of P1 and P2 (positions a and b in the diagram), but when P3 is at position c or d, I believe it returns the acute angle between the two vectors as opposed to the desired obtuse one. – Sebastian Echeverri Feb 26 '15 at 00:35

2 Answers2

4

To start with, an easier way to think about the angle between two 2D vectors with coordinates is to align an axis with your coordinate vectors and think about the relationship between two vectors. Using the drawing below, we can see that a relative angle can be found by subtracting one angle from the other.

Source: http://almaer.com/blog/uploads/atan2.png

It is not too hard to figure out looking at this graph that we can say

angle = atan2d(y2,x2) - atan2d(y1,x1)

However, since neither of your vectors are known to be aligned along the coordinate axis, the case can arise where the difference above is not in the range (-180, 180). This means we need to code in a check to add or subtract 360 degrees to get our desired angle:

if abs(angle) > 180
  angle = angle - 360*sign(angle)
end

Note, you are using a kind of reverse notation (CW positive) so the final code would look like:

v1 = p1 - p2;
x1 = v1(1);
y1 = v1(2);
v2 = p3 - p1;
x2 = v2(1);
y2 = v2(2);
angle = atan2d(y1,x1) - atan2d(y2,x2)

if abs(angle) > 180
    angle = angle - 360*sign(angle)
end

Where v1 and v2 have been changed to match your drawing.

whataberk
  • 56
  • 2
  • I think you need separate correction terms for `angle>180` (-360) and `angle<-180` (+360). – eigenchris Feb 26 '15 at 00:11
  • The sign() statement accounts for it. – whataberk Feb 26 '15 at 00:13
  • I just implemented this solution and it works. You are a hero for this! I've been bashing my head against this problem for days. This script is a key part of my current research. I'd love to credit you in the notes of the script - I'll use your username or something else if you would prefer. – Sebastian Echeverri Feb 26 '15 at 00:57
1

Like some of the comments mentioned, I'm a little confused about whether you want to use v2=p1-p3 or v2=p3-p1. Regardless, this method will work with any two vectors v and u where u is the reference vector (the vector we are measuring the angle against).

vx = v(1); vy= v(2); ux = u(1); uy = u(2);
va = -atan2d(vy,vx);         % angle of v relative to x-axis (clockwise = +ve)
ua = -atan2d(uy,ux);         % angle of u relative to x-axis (clockwise = +ve)
A = va - ua;                             % angle va relative to ua
A = A - 360*(A > 180) + 360*(A < -180)   % correction put in [-180,180]

This assumes you want the direction clockwise of u to be taken as the positive direction. Otherwise you just flip the sign of A.

eigenchris
  • 5,791
  • 2
  • 21
  • 30