/********************************************************
    © 2020 Continuum Graphics LLC. All Rights Reserved
 ********************************************************/

#if !defined _FILMIC_TONING_
#define _FILMIC_TONING_

//TODO: Move to ACES space (RRT Transform) and make more user friendly. Create more ODT profiles for support with nova.

struct CurveSegment {
	/*
		m_offsetX = 0.0f;
		m_offsetY = 0.0f;
		m_scaleX = 1.0f; // always 1 or -1
		m_lnA = 0.0f;
		m_B = 1.0f;
	*/

	float m_offsetX;
	float m_offsetY;
	float m_scaleX; // always 1 or -1
	float m_scaleY;
	float m_lnA;
	float m_B;
};

struct FullCurve {
	/*
		m_W = 1.0f;
		m_invW = 1.0f;

		m_x0 = .25f;
		m_y0 = .25f;
		m_x1 = .75f;
		m_y1 = .75f;
    */

	float m_W;
	float m_invW;

	float m_x0;
	float m_x1;
	float m_y0;
	float m_y1;

	CurveSegment m_segments[3];
	CurveSegment m_invSegments[3];
};

float EvalSegment(float x, CurveSegment seg) {
    float x0 = (x - seg.m_offsetX) * seg.m_scaleX;
    float y0 = 0.0;

    // log(0) is undefined but our function should evaluate to 0. There are better ways to handle this,
    // but it's doing it the slow way here for clarity.
    if (x0 > 0) {
        y0 = exp(seg.m_lnA + seg.m_B * log(x0));
    }

    return y0 * seg.m_scaleY + seg.m_offsetY;
}

float EvalSegmentInv(float y, CurveSegment seg) {
    float y0 = (y - seg.m_offsetY) / seg.m_scaleY;
    float x0 = 0.0;
    
    // watch out for log(0) again
    if (y0 > 0) {
        x0 = exp((log(y0) - seg.m_lnA) / seg.m_B);
    }

    return x0 / seg.m_scaleX + seg.m_offsetX;
}

float EvalFullCurve(float srcX, FullCurve curve) {
    float normX = srcX * curve.m_invW;
    int index = (normX < curve.m_x0) ? 0 : ((normX < curve.m_x1) ? 1 : 2);

    CurveSegment segment = curve.m_segments[index];

    return EvalSegment(normX, segment);
}

float EvalFullCurveInv(float y, FullCurve curve) {
    int index = (y < curve.m_y0) ? 0 : ((y < curve.m_y1) ? 1 : 2);
    CurveSegment segment = curve.m_segments[index];

    return EvalSegmentInv(y, segment) * curve.m_W;
}

#define TOE_STRENGTH 0.0 //[0-1]
#define TOE_LENGTH 0.15 //[0-1]
#define SHOULDER_STRENGTH 2.0 //[0-inf]
#define SHOULDER_LENGTH 0.5 //[1e-5 - inf]

#define SHOULDER_ANGLE 1.0 //[0-1]
#define GAMMA_TONEMAP 1.0

