tinyTVcolor.c

tinyTVcolor.c

// ttvColor.h
#ifndef TTVCOLOR_INCLUDED
#define TTVCOLOR_INCLUDED

typedef struct {
    float x, y;
} ttvCIEXY;

typedef struct {
    float r, g, b;
} ttvRGB;

typedef struct {
    float m[9];
} ttvMatrix3x3;

// calculate the from breakpoint as K0/phi
typedef struct {
    float gamma;      // gamma of log section
    float a;          // linear bias of log section
    float K0;         // to linear break point
    float phi;        // slope of to linear section, use 1/phi if from linear
} ttvTransferOp;

typedef struct {
    ttvTransferOp toLinear;
    ttvTransferOp fromLinear;
    ttvMatrix3x3  transform;
} ttvColorTransform;

typedef struct {
    const char* name;
    ttvCIEXY red, green, blue;
    ttvCIEXY white;
    ttvColorTransform ct;
} ttvColorSpace;

#ifdef __cplusplus
extern "C" {
#endif

ttvMatrix3x3      ttvGetRGBtoCIEXYZMatrix(ttvColorSpace* cs);
ttvMatrix3x3      ttvGetCIEXYZtoRGBMatrix(ttvColorSpace* cs);
ttvColorTransform ttvGetRGBtoRGBMatrix(ttvColorSpace* src, ttvColorSpace* dst);
ttvColorSpace     ttvGetNamedColorSpace(const char* name);
ttvRGB            ttvTransformColor(ttvColorTransform* ct, ttvRGB rgb);

#ifdef __cplusplus
}
#endif
#endif

// ttvColor.c
#include <stdbool.h>
#include <string.h>

/*
acescg
adobergb
g18_ap1
g22_ap1
g18_rec709
g22_rec709
lin_adobergb
lin_ap1 (alias for acesg)
lin_displayp3
lin_rec709
lin_rec2020
lin_srgb
srgb_displayp3
srgb_texture
*/

// White point chromaticities.
static const ttvCIEXY _WpD65 = { 0.3127, 0.3290 };
static const ttvCIEXY _WpACES = { 0.32168, 0.33767 };

static const ttvColorSpace _colorSpaces[] = {
    {
        "g22_rec709",
        { 0.640, 0.330 },
        { 0.300, 0.600 },
        { 0.150, 0.060 },
        _WpD65,
        {{ 1.0/2.4, 0.55, 0.003928571429, 12.92321018 }, // toLinear - Rec709 OETF
         { 2.4, 0, 0.0, 0.0 },            // fromLinear - BT.1886 EOTF, no gain or lift
         { 0,0,0, 0,0,0, 0,0,0 }}         // transform - zero must be computed
    },
    {
        "lin_ap0",
        { 0.7347, 0.2653 },
        { 0.0000, 1.0000 },
        { 0.0001, -0.0770 },
        _WpACES,
        {{ 1.0, 0.0, 0.0, 0.0 },  // toLinear
         { 1.0, 0.0, 0.0, 0.0 },  // fromLinear
         { 0,0,0, 0,0,0, 0,0,0 }} // transform - zero must be computed
    },
    {
        "lin_ap1",
        { 0.713, 0.293 },
        { 0.165, 0.830 },
        { 0.128, 0.044 },
        _WpACES,
        {{ 1.0, 0.0, 0.0, 0.0 },  // toLinear
         { 1.0, 0.0, 0.0, 0.0 },  // fromLinear
         { 0,0,0, 0,0,0, 0,0,0 }} // transform - zero must be computed
    },
    {
        "lin_rec709",
        { 0.640, 0.330 },
        { 0.300, 0.600 },
        { 0.150, 0.060 },
        _WpD65,
        {{ 1.0, 0.0, 0.0, 0.0 },  // toLinear
         { 1.0, 0.0, 0.0, 0.0 },  // fromLinear
         { 0,0,0, 0,0,0, 0,0,0 }} // transform - zero must be computed
    },
    {
        "lin_rec2020",
        { 0.708, 0.292 },
        { 0.170, 0.797 },
        { 0.131, 0.046 },
        _WpD65,
        {{ 1.0, 0.0, 0.0, 0.0 },  // toLinear
         { 1.0, 0.0, 0.0, 0.0 },  // fromLinear
         { 0,0,0, 0,0,0, 0,0,0 }} // transform - zero must be computed
    },
    {
        "lin_srgb",
        { 0.640, 0.330 },
        { 0.300, 0.600 },
        { 0.150, 0.060 },
        _WpD65,
        {{ 1.0, 0.0, 0.0, 0.0 },  // toLinear
         { 1.0, 0.0, 0.0, 0.0 },  // fromLinear
         { 0,0,0, 0,0,0, 0,0,0 }} // transform - zero must be computed
    },
};

static bool ttvMatrix3x3IsIdentity(const ttvMatrix3x3* m) {
    return m->m[0] == 1.0 && m->m[1] == 0.0 && m->m[2] == 0.0 &&
           m->m[3] == 0.0 && m->m[4] == 1.0 && m->m[5] == 0.0 &&
           m->m[6] == 0.0 && m->m[7] == 0.0 && m->m[8] == 1.0;
}

