mpv/video/csputils.c

/*
 * Common code related to colorspaces and conversion
 *
 * Copyleft (C) 2009 Reimar Döffinger <Reimar.Doeffinger@gmx.de>
 *
 * mp_invert_yuv2rgb based on DarkPlaces engine, original code (GPL2 or later)
 *
 * This file is part of MPlayer.
 *
 * MPlayer is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * MPlayer is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License along
 * with MPlayer; if not, write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 *
 * You can alternatively redistribute this file and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 */

#include "config.h"

#include <stdint.h>
#include <math.h>
#include <assert.h>
#include <libavutil/common.h>
#include <libavcodec/avcodec.h>

#include "mp_image.h"
#include "csputils.h"

const char *const mp_csp_names[MP_CSP_COUNT] = {
    "Autoselect",
    "BT.601 (SD)",
    "BT.709 (HD)",
    "SMPTE-240M",
    "BT.2020-NCL (UHD)",
    "BT.2020-CL (UHD)",
    "RGB",
    "XYZ",
    "YCgCo",
};

const char *const mp_csp_levels_names[MP_CSP_LEVELS_COUNT] = {
    "Autoselect",
    "TV",
    "PC",
};

const char *const mp_csp_prim_names[MP_CSP_PRIM_COUNT] = {
    "Autoselect",
    "BT.601 (525-line SD)",
    "BT.601 (625-line SD)",
    "BT.709 (HD)",
    "BT.2020 (UHD)",
};

const char *const mp_csp_equalizer_names[MP_CSP_EQ_COUNT] = {
    "brightness",
    "contrast",
    "hue",
    "saturation",
    "gamma",
};

const char *const mp_chroma_names[MP_CHROMA_COUNT] = {
    "unknown",
    "mpeg2/4/h264",
    "mpeg1/jpeg",
};

// The short name _must_ match with what vf_stereo3d accepts (if supported).
// The long name is closer to the Matroska spec (StereoMode element).
// The numeric index matches the Matroska StereoMode value. If you add entries
// that don't match Matroska, make sure demux_mkv.c rejects them properly.
// The long name is unused.
#define E(index, short, long) [index] = short
const char *const mp_stereo3d_names[MP_STEREO3D_COUNT] = {
    E(0,  "mono",   "mono"),                    // unsupported by vf_stereo3d
    E(1,  "sbs2l",  "side_by_side_left"),
    E(2,  "ab2r",   "top_bottom_right"),
    E(3,  "ab2l",   "top_bottom_left"),
    E(4,  "checkr", "checkboard_right"),        // unsupported by vf_stereo3d
    E(5,  "checkl", "checkboard_left"),         // unsupported by vf_stereo3d
    E(6,  "irr",    "row_interleaved_right"),
    E(7,  "irl",    "row_interleaved_left"),
    E(8,  "icr",    "column_interleaved_right"),// unsupported by vf_stereo3d
    E(9,  "icl",    "column_interleaved_left"), // unsupported by vf_stereo3d
    E(10, "arcc",   "anaglyph_cyan_red"),       // Matroska: unclear which mode
    E(11, "sbs2r",  "side_by_side_right"),
    E(12, "agmc",   "anaglyph_green_magenta"),  // Matroska: unclear which mode
};

enum mp_csp avcol_spc_to_mp_csp(int avcolorspace)
{
    switch (avcolorspace) {
    case AVCOL_SPC_BT709:       return MP_CSP_BT_709;
    case AVCOL_SPC_BT470BG:     return MP_CSP_BT_601;
#if HAVE_AVCOL_SPC_BT2020
    case AVCOL_SPC_BT2020_NCL:  return MP_CSP_BT_2020_NC;
    case AVCOL_SPC_BT2020_CL:   return MP_CSP_BT_2020_C;
#endif
    case AVCOL_SPC_SMPTE170M:   return MP_CSP_BT_601;
    case AVCOL_SPC_SMPTE240M:   return MP_CSP_SMPTE_240M;
    case AVCOL_SPC_RGB:         return MP_CSP_RGB;
    case AVCOL_SPC_YCOCG:       return MP_CSP_YCGCO;
    default:                    return MP_CSP_AUTO;
    }
}

