From d63eeb1ecc204be189d5958f1de5b4d89ee72db4 Mon Sep 17 00:00:00 2001
From: Niklas Haas <git@haasn.xyz>
Date: Fri, 19 Jun 2020 07:43:51 +0200
Subject: [PATCH] shaders/colorspace: introduce colorimetric gamut clipping

This change adds some smarter logic than the dumb per-channel clip. We
first get rid of negative (e.g. super-red) components by mixing in
enough gray until the result is no longer oversaturated, and then we
linearly scale all components down to get rid of too-bright values.

I experimented briefly with using e.g. BT.2390 tone-mapping instead of a
linear scale for the second step, but the results were not significantly
better enough to justify the performance hit of BT.2390.

I tested this by inspecting the resulting gamut-clipped values on a Yxy
chromaticity diagram and confirmed that the resulting color was indeed
at the same angle (relative to the white point) as the original, but
contained within the gamut (with Y preserved).

Subjective comparisons using a wide gamut monitor between the original
(unclipped) source and the clipped BT.709 versions of various test clips
seem visually pleasing, so I enabled it by default.

Fixes https://code.videolan.org/videolan/libplacebo/-/issues/101
---
 meson.build                                 |  2 +-
 src/colorspace.c                            | 30 +++++++++++++++++++++
 src/include/libplacebo/colorspace.h         |  5 ++++
 src/include/libplacebo/shaders/colorspace.h |  7 +++++
 src/shaders/colorspace.c                    | 29 +++++++++++++++++---
 5 files changed, 69 insertions(+), 4 deletions(-)

diff --git a/meson.build b/meson.build
index c700aedf7..dcfd58311 100644
--- a/meson.build
+++ b/meson.build
@@ -2,7 +2,7 @@ project('libplacebo', ['c', 'cpp'],
   license: 'LGPL2.1+',
   default_options: ['c_std=c99', 'cpp_std=c++11', 'warning_level=2'],
   meson_version: '>=0.49',
-  version: '2.79.0',
+  version: '2.80.0',
 )
 
 # Version number
diff --git a/src/colorspace.c b/src/colorspace.c
index e6bcce55d..25bf30278 100644
--- a/src/colorspace.c
+++ b/src/colorspace.c
@@ -782,6 +782,36 @@ struct pl_matrix3x3 pl_get_color_mapping_matrix(const struct pl_raw_primaries *s
     return xyz2rgb_d;
 }
 
