src/processing/Color.cpp

#include <cmath>

#include <FCam/Event.h>

namespace FCam {

    const float RGBtoXYZ[9] = {
    0.4124564, 0.3575761, 0.1804375,
    0.2126729, 0.7151522, 0.0721750,
    0.0193339, 0.1191920, 0.9503041
    };

    int xyToCCT(float x, float y) {
        const float x_e = 0.3366;
        const float y_e = 0.1735;
        const float A0 = -949.86315;
        const float A1 = 6253.80338;
        const float t1 = 0.92159;
        const float A2 = 28.70599;
        const float t2 = 0.20039;
        const float A3 = 0.00004;
        const float t3 = 0.07125;

        float n = (x - x_e)/(y - y_e);

        float CCT = A0 + A1 * std::exp(-n/t1) + A2 * std::exp(-n/t2) + A3 * std::exp(-n/t3);
        if (CCT < 3000 || CCT > 50000) {
            warning(Event::OutOfRange, "xyToCCT: Conversion only accurate within 3000 to 50000K, result was %d K.\n", CCT);
        }
        return std::floor(CCT+0.5f);
    }

    void kelvinToXY(int T, float *x, float *y) {
        if (!x || !y) return;

        if (T < 1667 || T > 25000) {
            warning(Event::OutOfRange, "KelvinToXY: Conversion only accurate within 1667K to 25000K, requested %d K\n", T);
        }
        // chromaticity x coefficients for T <= 4000K
        const float A_x00 = -0.2661239;
        const float A_x01 = -0.2343580;
        const float A_x02 = 0.8776956;
        const float A_x03 = 0.179910;
        // chromaticity x coefficients for T > 4000K
        const float A_x10 = -3.0258469;
        const float A_x11 = 2.1070379;
        const float A_x12 = 0.2226347;
        const float A_x13 = 0.24039;
        // chromaticity y coefficients for T <= 2222K
        const float A_y00 = -1.1063814;
        const float A_y01 = -1.34811020;
        const float A_y02 = 2.18555832;
        const float A_y03 = -0.20219683;
        // chromaticity y coefficients for 2222K < T <= 4000K
        const float A_y10 = -0.9549476;
        const float A_y11 = -1.37418593;
        const float A_y12 = 2.09137015;
        const float A_y13 = -0.16748867;
        // chromaticity y coefficients for T > 4000K
        const float A_y20 = 3.0817580;
        const float A_y21 = -5.87338670;
        const float A_y22 = 3.75112997;
        const float A_y23 = -0.37001483;

        float invKiloK = 1000.f/T;
        float xc, yc;
        if (T <= 4000) {
            xc = A_x00*invKiloK*invKiloK*invKiloK +
                A_x01*invKiloK*invKiloK +
                A_x02*invKiloK +
                A_x03;
        } else {
            xc = A_x10*invKiloK*invKiloK*invKiloK +
                A_x11*invKiloK*invKiloK +
                A_x12*invKiloK +
                A_x13;
        }

        if (T <= 2222) {
            yc = A_y00*xc*xc*xc +
                A_y01*xc*xc +
                A_y02*xc +
                A_y03;
        } else if (T <= 4000) {
            yc = A_y10*xc*xc*xc +
                A_y11*xc*xc +
                A_y12*xc +
                A_y13;
        } else {
            yc = A_y20*xc*xc*xc +
                A_y21*xc*xc +
                A_y22*xc +
                A_y23;
        }

        *x = xc;
        *y = yc;

    }

    void invert3x3(float *in, float *out) {
        out[0] = in[4]*in[8]-in[5]*in[7];
        out[3] = in[5]*in[6]-in[3]*in[8];
        out[6] = in[3]*in[7]-in[4]*in[6];

        float invdet = 1.0f/(in[0]*out[0] + in[1]*out[3] + in[2]*out[6]);

        out[0] *= invdet;
        out[3] *= invdet;
        out[6] *= invdet;

        out[1] = invdet*(in[7]*in[2]-in[8]*in[1]);
        out[4] = invdet*(in[8]*in[0]-in[6]*in[2]);
        out[7] = invdet*(in[6]*in[1]-in[7]*in[0]);
        
        out[2] = invdet*(in[1]*in[5]-in[2]*in[4]);
        out[5] = invdet*(in[2]*in[3]-in[0]*in[5]);
        out[8] = invdet*(in[0]*in[4]-in[1]*in[3]);

    }

    void invert3x3(double *in, double *out) {
        out[0] = in[4]*in[8]-in[5]*in[7];
        out[3] = in[5]*in[6]-in[3]*in[8];
        out[6] = in[3]*in[7]-in[4]*in[6];

        double invdet = 1.0/(in[0]*out[0] + in[1]*out[3] + in[2]*out[6]);

        out[0] *= invdet;
        out[3] *= invdet;
        out[6] *= invdet;

        out[1] = invdet*(in[7]*in[2]-in[8]*in[1]);
        out[4] = invdet*(in[8]*in[0]-in[6]*in[2]);
        out[7] = invdet*(in[6]*in[1]-in[7]*in[0]);
        
        out[2] = invdet*(in[1]*in[5]-in[2]*in[4]);
        out[5] = invdet*(in[2]*in[3]-in[0]*in[5]);
        out[8] = invdet*(in[0]*in[4]-in[1]*in[3]);

    }

    void colorMatrixInterpolate(const float *a, const float *b, float alpha, float *result) {
        // Check for identity interpolations
        if (alpha == 0.f) {
            for (int i=0; i < 12; i++)  result[i] = a[i];
            return;
        } 
        if (alpha == 1.f) {
            for (int i=0; i < 12; i++)  result[i] = b[i];
            return;
        }
        // just do a linear interpolation between two known settings
        for (int i = 0; i < 12; i++) {
            result[i] = (1-alpha)*a[i] + alpha*b[i];
        }

        // rescale the result so that the minimum row sum is one
        float rowSum[3] = {result[0] + result[1] + result[2],
                           result[4] + result[5] + result[6],
                           result[8] + result[9] + result[10]};
        float scale = 1.0f;
        if (rowSum[0] < rowSum[1] && rowSum[0] < rowSum[2]) {
            scale = 1.0f/rowSum[0];
        } else if (rowSum[1] < rowSum[2]) {
            scale = 1.0f/rowSum[1];
        } else {
            scale = 1.0f/rowSum[2];
        }

        for (int j = 0; j < 3; j++) {
            for (int i = 0; i < 3; i++) {
                result[j*4+i] *= scale;
            }
        }                           

    }

}