enum mp_csp_levels avcol_range_to_mp_csp_levels(int avrange)
{
    switch (avrange) {
    case AVCOL_RANGE_MPEG:      return MP_CSP_LEVELS_TV;
    case AVCOL_RANGE_JPEG:      return MP_CSP_LEVELS_PC;
    default:                    return MP_CSP_LEVELS_AUTO;
    }
}

enum mp_csp_prim avcol_pri_to_mp_csp_prim(int avpri)
{
    switch (avpri) {
    case AVCOL_PRI_SMPTE240M:   // Same as below
    case AVCOL_PRI_SMPTE170M:   return MP_CSP_PRIM_BT_601_525;
    case AVCOL_PRI_BT470BG:     return MP_CSP_PRIM_BT_601_625;
    case AVCOL_PRI_BT709:       return MP_CSP_PRIM_BT_709;
#if HAVE_AVCOL_SPC_BT2020
    case AVCOL_PRI_BT2020:      return MP_CSP_PRIM_BT_2020;
#endif
    default:                    return MP_CSP_PRIM_AUTO;
    }
}

int mp_csp_to_avcol_spc(enum mp_csp colorspace)
{
    switch (colorspace) {
    case MP_CSP_BT_709:         return AVCOL_SPC_BT709;
    case MP_CSP_BT_601:         return AVCOL_SPC_BT470BG;
#if HAVE_AVCOL_SPC_BT2020
    case MP_CSP_BT_2020_NC:     return AVCOL_SPC_BT2020_NCL;
    case MP_CSP_BT_2020_C:      return AVCOL_SPC_BT2020_CL;
#endif
    case MP_CSP_SMPTE_240M:     return AVCOL_SPC_SMPTE240M;
    case MP_CSP_RGB:            return AVCOL_SPC_RGB;
    case MP_CSP_YCGCO:          return AVCOL_SPC_YCOCG;
    default:                    return AVCOL_SPC_UNSPECIFIED;
    }
}

int mp_csp_levels_to_avcol_range(enum mp_csp_levels range)
{
    switch (range) {
    case MP_CSP_LEVELS_TV:      return AVCOL_RANGE_MPEG;
    case MP_CSP_LEVELS_PC:      return AVCOL_RANGE_JPEG;
    default:                    return AVCOL_RANGE_UNSPECIFIED;
    }
}

int mp_csp_prim_to_avcol_pri(enum mp_csp_prim prim)
{
    switch (prim) {
    case MP_CSP_PRIM_BT_601_525: return AVCOL_PRI_SMPTE170M;
    case MP_CSP_PRIM_BT_601_625: return AVCOL_PRI_BT470BG;
    case MP_CSP_PRIM_BT_709:     return AVCOL_PRI_BT709;
#if HAVE_AVCOL_SPC_BT2020
    case MP_CSP_PRIM_BT_2020:    return AVCOL_PRI_BT2020;
#endif
    default:                     return AVCOL_PRI_UNSPECIFIED;
    }
}

enum mp_csp mp_csp_guess_colorspace(int width, int height)
{
    return width >= 1280 || height > 576 ? MP_CSP_BT_709 : MP_CSP_BT_601;
}

enum mp_csp_prim mp_csp_guess_primaries(int width, int height)
{
    // HD content
    if (width >= 1280 || height > 576)
        return MP_CSP_PRIM_BT_709;

    switch (height) {
    case 576: // Typical PAL content, including anamorphic/squared
        return MP_CSP_PRIM_BT_601_625;

    case 480: // Typical NTSC content, including squared
    case 486: // NTSC Pro or anamorphic NTSC
        return MP_CSP_PRIM_BT_601_525;

    default: // No good metric, just pick BT.709 to minimize damage
        return MP_CSP_PRIM_BT_709;
    }
}

