/*

 * 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 "format.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 bounds

static inline ICh

desat_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 destination 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;

        I = I > 0 ? ((Qa * I + Qb) * I + Qc) * I : (Pa * I + Pb) * I;

        I += ctx->dst_knee;

    } else {

        I = ctx->I_scale * I + ctx->I_offset;

    }

    /**

     * Avoids raising saturation excessively when raising brightness, and

     * also desaturates when reducing brightness greatly to account for the

     * reduction in gamut volume.

     */

    desat = fminf(ipt.I / I, hull(I) / hull(ipt.I));

    return (IPT) {

        .I = I,

        .P = ipt.P * desat,

        .T = ipt.T * desat,

    };

}

static IPT perceptual(const CmsCtx *ctx, IPT ipt)

{

    ICh ich = ipt2ich(ipt);

    IPT mapped = rgb2ipt(ipt2rgb(ipt, ctx->tmp.lms2content), ctx->dst.content2lms);

    RGB rgb;

    float maxRGB;

    /* Protect in gamut region */

    const float maxC = fmaxf(ctx->tmp.peak.C, ctx->dst.peak.C);

    float k = smoothstepf(PERCEPTUAL_DEADZONE, 1.0f, ich.C / maxC);

    k *= PERCEPTUAL_STRENGTH;

    ipt.I = fmixf(ipt.I, mapped.I, k);

    ipt.P = fmixf(ipt.P, mapped.P, k);

    ipt.T = fmixf(ipt.T, mapped.T, k);

    rgb = ipt2rgb(ipt, ctx->dst.lms2content);

    maxRGB = fmaxf(rgb.R, fmaxf(rgb.G, rgb.B));

    rgb.R = fmaxf(softclip(rgb.R, maxRGB, ctx->dst.Lw), ctx->dst.Lb);

    rgb.G = fmaxf(softclip(rgb.G, maxRGB, ctx->dst.Lw), ctx->dst.Lb);

    rgb.B = fmaxf(softclip(rgb.B, maxRGB, ctx->dst.Lw), ctx->dst.Lb);

    return rgb2ipt(rgb, ctx->dst.content2lms);

}

static IPT relative(const CmsCtx *ctx, IPT ipt)

{

    return clip_gamma(ipt, COLORIMETRIC_GAMMA, ctx->dst);

}

static IPT absolute(const CmsCtx *ctx, IPT ipt)

{

    RGB rgb = ipt2rgb(ipt, ctx->dst.lms2encoding);

    float c[3] = { rgb.R, rgb.G, rgb.B };

    ff_sws_matrix3x3_apply(&ctx->adaptation, c);

    ipt = rgb2ipt((RGB) { c[0], c[1], c[2] }, ctx->dst.encoding2lms);

    return clip_gamma(ipt, COLORIMETRIC_GAMMA, ctx->dst);

}

static IPT saturation(const CmsCtx * ctx, IPT ipt)

{

    RGB rgb = ipt2rgb(ipt, ctx->tmp.lms2content);

    return rgb2ipt(rgb, ctx->dst.content2lms);

}

static av_always_inline av_const uint16_t av_round16f(float x)

{

    return av_clip_uint16(x * (UINT16_MAX - 1) + 0.5f);

}

/* Call this whenever the hue changes inside the loop body */

static av_always_inline void update_hue_peaks(CmsCtx *ctx, float P, float T)

{

    const float hue = atan2f(T, P);

    switch (ctx->map.intent) {

    case SWS_INTENT_PERCEPTUAL:

        ctx->tmp.peak = saturate(hue, ctx->tmp);

        /* fall through */

    case SWS_INTENT_RELATIVE_COLORIMETRIC:

    case SWS_INTENT_ABSOLUTE_COLORIMETRIC:

        ctx->dst.peak = saturate(hue, ctx->dst);

        return;

    default:

        return;

    }

}

