Compressed unit vectors

compressed_unit_vector.cpp

/*
    Reference:
    ----------

      A Survey of Efficient Representations for Independent Unit Vectors
*/

#include <cmath>
#include <cstdint>
#include <iostream>

#include "foundation/math/scalar.h"
#include "foundation/math/vector.h"
#include "foundation/utility/iostreamop.h"

using namespace std;
using namespace foundation;

class CompressedUnitVector32
{
  public:

    CompressedUnitVector32() {}

    explicit CompressedUnitVector32(const Vector3f& vec)
    {
        int16_t projected[2];

        const float inv_l1_norm = 1.0f / (fabs(vec[0]) + fabs(vec[1]) + fabs(vec[2]));

        if (vec[2] <= 0.0f)
        {
            projected[0] = floor16((1.0f - fabs(vec[1] * inv_l1_norm)) * copysign(1.0f, vec[0]));
            projected[1] = floor16((1.0f - fabs(vec[0] * inv_l1_norm)) * copysign(1.0f, vec[1]));
        }
        else
        {
            projected[0] = floor16(vec[0] * inv_l1_norm);
            projected[1] = floor16(vec[1] * inv_l1_norm);
        }

        int16_t best_projected[2];
        float error = 0.0f;

        unsigned int bits_x = projected[0];
        unsigned int bits_y = projected[1];

        for (int i = 0; i < 2; ++i)
        {
            for (int j = 0; j< 2; ++j)
            {
                projected[0] = bits_x + i;
                projected[1] = bits_y + j;

                Vector3f decoded = oct_decode(projected);

                const float alt_error = fabs(dot(vec, decoded));
                if (alt_error > error)
                {
                    error = alt_error;
                    best_projected[0] = projected[0];
                    best_projected[1] = projected[1];
                }
            }
        }

        m_bits[0] = best_projected[0];
        m_bits[1] = best_projected[1];
    }

    operator Vector3f() const
    {
        return oct_decode(m_bits);
    }

  private:
    static float clamp_minus_plus_one(const float x)
    {
        if (x < -1.0f) return -1.0f;
        if (x >  1.0f) return  1.0f;
        return x;
    }

    static int16_t round16(const float f)
    {
        return round(clamp_minus_plus_one(f) * 32767.0f);
    }

    static int16_t floor16(const float f)
    {
        return floor(clamp_minus_plus_one(f) * 32767.0f);
    }

    static Vector3f oct_decode(const int16_t bits[2])
    {
        Vector3f vec;

        const float rcp_32767 = 1.0f / 32767.0f;
        vec[0] = clamp_minus_plus_one(bits[0] * rcp_32767);
        vec[1] = clamp_minus_plus_one(bits[1] * rcp_32767);
        vec[2] = 1.0f - (fabs(vec[0]) + fabs(vec[1]));

        if (vec[2] < 0.0f)
        {
            const float tmp = vec[0];
            vec[0] = (1.0f - fabs(vec[1])) * copysign(1.0f, tmp);
            vec[1] = (1.0f - fabs(tmp))   * copysign(1.0f, vec[1]);
        }

        return normalize(vec);
    }

    int16_t m_bits[2];
};

int main(int argc, char** argv)
{
    cout << "Oct32 test program" << endl;

    Vector3f v(-0.1f, 2.11f, 7.0f);
    v = normalize(v);

    CompressedUnitVector32 oct_v(v);
    Vector3f v_rt(oct_v);

    cout << "V original = " << v    << endl;
    cout << "V decoded  = " << v_rt << endl;

    return 0;
}