enum mp_chroma_location avchroma_location_to_mp(int avloc)
{
    switch (avloc) {
    case AVCHROMA_LOC_LEFT:             return MP_CHROMA_LEFT;
    case AVCHROMA_LOC_CENTER:           return MP_CHROMA_CENTER;
    default:                            return MP_CHROMA_AUTO;
    }
}

int mp_chroma_location_to_av(enum mp_chroma_location mploc)
{
    switch (mploc) {
    case MP_CHROMA_LEFT:                return AVCHROMA_LOC_LEFT;
    case MP_CHROMA_CENTER:              return AVCHROMA_LOC_CENTER;
    default:                            return AVCHROMA_LOC_UNSPECIFIED;
    }
}

// Return location of chroma samples relative to luma samples. 0/0 means
// centered. Other possible values are -1 (top/left) and +1 (right/bottom).
void mp_get_chroma_location(enum mp_chroma_location loc, int *x, int *y)
{
    *x = 0;
    *y = 0;
    if (loc == MP_CHROMA_LEFT)
        *x = -1;
}

void mp_gen_gamma_map(uint8_t *map, int size, float gamma)
{
    if (gamma == 1.0) {
        for (int i = 0; i < size; i++)
            map[i] = 255 * i / (size - 1);
        return;
    }
    gamma = 1.0 / gamma;
    for (int i = 0; i < size; i++) {
        float tmp = (float)i / (size - 1.0);
        tmp = pow(tmp, gamma);
        if (tmp > 1.0)
            tmp = 1.0;
        if (tmp < 0.0)
            tmp = 0.0;
        map[i] = 255 * tmp;
    }
}

void mp_invert_matrix3x3(float m[3][3])
{
    float m00 = m[0][0], m01 = m[0][1], m02 = m[0][2],
          m10 = m[1][0], m11 = m[1][1], m12 = m[1][2],
          m20 = m[2][0], m21 = m[2][1], m22 = m[2][2];

    // calculate the adjoint
    m[0][0] =  (m11 * m22 - m21 * m12);
    m[0][1] = -(m01 * m22 - m21 * m02);
    m[0][2] =  (m01 * m12 - m11 * m02);
    m[1][0] = -(m10 * m22 - m20 * m12);
    m[1][1] =  (m00 * m22 - m20 * m02);
    m[1][2] = -(m00 * m12 - m10 * m02);
    m[2][0] =  (m10 * m21 - m20 * m11);
    m[2][1] = -(m00 * m21 - m20 * m01);
    m[2][2] =  (m00 * m11 - m10 * m01);

    // calculate the determinant (as inverse == 1/det * adjoint,
    // adjoint * m == identity * det, so this calculates the det)
    float det = m00 * m[0][0] + m10 * m[0][1] + m20 * m[0][2];
    det = 1.0f / det;

    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++)
            m[i][j] *= det;
    }
}

// A := A * B
void mp_mul_matrix3x3(float a[3][3], float b[3][3])
{
    float a00 = a[0][0], a01 = a[0][1], a02 = a[0][2],
          a10 = a[1][0], a11 = a[1][1], a12 = a[1][2],
          a20 = a[2][0], a21 = a[2][1], a22 = a[2][2];

    for (int i = 0; i < 3; i++) {
        a[0][i] = a00 * b[0][i] + a01 * b[1][i] + a02 * b[2][i];
        a[1][i] = a10 * b[0][i] + a11 * b[1][i] + a12 * b[2][i];
        a[2][i] = a20 * b[0][i] + a21 * b[1][i] + a22 * b[2][i];
    }
}

// return the primaries associated with a certain mp_csp_primaries val
struct mp_csp_primaries mp_get_csp_primaries(enum mp_csp_prim spc)
{
    /*
    Values from: ITU-R Recommendations BT.601-7, BT.709-5, BT.2020-0

    https://www.itu.int/dms_pubrec/itu-r/rec/bt/R-REC-BT.601-7-201103-I!!PDF-E.pdf
    https://www.itu.int/dms_pubrec/itu-r/rec/bt/R-REC-BT.709-5-200204-I!!PDF-E.pdf
    https://www.itu.int/dms_pubrec/itu-r/rec/bt/R-REC-BT.2020-0-201208-I!!PDF-E.pdf
    */

