
// CommonStrokeTessellatorEval.hpp - common GPU and CPU polar stroking evaluation routines

#ifndef __CommonStrokeTessellatorEval_hpp__
#define __CommonStrokeTessellatorEval_hpp__

#include "CommonStrokeTessellatorPrivate.hpp"

#include "CommonStrokeTessellatorTypes.hpp"

STATIC bool inRange(float t_min, float t, float t_max)
{
    return (t_min < t) && (t < t_max);
}

// Rotate v by angle theta, represented by theta's cosine and sine,
// and return rotated v's X component.
STATIC float rotate(vec2 v, float cos_theta, float sin_theta)
{
    return sin_theta*v.x - cos_theta*v.y;
}

STATIC vec2 evalQuadraticPoly(vec2 A, vec2 B, vec2 C, float t)
{
    return (A*t + B)*t + C;
}

STATIC float evalQuadraticPoly(float A, float B, float C, float t)
{
    return (A*t + B)*t + C;
}

STATIC ParametricDirection makeParametricDirection(float t, float theta)
{
    ParametricDirection rv;
    rv.t = t;
    rv.theta = theta;
    CONST vec2 gradient_direction = vec2(cos(theta), sin(theta));
    rv.normal = ortho(gradient_direction);
    return rv;
}

STATIC ParametricDirection makeParametricDirection(float t, float theta, vec2 normal)
{
    ParametricDirection rv;
    rv.t = t;
    rv.theta = theta;
    rv.normal = normal;
    return rv;
}

