Cheap algorithm to find measure of angle between vectors
Asked Answered
E

10

31

Finding the angle between two vectors is not hard using the cosine rule. However, because I am programming for a platform with very limited resources, I would like to avoid calculations such as sqrt and arccos. Even simple divisions should be limited as much as possible.

Fortunately, I do not need the angle per se, but only need some value that is proportional to said angle.

So I am looking for some computationally cheap algorithm to calculate a quantity that is related to the angle between two vectors. So far, I haven't found something that fits the bill, nor have I been able to come up with something myself.

Ephor answered 15/9, 2009 at 14:9 Comment(2)
hmm: important question: are the components of the vectors stored in fixed-point or floating-point format?Pelerine
Neither. Since the coordinates in question are pixel-coordinates, they are always integer values. No floating-point/fixed point is necessary. So I guess you could say they're fixed point with a multiplier of 1 :)Ephor
P
13

Have you tried a CORDIC algorithm? It's a general framework for solving polar ↔ rectangular problems with only add/subtract/bitshift + table, essentially doing rotation by angles of the form tan-1 (2-n). You can trade off accuracy with execution time by altering the number of iterations.

In your case, take one vector as a fixed reference, and copy the other to a temporary vector, which you rotate using the cordic angles towards the first vector (roughly bisection) until you reach a desired angular accuracy.

