LCOV - code coverage report
Current view: top level - src - tone_mapping.c (source / functions) Hit Total Coverage
Test: Code coverage Lines: 272 325 83.7 %
Date: 2025-03-29 09:04:10 Functions: 23 25 92.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 97 178 54.5 %

           Branch data     Line data    Source code
       1                 :            : /*
       2                 :            :  * This file is part of libplacebo.
       3                 :            :  *
       4                 :            :  * libplacebo is free software; you can redistribute it and/or
       5                 :            :  * modify it under the terms of the GNU Lesser General Public
       6                 :            :  * License as published by the Free Software Foundation; either
       7                 :            :  * version 2.1 of the License, or (at your option) any later version.
       8                 :            :  *
       9                 :            :  * libplacebo is distributed in the hope that it will be useful,
      10                 :            :  * but WITHOUT ANY WARRANTY; without even the implied warranty of
      11                 :            :  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      12                 :            :  * GNU Lesser General Public License for more details.
      13                 :            :  *
      14                 :            :  * You should have received a copy of the GNU Lesser General Public
      15                 :            :  * License along with libplacebo. If not, see <http://www.gnu.org/licenses/>.
      16                 :            :  */
      17                 :            : 
      18                 :            : #include <math.h>
      19                 :            : 
      20                 :            : #include "common.h"
      21                 :            : 
      22                 :            : #include <libplacebo/tone_mapping.h>
      23                 :            : 
      24                 :            : #define fclampf(x, lo, hi) fminf(fmaxf(x, lo), hi)
      25                 :         87 : static void fix_constants(struct pl_tone_map_constants *c)
      26                 :            : {
      27                 :            :     const float eps = 1e-6f;
      28                 :         87 :     c->knee_adaptation   = fclampf(c->knee_adaptation, 0.0f, 1.0f);
      29                 :         87 :     c->knee_minimum      = fclampf(c->knee_minimum, eps, 0.5f - eps);
      30                 :         87 :     c->knee_maximum      = fclampf(c->knee_maximum, 0.5f + eps, 1.0f - eps);
      31                 :         87 :     c->knee_default      = fclampf(c->knee_default, c->knee_minimum, c->knee_maximum);
      32                 :         87 :     c->knee_offset       = fclampf(c->knee_offset, 0.5f, 2.0f);
      33                 :         87 :     c->slope_tuning      = fclampf(c->slope_tuning, 0.0f, 10.0f);
      34                 :         87 :     c->slope_offset      = fclampf(c->slope_offset, 0.0f, 1.0f);
      35                 :         87 :     c->spline_contrast   = fclampf(c->spline_contrast, 0.0f, 1.5f);
      36                 :         87 :     c->reinhard_contrast = fclampf(c->reinhard_contrast, eps, 1.0f - eps);
      37                 :         87 :     c->linear_knee       = fclampf(c->linear_knee, eps, 1.0f - eps);
      38                 :         87 :     c->exposure          = fclampf(c->exposure, eps, 10.0f);
      39                 :         87 : }
      40                 :            : 
      41                 :            : static inline bool constants_equal(const struct pl_tone_map_constants *a,
      42                 :            :                                    const struct pl_tone_map_constants *b)
      43                 :            : {
      44                 :            :     pl_static_assert(sizeof(*a) % sizeof(float) == 0);
      45                 :         12 :     return !memcmp(a, b, sizeof(*a));
      46                 :            : }
      47                 :            : 
      48                 :         22 : bool pl_tone_map_params_equal(const struct pl_tone_map_params *a,
      49                 :            :                               const struct pl_tone_map_params *b)
      50                 :            : {
      51                 :         40 :     return a->function == b->function &&
      52         [ -  + ]:         18 :            a->param == b->param &&
      53         [ -  + ]:         18 :            a->input_scaling == b->input_scaling &&
      54         [ -  + ]:         18 :            a->output_scaling == b->output_scaling &&
      55         [ -  + ]:         18 :            a->lut_size == b->lut_size &&
      56         [ -  + ]:         18 :            a->input_min == b->input_min &&
      57         [ +  + ]:         18 :            a->input_max == b->input_max &&
      58         [ -  + ]:         16 :            a->input_avg == b->input_avg &&
      59         [ +  + ]:         16 :            a->output_min == b->output_min &&
      60         [ +  + ]:         14 :            a->output_max == b->output_max &&
      61   [ +  +  -  +  :         46 :            constants_equal(&a->constants, &b->constants) &&
                   -  + ]
      62                 :         12 :            pl_hdr_metadata_equal(&a->hdr, &b->hdr);
      63                 :            : }
      64                 :            : 
      65                 :         58 : bool pl_tone_map_params_noop(const struct pl_tone_map_params *p)
      66                 :            : {
      67                 :         58 :     float in_min = pl_hdr_rescale(p->input_scaling, PL_HDR_NITS, p->input_min);
      68                 :         58 :     float in_max = pl_hdr_rescale(p->input_scaling, PL_HDR_NITS, p->input_max);
      69                 :         58 :     float out_min = pl_hdr_rescale(p->output_scaling, PL_HDR_NITS, p->output_min);
      70                 :         58 :     float out_max = pl_hdr_rescale(p->output_scaling, PL_HDR_NITS, p->output_max);
      71                 :         58 :     bool can_inverse = p->function->map_inverse;
      72                 :            : 
      73                 :         94 :     return fabs(in_min - out_min) < 1e-4 && // no BPC
      74   [ +  +  +  - ]:         58 :            in_max < out_max + 1e-2 && // no range reduction
      75   [ -  +  -  - ]:         36 :            (out_max < in_max + 1e-2 || !can_inverse); // no inverse tone-mapping
      76                 :            : }
      77                 :            : 
      78                 :         87 : void pl_tone_map_params_infer(struct pl_tone_map_params *par)
      79                 :            : {
      80         [ -  + ]:         87 :     if (!par->function)
      81                 :          0 :         par->function = &pl_tone_map_clip;
      82                 :            : 
      83         [ -  + ]:         87 :     if (par->param) {
      84                 :            :         // Backwards compatibility for older API
      85   [ #  #  #  # ]:          0 :         if (par->function == &pl_tone_map_st2094_40 || par->function == &pl_tone_map_st2094_10)
      86                 :          0 :             par->constants.knee_adaptation = par->param;
      87         [ #  # ]:          0 :         if (par->function == &pl_tone_map_bt2390)
      88                 :          0 :             par->constants.knee_offset = par->param;
      89         [ #  # ]:          0 :         if (par->function == &pl_tone_map_spline)
      90                 :          0 :             par->constants.spline_contrast = par->param;
      91         [ #  # ]:          0 :         if (par->function == &pl_tone_map_reinhard)
      92                 :          0 :             par->constants.reinhard_contrast = par->param;
      93   [ #  #  #  # ]:          0 :         if (par->function == &pl_tone_map_mobius || par->function == &pl_tone_map_gamma)
      94                 :          0 :             par->constants.linear_knee = par->param;
      95   [ #  #  #  # ]:          0 :         if (par->function == &pl_tone_map_linear || par->function == &pl_tone_map_linear_light)
      96                 :          0 :             par->constants.exposure = par->param;
      97                 :            :     }
      98                 :            : 
      99                 :         87 :     fix_constants(&par->constants);
     100                 :            : 
     101                 :            :     // Constrain the input peak to be no less than target SDR white
     102                 :         87 :     float sdr = pl_hdr_rescale(par->output_scaling, par->input_scaling, par->output_max);
     103                 :         87 :     sdr = fminf(sdr, pl_hdr_rescale(PL_HDR_NITS, par->input_scaling, PL_COLOR_SDR_WHITE));
     104                 :         87 :     par->input_max = fmaxf(par->input_max, sdr);
     105                 :            : 
     106                 :            :     // Constrain the output peak if function does not support inverse mapping
     107         [ +  + ]:         87 :     if (!par->function->map_inverse)
     108                 :          8 :         par->output_max = fminf(par->output_max, par->input_max);
     109                 :         87 : }
     110                 :            : 
     111                 :            : // Infer params and rescale to function scaling
     112                 :         29 : static struct pl_tone_map_params fix_params(const struct pl_tone_map_params *params)
     113                 :            : {
     114                 :         29 :     struct pl_tone_map_params fixed = *params;
     115                 :         29 :     pl_tone_map_params_infer(&fixed);
     116                 :            : 
     117                 :         29 :     const struct pl_tone_map_function *fun = params->function;
     118                 :         29 :     fixed.input_scaling = fun->scaling;
     119                 :         29 :     fixed.output_scaling = fun->scaling;
     120                 :         29 :     fixed.input_min  = pl_hdr_rescale(params->input_scaling,  fun->scaling, fixed.input_min);
     121                 :         29 :     fixed.input_max  = pl_hdr_rescale(params->input_scaling,  fun->scaling, fixed.input_max);
     122                 :         29 :     fixed.input_avg  = pl_hdr_rescale(params->input_scaling,  fun->scaling, fixed.input_avg);
     123                 :         29 :     fixed.output_min = pl_hdr_rescale(params->output_scaling, fun->scaling, fixed.output_min);
     124                 :         29 :     fixed.output_max = pl_hdr_rescale(params->output_scaling, fun->scaling, fixed.output_max);
     125                 :            : 
     126                 :         29 :     return fixed;
     127                 :            : }
     128                 :            : 
     129                 :            : #define FOREACH_LUT(lut, V)                                                     \
     130                 :            :     for (float *_iter = lut, *_end = lut + params->lut_size, V;                 \
     131                 :            :          _iter < _end && ( V = *_iter, 1 ); *_iter++ = V)
     132                 :            : 
     133                 :         29 : static void map_lut(float *lut, const struct pl_tone_map_params *params)
     134                 :            : {
     135         [ +  + ]:         29 :     if (params->output_max > params->input_max + 1e-4) {
     136                 :            :         // Inverse tone-mapping
     137         [ -  + ]:          7 :         pl_assert(params->function->map_inverse);
     138                 :          7 :         params->function->map_inverse(lut, params);
     139                 :            :     } else {
     140                 :            :         // Forward tone-mapping
     141                 :         22 :         params->function->map(lut, params);
     142                 :            :     }
     143                 :         29 : }
     144                 :            : 
     145                 :         29 : void pl_tone_map_generate(float *out, const struct pl_tone_map_params *params)
     146                 :            : {
     147                 :         29 :     struct pl_tone_map_params fixed = fix_params(params);
     148                 :            : 
     149                 :            :     // Generate input values evenly spaced in `params->input_scaling`
     150         [ +  + ]:       5021 :     for (size_t i = 0; i < params->lut_size; i++) {
     151                 :       4992 :         float x = (float) i / (params->lut_size - 1);
     152                 :       4992 :         x = PL_MIX(params->input_min, params->input_max, x);
     153                 :       4992 :         out[i] = pl_hdr_rescale(params->input_scaling, fixed.function->scaling, x);
     154                 :            :     }
     155                 :            : 
     156                 :         29 :     map_lut(out, &fixed);
     157                 :            : 
     158                 :            :     // Sanitize outputs and adapt back to `params->scaling`
     159         [ +  + ]:       5021 :     for (size_t i = 0; i < params->lut_size; i++) {
     160   [ +  +  +  + ]:       4992 :         float x = PL_CLAMP(out[i], fixed.output_min, fixed.output_max);
     161                 :       4992 :         out[i] = pl_hdr_rescale(fixed.function->scaling, params->output_scaling, x);
     162                 :            :     }
     163                 :         29 : }
     164                 :            : 
     165                 :          0 : float pl_tone_map_sample(float x, const struct pl_tone_map_params *params)
     166                 :            : {
     167                 :          0 :     struct pl_tone_map_params fixed = fix_params(params);
     168                 :          0 :     fixed.lut_size = 1;
     169                 :            : 
     170   [ #  #  #  # ]:          0 :     x = PL_CLAMP(x, params->input_min, params->input_max);
     171                 :          0 :     x = pl_hdr_rescale(params->input_scaling, fixed.function->scaling, x);
     172                 :          0 :     map_lut(&x, &fixed);
     173   [ #  #  #  # ]:          0 :     x = PL_CLAMP(x, fixed.output_min, fixed.output_max);
     174                 :          0 :     x = pl_hdr_rescale(fixed.function->scaling, params->output_scaling, x);
     175                 :          0 :     return x;
     176                 :            : }
     177                 :            : 
     178                 :            : // Rescale from input-absolute to input-relative
     179                 :            : static inline float rescale_in(float x, const struct pl_tone_map_params *params)
     180                 :            : {
     181                 :        769 :     return (x - params->input_min) / (params->input_max - params->input_min);
     182                 :            : }
     183                 :            : 
     184                 :            : // Rescale from input-absolute to output-relative
     185                 :            : static inline float rescale(float x, const struct pl_tone_map_params *params)
     186                 :            : {
     187                 :        387 :     return (x - params->input_min) / (params->output_max - params->output_min);
     188                 :            : }
     189                 :            : 
     190                 :            : // Rescale from output-relative to output-absolute
     191                 :            : static inline float rescale_out(float x, const struct pl_tone_map_params *params)
     192                 :            : {
     193                 :       1024 :     return x * (params->output_max - params->output_min) + params->output_min;
     194                 :            : }
     195                 :            : 
     196                 :        896 : static inline float bt1886_eotf(float x, float min, float max)
     197                 :            : {
     198                 :        896 :     const float lb = powf(min, 1/2.4f);
     199                 :        896 :     const float lw = powf(max, 1/2.4f);
     200                 :        896 :     return powf((lw - lb) * x + lb, 2.4f);
     201                 :            : }
     202                 :            : 
     203                 :        896 : static inline float bt1886_oetf(float x, float min, float max)
     204                 :            : {
     205                 :        896 :     const float lb = powf(min, 1/2.4f);
     206                 :        896 :     const float lw = powf(max, 1/2.4f);
     207                 :        896 :     return (powf(x, 1/2.4f) - lb) / (lw - lb);
     208                 :            : }
     209                 :            : 
     210                 :          2 : static void noop(float *lut, const struct pl_tone_map_params *params)
     211                 :            : {
     212                 :          2 :     return;
     213                 :            : }
     214                 :            : 
     215                 :            : const struct pl_tone_map_function pl_tone_map_clip = {
     216                 :            :     .name = "clip",
     217                 :            :     .description = "No tone mapping (clip)",
     218                 :            :     .map = noop,
     219                 :            :     .map_inverse = noop,
     220                 :            : };
     221                 :            : 
     222                 :            : // Helper function to pick a knee point (for suitable methods) based on the
     223                 :            : // HDR10+ brightness metadata and scene brightness average matching.
     224                 :            : //
     225                 :            : // Inspired by SMPTE ST2094-10, with some modifications
     226                 :         16 : static void st2094_pick_knee(float *out_src_knee, float *out_dst_knee,
     227                 :            :                              const struct pl_tone_map_params *params)
     228                 :            : {
     229                 :         16 :     const float src_min = pl_hdr_rescale(params->input_scaling,  PL_HDR_PQ, params->input_min);
     230                 :         16 :     const float src_max = pl_hdr_rescale(params->input_scaling,  PL_HDR_PQ, params->input_max);
     231                 :         16 :     const float src_avg = pl_hdr_rescale(params->input_scaling,  PL_HDR_PQ, params->input_avg);
     232                 :         16 :     const float dst_min = pl_hdr_rescale(params->output_scaling, PL_HDR_PQ, params->output_min);
     233                 :         16 :     const float dst_max = pl_hdr_rescale(params->output_scaling, PL_HDR_PQ, params->output_max);
     234                 :            : 
     235                 :         16 :     const float min_knee = params->constants.knee_minimum;
     236                 :         16 :     const float max_knee = params->constants.knee_maximum;
     237                 :         16 :     const float def_knee = params->constants.knee_default;
     238                 :         16 :     const float src_knee_min = PL_MIX(src_min, src_max, min_knee);
     239                 :         16 :     const float src_knee_max = PL_MIX(src_min, src_max, max_knee);
     240                 :         16 :     const float dst_knee_min = PL_MIX(dst_min, dst_max, min_knee);
     241                 :         16 :     const float dst_knee_max = PL_MIX(dst_min, dst_max, max_knee);
     242                 :            : 
     243                 :            :     // Choose source knee based on source scene brightness
     244         [ +  + ]:         16 :     float src_knee = PL_DEF(src_avg, PL_MIX(src_min, src_max, def_knee));
     245                 :         16 :     src_knee = fclampf(src_knee, src_knee_min, src_knee_max);
     246                 :            : 
     247                 :            :     // Choose target adaptation point based on linearly re-scaling source knee
     248                 :         16 :     float target = (src_knee - src_min) / (src_max - src_min);
     249         [ -  + ]:         16 :     float adapted = PL_MIX(dst_min, dst_max, target);
     250                 :            : 
     251                 :            :     // Choose the destnation knee by picking the perceptual adaptation point
     252                 :            :     // between the source knee and the desired target. This moves the knee
     253                 :            :     // point, on the vertical axis, closer to the 1:1 (neutral) line.
     254                 :            :     //
     255                 :            :     // Adjust the adaptation strength towards 1 based on how close the knee
     256                 :            :     // point is to its extreme values (min/max knee)
     257                 :         16 :     float tuning = 1.0f - pl_smoothstep(max_knee, def_knee, target) *
     258                 :            :                           pl_smoothstep(min_knee, def_knee, target);
     259                 :         16 :     float adaptation = PL_MIX(params->constants.knee_adaptation, 1.0f, tuning);
     260                 :         16 :     float dst_knee = PL_MIX(src_knee, adapted, adaptation);
     261                 :         16 :     dst_knee = fclampf(dst_knee, dst_knee_min, dst_knee_max);
     262                 :            : 
     263                 :         16 :     *out_src_knee = pl_hdr_rescale(PL_HDR_PQ, params->input_scaling, src_knee);
     264                 :         16 :     *out_dst_knee = pl_hdr_rescale(PL_HDR_PQ, params->output_scaling, dst_knee);
     265                 :         16 : }
     266                 :            : 
     267                 :            : // Pascal's triangle
     268                 :            : static const uint16_t binom[17][17] = {
     269                 :            :     {1},
     270                 :            :     {1,1},
     271                 :            :     {1,2,1},
     272                 :            :     {1,3,3,1},
     273                 :            :     {1,4,6,4,1},
     274                 :            :     {1,5,10,10,5,1},
     275                 :            :     {1,6,15,20,15,6,1},
     276                 :            :     {1,7,21,35,35,21,7,1},
     277                 :            :     {1,8,28,56,70,56,28,8,1},
     278                 :            :     {1,9,36,84,126,126,84,36,9,1},
     279                 :            :     {1,10,45,120,210,252,210,120,45,10,1},
     280                 :            :     {1,11,55,165,330,462,462,330,165,55,11,1},
     281                 :            :     {1,12,66,220,495,792,924,792,495,220,66,12,1},
     282                 :            :     {1,13,78,286,715,1287,1716,1716,1287,715,286,78,13,1},
     283                 :            :     {1,14,91,364,1001,2002,3003,3432,3003,2002,1001,364,91,14,1},
     284                 :            :     {1,15,105,455,1365,3003,5005,6435,6435,5005,3003,1365,455,105,15,1},
     285                 :            :     {1,16,120,560,1820,4368,8008,11440,12870,11440,8008,4368,1820,560,120,16,1},
     286                 :            : };
     287                 :            : 
     288                 :          2 : static inline float st2094_intercept(uint8_t N, float Kx, float Ky)
     289                 :            : {
     290   [ +  -  -  + ]:          2 :     if (Kx <= 0 || Ky >= 1)
     291                 :          0 :         return 1.0f / N;
     292                 :            : 
     293                 :          2 :     const float slope = Ky / Kx * (1 - Kx) / (1 - Ky);
     294                 :          2 :     return fminf(slope / N, 1.0f);
     295                 :            : }
     296                 :            : 
     297                 :          2 : static void st2094_40(float *lut, const struct pl_tone_map_params *params)
     298                 :            : {
     299                 :          2 :     const float D = params->output_max;
     300                 :            : 
     301                 :            :     // Allocate space for the adjusted bezier control points, plus endpoints
     302                 :            :     float P[17], Kx, Ky, T;
     303                 :            :     uint8_t N;
     304                 :            : 
     305         [ -  + ]:          2 :     if (params->hdr.ootf.num_anchors) {
     306                 :            : 
     307                 :            :         // Use bezier curve from metadata
     308   [ #  #  #  # ]:          0 :         Kx = PL_CLAMP(params->hdr.ootf.knee_x, 0, 1);
     309   [ #  #  #  # ]:          0 :         Ky = PL_CLAMP(params->hdr.ootf.knee_y, 0, 1);
     310   [ #  #  #  # ]:          0 :         T = PL_CLAMP(params->hdr.ootf.target_luma, params->input_min, params->input_max);
     311                 :          0 :         N = params->hdr.ootf.num_anchors + 1;
     312         [ #  # ]:          0 :         pl_assert(N < PL_ARRAY_SIZE(P));
     313                 :          0 :         memcpy(P + 1, params->hdr.ootf.anchors, (N - 1) * sizeof(*P));
     314                 :          0 :         P[0] = 0.0f;
     315                 :          0 :         P[N] = 1.0f;
     316                 :            : 
     317                 :            :     } else {
     318                 :            : 
     319                 :            :         // Missing metadata, default to simple brightness matching
     320                 :            :         float src_knee, dst_knee;
     321                 :          2 :         st2094_pick_knee(&src_knee, &dst_knee, params);
     322                 :          2 :         Kx = src_knee / params->input_max;
     323                 :          2 :         Ky = dst_knee / params->output_max;
     324                 :            : 
     325                 :            :         // Solve spline to match slope at knee intercept
     326                 :          2 :         const float slope = Ky / Kx * (1 - Kx) / (1 - Ky);
     327   [ +  +  +  - ]:          2 :         N = PL_CLAMP((int) ceilf(slope), 2, PL_ARRAY_SIZE(P) - 1);
     328                 :          2 :         P[0] = 0.0f;
     329                 :          2 :         P[1] = st2094_intercept(N, Kx, Ky);
     330         [ +  + ]:          7 :         for (int i = 2; i <= N; i++)
     331                 :          5 :             P[i] = 1.0f;
     332                 :            :         T = D;
     333                 :            : 
     334                 :            :     }
     335                 :            : 
     336         [ -  + ]:          2 :     if (D < T) {
     337                 :            : 
     338                 :            :         // Output display darker than OOTF target, make brighter
     339                 :          0 :         const float Dmin = 0.0f, u = fmaxf(0.0f, (D - Dmin) / (T - Dmin));
     340                 :            : 
     341                 :            :         // Scale down the knee point to make more room for the OOTF
     342                 :          0 :         Kx *= u;
     343                 :          0 :         Ky *= u;
     344                 :            : 
     345                 :            :         // Make the slope of the knee more closely approximate a clip(),
     346                 :            :         // constrained to avoid exploding P[1]
     347                 :          0 :         const float beta = N * Kx / (1 - Kx);
     348                 :          0 :         const float Kxy = fminf(Kx * params->input_max / D, beta / (beta + 1));
     349                 :          0 :         Ky = PL_MIX(Kxy, Ky, u);
     350                 :            : 
     351         [ #  # ]:          0 :         for (int p = 2; p <= N; p++)
     352                 :          0 :             P[p] = PL_MIX(1.0f, P[p], u);
     353                 :            : 
     354                 :            :         // Make the OOTF intercept linear as D -> Dmin
     355                 :          0 :         P[1] = PL_MIX(st2094_intercept(N, Kx, Ky), P[1], u);
     356                 :            : 
     357         [ -  + ]:          2 :     } else if (D > T) {
     358                 :            : 
     359                 :            :         // Output display brighter than OOTF target, make more linear
     360         [ #  # ]:          0 :         pl_assert(params->input_max > T);
     361                 :          0 :         const float w = powf(1 - (D - T) / (params->input_max - T), 1.4f);
     362                 :            : 
     363                 :            :         // Constrain the slope of the input knee to prevent it from
     364                 :            :         // exploding and making the picture way too bright
     365                 :          0 :         Ky *= T / D;
     366                 :            : 
     367                 :            :         // Make the slope of the knee more linear by solving for f(Kx) = Kx
     368                 :          0 :         float Kxy = Kx * D / params->input_max;
     369                 :          0 :         Ky = PL_MIX(Kxy, Ky, w);
     370                 :            : 
     371         [ #  # ]:          0 :         for (int p = 2; p < N; p++) {
     372                 :          0 :             float anchor_lin = (float) p / N;
     373                 :          0 :             P[p] = PL_MIX(anchor_lin, P[p], w);
     374                 :            :         }
     375                 :            : 
     376                 :            :         // Make the OOTF intercept linear as D -> input_max
     377                 :          0 :         P[1] = PL_MIX(st2094_intercept(N, Kx, Ky), P[1], w);
     378                 :            : 
     379                 :            :     }
     380                 :            : 
     381   [ +  -  -  + ]:          2 :     pl_assert(Kx >= 0 && Kx <= 1);
     382   [ +  -  -  + ]:          2 :     pl_assert(Ky >= 0 && Ky <= 1);
     383                 :            : 
     384         [ +  + ]:        258 :     FOREACH_LUT(lut, x) {
     385                 :        256 :         x = bt1886_oetf(x, params->input_min, params->input_max);
     386                 :        256 :         x = bt1886_eotf(x, 0.0f, 1.0f);
     387                 :            : 
     388   [ +  -  +  + ]:        256 :         if (x <= Kx && Kx) {
     389                 :            :             // Linear section
     390                 :        110 :             x *= Ky / Kx;
     391                 :            :         } else {
     392                 :            :             // Bezier section
     393                 :        146 :             const float t = (x - Kx) / (1 - Kx);
     394                 :            : 
     395                 :            :             x = 0; // Bn
     396         [ +  + ]:        812 :             for (uint8_t p = 0; p <= N; p++)
     397                 :        666 :                 x += binom[N][p] * powf(t, p) * powf(1 - t, N - p) * P[p];
     398                 :            : 
     399                 :        146 :             x = Ky + (1 - Ky) * x;
     400                 :            :         }
     401                 :            : 
     402                 :        256 :         x = bt1886_oetf(x, 0.0f, 1.0f);
     403                 :        256 :         x = bt1886_eotf(x, params->output_min, params->output_max);
     404                 :            :     }
     405                 :          2 : }
     406                 :            : 
     407                 :            : const struct pl_tone_map_function pl_tone_map_st2094_40 = {
     408                 :            :     .name = "st2094-40",
     409                 :            :     .description = "SMPTE ST 2094-40 Annex B",
     410                 :            :     .param_desc = "Knee point target",
     411                 :            :     .param_min = 0.00f,
     412                 :            :     .param_def = 0.70f,
     413                 :            :     .param_max = 1.00f,
     414                 :            :     .scaling = PL_HDR_NITS,
     415                 :            :     .map = st2094_40,
     416                 :            : };
     417                 :            : 
     418                 :          1 : static void st2094_10(float *lut, const struct pl_tone_map_params *params)
     419                 :            : {
     420                 :            :     float src_knee, dst_knee;
     421                 :          1 :     st2094_pick_knee(&src_knee, &dst_knee, params);
     422                 :            : 
     423                 :          1 :     const float x1 = params->input_min;
     424                 :          1 :     const float x3 = params->input_max;
     425                 :          1 :     const float x2 = src_knee;
     426                 :            : 
     427                 :          1 :     const float y1 = params->output_min;
     428                 :          1 :     const float y3 = params->output_max;
     429                 :          1 :     const float y2 = dst_knee;
     430                 :            : 
     431                 :          1 :     const pl_matrix3x3 cmat = {{
     432                 :          1 :         { x2*x3*(y2 - y3), x1*x3*(y3 - y1), x1*x2*(y1 - y2) },
     433                 :          1 :         { x3*y3 - x2*y2,   x1*y1 - x3*y3,   x2*y2 - x1*y1   },
     434                 :          1 :         { x3 - x2,         x1 - x3,         x2 - x1         },
     435                 :            :     }};
     436                 :            : 
     437                 :          1 :     float coeffs[3] = { y1, y2, y3 };
     438                 :          1 :     pl_matrix3x3_apply(&cmat, coeffs);
     439                 :            : 
     440                 :          1 :     const float k = 1.0 / (x3*y3*(x1 - x2) + x2*y2*(x3 - x1) + x1*y1*(x2 - x3));
     441                 :          1 :     const float c1 = k * coeffs[0];
     442                 :          1 :     const float c2 = k * coeffs[1];
     443                 :          1 :     const float c3 = k * coeffs[2];
     444                 :            : 
     445         [ +  + ]:        129 :     FOREACH_LUT(lut, x)
     446                 :        128 :         x = (c1 + c2 * x) / (1 + c3 * x);
     447                 :          1 : }
     448                 :            : 
     449                 :            : const struct pl_tone_map_function pl_tone_map_st2094_10 = {
     450                 :            :     .name = "st2094-10",
     451                 :            :     .description = "SMPTE ST 2094-10 Annex B.2",
     452                 :            :     .param_desc = "Knee point target",
     453                 :            :     .param_min = 0.00f,
     454                 :            :     .param_def = 0.70f,
     455                 :            :     .param_max = 1.00f,
     456                 :            :     .scaling = PL_HDR_NITS,
     457                 :            :     .map = st2094_10,
     458                 :            : };
     459                 :            : 
     460                 :          1 : static void bt2390(float *lut, const struct pl_tone_map_params *params)
     461                 :            : {
     462                 :          1 :     const float minLum = rescale_in(params->output_min, params);
     463                 :          1 :     const float maxLum = rescale_in(params->output_max, params);
     464                 :          1 :     const float offset = params->constants.knee_offset;
     465                 :          1 :     const float ks = (1 + offset) * maxLum - offset;
     466         [ +  - ]:          1 :     const float bp = minLum > 0 ? fminf(1 / minLum, 4) : 4;
     467                 :          1 :     const float gain_inv = 1 + minLum / maxLum * powf(1 - maxLum, bp);
     468         [ +  - ]:          1 :     const float gain = maxLum < 1 ? 1 / gain_inv : 1;
     469                 :            : 
     470         [ +  + ]:        129 :     FOREACH_LUT(lut, x) {
     471                 :            :         x = rescale_in(x, params);
     472                 :            : 
     473                 :            :         // Piece-wise hermite spline
     474         [ +  - ]:        128 :         if (ks < 1) {
     475                 :        128 :             float tb = (x - ks) / (1 - ks);
     476                 :        128 :             float tb2 = tb * tb;
     477                 :        128 :             float tb3 = tb2 * tb;
     478                 :        128 :             float pb = (2 * tb3 - 3 * tb2 + 1) * ks +
     479                 :        128 :                        (tb3 - 2 * tb2 + tb) * (1 - ks) +
     480                 :        128 :                        (-2 * tb3 + 3 * tb2) * maxLum;
     481         [ +  + ]:        128 :             x = x < ks ? x : pb;
     482                 :            :         }
     483                 :            : 
     484                 :            :         // Black point adaptation
     485         [ +  - ]:        128 :         if (x < 1) {
     486                 :        128 :             x += minLum * powf(1 - x, bp);
     487                 :        128 :             x = gain * (x - minLum) + minLum;
     488                 :            :         }
     489                 :            : 
     490                 :        128 :         x = x * (params->input_max - params->input_min) + params->input_min;
     491                 :            :     }
     492                 :          1 : }
     493                 :            : 
     494                 :            : const struct pl_tone_map_function pl_tone_map_bt2390 = {
     495                 :            :     .name = "bt2390",
     496                 :            :     .description = "ITU-R BT.2390 EETF",
     497                 :            :     .scaling = PL_HDR_PQ,
     498                 :            :     .param_desc = "Knee offset",
     499                 :            :     .param_min = 0.50,
     500                 :            :     .param_def = 1.00,
     501                 :            :     .param_max = 2.00,
     502                 :            :     .map = bt2390,
     503                 :            : };
     504                 :            : 
     505                 :          1 : static void bt2446a(float *lut, const struct pl_tone_map_params *params)
     506                 :            : {
     507                 :          1 :     const float phdr = 1 + 32 * powf(params->input_max / 10000, 1/2.4f);
     508                 :          1 :     const float psdr = 1 + 32 * powf(params->output_max / 10000, 1/2.4f);
     509                 :            : 
     510         [ +  + ]:        129 :     FOREACH_LUT(lut, x) {
     511                 :        128 :         x = powf(rescale_in(x, params), 1/2.4f);
     512                 :        128 :         x = logf(1 + (phdr - 1) * x) / logf(phdr);
     513                 :            : 
     514         [ +  + ]:        128 :         if (x <= 0.7399f) {
     515                 :         94 :             x = 1.0770f * x;
     516         [ +  + ]:         34 :         } else if (x < 0.9909f) {
     517                 :         32 :             x = (-1.1510f * x + 2.7811f) * x - 0.6302f;
     518                 :            :         } else {
     519                 :          2 :             x = 0.5f * x + 0.5f;
     520                 :            :         }
     521                 :            : 
     522                 :        128 :         x = (powf(psdr, x) - 1) / (psdr - 1);
     523                 :        128 :         x = bt1886_eotf(x, params->output_min, params->output_max);
     524                 :            :     }
     525                 :          1 : }
     526                 :            : 
     527                 :          1 : static void bt2446a_inv(float *lut, const struct pl_tone_map_params *params)
     528                 :            : {
     529         [ +  + ]:        129 :     FOREACH_LUT(lut, x) {
     530                 :        128 :         x = bt1886_oetf(x, params->input_min, params->input_max);
     531                 :        128 :         x *= 255.0;
     532         [ +  + ]:        128 :         if (x > 70) {
     533                 :         67 :             x = powf(x, (2.8305e-6f * x - 7.4622e-4f) * x + 1.2528f);
     534                 :            :         } else {
     535                 :         61 :             x = powf(x, (1.8712e-5f * x - 2.7334e-3f) * x + 1.3141f);
     536                 :            :         }
     537                 :        128 :         x = powf(x / 1000, 2.4f);
     538                 :            :         x = rescale_out(x, params);
     539                 :            :     }
     540                 :          1 : }
     541                 :            : 
     542                 :            : const struct pl_tone_map_function pl_tone_map_bt2446a = {
     543                 :            :     .name = "bt2446a",
     544                 :            :     .description = "ITU-R BT.2446 Method A",
     545                 :            :     .scaling = PL_HDR_NITS,
     546                 :            :     .map = bt2446a,
     547                 :            :     .map_inverse = bt2446a_inv,
     548                 :            : };
     549                 :            : 
     550                 :         13 : static void spline(float *lut, const struct pl_tone_map_params *params)
     551                 :            : {
     552                 :            :     float src_pivot, dst_pivot;
     553                 :         13 :     st2094_pick_knee(&src_pivot, &dst_pivot, params);
     554                 :            : 
     555                 :            :     // Solve for linear knee (Pa = 0)
     556                 :         13 :     float slope = (dst_pivot - params->output_min) /
     557                 :         13 :                   (src_pivot - params->input_min);
     558                 :            : 
     559                 :            :     // Tune the slope at the knee point slightly: raise it to a user-provided
     560                 :            :     // gamma exponent, multiplied by an extra tuning coefficient designed to
     561                 :            :     // make the slope closer to 1.0 when the difference in peaks is low, and
     562                 :            :     // closer to linear when the difference between peaks is high.
     563                 :         13 :     float ratio = params->input_max / params->output_max - 1.0f;
     564                 :         13 :     ratio = fclampf(params->constants.slope_tuning * ratio,
     565                 :            :                     params->constants.slope_offset,
     566                 :            :                     1.0f + params->constants.slope_offset);
     567                 :         13 :     slope = powf(slope, (1.0f - params->constants.spline_contrast) * ratio);
     568                 :            : 
     569                 :            :     // Normalize everything the pivot to make the math easier
     570                 :         13 :     const float in_min = params->input_min - src_pivot;
     571                 :         13 :     const float in_max = params->input_max - src_pivot;
     572                 :         13 :     const float out_min = params->output_min - dst_pivot;
     573                 :         13 :     const float out_max = params->output_max - dst_pivot;
     574                 :            : 
     575                 :            :     // Solve P of order 2 for:
     576                 :            :     //  P(in_min) = out_min
     577                 :            :     //  P'(0.0) = slope
     578                 :            :     //  P(0.0) = 0.0
     579                 :         13 :     const float Pa = (out_min - slope * in_min) / (in_min * in_min);
     580                 :            :     const float Pb = slope;
     581                 :            : 
     582                 :            :     // Solve Q of order 3 for:
     583                 :            :     //  Q(in_max) = out_max
     584                 :            :     //  Q''(in_max) = 0.0
     585                 :            :     //  Q(0.0) = 0.0
     586                 :            :     //  Q'(0.0) = slope
     587                 :         13 :     const float t = 2 * in_max * in_max;
     588                 :         13 :     const float Qa = (slope * in_max - out_max) / (in_max * t);
     589                 :         13 :     const float Qb = -3 * (slope * in_max - out_max) / t;
     590                 :            :     const float Qc = slope;
     591                 :            : 
     592         [ +  + ]:       2957 :     FOREACH_LUT(lut, x) {
     593                 :       2944 :         x -= src_pivot;
     594         [ +  + ]:       2944 :         x = x > 0 ? ((Qa * x + Qb) * x + Qc) * x : (Pa * x + Pb) * x;
     595                 :       2944 :         x += dst_pivot;
     596                 :            :     }
     597                 :         13 : }
     598                 :            : 
     599                 :            : const struct pl_tone_map_function pl_tone_map_spline = {
     600                 :            :     .name = "spline",
     601                 :            :     .description = "Single-pivot polynomial spline",
     602                 :            :     .param_desc = "Contrast",
     603                 :            :     .param_min = 0.00f,
     604                 :            :     .param_def = 0.50f,
     605                 :            :     .param_max = 1.50f,
     606                 :            :     .scaling = PL_HDR_PQ,
     607                 :            :     .map = spline,
     608                 :            :     .map_inverse = spline,
     609                 :            : };
     610                 :            : 
     611                 :          1 : static void reinhard(float *lut, const struct pl_tone_map_params *params)
     612                 :            : {
     613                 :          1 :     const float peak = rescale(params->input_max, params),
     614                 :          1 :                 contrast = params->constants.reinhard_contrast,
     615                 :          1 :                 offset = (1.0 - contrast) / contrast,
     616                 :          1 :                 scale = (peak + offset) / peak;
     617                 :            : 
     618         [ +  + ]:        129 :     FOREACH_LUT(lut, x) {
     619                 :            :         x = rescale(x, params);
     620                 :        128 :         x = x / (x + offset);
     621                 :        128 :         x *= scale;
     622                 :            :         x = rescale_out(x, params);
     623                 :            :     }
     624                 :          1 : }
     625                 :            : 
     626                 :            : const struct pl_tone_map_function pl_tone_map_reinhard = {
     627                 :            :     .name = "reinhard",
     628                 :            :     .description = "Reinhard",
     629                 :            :     .param_desc = "Contrast",
     630                 :            :     .param_min = 0.001,
     631                 :            :     .param_def = 0.50,
     632                 :            :     .param_max = 0.99,
     633                 :            :     .map = reinhard,
     634                 :            : };
     635                 :            : 
     636                 :          1 : static void mobius(float *lut, const struct pl_tone_map_params *params)
     637                 :            : {
     638                 :          1 :     const float peak = rescale(params->input_max, params),
     639                 :          1 :                 j = params->constants.linear_knee;
     640                 :            : 
     641                 :            :     // Solve for M(j) = j; M(peak) = 1.0; M'(j) = 1.0
     642                 :            :     // where M(x) = scale * (x+a)/(x+b)
     643                 :          1 :     const float a = -j*j * (peak - 1.0f) / (j*j - 2.0f * j + peak);
     644                 :          1 :     const float b = (j*j - 2.0f * j * peak + peak) /
     645                 :          1 :                     fmaxf(1e-6f, peak - 1.0f);
     646                 :          1 :     const float scale = (b*b + 2.0f * b*j + j*j) / (b - a);
     647                 :            : 
     648         [ +  + ]:        129 :     FOREACH_LUT(lut, x) {
     649                 :            :         x = rescale(x, params);
     650         [ +  + ]:        128 :         x = x <= j ? x : scale * (x + a) / (x + b);
     651                 :            :         x = rescale_out(x, params);
     652                 :            :     }
     653                 :          1 : }
     654                 :            : 
     655                 :            : const struct pl_tone_map_function pl_tone_map_mobius = {
     656                 :            :     .name = "mobius",
     657                 :            :     .description = "Mobius",
     658                 :            :     .param_desc = "Knee point",
     659                 :            :     .param_min = 0.00,
     660                 :            :     .param_def = 0.30,
     661                 :            :     .param_max = 0.99,
     662                 :            :     .map = mobius,
     663                 :            : };
     664                 :            : 
     665                 :            : static inline float hable(float x)
     666                 :            : {
     667                 :            :     const float A = 0.15, B = 0.50, C = 0.10, D = 0.20, E = 0.02, F = 0.30;
     668                 :        129 :     return ((x * (A*x + C*B) + D*E) / (x * (A*x + B) + D*F)) - E/F;
     669                 :            : }
     670                 :            : 
     671                 :          1 : static void hable_map(float *lut, const struct pl_tone_map_params *params)
     672                 :            : {
     673                 :          1 :     const float peak = params->input_max / params->output_max,
     674                 :          1 :                 scale = 1.0f / hable(peak);
     675                 :            : 
     676         [ +  + ]:        129 :     FOREACH_LUT(lut, x) {
     677                 :        128 :         x = bt1886_oetf(x, params->input_min, params->input_max);
     678                 :        128 :         x = bt1886_eotf(x, 0, peak);
     679                 :        128 :         x = scale * hable(x);
     680                 :        128 :         x = bt1886_oetf(x, 0, 1);
     681                 :        128 :         x = bt1886_eotf(x, params->output_min, params->output_max);
     682                 :            :     }
     683                 :          1 : }
     684                 :            : 
     685                 :            : const struct pl_tone_map_function pl_tone_map_hable = {
     686                 :            :     .name = "hable",
     687                 :            :     .description = "Filmic tone-mapping (Hable)",
     688                 :            :     .map = hable_map,
     689                 :            : };
     690                 :            : 
     691                 :          1 : static void gamma_map(float *lut, const struct pl_tone_map_params *params)
     692                 :            : {
     693                 :          1 :     const float peak = rescale(params->input_max, params),
     694                 :          1 :                 cutoff = params->constants.linear_knee,
     695                 :          1 :                 gamma = logf(cutoff) / logf(cutoff / peak);
     696                 :            : 
     697         [ +  + ]:        129 :     FOREACH_LUT(lut, x) {
     698                 :            :         x = rescale(x, params);
     699         [ +  + ]:        128 :         x = x > cutoff ? powf(x / peak, gamma) : x;
     700                 :            :         x = rescale_out(x, params);
     701                 :            :     }
     702                 :          1 : }
     703                 :            : 
     704                 :            : const struct pl_tone_map_function pl_tone_map_gamma = {
     705                 :            :     .name = "gamma",
     706                 :            :     .description = "Gamma function with knee",
     707                 :            :     .param_desc = "Knee point",
     708                 :            :     .param_min = 0.001,
     709                 :            :     .param_def = 0.30,
     710                 :            :     .param_max = 1.00,
     711                 :            :     .map = gamma_map,
     712                 :            : };
     713                 :            : 
     714                 :          4 : static void linear(float *lut, const struct pl_tone_map_params *params)
     715                 :            : {
     716                 :          4 :     const float gain = params->constants.exposure;
     717                 :            : 
     718         [ +  + ]:        516 :     FOREACH_LUT(lut, x) {
     719                 :            :         x = rescale_in(x, params);
     720                 :        512 :         x *= gain;
     721                 :            :         x = rescale_out(x, params);
     722                 :            :     }
     723                 :          4 : }
     724                 :            : 
     725                 :            : const struct pl_tone_map_function pl_tone_map_linear = {
     726                 :            :     .name = "linear",
     727                 :            :     .description = "Perceptually linear stretch",
     728                 :            :     .param_desc = "Exposure",
     729                 :            :     .param_min = 0.001,
     730                 :            :     .param_def = 1.00,
     731                 :            :     .param_max = 10.0,
     732                 :            :     .scaling = PL_HDR_PQ,
     733                 :            :     .map = linear,
     734                 :            :     .map_inverse = linear,
     735                 :            : };
     736                 :            : 
     737                 :            : const struct pl_tone_map_function pl_tone_map_linear_light = {
     738                 :            :     .name = "linearlight",
     739                 :            :     .description = "Linear light stretch",
     740                 :            :     .param_desc = "Exposure",
     741                 :            :     .param_min = 0.001,
     742                 :            :     .param_def = 1.00,
     743                 :            :     .param_max = 10.0,
     744                 :            :     .scaling = PL_HDR_NORM,
     745                 :            :     .map = linear,
     746                 :            :     .map_inverse = linear,
     747                 :            : };
     748                 :            : 
     749                 :            : const struct pl_tone_map_function * const pl_tone_map_functions[] = {
     750                 :            :     &pl_tone_map_clip,
     751                 :            :     &pl_tone_map_st2094_40,
     752                 :            :     &pl_tone_map_st2094_10,
     753                 :            :     &pl_tone_map_bt2390,
     754                 :            :     &pl_tone_map_bt2446a,
     755                 :            :     &pl_tone_map_spline,
     756                 :            :     &pl_tone_map_reinhard,
     757                 :            :     &pl_tone_map_mobius,
     758                 :            :     &pl_tone_map_hable,
     759                 :            :     &pl_tone_map_gamma,
     760                 :            :     &pl_tone_map_linear,
     761                 :            :     &pl_tone_map_linear_light,
     762                 :            :     NULL
     763                 :            : };
     764                 :            : 
     765                 :            : const int pl_num_tone_map_functions = PL_ARRAY_SIZE(pl_tone_map_functions) - 1;
     766                 :            : 
     767                 :          0 : const struct pl_tone_map_function *pl_find_tone_map_function(const char *name)
     768                 :            : {
     769         [ #  # ]:          0 :     for (int i = 0; i < pl_num_tone_map_functions; i++) {
     770         [ #  # ]:          0 :         if (strcmp(name, pl_tone_map_functions[i]->name) == 0)
     771                 :          0 :             return pl_tone_map_functions[i];
     772                 :            :     }
     773                 :            : 
     774                 :            :     return NULL;
     775                 :            : }

Generated by: LCOV version 1.16