    static const struct mp_csp_col_xy d65 = {0.3127, 0.3290};

    switch (spc) {
    case MP_CSP_PRIM_BT_601_525:
        return (struct mp_csp_primaries) {
            .red   = {0.630, 0.340},
            .green = {0.310, 0.595},
            .blue  = {0.155, 0.070},
            .white = d65
        };
    case MP_CSP_PRIM_BT_601_625:
        return (struct mp_csp_primaries) {
            .red   = {0.640, 0.330},
            .green = {0.290, 0.600},
            .blue  = {0.150, 0.060},
            .white = d65
        };
    // This is the default assumption if no colorspace information could
    // be determined, eg. for files which have no video channel.
    case MP_CSP_PRIM_AUTO:
    case MP_CSP_PRIM_BT_709:
        return (struct mp_csp_primaries) {
            .red   = {0.640, 0.330},
            .green = {0.300, 0.600},
            .blue  = {0.150, 0.060},
            .white = d65
        };
    case MP_CSP_PRIM_BT_2020:
        return (struct mp_csp_primaries) {
            .red   = {0.708, 0.292},
            .green = {0.170, 0.797},
            .blue  = {0.131, 0.046},
            .white = d65
        };
    default:
        return (struct mp_csp_primaries) {{0}};
    }
}

// Compute the RGB/XYZ matrix as described here:
// http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
void mp_get_rgb2xyz_matrix(struct mp_csp_primaries space, float m[3][3])
{
    float S[3], X[4], Z[4];

    // Convert from CIE xyY to XYZ. Note that Y=1 holds true for all primaries
    X[0] = space.red.x   / space.red.y;
    X[1] = space.green.x / space.green.y;
    X[2] = space.blue.x  / space.blue.y;
    X[3] = space.white.x / space.white.y;

    Z[0] = (1 - space.red.x   - space.red.y)   / space.red.y;
    Z[1] = (1 - space.green.x - space.green.y) / space.green.y;
    Z[2] = (1 - space.blue.x  - space.blue.y)  / space.blue.y;
    Z[3] = (1 - space.white.x - space.white.y) / space.white.y;

    // S = XYZ^-1 * W
    for (int i = 0; i < 3; i++) {
        m[0][i] = X[i];
        m[1][i] = 1;
        m[2][i] = Z[i];
    }

    mp_invert_matrix3x3(m);

    for (int i = 0; i < 3; i++)
        S[i] = m[i][0] * X[3] + m[i][1] * 1 + m[i][2] * Z[3];

    // M = [Sc * XYZc]
    for (int i = 0; i < 3; i++) {
        m[0][i] = S[i] * X[i];
        m[1][i] = S[i] * 1;
        m[2][i] = S[i] * Z[i];
    }
}

// M := M * XYZd<-XYZs
void mp_apply_chromatic_adaptation(struct mp_csp_col_xy src,
                                   struct mp_csp_col_xy dest, float m[3][3])
{
    // If the white points are nearly identical, this is a wasteful identity
    // operation.
    if (fabs(src.x - dest.x) < 1e-6 && fabs(src.y - dest.y) < 1e-6)
        return;

    // XYZd<-XYZs = Ma^-1 * (I*[Cd/Cs]) * Ma
    // http://www.brucelindbloom.com/index.html?Eqn_ChromAdapt.html
    float C[3][2], tmp[3][3] = {{0}};

    // Ma = Bradford matrix, arguably most popular method in use today.
    // This is derived experimentally and thus hard-coded.
    float bradford[3][3] = {
        {  0.8951,  0.2664, -0.1614 },
        { -0.7502,  1.7135,  0.0367 },
        {  0.0389, -0.0685,  1.0296 },
    };

