video: Add support for non-BT.709 primaries

This add support for reading primary information from lavc, categorized
into BT.601-525, BT.601-625, BT.709 and BT.2020; and passes it on to the
vo. In vo_opengl, we always generate the 3dlut against the wider BT.2020
and transform our source into this colorspace in the shader.

9 files changed, 194 insertions, 18 deletions

diff --git a/DOCS/man/vo.rst b/DOCS/man/vo.rst
index 35df26e245..dc4f49097c 100644
--- a/DOCS/man/vo.rst
+++ b/DOCS/man/vo.rst
@@ -362,6 +362,12 @@ Available video output drivers are:
         LittleCMS 2. If both ``srgb`` and ``icc-profile`` are present, the
         latter takes precedence, as they are somewhat redundant.
 
+        Note: When playing back BT.2020 content with this option enabled, out
+        of gamut colors will be numerically clipped, which can potentially
+        change the hue and/or luminance. If this is not desired, it is
+        recommended to use ``icc-profile`` with an sRGB ICC profile instead,
+        when playing back wide-gamut BT.2020 content.
+
     ``pbo``
         Enable use of PBOs. This is slightly faster, but can sometimes lead to
         sporadic and temporary image corruption (in theory, because reupload
@@ -521,7 +527,7 @@ Available video output drivers are:
     ``3dlut-size=<r>x<g>x<b>``
         Size of the 3D LUT generated from the ICC profile in each dimension.
         Default is 128x256x64.
-        Sizes must be a power of two, and 256 at most.
+        Sizes must be a power of two, and 512 at most.
 
     ``alpha=<blend|yes|no>``
         Decides what to do if the input has an alpha component (default: blend).
diff --git a/video/csputils.c b/video/csputils.c
index 13398d6d10..32461697c0 100644
--- a/video/csputils.c
+++ b/video/csputils.c
@@ -42,7 +42,7 @@ const char *const mp_csp_names[MP_CSP_COUNT] = {
     "BT.601 (SD)",
     "BT.709 (HD)",
     "SMPTE-240M",
-    "BT.2020 (NC)",
+    "BT.2020-NC (UHD)",
     "RGB",
     "XYZ",
     "YCgCo",
@@ -54,6 +54,14 @@ const char *const mp_csp_levels_names[MP_CSP_LEVELS_COUNT] = {
     "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",
@@ -93,6 +101,20 @@ enum mp_csp_levels avcol_range_to_mp_csp_levels(int avrange)
     }
 }
 
+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) {
@@ -117,6 +139,19 @@ int mp_csp_levels_to_avcol_range(enum mp_csp_levels range)
     }
 }
 
