Color

/* -*-c++-*-
 */

#include "Math"
typedef unsigned char uint8_t;

float RGBMMaxRange = 8.0;

// https://gist.github.com/aras-p/1199797
void encodeRGBM(float rgb[3], uint8_t rgbm[4]) {
    // in our case,
    const float kRGBMMaxRange = RGBMMaxRange;
    const float kOneOverRGBMMaxRange = 1.0f / kRGBMMaxRange;

    // encode to RGBM, c = ARGB colors in 0..1 floats

    float r = rgb[0] * kOneOverRGBMMaxRange;
    float g = rgb[1] * kOneOverRGBMMaxRange;
    float b = rgb[2] * kOneOverRGBMMaxRange;

    float a = std::max(std::max(r, g), std::max(b, 1e-6f));
    a = ceilf(a * 255.0f) / 255.0f;

    rgbm[0] = uint8_t(std::min(r / a, 1.0f) * 255);
    rgbm[1] = uint8_t(std::min(g / a, 1.0f) * 255);
    rgbm[2] = uint8_t(std::min(b / a, 1.0f) * 255);
    rgbm[3] = uint8_t(std::min(a, 1.0f) * 255.0);
}

void decodeRGM(uint8_t rgbm[4], float rgb[3]) {
    rgb[0] = rgbm[0] / 255.0 * RGBMMaxRange * rgbm[3] / 255.0;
    rgb[1] = rgbm[1] / 255.0 * RGBMMaxRange * rgbm[3] / 255.0;
    rgb[2] = rgbm[2] / 255.0 * RGBMMaxRange * rgbm[3] / 255.0;
}

void encodeRGBE(float rgb[3], uint8_t rgbe[4]) {
    float maxRGB = std::max(rgb[0], std::max(rgb[1], rgb[2]));

    if (maxRGB < 1e-32) {
        rgbe[0] = rgbe[1] = rgbe[2] = rgbe[3] = 0;
    } else {
        int e;
        float v = frexp(maxRGB, &e) * 256.0 / maxRGB;
        rgbe[0] = (unsigned char)(rgb[0] * v);
        rgbe[1] = (unsigned char)(rgb[1] * v);
        rgbe[2] = (unsigned char)(rgb[2] * v);
        rgbe[3] = (unsigned char)(e + 128);
    }
}

void decodeRGBE(uint8_t rgba[4], float rgb[3]) {
    float f = pow(2.0, rgba[3] - (128.0 + 8.0));
    rgb[0] = rgba[0] * f;
    rgb[1] = rgba[1] * f;
    rgb[2] = rgba[2] * f;
}

Vec3f mul(const Vec3f& v, const float* M) {
    Vec3f result;

    result[0] = v[0] * M[0] + v[1] * M[3] + v[2] * M[6];
    result[1] = v[0] * M[1] + v[1] * M[4] + v[2] * M[7];
    result[2] = v[0] * M[2] + v[1] * M[5] + v[2] * M[8];

    return result;
}

// http://graphicrants.blogspot.fr/2009/04/rgbm-color-encoding.html
// M matrix, for encoding
const static float M[] = {0.2209, 0.3390, 0.4184, 0.1138, 0.6780,
                          0.7319, 0.0102, 0.1130, 0.2969};

// Inverse M matrix, for decoding
const static float InverseM[] = {6.0013,  -2.700, -1.7995, -1.332, 3.1029,
                                 -5.7720, 0.3007, -1.088,  5.6268};

Vec4f LogLuvEncode(const Vec3f& vRGB) {
    Vec4f vResult;
    Vec3f Xp_Y_XYZp;
    Vec3f test;

    Xp_Y_XYZp = mul(vRGB, M);

    Xp_Y_XYZp = Xp_Y_XYZp.max(1e-6f);

    vResult[0] = Xp_Y_XYZp[0] / Xp_Y_XYZp[2];
    vResult[1] = Xp_Y_XYZp[1] / Xp_Y_XYZp[2];

    float Le = 2.0 * log2(Xp_Y_XYZp[1]) + 127.0;
    vResult[3] = frac(Le);
    vResult[2] = (Le - (floor(vResult[3] * 255.0f)) / 255.0f) / 255.0f;

    return vResult;
}

void encodeLUV(float rgb[3], uint8_t luv[4]) {
    Vec4f result = LogLuvEncode(Vec3f(rgb[0], rgb[1], rgb[2]));
    luv[0] = uint8_t(result[0] * 255.0);
    luv[1] = uint8_t(result[1] * 255.0);
    luv[2] = uint8_t(result[2] * 255.0);
    luv[3] = uint8_t(result[3] * 255.0);
}

Vec3f LogLuvDecode(const Vec4f& vLogLuv) {
    float Le = vLogLuv[2] * 255.0 + vLogLuv[3];
    Vec3f Xp_Y_XYZp;
    Xp_Y_XYZp[1] = exp2((Le - 127.0) / 2.0);
    Xp_Y_XYZp[2] = Xp_Y_XYZp[1] / vLogLuv[1];
    Xp_Y_XYZp[0] = vLogLuv[0] * Xp_Y_XYZp[2];
    Vec3f vRGB;

    vRGB = mul(Xp_Y_XYZp, InverseM);

    return vRGB.max(0.0);
}

void decodeLUV(uint8_t luv[4], float rgb[3]) {
    Vec3f result =
        LogLuvDecode(Vec4f(luv[0] * 1.0 / 255.0, luv[1] * 1.0 / 255.0,
                           luv[2] * 1.0 / 255.0, luv[3] * 1.0 / 255.0));
    rgb[0] = result[0];
    rgb[1] = result[1];
    rgb[2] = result[2];
}