    for (int i = 0; i < 3; i++) {
        // source cone
        C[i][0] = bradford[i][0] * src.x / src.y
                + bradford[i][1] * 1
                + bradford[i][2] * (1 - src.x - src.y) / src.y;

        // dest cone
        C[i][1] = bradford[i][0] * dest.x / dest.y
                + bradford[i][1] * 1
                + bradford[i][2] * (1 - dest.x - dest.y) / dest.y;
    }

    // tmp := I * [Cd/Cs] * Ma
    for (int i = 0; i < 3; i++)
        tmp[i][i] = C[i][1] / C[i][0];

    mp_mul_matrix3x3(tmp, bradford);

    // M := M * Ma^-1 * tmp
    mp_invert_matrix3x3(bradford);
    mp_mul_matrix3x3(m, bradford);
    mp_mul_matrix3x3(m, tmp);
}

// get the coefficients of the source -> bt2020 cms matrix
void mp_get_cms_matrix(struct mp_csp_primaries src, struct mp_csp_primaries dest,
                       enum mp_render_intent intent, float m[3][3])
{
    float tmp[3][3];

    // In saturation mapping, we don't care about accuracy and just want
    // primaries to map to primaries, making this an identity transformation.
    if (intent == MP_INTENT_SATURATION) {
        for (int i = 0; i < 3; i++)
            m[i][i] = 1;
        return;
    }

    // RGBd<-RGBs = RGBd<-XYZd * XYZd<-XYZs * XYZs<-RGBs
    // Equations from: http://www.brucelindbloom.com/index.html?Math.html
    // Note: Perceptual is treated like relative colorimetric. There's no
    // definition for perceptual other than "make it look good".

    // RGBd<-XYZd, inverted from XYZd<-RGBd
    mp_get_rgb2xyz_matrix(dest, m);
    mp_invert_matrix3x3(m);

    // Chromatic adaptation, except in absolute colorimetric intent
    if (intent != MP_INTENT_ABSOLUTE_COLORIMETRIC)
        mp_apply_chromatic_adaptation(src.white, dest.white, m);

    // XYZs<-RGBs
    mp_get_rgb2xyz_matrix(src, tmp);
    mp_mul_matrix3x3(m, tmp);
}

// get the coefficients of an SMPTE 428-1 xyz -> rgb conversion matrix
// intent = the rendering intent used to convert to the target primaries
void mp_get_xyz2rgb_coeffs(struct mp_csp_params *params,
                           struct mp_csp_primaries prim,
                           enum mp_render_intent intent, struct mp_cmat *m)
{
    float brightness = params->brightness;
    mp_get_rgb2xyz_matrix(prim, m->m);
    mp_invert_matrix3x3(m->m);

    // All non-absolute mappings want to map source white to target white
    if (intent != MP_INTENT_ABSOLUTE_COLORIMETRIC) {
        // SMPTE 428-1 defines the calibration white point as CIE xy (0.314, 0.351)
        static const struct mp_csp_col_xy smpte428 = {0.314, 0.351};
        mp_apply_chromatic_adaptation(smpte428, prim.white, m->m);
    }

    // Since this outputs linear RGB rather than companded RGB, we
    // want to linearize any brightness additions. 2 is a reasonable
    // approximation for any sort of gamma function that could be in use.
    // As this is an aesthetic setting only, any exact values do not matter.
    if (brightness < 0) {
        brightness *= -brightness;
    } else {
        brightness *= brightness;
    }

    for (int i = 0; i < 3; i++)
        m->c[i] = brightness;
}

/* Fill in the Y, U, V vectors of a yuv2rgb conversion matrix
 * based on the given luma weights of the R, G and B components (lr, lg, lb).
 * lr+lg+lb is assumed to equal 1.
 * This function is meant for colorspaces satisfying the following
 * conditions (which are true for common YUV colorspaces):
 * - The mapping from input [Y, U, V] to output [R, G, B] is linear.
 * - Y is the vector [1, 1, 1].  (meaning input Y component maps to 1R+1G+1B)
 * - U maps to a value with zero R and positive B ([0, x, y], y > 0;
 *   i.e. blue and green only).
 * - V maps to a value with zero B and positive R ([x, y, 0], x > 0;
 *   i.e. red and green only).
 * - U and V are orthogonal to the luma vector [lr, lg, lb].
 * - The magnitudes of the vectors U and V are the minimal ones for which
 *   the image of the set Y=[0...1],U=[-0.5...0.5],V=[-0.5...0.5] under the
 *   conversion function will cover the set R=[0...1],G=[0...1],B=[0...1]
 *   (the resulting matrix can be converted for other input/output ranges
 *   outside this function).
 * Under these conditions the given parameters lr, lg, lb uniquely
 * determine the mapping of Y, U, V to R, G, B.
 */