void CreateCurve(inout FullCurve dstCurve) {
    // This is not actually the display gamma. It's just a UI space to avoid having to
    // enter small numbers for the input.
    const float perceptualGamma = 2.2;

    // constraints
    const float toeLength = pow(clamp01(TOE_LENGTH), perceptualGamma);
    const float toeStrength = clamp01(TOE_STRENGTH);
    const float shoulderAngle = SHOULDER_ANGLE;
    const float shoulderLength = max(1e-5, clamp01(SHOULDER_LENGTH));

    const float shoulderStrength = max(0.0, SHOULDER_STRENGTH);
    const float gamma = GAMMA_TONEMAP;

    // apply base params
    // toe goes from 0 to 0.5
    const float x0 = toeLength * 0.5;
    const float y0 = (1.0 - toeStrength) * x0; // lerp from 0 to x0

    const float remainingY = 1.0 - y0;

    const float y1_offset = (1.0 - shoulderLength) * remainingY;
    const float x1 = x0 + y1_offset;
    const float y1 = y0 + y1_offset;

    // filmic shoulder strength is in F stops
    const float extraW = exp2(shoulderStrength) - 1.0;
    const float W = x0 + remainingY + extraW;

    const float params_x0 = x0 / W; // srcParams.m_x0
    const float params_x1 = x1 / W; // srcParams.m_x1

    // normalize params to 1.0 range
    const float params_overshootX = ((W * 2.0) * shoulderAngle * shoulderStrength) / W;

    const float dx = (params_x1 - params_x0);
    const float m = dx == 0.0 ? 1.0 : (y1 - y0) / dx;
    const float b = y0 - params_x0 * m;
    
    // base function of linear section plus gamma is
    // y = (mx+b)^g

    // which we can rewrite as
    // y = exp(g*ln(m) + g*ln(x+b/m))

    // and our evaluation function is (skipping the if parts):
    /*
        float x0 = (x - m_offsetX)*m_scaleX;
        y0 = expf(m_lnA + m_B*logf(x0));
        return y0*m_scaleY + m_offsetY;
    */

    CurveSegment midSegment;
    midSegment.m_offsetX = -(b / m);
    midSegment.m_offsetY = 0.0;
    midSegment.m_scaleX = 1.0;
    midSegment.m_scaleY = 1.0;
    midSegment.m_lnA = gamma * log(m);
    midSegment.m_B = gamma;

    dstCurve.m_segments[1] = midSegment; // output

    const float toeM = gamma * m * pow(m * params_x0 + b, gamma - 1.0);
    const float shoulderM = gamma * m * pow(m * params_x1 + b, gamma - 1.0);

    // apply gamma to endpoints
    const float params_y0 = max(1e-5, pow(y0, gamma));
    const float params_y1 = max(1e-5, pow(y1, gamma));

    const float params_overshootY = pow(1.0 + 0.5 * shoulderAngle * shoulderStrength, gamma) - 1.0;

    // toe section
    CurveSegment toeSegment;
    toeSegment.m_offsetX = 0.0;
    toeSegment.m_offsetY = 0.0;
    toeSegment.m_scaleX = 1.0;
    toeSegment.m_scaleY = 1.0;

    toeSegment.m_B = (toeM * params_x0) / params_y0;
    toeSegment.m_lnA = log(params_y0) - toeSegment.m_B * log(params_x0);

    dstCurve.m_segments[0] = toeSegment;

    // shoulder section
    // use the simple version that is usually too flat
    const float shoulder_x0 = (1.0 + params_overshootX) - params_x1;
    const float shoulder_y0 = (1.0 + params_overshootY) - params_y1;

    const float B = (shoulderM * shoulder_x0) / shoulder_y0;
    const float lnA = log(shoulder_y0) - B * log(shoulder_x0);

    CurveSegment shoulderSegment;
    shoulderSegment.m_offsetX = (1.0 + params_overshootX);
    shoulderSegment.m_offsetY = (1.0 + params_overshootY);

    shoulderSegment.m_scaleX = -1.0f;
    shoulderSegment.m_scaleY = -1.0;
    shoulderSegment.m_lnA = lnA;
    shoulderSegment.m_B = B;

    dstCurve.m_segments[2] = shoulderSegment; // output

    dstCurve.m_W = W; // output
    dstCurve.m_invW = 1.0 / W; // output
    
    dstCurve.m_x0 = params_x0; // output
    dstCurve.m_x1 = params_x1; // output
    dstCurve.m_y0 = params_y0; // output
    dstCurve.m_y1 = params_y1; // output

    // Normalize so that we hit 1.0 at our white point. We wouldn't have do this if we
    // skipped the overshoot part.
    // evaluate shoulder at the end of the curve
    float invScale = 1.0 / EvalSegment(1.0, dstCurve.m_segments[2]);

    dstCurve.m_segments[0].m_offsetY *= invScale; // output
    dstCurve.m_segments[0].m_scaleY *= invScale; // output

    dstCurve.m_segments[1].m_offsetY *= invScale; // output
    dstCurve.m_segments[1].m_scaleY *= invScale; // output

    dstCurve.m_segments[2].m_offsetY *= invScale; // output
    dstCurve.m_segments[2].m_scaleY *= invScale; // output
}

vec3 FilmicTonemap(vec3 color) {
    FullCurve FinalCurve;
    CreateCurve(FinalCurve);

    color.r = EvalFullCurve(color.r, FinalCurve);
    color.g = EvalFullCurve(color.g, FinalCurve);
    color.b = EvalFullCurve(color.b, FinalCurve);

    return color;
}


#endif
