Freesteel Blog

Encoding 2D angles without trigonometry

Friday, June 5th, 2009 at 3:03 pm Written by: Julian

To prove I can still do work, here’s for the machining category.

There are times when you need to encode the direction of a 2D vector. For example, the sequence of arc segments on the bottom rim of the cutter that are in material for the purposes of running the Adaptive Clearing algorithm.

These vector directions from the centre point to different points on the cutter need to be in angular order and analysed in relation to a reference direction.

Often when programmers see a problem to do with angles they use school mathematics and apply trigonometry and convert the direction into its representation in radians whose values they can then compare:

double RadianAngle(x, y)
{
    if (y > 0)
        return atan(x / y) + 1.5707963267948966;
    else if (y < 0)
        return atan(x / y) + 4.7123889803846897; 
    else if (x > 0)
        return 0.0;
    else
        return 3.1415926535897931; 
}

I really hate using radians. How can anyone justify a numbering system that cannot unambiguously represent the important perfect 90 degree right angle in floating point notation without needing to compare to pi/2 to epsilon accuracy all the time? But that’s another story.

The above function, which has the simple inverse (cos(a), sin(a)), has a lot of instability problems for small values of y that can be avoided by introducing more cases:

   if ((y > 0) && (y < x))
      return atan(y / x); 
   else if ((y < 0) && (-y < x))
      return atan(y / x) + 6.2831853071795862; 
   else if ((x < 0) && (x < fabs(y))
      return atan(y / x) + 3.1415926535897931

There is a standard C library function that encapsulates all of this mess into the single function atan2(x, y), which programmers may be happy to use because it hides all this complicated and slow calculation.

But I'm not happy with it because I don't like trigonometry and I need the speed. So I use the following function (subject to the special cases):

double DiamondAngle(x, y)
{
    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 result ranges from 0 to 4 with right angles being integers. Note how the calculation requires a single division and some sums, and you get a number you can use to sort a list of vectors by their angles. The inverse function is:

   return P2((a < 2 ? 1-a : a-3), 
               (a < 3 ? ((a > 1) ? 2-a : a) : a-4); 

This doesn't result in a unit vector, so if you need it to be length 1 you need to normalize the result at the cost of 4 multiplications and one square root, which is still a lot better than any sine or cosine combined.

Check the calculation is correct for x,y>0:

    DiamondAngle(x,y) = y/(x+y) 

which is between 0 and 1, so the inverse comes out as:

  InvDiamondArg(DiamondArg(x,y)) 
                = P2(1-y/(x+y), y/(x+y))

which is inline with (x,y) because both coordinates are positive and if we dot it with the perpendicular:

  Dot(P2(1-y/(x+y), y/(x+y)), P2(y,-x)) 
                = y - y^2/(x+y) - x y /(x+y) 
                = (x y + y^2 - y^2 - x y) / (x+y)

we get zero.

I guess the reason this way of doing it is unpopular is that it looks more complicated.

But if you think about it as projecting the vector down to the piecewise linear diamond shape and then parametrizing by length, it's easy. There's a bit of mess to compress the calculations down and merge the different cases to shorten the code, but once it's there, you don't need to look at it again.

Next week: Encoding 3D angles using the Octohedron