/* * Copyright (C) 2024 Niklas Haas * * This file is part of FFmpeg. * * FFmpeg is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public * License as published by the Free Software Foundation; either * version 2.1 of the License, or (at your option) any later version. * * FFmpeg is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU * Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with FFmpeg; if not, write to the Free Software * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ #include <math.h>#include <string.h> #include "libavutil/attributes.h"#include "libavutil/avassert.h"#include "libavutil/csp.h"#include "libavutil/slicethread.h" #include "cms.h"#include "csputils.h"#include "libswscale/swscale.h"#include "utils.h" bool ff_sws_color_map_noop(const SwsColorMap *map){    /* If the encoding space is different, we must go through a conversion */    if (map->src.prim != map->dst.prim || map->src.trc != map->dst.trc)        return false;     /* If the black point changes, we have to perform black point compensation */    if (av_cmp_q(map->src.min_luma, map->dst.min_luma))        return false;     switch (map->intent) {    case SWS_INTENT_ABSOLUTE_COLORIMETRIC:    case SWS_INTENT_RELATIVE_COLORIMETRIC:        return ff_prim_superset(&map->dst.gamut, &map->src.gamut) &&               av_cmp_q(map->src.max_luma, map->dst.max_luma) <= 0;    case SWS_INTENT_PERCEPTUAL:    case SWS_INTENT_SATURATION:        return ff_prim_equal(&map->dst.gamut, &map->src.gamut) &&               !av_cmp_q(map->src.max_luma, map->dst.max_luma);    default:        av_assert0(!"Invalid gamut mapping intent?");        return true;    }} /* Approximation of gamut hull at a given intensity level */static const float hull(float I){    return ((I - 6.0f) * I + 9.0f) * I;} /* For some minimal type safety, and code cleanliness */typedef struct RGB {    float R, G, B; /* nits */} RGB; typedef struct IPT {    float I, P, T;} IPT; typedef struct ICh {    float I, C, h;} ICh; static av_always_inline ICh ipt2ich(IPT c){    return (ICh) {        .I = c.I,        .C = sqrtf(c.P * c.P + c.T * c.T),        .h = atan2f(c.T, c.P),    };} static av_always_inline IPT ich2ipt(ICh c){    return (IPT) {        .I = c.I,        .P = c.C * cosf(c.h),        .T = c.C * sinf(c.h),    };} /* Helper struct containing pre-computed cached values describing a gamut */typedef struct Gamut {    SwsMatrix3x3 encoding2lms;    SwsMatrix3x3 lms2encoding;    SwsMatrix3x3 lms2content;    SwsMatrix3x3 content2lms;    av_csp_eotf_function eotf;    av_csp_eotf_function eotf_inv;    float Iavg_frame;    float Imax_frame;    float Imin, Imax;    float Lb, Lw;    AVCIExy wp;    ICh peak; /* updated as needed in loop body when hue changes */} Gamut; static Gamut gamut_from_colorspace(SwsColor fmt){    const AVColorPrimariesDesc *encoding = av_csp_primaries_desc_from_id(fmt.prim);    const AVColorPrimariesDesc content = {        .prim = fmt.gamut,        .wp   = encoding->wp,    };     const float Lw = av_q2d(fmt.max_luma), Lb = av_q2d(fmt.min_luma);    const float Imax = pq_oetf(Lw);     return (Gamut) {        .encoding2lms = ff_sws_ipt_rgb2lms(encoding),        .lms2encoding = ff_sws_ipt_lms2rgb(encoding),        .lms2content  = ff_sws_ipt_lms2rgb(&content),        .content2lms  = ff_sws_ipt_rgb2lms(&content),        .eotf         = av_csp_itu_eotf(fmt.trc),        .eotf_inv     = av_csp_itu_eotf_inv(fmt.trc),        .wp           = encoding->wp,        .Imin         = pq_oetf(Lb),        .Imax         = Imax,        .Imax_frame   = fmt.frame_peak.den ? pq_oetf(av_q2d(fmt.frame_peak)) : Imax,        .Iavg_frame   = fmt.frame_avg.den  ? pq_oetf(av_q2d(fmt.frame_avg))  : 0.0f,        .Lb           = Lb,        .Lw           = Lw,    };} static av_always_inline IPT rgb2ipt(RGB c, const SwsMatrix3x3 rgb2lms){    const float L = rgb2lms.m[0][0] * c.R +                    rgb2lms.m[0][1] * c.G +                    rgb2lms.m[0][2] * c.B;    const float M = rgb2lms.m[1][0] * c.R +                    rgb2lms.m[1][1] * c.G +                    rgb2lms.m[1][2] * c.B;    const float S = rgb2lms.m[2][0] * c.R +                    rgb2lms.m[2][1] * c.G +                    rgb2lms.m[2][2] * c.B;    const float Lp = pq_oetf(L);    const float Mp = pq_oetf(M);    const float Sp = pq_oetf(S);    return (IPT) {        .I = 0.4000f * Lp + 0.4000f * Mp + 0.2000f * Sp,        .P = 4.4550f * Lp - 4.8510f * Mp + 0.3960f * Sp,        .T = 0.8056f * Lp + 0.3572f * Mp - 1.1628f * Sp,    };} static av_always_inline RGB ipt2rgb(IPT c, const SwsMatrix3x3 lms2rgb){    const float Lp = c.I + 0.0975689f * c.P + 0.205226f * c.T;    const float Mp = c.I - 0.1138760f * c.P + 0.133217f * c.T;    const float Sp = c.I + 0.0326151f * c.P - 0.676887f * c.T;    const float L = pq_eotf(Lp);    const float M = pq_eotf(Mp);    const float S = pq_eotf(Sp);    return (RGB) {        .R = lms2rgb.m[0][0] * L +             lms2rgb.m[0][1] * M +             lms2rgb.m[0][2] * S,        .G = lms2rgb.m[1][0] * L +             lms2rgb.m[1][1] * M +             lms2rgb.m[1][2] * S,        .B = lms2rgb.m[2][0] * L +             lms2rgb.m[2][1] * M +             lms2rgb.m[2][2] * S,    };} static inline bool ingamut(IPT c, Gamut gamut){    const float min_rgb = gamut.Lb - 1e-4f;    const float max_rgb = gamut.Lw + 1e-2f;    const float Lp = c.I + 0.0975689f * c.P + 0.205226f * c.T;    const float Mp = c.I - 0.1138760f * c.P + 0.133217f * c.T;    const float Sp = c.I + 0.0326151f * c.P - 0.676887f * c.T;    if (Lp < gamut.Imin || Lp > gamut.Imax ||        Mp < gamut.Imin || Mp > gamut.Imax ||        Sp < gamut.Imin || Sp > gamut.Imax)    {        /* Values outside legal LMS range */        return false;    } else {        const float L = pq_eotf(Lp);        const float M = pq_eotf(Mp);        const float S = pq_eotf(Sp);        RGB rgb = {            .R = gamut.lms2content.m[0][0] * L +                 gamut.lms2content.m[0][1] * M +                 gamut.lms2content.m[0][2] * S,            .G = gamut.lms2content.m[1][0] * L +                 gamut.lms2content.m[1][1] * M +                 gamut.lms2content.m[1][2] * S,            .B = gamut.lms2content.m[2][0] * L +                 gamut.lms2content.m[2][1] * M +                 gamut.lms2content.m[2][2] * S,        };        return rgb.R >= min_rgb && rgb.R <= max_rgb &&               rgb.G >= min_rgb && rgb.G <= max_rgb &&               rgb.B >= min_rgb && rgb.B <= max_rgb;    }} static const float maxDelta = 5e-5f; // Find gamut intersection using specified boundsstatic inline IChdesat_bounded(float I, float h, float Cmin, float Cmax, Gamut gamut){    if (I <= gamut.Imin)        return (ICh) { .I = gamut.Imin, .C = 0, .h = h };    else if (I >= gamut.Imax)        return (ICh) { .I = gamut.Imax, .C = 0, .h = h };    else {        const float maxDI = I * maxDelta;        ICh res = { .I = I, .C = (Cmin + Cmax) / 2, .h = h };        do {            if (ingamut(ich2ipt(res), gamut)) {                Cmin = res.C;            } else {                Cmax = res.C;            }            res.C = (Cmin + Cmax) / 2;        } while (Cmax - Cmin > maxDI);         return res;    }} // Finds maximally saturated in-gamut color (for given hue)static inline ICh saturate(float hue, Gamut gamut){    static const float invphi = 0.6180339887498948f;    static const float invphi2 = 0.38196601125010515f;     ICh lo = { .I = gamut.Imin, .h = hue };    ICh hi = { .I = gamut.Imax, .h = hue };    float de = hi.I - lo.I;    ICh a = { .I = lo.I + invphi2 * de };    ICh b = { .I = lo.I + invphi  * de };    a = desat_bounded(a.I, hue, 0.0f, 0.5f, gamut);    b = desat_bounded(b.I, hue, 0.0f, 0.5f, gamut);     while (de > maxDelta) {        de *= invphi;        if (a.C > b.C) {            hi = b;            b = a;            a.I = lo.I + invphi2 * de;            a = desat_bounded(a.I, hue, lo.C - maxDelta, 0.5f, gamut);        } else {            lo = a;            a = b;            b.I = lo.I + invphi * de;            b = desat_bounded(b.I, hue, hi.C - maxDelta, 0.5f, gamut);        }    }     return a.C > b.C ? a : b;} static float softclip(float value, float source, float target){    const float j = SOFTCLIP_KNEE;    float peak, x, a, b, scale;    if (!target)        return 0.0f;     peak = source / target;    x = fminf(value / target, peak);    if (x <= j || peak <= 1.0)        return value;     /* Apply simple mobius function */    a = -j*j * (peak - 1.0f) / (j*j - 2.0f * j + peak);    b = (j*j - 2.0f * j * peak + peak) / fmaxf(1e-6f, peak - 1.0f);    scale = (b*b + 2.0f * b*j + j*j) / (b - a);     return scale * (x + a) / (x + b) * target;} /** * Something like fmixf(base, c, x) but follows an exponential curve, note * that this can be used to extend 'c' outwards for x > 1 */static inline ICh mix_exp(ICh c, float x, float gamma, float base){    return (ICh) {        .I = base + (c.I - base) * powf(x, gamma),        .C = c.C * x,        .h = c.h,    };} /** * Drop gamma for colors approaching black and achromatic to avoid numerical * instabilities, and excessive brightness boosting of grain, while also * strongly boosting gamma for values exceeding the target peak */static inline float scale_gamma(float gamma, ICh ich, Gamut gamut){    const float Imin = gamut.Imin;    const float Irel = fmaxf((ich.I - Imin) / (gamut.peak.I - Imin), 0.0f);    return gamma * powf(Irel, 3) * fminf(ich.C / gamut.peak.C, 1.0f);} /* Clip a color along the exponential curve given by `gamma` */static inline IPT clip_gamma(IPT ipt, float gamma, Gamut gamut){    float lo = 0.0f, hi = 1.0f, x = 0.5f;    const float maxDI = fmaxf(ipt.I * maxDelta, 1e-7f);    ICh ich;     if (ipt.I <= gamut.Imin)        return (IPT) { .I = gamut.Imin };    if (ingamut(ipt, gamut))        return ipt;     ich = ipt2ich(ipt);    if (!gamma)        return ich2ipt(desat_bounded(ich.I, ich.h, 0.0f, ich.C, gamut));     gamma = scale_gamma(gamma, ich, gamut);    do {        ICh test = mix_exp(ich, x, gamma, gamut.peak.I);        if (ingamut(ich2ipt(test), gamut)) {            lo = x;        } else {            hi = x;        }        x = (lo + hi) / 2.0f;    } while (hi - lo > maxDI);     return ich2ipt(mix_exp(ich, x, gamma, gamut.peak.I));} typedef struct CmsCtx CmsCtx;struct CmsCtx {    /* Tone mapping parameters */    float Qa, Qb, Qc, Pa, Pb, src_knee, dst_knee; /* perceptual */    float I_scale, I_offset; /* linear methods */     /* Colorspace parameters */    Gamut src;    Gamut tmp; /* after tone mapping */    Gamut dst;    SwsMatrix3x3 adaptation; /* for absolute intent */     /* Invocation parameters */    SwsColorMap map;    float (*tone_map)(const CmsCtx *ctx, float I);    IPT (*adapt_colors)(const CmsCtx *ctx, IPT ipt);    v3u16_t *input;    v3u16_t *output;     /* Threading parameters */    int slice_size;    int size_input;    int size_output_I;    int size_output_PT;}; /** * Helper function to pick a knee point based on the * HDR10+ brightness * metadata and scene brightness average matching. * * Inspired by SMPTE ST2094-10, with some modifications */static void st2094_pick_knee(float src_max, float src_min, float src_avg,                             float dst_max, float dst_min,                             float *out_src_knee, float *out_dst_knee){    const float min_knee = PERCEPTUAL_KNEE_MIN;    const float max_knee = PERCEPTUAL_KNEE_MAX;    const float def_knee = PERCEPTUAL_KNEE_DEF;    const float src_knee_min = fmixf(src_min, src_max, min_knee);    const float src_knee_max = fmixf(src_min, src_max, max_knee);    const float dst_knee_min = fmixf(dst_min, dst_max, min_knee);    const float dst_knee_max = fmixf(dst_min, dst_max, max_knee);    float src_knee, target, adapted, tuning, adaptation, dst_knee;     /* Choose source knee based on dynamic source scene brightness */    src_knee = src_avg ? src_avg : fmixf(src_min, src_max, def_knee);    src_knee = av_clipf(src_knee, src_knee_min, src_knee_max);     /* Choose target adaptation point based on linearly re-scaling source knee */    target = (src_knee - src_min) / (src_max - src_min);    adapted = fmixf(dst_min, dst_max, target);     /**     * Choose the destnation knee by picking the perceptual adaptation point     * between the source knee and the desired target. This moves the knee     * point, on the vertical axis, closer to the 1:1 (neutral) line.     *     * Adjust the adaptation strength towards 1 based on how close the knee     * point is to its extreme values (min/max knee)     */    tuning = smoothstepf(max_knee, def_knee, target) *             smoothstepf(min_knee, def_knee, target);    adaptation = fmixf(1.0f, PERCEPTUAL_ADAPTATION, tuning);    dst_knee = fmixf(src_knee, adapted, adaptation);    dst_knee = av_clipf(dst_knee, dst_knee_min, dst_knee_max);     *out_src_knee = src_knee;    *out_dst_knee = dst_knee;} static void tone_map_setup(CmsCtx *ctx, bool dynamic){    const float dst_min = ctx->dst.Imin;    const float dst_max = ctx->dst.Imax;    const float src_min = ctx->src.Imin;    const float src_max = dynamic ? ctx->src.Imax_frame : ctx->src.Imax;    const float src_avg = dynamic ? ctx->src.Iavg_frame : 0.0f;    float slope, ratio, in_min, in_max, out_min, out_max, t;     switch (ctx->map.intent) {    case SWS_INTENT_PERCEPTUAL:        st2094_pick_knee(src_max, src_min, src_avg, dst_max, dst_min,                         &ctx->src_knee, &ctx->dst_knee);         /* Solve for linear knee (Pa = 0) */        slope = (ctx->dst_knee - dst_min) / (ctx->src_knee - src_min);         /**         * Tune the slope at the knee point slightly: raise it to a user-provided         * gamma exponent, multiplied by an extra tuning coefficient designed to         * make the slope closer to 1.0 when the difference in peaks is low, and         * closer to linear when the difference between peaks is high.         */        ratio = src_max / dst_max - 1.0f;        ratio = av_clipf(SLOPE_TUNING * ratio, SLOPE_OFFSET, 1.0f + SLOPE_OFFSET);        slope = powf(slope, (1.0f - PERCEPTUAL_CONTRAST) * ratio);         /* Normalize everything the pivot to make the math easier */        in_min  = src_min - ctx->src_knee;        in_max  = src_max - ctx->src_knee;        out_min = dst_min - ctx->dst_knee;        out_max = dst_max - ctx->dst_knee;         /**         * Solve P of order 2 for:         *  P(in_min) = out_min         *  P'(0.0) = slope         *  P(0.0) = 0.0         */        ctx->Pa = (out_min - slope * in_min) / (in_min * in_min);        ctx->Pb = slope;         /**         * Solve Q of order 3 for:         *  Q(in_max) = out_max         *  Q''(in_max) = 0.0         *  Q(0.0) = 0.0         *  Q'(0.0) = slope         */        t = 2 * in_max * in_max;        ctx->Qa = (slope * in_max - out_max) / (in_max * t);        ctx->Qb = -3 * (slope * in_max - out_max) / t;        ctx->Qc = slope;        break;    case SWS_INTENT_SATURATION:        /* Linear stretch */        ctx->I_scale = (dst_max - dst_min) / (src_max - src_min);        ctx->I_offset = dst_min - src_min * ctx->I_scale;        break;    case SWS_INTENT_RELATIVE_COLORIMETRIC:        /* Pure black point adaptation */        ctx->I_scale = src_max / (src_max - src_min) /                      (dst_max / (dst_max - dst_min));        ctx->I_offset = dst_min - src_min * ctx->I_scale;        break;    case SWS_INTENT_ABSOLUTE_COLORIMETRIC:        /* Hard clip */        ctx->I_scale = 1.0f;        ctx->I_offset = 0.0f;        break;    }} static av_always_inline IPT tone_map_apply(const CmsCtx *ctx, IPT ipt){    float I = ipt.I, desat;     if (ctx->map.intent == SWS_INTENT_PERCEPTUAL) {        const float Pa = ctx->Pa, Pb = ctx->Pb;        const float Qa = ctx->Qa, Qb = ctx->Qb, Qc = ctx->Qc;        I -= ctx->src_knee;

Lines 1–500 of 766. Load more