+// Test the sign of 'p' relative to the line 'ab' (barycentric coordinates)
+static float test_point_line(const struct pl_cie_xy p,
+                             const struct pl_cie_xy a,
+                             const struct pl_cie_xy b)
+{
+    return (p.x - b.x) * (a.y - b.y) - (a.x - b.x) * (p.y - b.y);
+}
+
+// Test if a point is entirely inside a gamut
+static float test_point_gamut(struct pl_cie_xy point,
+                              const struct pl_raw_primaries *prim)
+{
+    float d1 = test_point_line(point, prim->red, prim->green),
+          d2 = test_point_line(point, prim->green, prim->blue),
+          d3 = test_point_line(point, prim->blue, prim->red);
+
+    bool has_neg = d1 < 0 || d2 < 0 || d3 < 0,
+         has_pos = d1 > 0 || d2 > 0 || d3 > 0;
+
+    return !(has_neg && has_pos);
+}
+
+bool pl_primaries_superset(const struct pl_raw_primaries *a,
+                           const struct pl_raw_primaries *b)
+{
+    return test_point_gamut(b->red, a) &&
+           test_point_gamut(b->green, a) &&
+           test_point_gamut(b->blue, a);
+}
+
 /* Fill in the Y, U, V vectors of a yuv-to-rgb 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.
diff --git a/src/include/libplacebo/colorspace.h b/src/include/libplacebo/colorspace.h
index 49671f577..cf8e5a1f2 100644
--- a/src/include/libplacebo/colorspace.h
+++ b/src/include/libplacebo/colorspace.h
@@ -385,6 +385,11 @@ struct pl_matrix3x3 pl_get_color_mapping_matrix(const struct pl_raw_primaries *s
                                                 const struct pl_raw_primaries *dst,
                                                 enum pl_rendering_intent intent);
 
+// Returns true if 'b' is entirely contained in 'a'. Useful for figuring out if
+// colorimetric clipping will occur or not.
+bool pl_primaries_superset(const struct pl_raw_primaries *a,
+                           const struct pl_raw_primaries *b);
+
 // Cone types involved in human vision
 enum pl_cone {
     PL_CONE_L = 1 << 0,
diff --git a/src/include/libplacebo/shaders/colorspace.h b/src/include/libplacebo/shaders/colorspace.h
index 82dd29c25..bc43134fb 100644
--- a/src/include/libplacebo/shaders/colorspace.h
+++ b/src/include/libplacebo/shaders/colorspace.h
@@ -250,6 +250,13 @@ struct pl_color_map_params {
     // all out-of-gamut colors (by inverting them), if they would have been
     // clipped as a result of gamut or tone mapping.
     bool gamut_warning;
+
+    // If true, enables colorimetric clipping. This will colorimetrically clip
+    // out-of-gamut colors by desaturating them until they hit the boundary of
+    // the permissible color volume, rather than by hard-clipping. This mode of
+    // clipping preserves luminance between the source and the destination, at
+    // the cost of introducing some color distortion in the opposite direction.
+    bool gamut_clipping;
 };
 
 extern const struct pl_color_map_params pl_color_map_default_params;
diff --git a/src/shaders/colorspace.c b/src/shaders/colorspace.c
index 4b8f67ad8..d592533b7 100644
--- a/src/shaders/colorspace.c
+++ b/src/shaders/colorspace.c
@@ -813,6 +813,7 @@ const struct pl_color_map_params pl_color_map_default_params = {
     .desaturation_strength  = 0.75,
     .desaturation_exponent  = 1.50,
     .desaturation_base      = 0.18,
+    .gamut_clipping         = true,
 };
 
 static void pl_shader_tone_map(struct pl_shader *sh, struct pl_color_space src,
@@ -1061,7 +1062,7 @@ void pl_shader_color_map(struct pl_shader *sh,
                        src.sig_avg != dst.sig_avg ||
                        src.sig_scale != dst.sig_scale ||
                        src.light != dst.light;
-
+    bool need_gamut_warn = false;
     bool is_linear = prelinearized;
     if (need_linear && !is_linear) {
         pl_shader_linearize(sh, src.transfer);
@@ -1072,8 +1073,10 @@ void pl_shader_color_map(struct pl_shader *sh,
         pl_shader_ootf(sh, src);
 
     // Tone map to rescale the signal average/peak if needed
-    if (src.sig_peak * src.sig_scale > dst.sig_peak * dst.sig_scale + 1e-6)
+    if (src.sig_peak * src.sig_scale > dst.sig_peak * dst.sig_scale + 1e-6) {
         pl_shader_tone_map(sh, src, dst, peak_detect_state, params);
+        need_gamut_warn = true;
+    }
 
     // Adapt to the right colorspace (primaries) if necessary
     if (src.primaries != dst.primaries) {
@@ -1082,14 +1085,34 @@ void pl_shader_color_map(struct pl_shader *sh,
         csp_dst = pl_raw_primaries_get(dst.primaries);
         struct pl_matrix3x3 cms_mat;
         cms_mat = pl_get_color_mapping_matrix(csp_src, csp_dst, params->intent);
+
         GLSL("color.rgb = %s * color.rgb;\n", sh_var(sh, (struct pl_shader_var) {
             .var = pl_var_mat3("cms_matrix"),
             .data = PL_TRANSPOSE_3X3(cms_mat.m),
         }));
+
+        if (!pl_primaries_superset(csp_dst, csp_src)) {
+            if (params->gamut_clipping) {
+                GLSL("float cmin = min(min(color.r, color.g), color.b);     \n"
+                     "if (cmin < 0.0) {                                     \n"
+                     "    float luma = dot(%s, color.rgb);                  \n"
+                     "    float coeff = cmin / (cmin - luma);               \n"
+                     "    color.rgb = mix(color.rgb, vec3(luma), coeff);    \n"
+                     "}                                                     \n"
+                     "float cmax = max(max(color.r, color.g), color.b);     \n"
+                     "if (cmax > 1.0)                                       \n"
+                     "    color.rgb /= cmax;                                \n"
+                     ,
+                     sh_luma_coeffs(sh, dst.primaries));
+
+            } else {
+                need_gamut_warn = true;
+            }
+        }
     }
 
     // Warn for remaining out-of-gamut colors if enabled
-    if (params->gamut_warning) {
+    if (params->gamut_warning && need_gamut_warn) {
         GLSL("if (any(greaterThan(color.rgb, vec3(%f + 0.005))) ||\n"
              "    any(lessThan(color.rgb, vec3(-0.005))))\n"
              "    color.rgb = vec3(%f) - color.rgb; // invert\n",
-- 
GitLab