static void luma_coeffs(struct mp_cmat *mat, float lr, float lg, float lb)
{
    assert(fabs(lr+lg+lb - 1) < 1e-6);
    *mat = (struct mp_cmat) {
        { {1, 0,                    2 * (1-lr)          },
          {1, -2 * (1-lb) * lb/lg, -2 * (1-lr) * lr/lg  },
          {1,  2 * (1-lb),          0                   } },
        // Constant coefficients (mat->c) not set here
    };
}

// get the coefficients of the yuv -> rgb conversion matrix
void mp_get_yuv2rgb_coeffs(struct mp_csp_params *params, struct mp_cmat *m)
{
    int colorspace = params->colorspace;
    if (colorspace <= MP_CSP_AUTO || colorspace >= MP_CSP_COUNT)
        colorspace = MP_CSP_BT_601;
    int levels_in = params->levels_in;
    if (levels_in <= MP_CSP_LEVELS_AUTO || levels_in >= MP_CSP_LEVELS_COUNT)
        levels_in = MP_CSP_LEVELS_TV;

    switch (colorspace) {
    case MP_CSP_BT_601:     luma_coeffs(m, 0.299,  0.587,  0.114 ); break;
    case MP_CSP_BT_709:     luma_coeffs(m, 0.2126, 0.7152, 0.0722); break;
    case MP_CSP_SMPTE_240M: luma_coeffs(m, 0.2122, 0.7013, 0.0865); break;
    case MP_CSP_BT_2020_NC: luma_coeffs(m, 0.2627, 0.6780, 0.0593); break;
    case MP_CSP_BT_2020_C: {
        // Note: This outputs into the [-0.5,0.5] range for chroma information.
        // If this clips on any VO, a constant 0.5 coefficient can be added
        // to the chroma channels to normalize them into [0,1]. This is not
        // currently needed by anything, though.
        *m = (struct mp_cmat){{{0, 0, 1}, {1, 0, 0}, {0, 1, 0}}};
        break;
    }
    case MP_CSP_RGB: {
        *m = (struct mp_cmat){{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}};
        levels_in = -1;
        break;
    }
    case MP_CSP_XYZ: {
        // The vo should probably not be using a matrix generated by this
        // function for XYZ sources, but if it does, let's just assume it
        // wants BT.709 with D65 white point (virtually all other content).
        mp_get_xyz2rgb_coeffs(params, mp_get_csp_primaries(MP_CSP_PRIM_BT_709),
                              MP_INTENT_RELATIVE_COLORIMETRIC, m);
        levels_in = -1;
        break;
    }
    case MP_CSP_YCGCO: {
        *m = (struct mp_cmat) {
            {{1,  -1,  1},
             {1,   1,  0},
             {1,  -1, -1}},
        };
        break;
    }
    default:
        abort();
    };

    // Hue is equivalent to rotating input [U, V] subvector around the origin.
    // Saturation scales [U, V].
    float huecos = params->saturation * cos(params->hue);
    float huesin = params->saturation * sin(params->hue);
    for (int i = 0; i < 3; i++) {
        float u = m->m[i][1], v = m->m[i][2];
        m->m[i][1] = huecos * u - huesin * v;
        m->m[i][2] = huesin * u + huecos * v;
    }

    assert(params->input_bits >= 8);
    assert(params->texture_bits >= params->input_bits);
    double s = (1 << (params->input_bits-8)) / ((1<<params->texture_bits)-1.);
    // The values below are written in 0-255 scale
    struct yuvlevels { double ymin, ymax, cmin, cmid; }
        yuvlim =  { 16*s, 235*s, 16*s, 128*s },
        yuvfull = {  0*s, 255*s,  1*s, 128*s },  // '1' for symmetry around 128
        anyfull = {  0*s, 255*s, -255*s/2, 0 },
        yuvlev;
    switch (levels_in) {
    case MP_CSP_LEVELS_TV: yuvlev = yuvlim; break;
    case MP_CSP_LEVELS_PC: yuvlev = yuvfull; break;
    case -1: yuvlev = anyfull; break;
    default:
        abort();
    }

    int levels_out = params->levels_out;
    if (levels_out <= MP_CSP_LEVELS_AUTO || levels_out >= MP_CSP_LEVELS_COUNT)
        levels_out = MP_CSP_LEVELS_PC;
    struct rgblevels { double min, max; }
        rgblim =  { 16/255., 235/255. },
        rgbfull = {      0,        1  },
        rgblev;
    switch (levels_out) {
    case MP_CSP_LEVELS_TV: rgblev = rgblim; break;
    case MP_CSP_LEVELS_PC: rgblev = rgbfull; break;
    default:
        abort();
    }

    double ymul = (rgblev.max - rgblev.min) / (yuvlev.ymax - yuvlev.ymin);
    double cmul = (rgblev.max - rgblev.min) / (yuvlev.cmid - yuvlev.cmin) / 2;
    for (int i = 0; i < 3; i++) {
        m->m[i][0] *= ymul;
        m->m[i][1] *= cmul;
        m->m[i][2] *= cmul;
        // Set c so that Y=umin,UV=cmid maps to RGB=min (black to black)
        m->c[i] = rgblev.min - m->m[i][0] * yuvlev.ymin
                  -(m->m[i][1] + m->m[i][2]) * yuvlev.cmid;
    }

    // Brightness adds a constant to output R,G,B.
    // Contrast scales Y around 1/2 (not 0 in this implementation).
    for (int i = 0; i < 3; i++) {
        m->c[i] += params->brightness;
        m->m[i][0] *= params->contrast;
        m->c[i] += (rgblev.max-rgblev.min) * (1 - params->contrast)/2;
    }

    int in_bits = FFMAX(params->int_bits_in, 1);
    int out_bits = FFMAX(params->int_bits_out, 1);
    double in_scale = (1 << in_bits) - 1.0;
    double out_scale = (1 << out_bits) - 1.0;
    for (int i = 0; i < 3; i++) {
        m->c[i] *= out_scale; // constant is 1.0
        for (int x = 0; x < 3; x++)
            m->m[i][x] *= out_scale / in_scale;
    }
}