static ttvMatrix3x3 ttvMatrix3x3Invert(ttvMatrix3x3 m) {
    ttvMatrix3x3 inv;
    float det = m.m[0] * (m.m[4] * m.m[8] - m.m[5] * m.m[7]) -
                m.m[1] * (m.m[3] * m.m[8] - m.m[5] * m.m[6]) +
                m.m[2] * (m.m[3] * m.m[7] - m.m[4] * m.m[6]);
    float invdet = 1.0 / det;
    inv.m[0] = (m.m[4] * m.m[8] - m.m[5] * m.m[7]) * invdet;
    inv.m[1] = (m.m[2] * m.m[7] - m.m[1] * m.m[8]) * invdet;
    inv.m[2] = (m.m[1] * m.m[5] - m.m[2] * m.m[4]) * invdet;
    inv.m[3] = (m.m[5] * m.m[6] - m.m[3] * m.m[8]) * invdet;
    inv.m[4] = (m.m[0] * m.m[8] - m.m[2] * m.m[6]) * invdet;
    inv.m[5] = (m.m[2] * m.m[3] - m.m[0] * m.m[5]) * invdet;
    inv.m[6] = (m.m[3] * m.m[7] - m.m[4] * m.m[6]) * invdet;
    inv.m[7] = (m.m[1] * m.m[6] - m.m[0] * m.m[7]) * invdet;
    inv.m[8] = (m.m[0] * m.m[4] - m.m[1] * m.m[3]) * invdet;
    return inv;
}

static ttvMatrix3x3 ttvMatrix3x3Multiply(ttvMatrix3x3 lh, ttvMatrix3x3 rh) {
    ttvMatrix3x3 m;
    m.m[0] = lh.m[0] * rh.m[0] + lh.m[1] * rh.m[3] + lh.m[2] * rh.m[6];
    m.m[1] = lh.m[0] * rh.m[1] + lh.m[1] * rh.m[4] + lh.m[2] * rh.m[7];
    m.m[2] = lh.m[0] * rh.m[2] + lh.m[1] * rh.m[5] + lh.m[2] * rh.m[8];
    m.m[3] = lh.m[3] * rh.m[0] + lh.m[4] * rh.m[3] + lh.m[5] * rh.m[6];
    m.m[4] = lh.m[3] * rh.m[1] + lh.m[4] * rh.m[4] + lh.m[5] * rh.m[7];
    m.m[5] = lh.m[3] * rh.m[2] + lh.m[4] * rh.m[5] + lh.m[5] * rh.m[8];
    m.m[6] = lh.m[6] * rh.m[0] + lh.m[7] * rh.m[3] + lh.m[8] * rh.m[6];
    m.m[7] = lh.m[6] * rh.m[1] + lh.m[7] * rh.m[4] + lh.m[8] * rh.m[7];
    m.m[8] = lh.m[6] * rh.m[2] + lh.m[7] * rh.m[5] + lh.m[8] * rh.m[8];
    return m;
}

ttvMatrix3x3 ttvGetRGBtoCIEXYZMatrix(ttvColorSpace* cs) {
    if (cs->ct.transform.m[8] != 0.0)
        return cs->ct.transform;

    ttvMatrix3x3 m;
    // To be consistent, let us simply use SMPTE RP 177-1993
    // compute xyz [little xyz]
    float red[3]   = { cs->red.x,   cs->red.y,   1.f - cs->red.x   - cs->red.y };
    float green[3] = { cs->green.x, cs->green.y, 1.f - cs->green.x - cs->green.y };
    float blue[3]  = { cs->blue.x,  cs->blue.y,  1.f - cs->blue.x  - cs->blue.y };
    float white[3] = { cs->white.x, cs->white.y, 1.f - cs->white.x - cs->white.y };

    // Build the P matrix by column binding red, green, and blue:
    m.m[0] = red[0];
    m.m[1] = green[0];
    m.m[2] = blue[0];
    m.m[3] = red[1];
    m.m[4] = green[1];
    m.m[5] = blue[1];
    m.m[6] = red[2];
    m.m[7] = green[2];
    m.m[8] = blue[2];

    // and W
    // white has luminance factor of 1.0, ie Y = 1
    float W[3] = { white[0] / white[1], white[1] / white[1], white[2] / white[1] };

    // compute the coefficients to scale primaries
    ttvMatrix3x3 mInv = ttvMatrix3x3Invert(m);
    float C[3] = {
        mInv.m[0] * W[0] + mInv.m[1] * W[1] + mInv.m[2] * W[2],
        mInv.m[3] * W[0] + mInv.m[4] * W[1] + mInv.m[5] * W[2],
        mInv.m[6] * W[0] + mInv.m[7] * W[1] + mInv.m[8] * W[2]
    };

    // multiply the P matrix by the diagonal matrix of coefficients
    m.m[0] *= C[0];
    m.m[1] *= C[1];
    m.m[2] *= C[2];
    m.m[3] *= C[0];
    m.m[4] *= C[1];
    m.m[5] *= C[2];
    m.m[6] *= C[0];
    m.m[7] *= C[1];
    m.m[8] *= C[2];

    // overwrite cs's transform. It's fine if two threads do it simultaneously
    // because they will both write the same value.
    cs->ct.transform = m;
    return m;
}

