2

I've written a function to calculate the cosine, sine and degrees of the angle between three points of which I have the x and y coordinates - point 1 (x1, y1), point 2 (x2, y2) and point 3 (x3, y3). I've written the function and been trying to test it out, but I'm not completely confident in how accurate it is. Does anyone know if I've made a mistake in my calculations?

def path_angle_degree(x1, y1, x2, y2, x3, y3):
    u = (x2 - x1, y2 - y1)
    v = (x3 - x2, y3 - y2)
    norm_u = math.sqrt(u[0] * u[0] + u[1] * u[1])
    norm_v = math.sqrt(v[0] * v[0] + v[1] * v[1])

    # this conditional is to check there has been movement between the points
    if norm_u < 0.001 or norm_v < 0.001:
        return (None, None, None)
    prod_n = norm_u * norm_v
    dot_uv = u[0] * v[0] + u[1] * v[1]
    cos_uv = dot_uv / prod_n

    # fixes floating point rounding
    if cos_uv > 1.0 or cos_uv < -1.0:
        cos_uv = round(cos_uv)
    radians = math.acos(cos_uv)
    sin_uv = math.sin(radians)
    degree = math.degrees(radians)
    return (cos_uv, sin_uv, degree)

An example of this function being called on a straight path would be:

print(path_angle_degree(6,6,7,6,8,6))

Thanks so much!

Jack Simpson
  • 1,681
  • 3
  • 30
  • 54
  • excuse me, how do you locate an angle between 3 collinear points? in planar geometry angle is a figure formed by two rays. and you cant pass 2 rays from this points. – marmeladze Apr 16 '15 at 07:32
  • 2
    @marmeladze without knowing what the OP intends, you could simply construct two vectors from the three points, assuming one point is common to both vectors -- for example, the points (a,b,c) could be used to contruct one vector (a->b) and a second (b->c). (or (a->b) and (a->c), etc.) – jedwards Apr 16 '15 at 07:36
  • 2
    *"I'm not completely confident in how accurate it is"* - why not? Why don't you calculate a few examples yourself, and compare outputs? – jonrsharpe Apr 16 '15 at 07:39
  • 1
    Do you care about the sign of your angle? E.g., would you want `path_angle_degree(6,6,7,6,8,7)` and `path_angle_degree(6,6,7,6,8,5)` to give different results? (Or equivalently, do you want to be able to distinguish between a turn to the left and a turn of the same magnitude to the right?) – Mark Dickinson Apr 16 '15 at 07:58
  • FWIW, you _could_ simplify your norm calcs: `norm_u = math.hypot(*u)`. You might like to round values of `cos_uv` when `abs(cos_uv)>1-EPSILON`, where `EPSILON` is something small, like 5.0E-12. – PM 2Ring Apr 16 '15 at 07:59
  • 1
    For the angle, all you need is a single call to `math.atan2`: no square roots or other trig calls necessary. `math.atan2(u[0]*v[1]-u[1]*v[0], u[0]*v[0]+u[1]*v[1])` should do it. – Mark Dickinson Apr 16 '15 at 08:05
  • Thanks so much for the responses guys, so I've been tracking a tagged insect for an experiment and so a typical path (filmed at 25 fps) will have several thousand x and y coordinates where the insect was located. I'm trying to see how frequently it changed direction, because then I can separate out extensive from intensive activity. So ever 3 points (or more if I try to smooth by subsampling) I'll calculate the change in angle, cosine and sine. Hope that clears up the reason I'm doing this :) – Jack Simpson Apr 16 '15 at 13:14
  • Yes, I've been trying to work out how to distinguish between cases when the insect completely reverses direction like you showed Mark - I was getting frustrated when I couldn't distinguish between them. – Jack Simpson Apr 16 '15 at 13:22
  • 1
    @JackSimpson: You're safe for direction reversals with just the `acos` method: it'll give an unsigned angle between 0 and pi. It's only if you want to be able to distinguish clockwise from anticlockwise turns that you need more than the dot product. And can I just say: Insect Math! Cool! – Mark Dickinson Apr 16 '15 at 15:20

2 Answers2

1

In addition to being expensive, sin_uv is unstable and doesn't cope with numerous edge cases. Use the cross product instead (you only need the z-component).

Also, you'll find it simpler and cheaper to normalize u and v before computing the products.

Once you have cos_uv and sin_uv, use atan2 to get the angle.

Marcelo Cantos
  • 181,030
  • 38
  • 327
  • 365
  • As a clarification, the analogue of the cross product in 2D is generally known as the "perpendicular dot product" or "perp dot product". – Sneftel Apr 16 '15 at 07:59
  • The cross product by itself doesn't provide enough information, though: it can't distinguish between an angle of 45° and 135°, for example. – Mark Dickinson Apr 16 '15 at 12:07
  • @MarkDickinson: Together with the dot product it does, but I should have spelled it out. I've amended my answer. – Marcelo Cantos Apr 16 '15 at 21:52
0

You might want to look at the arctan2() function.

user50619
  • 325
  • 5
  • 14