// size of gamma map use to avoid slow exp function in gen_yuv2rgb_map
#define GMAP_SIZE (1024)
// generate a 3D YUV -> RGB map
// map must provide space for (size + 2)^3 elements
void mp_gen_yuv2rgb_map(struct mp_csp_params *params, unsigned char *map, int size)
{
    int i, j, k, l;
    float step = 1.0 / size;
    float y, u, v;
    struct mp_cmat yuv2rgb;
    unsigned char gmaps[3][GMAP_SIZE];
    mp_gen_gamma_map(gmaps[0], GMAP_SIZE, params->rgamma);
    mp_gen_gamma_map(gmaps[1], GMAP_SIZE, params->ggamma);
    mp_gen_gamma_map(gmaps[2], GMAP_SIZE, params->bgamma);
    mp_get_yuv2rgb_coeffs(params, &yuv2rgb);
    for (i = 0; i < 3; i++) {
        for (j = 0; j < 3; j++)
            yuv2rgb.m[i][j] *= GMAP_SIZE - 1;
        yuv2rgb.c[i] *= GMAP_SIZE - 1;
    }
    v = 0;
    for (i = -1; i <= size; i++) {
        u = 0;
        for (j = -1; j <= size; j++) {
            y = 0;
            for (k = -1; k <= size; k++) {
                for (l = 0; l < 3; l++) {
                    float rgb = yuv2rgb.m[l][0] * y + yuv2rgb.m[l][1] * u +
                                yuv2rgb.m[l][2] * v + yuv2rgb.c[l];
                    *map++ = gmaps[l][av_clip(rgb, 0, GMAP_SIZE - 1)];
                }
                y += (k == -1 || k == size - 1) ? step / 2 : step;
            }
            u += (j == -1 || j == size - 1) ? step / 2 : step;
        }
        v += (i == -1 || i == size - 1) ? step / 2 : step;
    }
}