// From polar angle theta_s, generate a parametric t + gradient direction.
ParametricDirection StrokeTessellatorMETHOD evalCurveSampleInRange(float theta,
                                                                   float t_min,
                                                                   float t_max) CONST
{
    // Output cosine and sine of theta.
    CONST float cos_theta = cos(theta);
    CONST float sin_theta = sin(theta);
    // Returns (normalized & never degenerate) normal to caller.
    CONST vec2 gradient_direction = vec2(cos_theta, sin_theta);
    CONST vec2 normal = ortho(gradient_direction);
    
    // Continue to compute proper t...

    // Rotate the cubic segment's quadratic gradient function by the angle theta.
    CONST vec2 A = vec2(A_C[0], A_C[2]);
    CONST vec2 C = vec2(A_C[1], A_C[3]);
    CONST vec2 B = vec2(B_bogus[0],B_bogus[1]);
    CONST float N2 = rotate(A, cos_theta, sin_theta);
    CONST float N1 = rotate(B, cos_theta, sin_theta);
    CONST float N0 = rotate(C, cos_theta, sin_theta);

    float ts[2];  // For two possible t values from solving a quadratic equation.
#if !defined(NDEBUG)
    // Compute quadratic solutions at double precision for reference.
    CONST double N2d = N2;
    CONST double N1d = N1;
    CONST double N0d = N0;
    CONST double discriminant_fp64 = N1d*N1d - 4*N2d*N0d;
    CONST double sqrt_discr_fp64 = sqrt(discriminant_fp64);
    UNUSED CONST double t0d = (-N1 - sqrt_discr_fp64)/(2*N2);
    UNUSED CONST double t1d = (-N1 + sqrt_discr_fp64)/(2*N2);
#endif
    // Solve as a proper quadratic...
    CONST float discriminant = N1*N1 - 4*N2*N0;
    if (discriminant >= 0) {
        // Solve quadratic equation
        //   0 = N2^2*t + N1*t + N0
        // for t.

        // Use unconventional "citardauq" quadratic equation solver for more precision around zero.
        // This helps make up for the GPU's lack of precision.
        CONST float sqrt_discr = sqrt(discriminant);
        CONST float q = -0.5f * (N1 + (N1>=0 ? +sqrt_discr : -sqrt_discr));
        CONST float alt_t0 = q/N2;
        CONST float alt_t1 = N0/q;

        // Evaluate the conventonal quadratic equation solver for comparison.
        UNUSED CONST float t0 = (-N1 - sqrt_discr)/(2*N2);
        UNUSED CONST float t1 = (-N1 + sqrt_discr)/(2*N2);

        // Compute residual of different solutions.
        UNUSED CONST float s1alt = evalQuadraticPoly(N2, N1, N0, alt_t0);
        UNUSED CONST float s2alt = evalQuadraticPoly(N2, N1, N0, alt_t1);
        UNUSED CONST float s1 = evalQuadraticPoly(N2, N1, N0, t0);
        UNUSED CONST float s2 = evalQuadraticPoly(N2, N1, N0, t1);

        ts[0] = alt_t0;
        ts[1] = alt_t1;
    } else {
        ts[0] = t_min;
        ts[1] = t_max;
    }

    // Now pick the better of the two (quantized) solutions.
    CONST float t0 = ts[0];
    CONST float t1 = ts[1];
    // Are both t0 & t1 in the range of interest [t_min,t_max] exclusive?
    if (inRange(t_min, t0, t_max) && inRange(t_min, t1, t_max)) {
        // Yes, so take care to pick one of the two solutions...
        
        // Evaluate (one-third) gradient direction in Horner form for numerical accuracy
        CONST vec2 Q0 = evalQuadraticPoly(A, B, C, t0);
        CONST vec2 Q1 = evalQuadraticPoly(A, B, C, t1);

        // Rotate each gradient by theta.
        CONST float r0 = cos_theta*Q0[0] + sin_theta*Q0[1];
        CONST float r1 = cos_theta*Q1[0] + sin_theta*Q1[1];
        // Which rotate gradient has the larger magnitude?
        if (abs(r0) > abs(r1)) {
            // Pick the more positive solution.
            CONST float t = (r0 > 0) ? t0 : t1;
            return makeParametricDirection(t, theta, normal);
        } else {
            // Pick the more positive solution.
            CONST float t = (r1 > 0) ? t1 : t0;
            return makeParametricDirection(t, theta, normal);
        }
    } else {
        // None or one solution is in range...
        if (inRange(t_min, t0, t_max)) {
            // If t0 is in range, pick it!
            return makeParametricDirection(t0, theta, normal);
        } else if (inRange(t_min, t1, t_max)) {
            // If t1 is in range, pick it!
            return makeParametricDirection(t1, theta, normal);
        } else {
            // No solutions in valid ranges, likely due to unstable solutions.
            // Unstable solutions might be way outside the parameteric
            // range of Bezier interest to us so saturate to [0,1] range.
            CONST float saturated_t0 = saturate(t0);
            CONST float saturated_t1 = saturate(t1);

            // Now Horner evaluate where each (saturated solution) would be in gradient space.
            CONST vec2 Q0 = evalQuadraticPoly(A, B, C, saturated_t0);
            CONST vec2 Q1 = evalQuadraticPoly(A, B, C, saturated_t1);
            // Also Horner evaluate the bounds of the range of interest [t_min,t_max] are in gradient space.
            CONST vec2 Qmin = evalQuadraticPoly(A, B, C, t_min);
            CONST vec2 Qmax = evalQuadraticPoly(A, B, C, t_max);
     
            // Compute each square distance paring Q0 & Q1 to Qmin & Qmax.
            CONST float d_0min = sqdist(Q0-Qmin);
            CONST float d_1min = sqdist(Q1-Qmin);
            CONST float d_0max = sqdist(Q0-Qmax);
            CONST float d_1max = sqdist(Q1-Qmax);

#if !defined(NDEBUG)
            // For debugging, evaluate half-way point of range.
            CONST float t_avg = 0.5f*(t_min + t_max);
            CONST vec2 Qavg = evalQuadraticPoly(A, B, C, t_avg);
            UNUSED CONST float d_0avg = sqdist(Q0-Qavg);
            UNUSED CONST float d_1avg = sqdist(Q1-Qavg);
#endif
            
            // Compute the minimum square distance.
            CONST float d_min = min(min(d_0min, d_0max), min(d_1min, d_1max));
            // If d_0max or d_1max is closer...
            if (d_min == d_0max || d_min == d_1max) {
                return makeParametricDirection(t_max, theta, normal);
            } else {
                return makeParametricDirection(t_min, theta, normal);
            }
        }
    }
}

// From polar angle theta_s, generate a parametric t + gradient direction.
ParametricDirection StrokeTessellatorMETHOD evalCurveSample(uint ndx,
                                                            float delta_theta) CONST
{
    CONST int orientationBit = (1<<ndx);  // convert index to a single-bit mask
    CONST uint orientationBitMask = OrientationBitMask();
    CONST bool negate = (orientationBitMask & orientationBit) != 0;
    CONST float signed_delta_theta = negate ? -delta_theta : delta_theta;  // conditional negate
    CONST float theta_t = theta[ndx] + signed_delta_theta;
    CONST float t_min = tbounds[ndx];
    CONST float t_max = tbounds[ndx+1];

    return evalCurveSampleInRange(theta_t, t_min, t_max);
}