+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;
@@ -209,6 +244,53 @@ static void luma_coeffs(float m[3][4], float lr, float lg, float lb)
 }
 
 /**
+ * Get the coefficients of the source -> bt2020 gamut mapping matrix,
+ * given an identifier describing the source gamut primaries.
+ */
+void mp_get_cms_matrix(enum mp_csp_prim source, float m[3][3])
+{
+    // Conversion matrices to BT.2020 primaries
+    // These were computed using: http://lpaste.net/101796
+    switch (source) {
+        case MP_CSP_PRIM_BT_601_525: {
+            static const float from_525[3][3] = {
+                {0.5952542, 0.3493139, 0.0554319},
+                {0.0812437, 0.8915033, 0.0272530},
+                {0.0155123, 0.0819116, 0.9025760},
+            };
+            memcpy(m, from_525, sizeof(from_525));
+            break;
+        }
+        case MP_CSP_PRIM_BT_601_625: {
+            static const float from_625[3][3] = {
+                {0.6550368, 0.3021610, 0.0428023},
+                {0.0721406, 0.9166311, 0.0112283},
+                {0.0171134, 0.0978535, 0.8850332},
+            };
+            memcpy(m, from_625, sizeof(from_625));
+            break;
+        }
+        case MP_CSP_PRIM_BT_709: {
+            static const float from_709[3][3] = {
+                {0.6274039, 0.3292830, 0.0433131},
+                {0.0690973, 0.9195404, 0.0113623},
+                {0.0163914, 0.0880133, 0.8955953},
+            };
+            memcpy(m, from_709, sizeof(from_709));
+            break;
+        }
+        case MP_CSP_PRIM_BT_2020:
+        default: {
+            // No conversion necessary. This matrix should not even be used
+            // in this case, but assign it to the identity just to be safe.
+            static const float ident[3][3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
+            memcpy(m, ident, sizeof(ident));
+            break;
+        }
+    }
+}
+
+/**
  * \brief get the coefficients of the yuv -> rgb conversion matrix
  * \param params struct specifying the properties of the conversion like
  *  brightness, ...
diff --git a/video/csputils.h b/video/csputils.h
index 57e10b3c27..d7d052a753 100644
--- a/video/csputils.h
+++ b/video/csputils.h
@@ -57,6 +57,18 @@ enum mp_csp_levels {
 // Any enum mp_csp_levels value is a valid index (except MP_CSP_LEVELS_COUNT)
 extern const char *const mp_csp_levels_names[MP_CSP_LEVELS_COUNT];
 
+enum mp_csp_prim {
+    MP_CSP_PRIM_AUTO,
+    MP_CSP_PRIM_BT_601_525,
+    MP_CSP_PRIM_BT_601_625,
+    MP_CSP_PRIM_BT_709,
+    MP_CSP_PRIM_BT_2020,
+    MP_CSP_PRIM_COUNT
+};
+
+// Any enum mp_csp_prim value is a valid index (except MP_CSP_PRIM_COUNT)
+extern const char *const mp_csp_prim_names[MP_CSP_PRIM_COUNT];
+
 struct mp_csp_details {
     enum mp_csp format;
     enum mp_csp_levels levels_in;      // encoded video
@@ -141,10 +153,14 @@ enum mp_csp avcol_spc_to_mp_csp(int avcolorspace);
 
 enum mp_csp_levels avcol_range_to_mp_csp_levels(int avrange);
 
+enum mp_csp_prim avcol_pri_to_mp_csp_prim(int avpri);
+
 int mp_csp_to_avcol_spc(enum mp_csp colorspace);
 
 int mp_csp_levels_to_avcol_range(enum mp_csp_levels range);
 
+int mp_csp_prim_to_avcol_pri(enum mp_csp_prim prim);
+
 enum mp_csp mp_csp_guess_colorspace(int width, int height);
 
 enum mp_chroma_location avchroma_location_to_mp(int avloc);
@@ -160,6 +176,8 @@ void mp_gen_gamma_map(unsigned char *map, int size, float gamma);
 #define COL_U 1
 #define COL_V 2
 #define COL_C 3
+void mp_get_cms_matrix(enum mp_csp_prim source, float m[3][3]);
+
 void mp_get_yuv2rgb_coeffs(struct mp_csp_params *params, float yuv2rgb[3][4]);
 void mp_gen_yuv2rgb_map(struct mp_csp_params *params, uint8_t *map, int size);
 
diff --git a/video/decode/vd_lavc.c b/video/decode/vd_lavc.c
index 721461fdb7..a59fffa790 100644
--- a/video/decode/vd_lavc.c
+++ b/video/decode/vd_lavc.c
@@ -477,6 +477,7 @@ static void update_image_params(struct dec_video *vd, AVFrame *frame,
         .d_h = d_h,
         .colorspace = avcol_spc_to_mp_csp(ctx->avctx->colorspace),
         .colorlevels = avcol_range_to_mp_csp_levels(ctx->avctx->color_range),
+        .primaries = avcol_pri_to_mp_csp_prim(ctx->avctx->color_primaries),
         .chroma_location =
             avchroma_location_to_mp(ctx->avctx->chroma_sample_location),
         .rotate = vd->header->video->rotate,
diff --git a/video/mp_image.c b/video/mp_image.c
index 64c310aaed..1c686e5193 100644
--- a/video/mp_image.c
+++ b/video/mp_image.c
@@ -371,6 +371,7 @@ void mp_image_copy_attributes(struct mp_image *dst, struct mp_image *src)
     if ((dst->flags & MP_IMGFLAG_YUV) == (src->flags & MP_IMGFLAG_YUV)) {
         dst->params.colorspace = src->params.colorspace;
         dst->params.colorlevels = src->params.colorlevels;
+        dst->params.primaries = src->params.primaries;
         dst->params.chroma_location = src->params.chroma_location;
     }
     if ((dst->fmt.flags & MP_IMGFLAG_PAL) && (src->fmt.flags & MP_IMGFLAG_PAL)) {
@@ -486,6 +487,7 @@ bool mp_image_params_equal(const struct mp_image_params *p1,
            p1->colorspace == p2->colorspace &&
            p1->colorlevels == p2->colorlevels &&
            p1->outputlevels == p2->outputlevels &&
+           p1->primaries == p2->primaries &&
            p1->chroma_location == p2->chroma_location &&
            p1->rotate == p2->rotate;
 }
@@ -533,16 +535,43 @@ void mp_image_params_guess_csp(struct mp_image_params *params)
             params->colorspace = mp_csp_guess_colorspace(params->w, params->h);
         if (params->colorlevels == MP_CSP_LEVELS_AUTO)
             params->colorlevels = MP_CSP_LEVELS_TV;
+        if (params->primaries == MP_CSP_PRIM_AUTO) {
+            // We assume BT.709 primaries for all untagged BT.609/BT.709
+            // content, because it offers the minimal deviation from all three,
+            // including both NTSC and PAL/SECAM.
+            if (params->colorspace == MP_CSP_BT_2020_NC) {
+                params->primaries = MP_CSP_PRIM_BT_2020;
+            } else {
+                params->primaries = MP_CSP_PRIM_BT_709;
+            }
+        }
     } else if (fmt.flags & MP_IMGFLAG_RGB) {
         params->colorspace = MP_CSP_RGB;
         params->colorlevels = MP_CSP_LEVELS_PC;
+
+        // The majority of RGB content is either sRGB or (rarely) some other
+        // color space which we don't even handle, like AdobeRGB or
+        // ProPhotoRGB. The only reasonable thing we can do is assume it's
+        // sRGB and hope for the best, which should usually just work out fine.
+        // Note: sRGB primaries = BT.709 primaries
+        if (params->primaries == MP_CSP_PRIM_AUTO)
+            params->primaries = MP_CSP_PRIM_BT_709;
     } else if (fmt.flags & MP_IMGFLAG_XYZ) {
         params->colorspace = MP_CSP_XYZ;
         params->colorlevels = MP_CSP_LEVELS_PC;
+
+        // The default XYZ matrix converts it to BT.709 color space
+        // since that's the most likely scenario. Proper VOs should ignore
+        // this field as well as the matrix and treat XYZ input as absolute,
+        // but for VOs which use the matrix (and hence, consult this field)
+        // this is the correct parameter.
+        if (params->primaries == MP_CSP_PRIM_AUTO)
+            params->primaries = MP_CSP_PRIM_BT_709;
     } else {
         // We have no clue.
         params->colorspace = MP_CSP_AUTO;
         params->colorlevels = MP_CSP_LEVELS_AUTO;
+        params->primaries = MP_CSP_PRIM_AUTO;
     }
 }
 
diff --git a/video/mp_image.h b/video/mp_image.h
index 61a3d50f53..5ab12ae3d8 100644
--- a/video/mp_image.h
+++ b/video/mp_image.h
@@ -47,6 +47,7 @@ struct mp_image_params {
     int d_w, d_h;               // define display aspect ratio (never 0/0)
     enum mp_csp colorspace;
     enum mp_csp_levels colorlevels;
+    enum mp_csp_prim primaries;
     enum mp_chroma_location chroma_location;
     // The image should be converted to these levels. Unlike colorlevels, it
     // does not describe the current state of the image. (Somewhat similar to
diff --git a/video/out/gl_lcms.c b/video/out/gl_lcms.c
index 75388ada05..b1967ebf0e 100644
--- a/video/out/gl_lcms.c
+++ b/video/out/gl_lcms.c
@@ -46,7 +46,7 @@ static bool parse_3dlut_size(const char *arg, int *p1, int *p2, int *p3)
         return false;
     for (int n = 0; n < 3; n++) {
         int s = ((int[]) { *p1, *p2, *p3 })[n];
-        if (s < 2 || s > 256 || ((s - 1) & s))
+        if (s < 2 || s > 512 || ((s - 1) & s))
             return false;
     }
     return true;
@@ -135,8 +135,8 @@ struct lut3d *mp_load_icc(struct mp_icc_opts *opts, struct mp_log *log,
     char *cache_info =
         // Gamma is included in the header to help uniquely identify it,
         // because we may change the parameter in the future or make it
-        // customizable.
-        talloc_asprintf(tmp, "intent=%d, size=%dx%dx%d, gamma=2.4",
+        // customizable, same for the primaries.
+        talloc_asprintf(tmp, "intent=%d, size=%dx%dx%d, gamma=2.4, prim=bt2020\n",
                         opts->intent, s_r, s_g, s_b);
 
     // check cache
@@ -165,17 +165,19 @@ struct lut3d *mp_load_icc(struct mp_icc_opts *opts, struct mp_log *log,
     if (!profile)
         goto error_exit;
 
-    cmsCIExyY d65 = {0.3127, 0.3290, 1.0};
-    static const cmsCIExyYTRIPLE bt709prim = {
-        .Red   = {0.64, 0.33, 1.0},
-        .Green = {0.30, 0.60, 1.0},
-        .Blue  = {0.15, 0.06, 1.0},
+    // We always generate the 3DLUT against BT.2020, and transform into this
+    // space inside the shader if the source differs.
+    static const cmsCIExyY d65 = {0.3127, 0.3290, 1.0};
+    static const cmsCIExyYTRIPLE bt2020prim = {
+        .Red   = {0.708, 0.292, 1.0},
+        .Green = {0.170, 0.797, 1.0},
+        .Blue  = {0.131, 0.046, 1.0},
     };
 
     // 2.4 is arbitrarily used as a gamma compression factor for the 3DLUT,
     // reducing artifacts due to rounding errors on wide gamut profiles
     cmsToneCurve *tonecurve = cmsBuildGamma(cms, 2.4);
-    cmsHPROFILE vid_profile = cmsCreateRGBProfileTHR(cms, &d65, &bt709prim,
+    cmsHPROFILE vid_profile = cmsCreateRGBProfileTHR(cms, &d65, &bt2020prim,
                         (cmsToneCurve*[3]){tonecurve, tonecurve, tonecurve});
     cmsFreeToneCurve(tonecurve);
     cmsHTRANSFORM trafo = cmsCreateTransformTHR(cms, vid_profile, TYPE_RGB_16,
diff --git a/video/out/gl_video.c b/video/out/gl_video.c
index 9ba5c54284..e50cef3686 100644
--- a/video/out/gl_video.c
+++ b/video/out/gl_video.c
@@ -639,6 +639,13 @@ static void update_uniforms(struct gl_video *p, GLuint program)
 
     gl->Uniform1i(gl->GetUniformLocation(program, "lut_3d"), TEXUNIT_3DLUT);
 
+    loc = gl->GetUniformLocation(program, "cms_matrix");
+    if (loc >= 0) {
+        float cms_matrix[3][3] = {{0}};
+        mp_get_cms_matrix(p->image_params.primaries, cms_matrix);
+        gl->UniformMatrix3fv(loc, 1, GL_TRUE, &cms_matrix[0][0]);
+    }
+
     for (int n = 0; n < 2; n++) {
         const char *lut = p->scalers[n].lut_name;
         if (lut)
@@ -840,6 +847,9 @@ static void compile_shaders(struct gl_video *p)
     char *header = talloc_asprintf(tmp, "#version %d\n%s%s", gl->glsl_version,
                                    shader_prelude, PRELUDE_END);
 
+    bool use_cms = p->opts.srgb || p->use_lut_3d;
+    bool use_cms_matrix = use_cms && (p->image_params.primaries != MP_CSP_PRIM_BT_2020);
+
     if (p->gl_target == GL_TEXTURE_RECTANGLE) {
         shader_def(&header, "VIDEO_SAMPLER", "sampler2DRect");
         shader_def_opt(&header, "USE_RECTANGLE", true);
@@ -852,8 +862,8 @@ static void compile_shaders(struct gl_video *p)
         shader_def_opt(&header, "USE_ALPHA", p->has_alpha);
 
     char *header_osd = talloc_strdup(tmp, header);
-    shader_def_opt(&header_osd, "USE_OSD_LINEAR_CONV", p->opts.srgb ||
-                                                       p->use_lut_3d);
+    shader_def_opt(&header_osd, "USE_OSD_LINEAR_CONV", use_cms);
+    shader_def_opt(&header_osd, "USE_OSD_CMS_MATRIX", use_cms_matrix);
     shader_def_opt(&header_osd, "USE_OSD_3DLUT", p->use_lut_3d);
     // 3DLUT overrides SRGB
     shader_def_opt(&header_osd, "USE_OSD_SRGB", !p->use_lut_3d && p->opts.srgb);
@@ -888,7 +898,7 @@ static void compile_shaders(struct gl_video *p)
     // Linear light scaling is only enabled when either color correction
     // option (3dlut or srgb) is enabled, otherwise scaling is done in the
     // source space.
-    bool convert_to_linear_gamma = !p->is_linear_rgb && (p->opts.srgb || p->use_lut_3d);
+    bool convert_to_linear_gamma = !p->is_linear_rgb && use_cms;
 
     if (p->image_format == IMGFMT_NV12 || p->image_format == IMGFMT_NV21) {
         shader_def(&header_conv, "USE_CONV", "CONV_NV12");
@@ -912,6 +922,7 @@ static void compile_shaders(struct gl_video *p)
         shader_def(&header_conv, "USE_ALPHA_BLEND", "1");
 
     shader_def_opt(&header_final, "USE_GAMMA_POW", p->opts.gamma > 0);
+    shader_def_opt(&header_final, "USE_CMS_MATRIX", use_cms_matrix);
     shader_def_opt(&header_final, "USE_3DLUT", p->use_lut_3d);
     // 3DLUT overrides SRGB
     shader_def_opt(&header_final, "USE_SRGB", p->opts.srgb && !p->use_lut_3d);
diff --git a/video/out/gl_video_shaders.glsl b/video/out/gl_video_shaders.glsl
index b0058f8ca6..a6ecc2eeeb 100644
--- a/video/out/gl_video_shaders.glsl
+++ b/video/out/gl_video_shaders.glsl
@@ -56,6 +56,13 @@ vec3 bt2020_expand(vec3 v)
 }
 #endif
 
+// Constant matrix for conversion from BT.2020 to sRGB
+const mat3 srgb_matrix = mat3(
+     1.6604910, -0.1245505, -0.0181508,
+    -0.5876411,  1.1328999, -0.1005789,
+    -0.0728499, -0.0083494,  1.1187297
+);
+
 #!section vertex_all
 
 #if __VERSION__ < 130
@@ -66,6 +73,7 @@ vec3 bt2020_expand(vec3 v)
 
 uniform mat3 transform;
 uniform sampler3D lut_3d;
+uniform mat3 cms_matrix; // transformation from file's gamut to bt.2020
 
 in vec2 vertex_position;
 in vec4 vertex_color;
@@ -89,12 +97,17 @@ void main() {
     // NOTE: This always applies the true BT2020, maybe we need to use
     // approx-gamma here too?
 #endif
+#ifdef USE_OSD_CMS_MATRIX
+    // Convert to the right target gamut first (to BT.709 for sRGB,
+    // and to BT.2020 for 3DLUT).
+    color.rgb = clamp(cms_matrix * color.rgb, 0, 1);
+#endif
 #ifdef USE_OSD_3DLUT
     color.rgb = pow(color.rgb, vec3(1/2.4)); // linear -> 2.4 3DLUT space
     color = vec4(texture3D(lut_3d, color.rgb).rgb, color.a);
 #endif
 #ifdef USE_OSD_SRGB
-    color.rgb = srgb_compand(color.rgb);
+    color.rgb = srgb_compand(clamp(srgb_matrix * color.rgb, 0, 1));
 #endif
 
     texcoord = vertex_texcoord;
@@ -136,6 +149,7 @@ uniform sampler2D lut_l_2d;
 uniform sampler3D lut_3d;
 uniform sampler2D dither;
 uniform mat4x3 colormatrix;
+uniform mat3 cms_matrix;
 uniform mat2 dither_trafo;
 uniform vec3 inv_gamma;
 uniform float input_gamma;
@@ -406,16 +420,28 @@ void main() {
     // User-defined gamma correction factor (via the gamma sub-option)
     color = pow(color, inv_gamma);
 #endif
+#ifdef USE_CMS_MATRIX
+    // Convert to the right target gamut first (to BT.709 for sRGB,
+    // and to BT.2020 for 3DLUT).
+    color = cms_matrix * color;
+#endif
 #ifdef USE_3DLUT
     // For the 3DLUT we are arbitrarily using 2.4 as input gamma to reduce
     // the amount of rounding errors, so we pull up to that space first and
     // then pass it through the 3D texture.
-    color = pow(color, vec3(1/2.4));
+    //
+    // The value is clamped to [0,1] first because the gamma function is not
+    // well-defined outside it. This should not be a problem because the 3dlut
+    // is not defined for values outside its boundaries either way, and no
+    // media can possibly exceed its BT.2020 source gamut either way due to
+    // that being the biggest taggable color space. This is just to avoid
+    // numerical quirks like -1e-30 turning into NaN.
+    color = pow(clamp(color, 0, 1), vec3(1/2.4));
     color = texture3D(lut_3d, color).rgb;
 #endif
 #ifdef USE_SRGB
-    // Compand from the linear scaling gamma to the sRGB output gamma
-    color = srgb_compand(color.rgb);
+    // Adapt and compand from the linear BT2020 source to the sRGB output
+    color = srgb_compand(clamp(srgb_matrix * color, 0, 1));
 #endif
 #ifdef USE_DITHER
     vec2 dither_pos = gl_FragCoord.xy / dither_size;