(edit: use sign of dot product to determine at each step whether to rotate forward or backward. Although if multiplies are cheap enough to allow using dot product, then don't bother with CORDIC, perhaps use a table of sin/cos pairs for rotation matrices of angles π/2n to solve the problem with bisection.)

(edit: I like Eric Bainville's suggestion in the comments: rotate both vectors towards zero and keep track of the angle difference.)

Pelerine answered 15/9, 2009 at 14:34 Comment(3)
+1. CORDIC is what you want. It may be faster to rotate the two vectors simultaneously towards the X axis, keeping track of the angle difference.3d
... and CORDIC will give you the norm of each vector for free.3d
Nice. CORDIC seems like it will do what I need. I haven't implemented this yet, but will mark your answer as "accepted" anyway. Thanks!Ephor
A
26

If you don't need the actual euclidean angle, but something that you can use as a base for angle comparisons, then changing to taxicab geometry may be a choice, because you can drop trigonometry and it's slowness while MAINTAINING THE ACCURACY (or at least with really minor loosing of accuracy, see below).

In main modern browser engines the speedup factor is between 1.44 - 15.2 and the accuracy is nearly the same as in atan2. Calculating diamond angle is average 5.01 times faster than atan2 and using inline code in Firefox 18 the speedup reaches factor 15.2. Speed comparison: http://jsperf.com/diamond-angle-vs-atan2/2.

The code is very simple:

function DiamondAngle(y, x)
{
    if (y >= 0)
        return (x >= 0 ? y/(x+y) : 1-x/(-x+y)); 
    else
        return (x < 0 ? 2-y/(-x-y) : 3+x/(x-y)); 
}

The above code gives you angle between 0 and 4, while atan2 gives you angle between -PI and PI as the following table shows:

enter image description here

Note that diamond angle is always positive and in the range of 0-4, while atan2 gives also negative radians. So diamond angle is sort of more normalized. And another note is that atan2 gives a little more precise result, because the range length is 2*pi (ie 6.283185307179586), while in diamond angles it is 4. In practice this is not very significant, eg. rad 2.3000000000000001 and 2.3000000000000002 are both in diamond angles 1.4718731421442295, but if we lower the precision by dropping one zero, rad 2.300000000000001 and 2.300000000000002 gives both different diamond angle. This "precision loosing" in diamond angles is so small, that it has some significant influence only if the distances are huge. You can play with conversions in http://jsbin.com/bewodonase/1/edit?output (Old version: http://jsbin.com/idoyon/1):

enter image description here

The above code is enough for fast angle comparisons, but in many cases there is need to convert diamond angle to radians and vice verca. If you eg. have some tolerance as radian angles, and then you have a 100,000 times loop where this tolerance is compared to other angles, it's not wise to make comparisons using atan2. Instead, before looping, you change the radian tolerance to taxicab (diamond angles) tolerance and make in-loop comparisons using diamond tolerance and this way you don't have to use slow trigonometric functions in speed-critical parts of the code ( = in loops).

The code that makes this conversion is this:

function RadiansToDiamondAngle(rad)
{
  var P = {"x": Math.cos(rad), "y": Math.sin(rad) };
  return DiamondAngle(P.y, P.x);
}  

As you notice there is cos and sin. As you know, they are slow, but you don't have to make the conversion in-loop, but before looping and the speedup is huge.

And if for some reason, you have to convert diamond angle to radians, eg. after looping and making angle comparisons to return eg. the minimum angle of comparisons or whatever as radians, the code is as follows:

function DiamondAngleToRadians(dia)
{
  var P = DiamondAngleToPoint(dia);
  return Math.atan2(P.y,P.x);
}

function DiamondAngleToPoint(dia)
{
  return {"x": (dia < 2 ? 1-dia : dia-3), 
  "y": (dia < 3 ? ((dia > 1) ? 2-dia : dia) : dia-4)};
}

Here you use atan2, which is slow, but idea is to use this outside any loops. You cannot convert diamond angle to radians by simply multiplying by some factor, but instead finding a point in taxicab geometry of which diamond angle between that point and the positive X axis is the diamond angle in question and converting this point to radians using atan2.

This should be enough for fast angle comparisons.

Of course there is other atan2 speedup techniques (eg. CORDIC and lookup tables), but AFAIK they all loose accuracy and still may be slower than atan2.

BACKGROUND: I have tested several techniques: dot products, inner products, law of cosine, unit circles, lookup tables etc. but nothing was sufficient in case where both speed and accuracy are important. Finally I found a page in http://www.freesteel.co.uk/wpblog/2009/06/encoding-2d-angles-without-trigonometry/ which had the desired functions and principles.

I assumed first that also taxicab distances could be used for accurate and fast distance comparisons, because the bigger distance in euclidean is bigger also in taxicab. I realized that contrary to euclidean distances, the angle between start and end point has affect to the taxicab distance. Only lengths of vertical and horizontal vectors can be converted easily and fast between euclidean and taxicab, but in every other case you have to take the angle into account and then the process is too slow (?).

So as a conclusion I am thinking that in speed critical applications where is a loop or recursion of several comparisons of angles and/or distances, angles are faster to compare in taxicab space and distances in euclidean (squared, without using sqrt) space.

Apiece answered 3/2, 2013 at 18:51 Comment(3)
This is exactly what I needed—thank you. But I named the parameters dX and dY, since they're really axis deltas, not point ordinates. I think it's clearer that way.Jiva
Very clever solution! Some people might think this is only useful for 2D vectors, but actually it can be easily extended to 3D, 4D, etc. Given any two vectors A, B, you set x = dot(A,B) and y = cross(A,B) and then d = DiamondAngle(y,x) gives the angle between A and B in any dimension.Benedix
from what I see, the ifology case can be simplified to e.g. return (y >= 0) ? (x >= 0 ? y/(x+y) : 1+x/(x-y)) : (x < 0 ? 2+y/(x+y) : 3+x/(x-y)); - which highlights other possible simplifications here IMVHOAquarius
E
13

Back in the day of a few K of RAM and machines with limited mathematical capabilities I used lookup tables and linear interpolation. The basic idea is simple: create an array with as much resolution as you need (more elements reduce the error created by interpolation). Then interpolate between lookup values.

Here is an example in processing (original dead link).

You can do this with your other trig functions as well. On the 6502 processor this allowed for full 3D wire frame graphics to be computed with an order of magnitude speed increase.

Emmet answered 15/9, 2009 at 14:16 Comment(5)
+1: interpolation into tables -- fast and small. Not so accurate, but often accurate enough. 256 distinct cosine values is often enough.Eruct
Thanks for the suggestion, but that would still only get rid of the arccos. determining the remaining elements of the cosine rule in this case (the lengths of all sides of the triangle) would still require several sqrt's.Ephor
Less than that can be enough. If you're willing to do a little extra calculation, you only need values between zero and half of pi/90 degrees, and of course you need to decide the granularity and accuracy you need. It's hard to say anything definite about details without considering the individual case.Karwan
SQRT is trivial to replace with an approximation loop that converges rapidly with plenty of accuracy. And it involves no divides, just multiplies.Eruct
CORDIC is a great solution if you don't have a multiplier, but if you do have one on your system I suspect that the interpolation will be faster (realizing that you can build such a table for SQRT as well over your domain). CORDIC is iterative while lookup/interpolate is fixed cost. My domain of experience is 3D on low end hardware though, so maybe your hardware will find CORDIC amiable.Emmet
P
13

Have you tried a CORDIC algorithm? It's a general framework for solving polar ↔ rectangular problems with only add/subtract/bitshift + table, essentially doing rotation by angles of the form tan-1 (2-n). You can trade off accuracy with execution time by altering the number of iterations.

In your case, take one vector as a fixed reference, and copy the other to a temporary vector, which you rotate using the cordic angles towards the first vector (roughly bisection) until you reach a desired angular accuracy.

(edit: use sign of dot product to determine at each step whether to rotate forward or backward. Although if multiplies are cheap enough to allow using dot product, then don't bother with CORDIC, perhaps use a table of sin/cos pairs for rotation matrices of angles π/2n to solve the problem with bisection.)

(edit: I like Eric Bainville's suggestion in the comments: rotate both vectors towards zero and keep track of the angle difference.)

Pelerine answered 15/9, 2009 at 14:34 Comment(3)
+1. CORDIC is what you want. It may be faster to rotate the two vectors simultaneously towards the X axis, keeping track of the angle difference.3d
... and CORDIC will give you the norm of each vector for free.3d
Nice. CORDIC seems like it will do what I need. I haven't implemented this yet, but will mark your answer as "accepted" anyway. Thanks!Ephor
T
5

Here on SO I still don't have the privilege to comment (though I have at math.se) so this is actually a reply to Timo's post on diamond angles.

The whole concept of diamond angles based on the L1 norm is most interesting and if it were merely a comparison of which vector has a greater/lesser w.r.t. the positive X axis it would be sufficient. However, the OP did mention angle between two generic vectors, and I presume the OP wants to compare it to some tolerance for finding smoothness/corner status or sth like that, but unfortunately, it seems that with only the formulae provided on jsperf.com or freesteel.co.uk (links above) it seems it is not possible to do this using diamond angles.

Observe the following output from my Asymptote implementation of the formulae:

Vectors : 50,20 and -40,40
Angle diff found by acos      : 113.199
Diff of angles found by atan2 : 113.199
Diamond minus diamond         : 1.21429
Convert that to degrees       : 105.255
Rotate same vectors by 30 deg.
Vectors : 33.3013,42.3205 and -54.641,14.641
Angle diff found by acos      : 113.199
Diff of angles found by atan2 : 113.199
Diamond minus diamond         : 1.22904
Convert that to degrees       : 106.546
Rotate same vectors by 30 deg.
Vectors : 7.67949,53.3013 and -54.641,-14.641
Angle diff found by acos      : 113.199
Diff of angles found by atan2 : 113.199
Diamond minus diamond         : 1.33726
Convert that to degrees       : 116.971

So the point is you can't do diamond(alpha)-diamond(beta) and compare it to some tolerance unlike you can do with the output of atan2. If all you want to do is diamond(alpha)>diamond(beta) then I suppose diamond is fine.

Truffle answered 3/6, 2013 at 14:2 Comment(1)
Yes, you are right. Diamond angles can only be compared but not added, subtracted, multiplied or divided. The only exceptions are degrees 0, 45, 90, 135, 180, 225, 270, 315, 360. You can safely transfer these angles between taxicab and euclidean spaces using divide/multiply. Transferring any other angle causes erroneous and this error varies from 0 to 4.074568596222775 degrees. I have tried to figure out the logic in this variance and if it can be take into account to increase the space-to-space-transferring accuracy.Marriott
C
2

The cross product is proportional to the angle between two vectors, and when the vectors are normalized and the angle is small the cross product is very close to the actual angle in radians due to the small angle approximation.

specifically:

I1Q2-I2Q1 is proportional to the angle between I1Q1 and I2Q2.

Compton answered 18/12, 2010 at 23:9 Comment(0)
O
1

The solution would be trivial if the vectors were defined/stored using polar coordinates instead of cartesian coordinates (or, 'as well as' using cartesian coordinates).

Officinal answered 15/9, 2009 at 14:13 Comment(1)
Sadly, this is not the case. The only way would be to convert the coordinates to polar in my code, but that runs into the same problem that I'm trying to solve in the first place.Ephor
R
1

dot product of two vectors (x1, y1) and (x2, y2) is

x1 * x2 + y1 * y2 

and is equivilent to the product of the lengths of the two vectors times the cosine of the angle between them.

So if you normalize the two vectors first (divide the coordinates by the length)

Where length of V1 L1 = sqrt(x1^2 + y1^2),  
  and length of V2 L2 = sqrt(x2^2 + y2^2),

Then normalized vectors are

(x1/L1, y1/L1),  and (x2/L2, y2/L2),  

And dot product of normalized vectors (which is the same as the cosine of angle between the vectors) would be

 (x1*x2 + y1*y2)
 -----------------
     (L1*L2)

of course this may be just as computationally difficult as calculating the cosine

Roswell answered 15/9, 2009 at 14:22 Comment(2)
OP knows that, dot product requires an arccos to yield an anglePelerine
@Jason, yes, but dot product by itself is a "quantity that is related to the angle" you don;t necessarily have to convert it to radians or degrees to use it...Roswell
T
1

if you need to compute the square root, then consider using the invsqrt hack.

acos((x1*x2 + y1*y2) * invsqrt((x1*x1+y1*y1)*(x2*x2+y2*y2)));
Thump answered 15/9, 2009 at 14:33 Comment(1)
Unless there are overflow issues, I'd combine the invsqrt's into one call invsqrt((x1*x1+y1*y1)*(x2*x2+y2*y2)); better to run a complex algorithm once rather than twice.Pelerine
G
0

The dot product might work in your case. It's not proportional to the angle, but "related".

Gatling answered 15/9, 2009 at 14:16 Comment(2)
OP knows that, dot product requires an arccos to yield an anglePelerine
Or square roots to normalize lengths. Both of which the OP is trying to avoid.Karwan
L
0

To make more cheap to compute metric of vectors difference angle I make simplified version of 'full radian' equation and it uses only 1 division and many multiplications (can be put to SIMD easily enough to process many pairs of verctors in parallel):

Also only one condition (sign check only) - it can be also put to SIMD with sign bit extraction and use as mask (or compare with zero).

No sqrt() and no arccos() used.

It returns vectors difference as positive only metric in fixed range of 0..2 (float). Where 0 is full coherent vectors (zero angle difference) and 2 is full counter-directed (max possible angle difference).

Current C-function is:

float fDiffAngleVect(int x1, int y1, int x2, int y2)
{
  float fResult = 0.0f;
  // check if any of 2 input vectors is zero vector - return 0
  if ((x1 == 0) && (y1 == 0) || (x2 == 0) && (y2 == 0))
    return 0.0f;

  int iUpper = x1 * x2 + y1 * y2;

  if (iUpper > 0)
  {
    fResult = 1.0f - ((float)(iUpper * iUpper) / (float((x1 * x1 + y1 * y1) * (x2 * x2 + y2 * y2))));
  }
  else
  {
    fResult = 1.0f + ((float)(iUpper * iUpper) / (float((x1 * x1 + y1 * y1) * (x2 * x2 + y2 * y2))));
  }

  return fResult;

}

My use case is integer only x,y vectors to compare so you can rewrite it to all-float input and output.

Lindsaylindsey answered 9/4, 2024 at 14:19 Comment(0)

© 2022 - 2025 — McMap. All rights reserved.