static void generate_slice(void *priv, int jobnr, int threadnr, int nb_jobs,

                           int nb_threads)

{

    CmsCtx ctx = *(const CmsCtx *) priv;

    const int slice_start  = jobnr * ctx.slice_size;

    const int slice_stride = ctx.size_input * ctx.size_input;

    const int slice_end    = FFMIN((jobnr + 1) * ctx.slice_size, ctx.size_input);

    v3u16_t *input = &ctx.input[slice_start * slice_stride];

    const int output_slice_h = (ctx.size_output_PT + nb_jobs - 1) / nb_jobs;

    const int output_start   = jobnr * output_slice_h;

    const int output_stride  = ctx.size_output_PT * ctx.size_output_I;

    const int output_end     = FFMIN((jobnr + 1) * output_slice_h, ctx.size_output_PT);

    v3u16_t *output = ctx.output ? &ctx.output[output_start * output_stride] : NULL;

    const float I_scale   = 1.0f / (ctx.src.Imax - ctx.src.Imin);

    const float I_offset  = -ctx.src.Imin * I_scale;

    const float PT_offset = (float) (1 << 15) / (UINT16_MAX - 1);

    const float input_scale     = 1.0f / (ctx.size_input - 1);

    const float output_scale_PT = 1.0f / (ctx.size_output_PT - 1);

    const float output_scale_I  = (ctx.tmp.Imax - ctx.tmp.Imin) /

                                  (ctx.size_output_I - 1);

    for (int Bx = slice_start; Bx < slice_end; Bx++) {

        const float B = input_scale * Bx;

        for (int Gx = 0; Gx < ctx.size_input; Gx++) {

            const float G = input_scale * Gx;

            for (int Rx = 0; Rx < ctx.size_input; Rx++) {

                double c[3] = { input_scale * Rx, G, B };

                RGB rgb;

                IPT ipt;

                ctx.src.eotf(ctx.src.Lw, ctx.src.Lb, c);

                rgb = (RGB) { c[0], c[1], c[2] };

                ipt = rgb2ipt(rgb, ctx.src.encoding2lms);

                if (output) {

                    /* Save intermediate value to 3DLUT */

                    *input++ = (v3u16_t) {

                        av_round16f(I_scale * ipt.I + I_offset),

                        av_round16f(ipt.P + PT_offset),

                        av_round16f(ipt.T + PT_offset),

                    };

                } else {

                    update_hue_peaks(&ctx, ipt.P, ipt.T);

                    ipt = tone_map_apply(&ctx, ipt);

                    ipt = ctx.adapt_colors(&ctx, ipt);

                    rgb = ipt2rgb(ipt, ctx.dst.lms2encoding);

                    c[0] = rgb.R;

                    c[1] = rgb.G;

                    c[2] = rgb.B;

                    ctx.dst.eotf_inv(ctx.dst.Lw, ctx.dst.Lb, c);

                    *input++ = (v3u16_t) {

                        av_round16f(c[0]),

                        av_round16f(c[1]),

                        av_round16f(c[2]),

                    };

                }

            }

        }

    }

    if (!output)

        return;

    /* Generate split gamut mapping LUT */

    for (int Tx = output_start; Tx < output_end; Tx++) {

        const float T = output_scale_PT * Tx - PT_offset;

        for (int Px = 0; Px < ctx.size_output_PT; Px++) {

            const float P = output_scale_PT * Px - PT_offset;

            update_hue_peaks(&ctx, P, T);

            for (int Ix = 0; Ix < ctx.size_output_I; Ix++) {

                const float I = output_scale_I * Ix + ctx.tmp.Imin;

                IPT ipt = ctx.adapt_colors(&ctx, (IPT) { I, P, T });

                RGB rgb = ipt2rgb(ipt, ctx.dst.lms2encoding);

                double c[3] = { rgb.R, rgb.G, rgb.B };

                ctx.dst.eotf_inv(ctx.dst.Lw, ctx.dst.Lb, c);

                *output++ = (v3u16_t) {

                    av_round16f(c[0]),

                    av_round16f(c[1]),

                    av_round16f(c[2]),

                };

            }

        }

    }

}