// Set colorspace related fields in p from f. Don't touch other fields.
void mp_csp_set_image_params(struct mp_csp_params *params,
                             const struct mp_image_params *imgparams)
{
    struct mp_image_params p = *imgparams;
    mp_image_params_guess_csp(&p); // ensure consistency
    params->colorspace = p.colorspace;
    params->levels_in = p.colorlevels;
    params->levels_out = p.outputlevels;
}

// Copy settings from eq into params.
void mp_csp_copy_equalizer_values(struct mp_csp_params *params,
                                  const struct mp_csp_equalizer *eq)
{
    params->brightness = eq->values[MP_CSP_EQ_BRIGHTNESS] / 100.0;
    params->contrast = (eq->values[MP_CSP_EQ_CONTRAST] + 100) / 100.0;
    params->hue = eq->values[MP_CSP_EQ_HUE] / 100.0 * 3.1415927;
    params->saturation = (eq->values[MP_CSP_EQ_SATURATION] + 100) / 100.0;
    float gamma = exp(log(8.0) * eq->values[MP_CSP_EQ_GAMMA] / 100.0);
    params->rgamma = gamma;
    params->ggamma = gamma;
    params->bgamma = gamma;
}

static int find_eq(int capabilities, const char *name)
{
    for (int i = 0; i < MP_CSP_EQ_COUNT; i++) {
        if (strcmp(name, mp_csp_equalizer_names[i]) == 0)
            return ((1 << i) & capabilities) ? i : -1;
    }
    return -1;
}

int mp_csp_equalizer_get(struct mp_csp_equalizer *eq, const char *property,
                         int *out_value)
{
    int index = find_eq(eq->capabilities, property);
    if (index < 0)
        return -1;

    *out_value = eq->values[index];

    return 0;
}

int mp_csp_equalizer_set(struct mp_csp_equalizer *eq, const char *property,
                         int value)
{
    int index = find_eq(eq->capabilities, property);
    if (index < 0)
        return 0;

    eq->values[index] = value;

    return 1;
}

void mp_invert_yuv2rgb(struct mp_cmat *out, struct mp_cmat *in)
{
    *out = *in;
    mp_invert_matrix3x3(out->m);

    // fix the constant coefficient
    // rgb = M * yuv + C
    // M^-1 * rgb = yuv + M^-1 * C
    // yuv = M^-1 * rgb - M^-1 * C
    //                  ^^^^^^^^^^
    out->c[0] = -(out->m[0][0] * in->c[0] + out->m[0][1] * in->c[1] + out->m[0][2] * in->c[2]);
    out->c[1] = -(out->m[1][0] * in->c[0] + out->m[1][1] * in->c[1] + out->m[1][2] * in->c[2]);
    out->c[2] = -(out->m[2][0] * in->c[0] + out->m[2][1] * in->c[1] + out->m[2][2] * in->c[2]);
}

// Multiply the color in c with the given matrix.
// c is {R, G, B} or {Y, U, V} (depending on input/output and matrix).
// Output is clipped to the given number of bits.
void mp_map_int_color(struct mp_cmat *matrix, int clip_bits, int c[3])
{
    int in[3] = {c[0], c[1], c[2]};
    for (int i = 0; i < 3; i++) {
        double val = matrix->c[i];
        for (int x = 0; x < 3; x++)
            val += matrix->m[i][x] * in[x];
        int ival = lrint(val);
        c[i] = av_clip(ival, 0, (1 << clip_bits) - 1);
    }
}