#if !defined INCLUDE_ACES_TONESCALES
#define INCLUDE_ACES_TONESCALES

#include "utility.glsl"

struct SegmentedSplineParamsC5 {
    float slope_low; // log-log slope of low linear extension
    float slope_high; // log-log slope of high linear extension
    vec2 log_min_point; // {log luminance, log luminance} linear extension below
                        // this
    vec2 log_mid_point; // {log luminance, log luminance}
    vec2 log_max_point; // {log luminance, log luminance} linear extension above
                        // this
    float[6] coeff_low; // Coefficients for B-spline between mid point and max
                        // point (units of log luminance)
    float[6] coeff_high; // Coefficients for B-spline between min point and mid
                         // point (units of log luminance)
};

struct SegmentedSplineParamsC9 {
    float slope_low; // log-log slope of low linear extension
    float slope_high; // log-log slope of high linear extension
    vec2 log_min_point; // {log luminance, log luminance} linear extension below
                        // this
    vec2 log_mid_point; // {log luminance, log luminance}
    vec2 log_max_point; // {log luminance, log luminance} linear extension above
                        // this
    float[10] coeff_low; // Coefficients for B-spline between mid point and max
                         // point (units of log luminance)
    float[10] coeff_high; // Coefficients for B-spline between min point and mid
                          // point (units of log luminance)
};

// Textbook monomial to basis-function conversion matrix
const mat3 M = mat3(0.5, -1.0, 0.5, -1.0, 1.0, 0.5, 0.5, 0.0, 0.0);

float segmented_spline_c5_fwd(float x) {
    // RRT parameters
    const SegmentedSplineParamsC5 params = SegmentedSplineParamsC5(
        0.0, // slope_low
        0.0, // slope_high
        log(vec2(0.18 * exp2(-15.0), 0.0001)) * rcp(log(10.0)), // log_min_point
        log(vec2(0.18, 4.8)) * rcp(log(10.0)), // log_mid_point
        log(vec2(0.18 * exp2(18.0), 1000.0)) * rcp(log(10.0)), // log_max_point
        // coeff_low
        float[6](
            -4.0000000000,
            -4.0000000000,
            -3.1573765773,
            -0.4852499958,
            1.8477324706,
            1.8477324706
        ),
        // coeff_high
        float[6](
            -0.7185482425,
            2.0810307172,
            3.6681241237,
            4.0000000000,
            4.0000000000,
            4.0000000000
        )
    );

    // Check for negatives or zero before taking the log. If negative or zero,
    // set to 1e-6
    float log_x = log10(max(x, eps));
    float log_y;

    if (log_x <= params.log_min_point.x) {
        log_y = log_x * params.slope_low +
            (params.log_min_point.y -
             params.slope_low * params.log_min_point.x);
    } else if (log_x > params.log_min_point.x &&
               log_x < params.log_mid_point.x) {
        float knot_coord = 3.0 * (log_x - params.log_min_point.x) /
            (params.log_mid_point.x - params.log_min_point.x);
        uint i = uint(knot_coord);
        float f = fract(knot_coord);

        vec3 cf = vec3(
            params.coeff_low[i],
            params.coeff_low[i + 1],
            params.coeff_low[i + 2]
        );

        vec3 monomials = vec3(f * f, f, 1.0);

        log_y = dot(monomials, M * cf);
    } else if (log_x >= params.log_mid_point.x &&
               log_x <= params.log_max_point.x) {
        float knot_coord = 3.0 * (log_x - params.log_mid_point.x) /
            (params.log_max_point.x - params.log_mid_point.x);
        uint i = uint(knot_coord);
        float f = fract(knot_coord);

        vec3 cf = vec3(
            params.coeff_high[i],
            params.coeff_high[i + 1],
            params.coeff_high[i + 2]
        );

        vec3 monomials = vec3(f * f, f, 1.0);

        log_y = dot(monomials, M * cf);
    } else {
        log_y = log_x * params.slope_high +
            (params.log_max_point.y -
             params.slope_high * params.log_max_point.x);
    }

    return pow(10.0, log_y);
}

float segmented_spline_c9_fwd(float x) {
    // 48nit ODT parameters
    const SegmentedSplineParamsC9 params = SegmentedSplineParamsC9(
        0.0, // slope_low
        0.04, // slope_high
        // ODT min/max/mid points precomputed using a small C++ program
        // This avoids 3 unnecessary SegmentedSplineC5 calls per fragment
        vec2(-2.5406231880, -1.6989699602),
        vec2(0.6812411547, 0.6812412143),
        vec2(3.0024764538, 1.6812412739),
        // coeff_low
        float[10](
            -1.6989700043,
            -1.6989700043,
            -1.4779000000,
            -1.2291000000,
            -0.8648000000,
            -0.4480000000,
            0.0051800000,
            0.4511080334,
            0.9113744414,
            0.9113744414
        ),
        // coeff_high
        float[10](
            0.5154386965,
            0.8470437783,
            1.1358000000,
            1.3802000000,
            1.5197000000,
            1.5985000000,
            1.6467000000,
            1.6746091357,
            1.6878733390,
            1.6878733390
        )
    );

    // Check for negatives or zero before taking the log. If negative or zero,
    // set to 1e-6
    float log_x = log10(max(x, eps));
    float log_y;

    if (log_x <= params.log_min_point.x) {
        log_y = log_x * params.slope_low +
            (params.log_min_point.y -
             params.slope_low * params.log_min_point.x);
    } else if ((log_x > params.log_min_point.x) &&
               (log_x < params.log_mid_point.x)) {
        float knot_coord = 7.0 * (log_x - params.log_min_point.x) /
            (params.log_mid_point.x - params.log_min_point.x);
        uint i = uint(knot_coord);
        float f = fract(knot_coord);

        vec3 cf = vec3(
            params.coeff_low[i],
            params.coeff_low[i + 1],
            params.coeff_low[i + 2]
        );

        vec3 monomials = vec3(f * f, f, 1.0);

        log_y = dot(monomials, M * cf);
    } else if ((log_x >= params.log_mid_point.x) &&
               (log_x <= params.log_max_point.x)) {
        float knot_coord = 7.0 * (log_x - params.log_mid_point.x) /
            (params.log_max_point.x - params.log_mid_point.x);
        uint i = uint(knot_coord);
        float f = fract(knot_coord);

        vec3 cf = vec3(
            params.coeff_high[i],
            params.coeff_high[i + 1],
            params.coeff_high[i + 2]
        );

        vec3 monomials = vec3(f * f, f, 1.0);

        log_y = dot(monomials, M * cf);
    } else {
        log_y = log_x * params.slope_high +
            (params.log_max_point.y -
             params.slope_high * params.log_max_point.x);
    }

    return pow(10.0, log_y);
}

#endif // INCLUDE_ACES_TONESCALES