int ff_sws_color_map_generate_static(v3u16_t *lut, int size, const SwsColorMap *map)

{

    return ff_sws_color_map_generate_dynamic(lut, NULL, size, 1, 1, map);

}

int ff_sws_color_map_generate_dynamic(v3u16_t *input, v3u16_t *output,

                                      int size_input, int size_I, int size_PT,

                                      const SwsColorMap *map)

{

    AVSliceThread *slicethread;

    int ret, num_slices;

    CmsCtx ctx = {

        .map            = *map,

        .input          = input,

        .output         = output,

        .size_input     = size_input,

        .size_output_I  = size_I,

        .size_output_PT = size_PT,

        .src            = gamut_from_colorspace(map->src),

        .dst            = gamut_from_colorspace(map->dst),

    };

    switch (ctx.map.intent) {

    case SWS_INTENT_PERCEPTUAL:            ctx.adapt_colors = perceptual; break;

    case SWS_INTENT_RELATIVE_COLORIMETRIC: ctx.adapt_colors = relative;   break;

    case SWS_INTENT_SATURATION:            ctx.adapt_colors = saturation; break;

    case SWS_INTENT_ABSOLUTE_COLORIMETRIC: ctx.adapt_colors = absolute;   break;

    default: return AVERROR(EINVAL);

    }

    if (!output) {

        /* Tone mapping is handled in a separate step when using dynamic TM */

        tone_map_setup(&ctx, false);

    }

    /* Intermediate color space after tone mapping */

    ctx.tmp      = ctx.src;

    ctx.tmp.Lb   = ctx.dst.Lb;

    ctx.tmp.Lw   = ctx.dst.Lw;

    ctx.tmp.Imin = ctx.dst.Imin;

    ctx.tmp.Imax = ctx.dst.Imax;

    if (ctx.map.intent == SWS_INTENT_ABSOLUTE_COLORIMETRIC) {

        /**

         * The IPT transform already implies an explicit white point adaptation

         * from src to dst, so to get absolute colorimetric semantics we have

         * to explicitly undo this adaptation with a * corresponding inverse.

         */

        ctx.adaptation = ff_sws_get_adaptation(&ctx.map.dst.gamut,

                                               ctx.dst.wp, ctx.src.wp);

    }

    ret = avpriv_slicethread_create(&slicethread, &ctx, generate_slice, NULL, 0);

    if (ret < 0)

        return ret;

    ctx.slice_size = (ctx.size_input + ret - 1) / ret;

    num_slices = (ctx.size_input + ctx.slice_size - 1) / ctx.slice_size;

    avpriv_slicethread_execute(slicethread, num_slices, 0);

    avpriv_slicethread_free(&slicethread);

    return 0;

}

void ff_sws_tone_map_generate(v2u16_t *lut, int size, const SwsColorMap *map)

{

    CmsCtx ctx = {

        .map = *map,

        .src = gamut_from_colorspace(map->src),

        .dst = gamut_from_colorspace(map->dst),

    };

    const float src_scale  = (ctx.src.Imax - ctx.src.Imin) / (size - 1);

    const float src_offset = ctx.src.Imin;

    const float dst_scale  = 1.0f / (ctx.dst.Imax - ctx.dst.Imin);

    const float dst_offset = -ctx.dst.Imin * dst_scale;

    tone_map_setup(&ctx, true);

    for (int i = 0; i < size; i++) {

        const float I = src_scale * i + src_offset;

        IPT ipt = tone_map_apply(&ctx, (IPT) { I, 1.0f });

        lut[i] = (v2u16_t) {

            av_round16f(dst_scale * ipt.I + dst_offset),

            av_clip_uint16(ipt.P * (1 << 15) + 0.5f),

        };

    }

}