ttvMatrix3x3 ttvGetCIEXYZtoRGBMatrix(ttvColorSpace* cs) {
    return ttvMatrix3x3Invert(ttvGetRGBtoCIEXYZMatrix(cs));
}

ttvColorTransform ttvGetRGBtoRGBMatrix(ttvColorSpace* src, ttvColorSpace* dst) {
    ttvColorTransform t;
/*    ttvMatrix3x3 s = ttvGetRGBtoCIEXYZMatrix(src);
    ttvMatrix3x3 d = ttvGetRGBtoCIEXYZMatrix(dst);
    ttvMatrix3x3 di = ttvMatrix3x3Invert(d);
    ttvMatrix3x3 sd = ttvMatrix3x3Multiply(s, di); */
    
    t.transform = ttvMatrix3x3Multiply(ttvGetRGBtoCIEXYZMatrix(src), ttvGetCIEXYZtoRGBMatrix(dst));
    t.toLinear = dst->ct.toLinear;
    t.fromLinear = src->ct.fromLinear;
    return t;
}

void testme() {
    // this gives the correct result per Annex C
    ttvMatrix3x3 s = { 0.56711181859, 0.2793268677, 0, 0.1903210663, 0.6434664624, 0.0725032634, 0.1930166748, 0.0772066699, 1.0165544874 };
    ttvMatrix3x3 d = { 0.4123907993, 0.2126390059, 0.0193308187, 0.3575843394, 0.7151686788, 0.1191947798, 0.1804807884, 0.0721923154, 0.9505321522 };
    ttvMatrix3x3 di = ttvMatrix3x3Invert(d);
    ttvMatrix3x3 sd = ttvMatrix3x3Multiply(s, di);
}

ttvRGB ttvTransformColor(ttvColorTransform* ct, ttvRGB rgb) {
    if (!ct)
        return rgb;
    
    bool fromLinear = ct->fromLinear.gamma != 1.0;
    bool toLinear = ct->toLinear.gamma != 1.0;
    
    ttvRGB in;
    if (toLinear) {
        in.r = in.r >= ct->toLinear.K0 ?
               (1.0 + ct->toLinear.a) * powf(in.r, ct->toLinear.gamma) - ct->toLinear.a :
               in.r * ct->toLinear.phi;
        in.g = in.g >= ct->toLinear.K0 ?
               (1.0 + ct->toLinear.a) * powf(in.r, ct->toLinear.gamma) - ct->toLinear.a :
               in.g * ct->toLinear.phi;
        in.b = in.b >= ct->toLinear.K0 ?
               (1.0 + ct->toLinear.a) * powf(in.r, ct->toLinear.gamma) - ct->toLinear.a :
               in.b * ct->toLinear.phi;
    }
    else {
        in = rgb;
    }
    
    ttvRGB out;
    out.r = ct->transform.m[0] * rgb.r + ct->transform.m[1] * rgb.g + ct->transform.m[2] * rgb.b;
    out.g = ct->transform.m[3] * rgb.r + ct->transform.m[4] * rgb.g + ct->transform.m[5] * rgb.b;
    out.b = ct->transform.m[6] * rgb.r + ct->transform.m[7] * rgb.g + ct->transform.m[8] * rgb.b;
    
    if (fromLinear) {
        out.r = out.r < ct->fromLinear.K0 ?
                out.r / ct->fromLinear.phi :
                powf((out.r + ct->fromLinear.a) / (1.0 + ct->fromLinear.a), ct->fromLinear.gamma);
        out.g = out.g < ct->fromLinear.K0 ?
                out.g / ct->fromLinear.phi :
                powf((out.g + ct->fromLinear.a) / (1.0 + ct->fromLinear.a), ct->fromLinear.gamma);
        out.b = out.b < ct->fromLinear.K0 ?
                out.b / ct->fromLinear.phi :
                powf((out.b + ct->fromLinear.a) / (1.0 + ct->fromLinear.a), ct->fromLinear.gamma);
    }
    return out;
}

static ttvColorSpace ttvColorSpaceMakeEmpty() {
    ttvColorSpace cs;
    memset(&cs, 0, sizeof(cs));
    return cs;
}

ttvColorSpace ttvGetNamedColorSpace(const char* name) {
    if (!name)
        return ttvColorSpaceMakeEmpty();

    // is this a known color space?
    for (int i = 0; i < sizeof(_colorSpaces) / sizeof(_colorSpaces[0]); i++) {
        if (strcmp(name, _colorSpaces[i].name) == 0) {
            return _colorSpaces[i];
        }
    }

    return ttvColorSpaceMakeEmpty();
}