int StrokeTessellatorMETHOD thetaToIndex(float theta_s) CONST
{
    // With indexing (rather than branching) determine the range of interest
    // [t_min,t_max]] and theta_s within that range.
    if (theta_s <= delta[1]) {  // delta[1] holds D0A
        return 0;
    } else if (theta_s < delta[2]) {  // delta[2] holds D0AB
        return 1;
    } else {
        return 2;
    }
}

float StrokeTessellatorMETHOD indexedRibToUniformSpacedDeltaTheta(uint ndx, uint i) CONST
{
    CONST float stepping = i - delta_steps[ndx];
    return stepping*dtheta_dsteps[ndx];
}

#if HALF_STEPS
float StrokeTessellatorMETHOD indexedRibToDeltaThetaWithHalfSteps(uint ndx, uint i) CONST
{
    assert(useHalfSteps());
    assert(stashed_link_type == LT_CUBIC_SEGMENT || stashed_link_type == LT_QUADRATIC_SEGMENT);
    // Is this i at the beginning of a sequence of delta steps?
    if (i == delta_steps[ndx]+1) {
        // Yes, make it an initial half segment.
        CONST float half_step = 0.5f;
        CONST float stepping = half_step;
        return stepping*dtheta_dsteps[ndx];
    }
    // Is this i at the end of a sequence of delta steps?
    else if(i == delta_steps[ndx+1]-1) {
        // Yes, make it a final half segment.
        CONST float half_step = 0.5f;
        CONST int steps_for_half_step = 2;  // Subtract out both initial and final half-step.
        CONST float stepping = i - delta_steps[ndx] - steps_for_half_step + half_step;
        return stepping*dtheta_dsteps[ndx];
    } 
    // Otherwise in the middle...
    else {
        // Neither of the above so is an interior integral point in the segment strip.
        CONST int segment_for_half_step = 1;  // Subtract out the initial half-step.
        CONST float stepping = i - delta_steps[ndx] - segment_for_half_step;
        return stepping*dtheta_dsteps[ndx];
    }
}
#endif

float StrokeTessellatorMETHOD indexedRibToDeltaTheta(uint ndx, uint i) CONST
{
    assert(ndx < 3);

    return indexedRibToUniformSpacedDeltaTheta(ndx, i);
}

// Evaluate at a particular delta step point.  Such points are either
// end points (t=0 or 1) or inflection points.
ParametricDirection StrokeTessellatorMETHOD evalDeltaStepRibIndex(uint ndx) CONST
{
    CONST float t = tbounds[ndx];
    CONST float endPointTheta = theta[ndx];
    CONST vec2 gradient_direction = vec2(cos(endPointTheta),sin(endPointTheta));
    CONST vec2 normal = ortho(gradient_direction);
    return makeParametricDirection(t, endPointTheta, normal);
}

// Convert index i into a (t,theta,normal) assuming a conformal transform
// (so no non-uniform scaling or projection) so ignoring any such object-to-window transform.
// NOT using theta map.
ParametricDirection StrokeTessellatorMETHOD evalConformalRibIndex(uint i) CONST
{
    CONST uint ndx = i<=0 ? 0 : int(delta_steps[1] <= i) + int(delta_steps[2] <= i) + int(delta_steps[3] <= i);
    assert(ndx < 4);

    // Is this segment index an interval boundary?
    if (i == delta_steps[ndx]) {
        // Yes, use a tbounds value and skip evaluation.
        return evalDeltaStepRibIndex(ndx);
    } else {
        // Accounting for initial and final "half steps", find delta_theta within step index.
        CONST float delta_theta = indexedRibToDeltaTheta(ndx, i);
        return evalCurveSample(ndx, delta_theta);
    }
}

// Evaluate segment gradient direction & t at a particular rib index.
ParametricDirection StrokeTessellatorMETHOD evalRibIndex(uint i) CONST
{
    return evalConformalRibIndex(i);
}

vec2 StrokeTessellatorMETHOD startDirection() CONST
{
    CONST float startTheta = theta[0];  // theta0
    CONST vec2 gradient_direction = vec2(cos(startTheta),sin(startTheta));
    return gradient_direction;
}

vec2 StrokeTessellatorMETHOD endDirection() CONST
{
    CONST float endTheta = theta[3];  // theta1
    CONST vec2 gradient_direction = vec2(cos(endTheta),sin(endTheta));
    return gradient_direction;
}

#endif // __CommonStrokeTessellatorEval_hpp__
