vo_gpu: add BT.2390 tone-mapping

Implementation copy/pasted from:
https://code.videolan.org/videolan/libplacebo/-/commit/f793fc0480f

This brings mpv's tone mapping more in line with industry standard
practices, for a hopefully more consistent result across the board.

Note that we ignore the black point adjustment of the tone mapping
entirely. In theory we could revisit this, if we ever make black point
compensation part of the mpv rendering pipeline.

DOCS/man/options.rst

@@ -6058,6 +6058,9 @@ The following video options are currently all specific to ``--vo=gpu`` and
         color/brightness accuracy. This is roughly equivalent to
         ``--tone-mapping=reinhard --tone-mapping-param=0.24``. If possible,
         you should also enable ``--hdr-compute-peak`` for the best results.
+    bt.2390
+        Perceptual tone mapping curve (EETF) specified in ITU-R Report BT.2390.
+        This is the recommended curve to use for typical HDR-mastered content.
         (Default)
     gamma
         Fits a logarithmic transfer between the tone curves.

video/out/gpu/video.c

@@ -322,7 +322,7 @@ static const struct gl_video_opts gl_video_opts_def = {
     .background = {0, 0, 0, 255},
     .gamma = 1.0f,
     .tone_map = {
-        .curve = TONE_MAPPING_HABLE,
+        .curve = TONE_MAPPING_BT_2390,
         .curve_param = NAN,
         .max_boost = 1.0,
         .decay_rate = 100.0,
@@ -382,7 +382,8 @@ const struct m_sub_options gl_video_conf = {
             {"reinhard", TONE_MAPPING_REINHARD},
             {"hable",    TONE_MAPPING_HABLE},
             {"gamma",    TONE_MAPPING_GAMMA},
-            {"linear",   TONE_MAPPING_LINEAR})},
+            {"linear",   TONE_MAPPING_LINEAR},
+            {"bt.2390",  TONE_MAPPING_BT_2390})},
         {"hdr-compute-peak", OPT_CHOICE(tone_map.compute_peak,
             {"auto", 0},
             {"yes", 1},

video/out/gpu/video.h

@@ -94,6 +94,7 @@ enum tone_mapping {
     TONE_MAPPING_HABLE,
     TONE_MAPPING_GAMMA,
     TONE_MAPPING_LINEAR,
+    TONE_MAPPING_BT_2390,
 };
 
 struct gl_tone_map_opts {

video/out/gpu/video_shaders.c

@@ -649,6 +649,15 @@ static void hdr_update_peak(struct gl_shader_cache *sc,
     GLSL(})
 }
 
+static inline float pq_delinearize(float x)
+{
+    x *= MP_REF_WHITE / 10000.0;
+    x = powf(x, PQ_M1);
+    x = (PQ_C1 + PQ_C2 * x) / (1.0 + PQ_C3 * x);
+    x = pow(x, PQ_M2);
+    return x;
+}
+
 // Tone map from a known peak brightness to the range [0,1]. If ref_peak
 // is 0, we will use peak detection instead
 static void pass_tone_map(struct gl_shader_cache *sc,
@@ -672,12 +681,18 @@ static void pass_tone_map(struct gl_shader_cache *sc,
 
     GLSLF("vec3 sig = color.rgb;\n");
 
+    // This function always operates on an absolute scale, so ignore the
+    // dst_peak normalization for it
+    float dst_scale = dst_peak;
+    if (opts->curve == TONE_MAPPING_BT_2390)
+        dst_scale = 1.0;
+
     // Rescale the variables in order to bring it into a representation where
     // 1.0 represents the dst_peak. This is because all of the tone mapping
     // algorithms are defined in such a way that they map to the range [0.0, 1.0].
-    if (dst_peak > 1.0) {
+    if (dst_scale > 1.0) {
-        GLSLF("sig *= 1.0/%f;\n", dst_peak);
+        GLSLF("sig *= 1.0/%f;\n", dst_scale);
-        GLSLF("sig_peak *= 1.0/%f;\n", dst_peak);
+        GLSLF("sig_peak *= 1.0/%f;\n", dst_scale);
     }
 
     GLSL(float sig_orig = sig[sig_idx];)
@@ -744,6 +759,40 @@ static void pass_tone_map(struct gl_shader_cache *sc,
         break;
     }
 
+    case TONE_MAPPING_BT_2390:
+        // We first need to encode both sig and sig_peak into PQ space
+        GLSLF("vec4 sig_pq = vec4(sig.rgb, sig_peak);                           \n"
+              "sig_pq *= vec4(1.0/%f);                                          \n"
+              "sig_pq = pow(sig_pq, vec4(%f));                                  \n"
+              "sig_pq = (vec4(%f) + vec4(%f) * sig_pq)                          \n"
+              "          / (vec4(1.0) + vec4(%f) * sig_pq);                     \n"
+              "sig_pq = pow(sig_pq, vec4(%f));                                  \n",
+              10000.0 / MP_REF_WHITE, PQ_M1, PQ_C1, PQ_C2, PQ_C3, PQ_M2);
+        // Encode both the signal and the target brightness to be relative to
+        // the source peak brightness, and figure out the target peak in this space
+        GLSLF("float scale = 1.0 / sig_pq.a;                                    \n"
+              "sig_pq.rgb *= vec3(scale);                                       \n"
+              "float maxLum = %f * scale;                                       \n",
+              pq_delinearize(dst_peak));
+        // Apply piece-wise hermite spline
+        GLSLF("float ks = 1.5 * maxLum - 0.5;                                   \n"
+              "vec3 tb = (sig_pq.rgb - vec3(ks)) / vec3(1.0 - ks);              \n"
+              "vec3 tb2 = tb * tb;                                              \n"
+              "vec3 tb3 = tb2 * tb;                                             \n"
+              "vec3 pb = (2.0 * tb3 - 3.0 * tb2 + vec3(1.0)) * vec3(ks) +       \n"
+              "          (tb3 - 2.0 * tb2 + tb) * vec3(1.0 - ks) +              \n"
+              "          (-2.0 * tb3 + 3.0 * tb2) * vec3(maxLum);               \n"
+              "sig = mix(pb, sig_pq.rgb, lessThan(sig_pq.rgb, vec3(ks)));       \n");
+        // Convert back from PQ space to linear light
+        GLSLF("sig *= vec3(sig_pq.a);                                           \n"
+              "sig = pow(sig, vec3(1.0/%f));                                    \n"
+              "sig = max(sig - vec3(%f), 0.0) /                                 \n"
+              "          (vec3(%f) - vec3(%f) * sig);                           \n"
+              "sig = pow(sig, vec3(1.0/%f));                                    \n"
+              "sig *= vec3(%f);                                                 \n",
+              PQ_M2, PQ_C1, PQ_C2, PQ_C3, PQ_M1, 10000.0 / MP_REF_WHITE);
+        break;
+
     default:
         abort();
     }
@@ -754,11 +803,11 @@ static void pass_tone_map(struct gl_shader_cache *sc,
     // Mix between the per-channel tone mapped and the linear tone mapped
     // signal based on the desaturation strength
     if (opts->desat > 0) {
-        float base = 0.18 * dst_peak;
+        float base = 0.18 * dst_scale;
         GLSLF("float coeff = max(sig[sig_idx] - %f, 1e-6) / "
               "              max(sig[sig_idx], 1.0);\n", base);
         GLSLF("coeff = %f * pow(coeff, %f);\n", opts->desat, opts->desat_exp);
-        GLSLF("color.rgb = mix(sig_lin, %f * sig, coeff);\n", dst_peak);
+        GLSLF("color.rgb = mix(sig_lin, %f * sig, coeff);\n", dst_scale);
     } else {
         GLSL(color.rgb = sig_lin;)
     }