Fast computation of the direction of a vector

Question

I need to compute the direction of 2D vectors that are provided as a pair of 16 bit signed integers, and return an integer in range $[0,256)$.

This can be achieved by the expression

$$\frac{128}\pi\text{atan}_2(y, x)+ 128$$

or alternatively

$$\frac{128}\pi\arctan\frac yx\stackrel?+128$$ with extra logics to adjust the quadrant.

Given the relatively low accuracy requirement, I hope to find fast approximation. The goal is of course to outperform the above expressions when computed by the FPU (assuming a PC architecture).

Answer 1

You could make use of an existing implementation and map it to your needs:

//  map 16-bit to 8-bit
static byte atan2(short y, short x)
{
    uint at = fxpt_atan2(y, x);
    byte ret = (byte)((at - 0x7FFF) >> 8);

    return ret;
}

/**
 * adapted from https://geekshavefeelings.com/posts/fixed-point-atan2#
 * 
 * 16-bit fixed point four-quadrant arctangent. Given some Cartesian vector
 * (x, y), find the angle subtended by the vector and the positive x-axis.
 *
 * The value returned is in units of 1/65536ths of one turn. This allows the use
 * of the full 16-bit unsigned range to represent a turn. e.g. 0x0000 is 0
 * radians, 0x8000 is pi radians, and 0xFFFF is (65536 / 128) * pi radians.
 *
 * Because the magnitude of the input vector does not change the angle it
 * represents, the inputs can be in any signed 16-bit fixed-point format.
 *
 * @param y y-coordinate in signed 16-bit
 * @param x x-coordinate in signed 16-bit
 * @return angle in (val / 32768) * pi radian increments from 0x0000 to 0xFFFF
 */
const double M_1_PI = 0.31830988618379067154; // The constant 1/pi.

static ushort fxpt_atan2(short y, short x) {
    if (x == y) { // x/y or y/x would return -1 since 1 isn't representable
        if (y > 0) { // 1/8
            return 8192;
        } else if (y< 0) { // 5/8
            return 40960;
        } else { // x = y = 0
            return 0;
        }
    }
    short nabs_y = s16_nabs(y);
    short nabs_x = s16_nabs(x);

    if (nabs_x<nabs_y) { // octants 1, 4, 5, 8
        short y_over_x = q15_div(y, x);
        short correction = q15_mul(q15_from_double(0.273 * M_1_PI), s16_nabs(y_over_x));
        short unrotated = q15_mul((short)(q15_from_double(0.25 + 0.273 * M_1_PI) + correction), y_over_x);

        if (x > 0) { // octants 1, 8   ***  not sure about the '>', could be '<' ***
            return (ushort)unrotated;
        } else { // octants 4, 5
            return (ushort)(32768 + unrotated);
        }
    } else { // octants 2, 3, 6, 7
        short x_over_y = q15_div(x, y);
        short correction = q15_mul(q15_from_double(0.273 * M_1_PI), s16_nabs(x_over_y));
        short unrotated = q15_mul((short)(q15_from_double(0.25 + 0.273 * M_1_PI) + correction), x_over_y);

        if (y > 0) { // octants 2, 3
            return (ushort)(16384 - unrotated);
        } else { // octants 6, 7
            return (ushort)(49152 - unrotated);
        }
    }
}

static short q15_from_double(double d) 
{
    return (short)Math.Round(d * 32768);
}

//  return negative absolute value
static short s16_nabs(short i)
{
    return (i < 0) ? i : (short)-i;
}

static short q15_mul(short j, short k) 
{
    int intermediate = j * (int)k;

    return (short)((intermediate + ((intermediate & 0x7FFF) == 0x4000 ? 0 : 0x4000)) >> 15);
}

static short q15_div(short numer, short denom) 
{
    return (short)(((int) numer << 15) / denom);
}