
// CommonStrokeTessellatorQuadratic.hpp - common GPU and CPU polar stroking for quadratic Bezier segments

#ifndef __CommonStrokeTessellatorQuadratic_hpp__
#define __CommonStrokeTessellatorQuadratic_hpp__

#include "CommonStrokeTessellatorPrivate.hpp"
#include "CommonStrokeTessellatorMath.hpp"

#include <Cg/all.hpp>
#include <Cg/degrees.hpp>
#include <Cg/isfinite.hpp>

void StrokeTessellatorMETHOD passVectorCoefficients(vec2 A, vec2 B, vec2 C)
{
    // A_C[0] = A.x
    // A_C[1] = C.x
    // A_C[2] = A.y
    // A_C[3] = C.y
    PER_THREAD if ((TID & 1) == 0) A_C[TID] = A[TID>>1];
    PER_THREAD if ((TID & 1) != 0) A_C[TID] = C[TID>>1];

    // B_bogus[0] = B.x
    // B_bogus[1] = B.y
#if PEDANTIC
    FOR_THREAD(0) B_bogus[TID] = B.x;
    FOR_THREAD(1) B_bogus[TID] = B.y;
#else
    PER_THREAD B_bogus[TID] = B[TID];
#endif
}

STATIC bool essentiallyParabolic(float k)
{
    CONST float epsilon = 1/float(1048576); // 1 over about a million, or 0.00000095367431640625
    CONST bool essentially_parabolic = (abs(k-1) <= epsilon);
    if (essentially_parabolic) {
        if (k != 1) {
            printf("essentiallyParabolic but not actually parabolic\n");
        }
    }
    return essentially_parabolic;
}

STATIC bool isNotVeryNearlyOne(float x)
{
    const float epsilon = 0.000001f;
    bool result = (abs(x-1) >= epsilon);
    return result;
}

#if !AFFINE
void StrokeTessellatorMETHOD computeRationalQuadraticVectorCoefficients()
{
    //> Px:=((1-t)^2*P0x + 2*(1-t)*t*P1x + t^2*P2x)/
    //>      ((1-t)^2*w0 + 2*(1-t)*t*w1 + t^2*w2);
    //> Py:=((1-t)^2*P0y + 2*(1-t)*t*P1y + t^2*P2y)/
    //> ((1-t)^2*w0 + 2*(1-t)*t*w1 + t^2*w2);
    //
    //                         2                          2
    //                  (1 - t)  P0x + 2 (1 - t) t P1x + t  P2x
    //            Px := ---------------------------------------
    //                          2                        2
    //                   (1 - t)  w0 + 2 (1 - t) t w1 + t  w2
    //
    //
    //                         2                          2
    //                  (1 - t)  P0y + 2 (1 - t) t P1y + t  P2y
    //            Py := ---------------------------------------
    //                          2                        2
    //                   (1 - t)  w0 + 2 (1 - t) t w1 + t  w2
    //
    //> dPx:=radsimp(diff(Px,t));
    //> dPy:=radsimp(diff(Py,t));
    //
    //                                     2                  2
    //  dPx := -2 (-P1x w0 + P0x w1 + P1x t  w2 - t P2x w0 - t  P1x w0
    //
    //            2           2                            2           2
    //         + t  P2x w0 - t  P2x w1 - 2 P0x t w1 + P0x t  w1 - P0x t  w2
    //
    //                                    /
    //         + 2 t P1x w0 + P0x t w2)  /
    //                                  /
    //
    //                           2               2       2    2
    //        (w0 - 2 w0 t + w0 t  + 2 t w1 - 2 t  w1 + t  w2)
    //
    //
    //             2           2           2
    //  dPy := 2 (t  P1y w0 - t  P2y w0 + t  P2y w1 + 2 P0y t w1
    //
    //                2           2                        2
    //         - P0y t  w1 + P0y t  w2 - 2 t P1y w0 - P1y t  w2 + t P2y w0
    //
    //                                         /
    //         - P0y t w2 + P1y w0 - P0y w1)  /
    //                                       /
    //
    //                           2               2       2    2
    //        (w0 - 2 w0 t + w0 t  + 2 t w1 - 2 t  w1 + t  w2)
    //
    //> eq := sinT*dPx - cosT*dPy = 0;
    //
    //                                         2                  2
    //  eq := -2 sinT (-P1x w0 + P0x w1 + P1x t  w2 - t P2x w0 - t  P1x w0
    //
    //            2           2                            2           2
    //         + t  P2x w0 - t  P2x w1 - 2 P0x t w1 + P0x t  w1 - P0x t  w2
    //
    //                                    /
    //         + 2 t P1x w0 + P0x t w2)  /
    //                                  /
    //
    //                           2               2       2    2
    //        (w0 - 2 w0 t + w0 t  + 2 t w1 - 2 t  w1 + t  w2)  - 2 cosT (
    //
    //         2           2           2                            2
    //        t  P1y w0 - t  P2y w0 + t  P2y w1 + 2 P0y t w1 - P0y t  w1
    //
    //                2                        2
    //         + P0y t  w2 - 2 t P1y w0 - P1y t  w2 + t P2y w0 - P0y t w2
    //
    //                              /
    //         + P1y w0 - P0y w1)  /
    //                            /
    //
    //                           2               2       2    2
    //        (w0 - 2 w0 t + w0 t  + 2 t w1 - 2 t  w1 + t  w2)  = 0
    //
    //> poly := collect(numer(radsimp(lhs(eq)))/(2),t);
    //
    //  poly := (-sinT P1x w2 - cosT P1y w0 + cosT P2y w0 - sinT P2x w0
    //
    //         + cosT P1y w2 + sinT P2x w1 - sinT P0x w1 + sinT P1x w0
    //
    //                                                                   2
    //         + sinT P0x w2 + cosT P0y w1 - cosT P0y w2 - cosT P2y w1) t
    //
    //         + (-2 cosT P0y w1 - cosT P2y w0 + 2 sinT P0x w1
    //
    //         - sinT P0x w2 - 2 sinT P1x w0 + cosT P0y w2 + sinT P2x w0
    //
    //         + 2 cosT P1y w0) t + sinT P1x w0 - sinT P0x w1 + cosT P0y w1
    //
    //         - cosT P1y w0
    //
    //> A := collect(coeff(poly,t,2),[sinT,cosT]);
    //
    //  A := (-P1x w2 - P2x w0 + P0x w2 + P2x w1 - P0x w1 + P1x w0) sinT
    //
    //         + (P1y w2 - P1y w0 + P2y w0 - P2y w1 + P0y w1 - P0y w2) cosT
    //
    //> B := collect(coeff(poly,t,1),[sinT,cosT]);
    //
    //  B := (2 P0x w1 - 2 P1x w0 + P2x w0 - P0x w2) sinT
    //
    //         + (-2 P0y w1 - P2y w0 + 2 P1y w0 + P0y w2) cosT
    //
    //> C := collect(coeff(poly,t,0),[sinT,cosT]);
    //
    //         C := (P1x w0 - P0x w1) sinT + (P0y w1 - P1y w0) cosT
    //
    //> Ax := factor(op(1,A)/(sinT));
    //> Ay := factor(op(2,A)/(-cosT));
    //
    //      Ax := -P1x w2 - P2x w0 + P0x w2 + P2x w1 - P0x w1 + P1x w0
    //
    //
    //      Ay := P1y w0 - P2y w0 - P1y w2 - P0y w1 + P0y w2 + P2y w1
    //
    //> Bx := collect(op(1,B)/sinT,[P0x]);
    //> By := collect(op(2,B)/(-cosT),[P0y]);
    //
    //              Bx := (2 w1 - w2) P0x - 2 P1x w0 + P2x w0
    //
    //
    //              By := (2 w1 - w2) P0y + P2y w0 - 2 P1y w0
    //
    //> Cx := op(1,C)/(sinT);
    //> Cy := op(2,C)/(-cosT);
    //
    //                        Cx := P1x w0 - P0x w1
    //
    //
    //                        Cy := -P0y w1 + P1y w0
    //
    CONST vec2 P0 = cp[0].xy;
    CONST vec2 P1 = cp[1].xy;
    assert(all(P1 == cp[2].xy));
    CONST vec2 P2 = cp[3].xy;

    CONST float w0 = cp[0].z;
    CONST float w1 = cp[1].z;
    assert(w1 == cp[2].z);
    CONST float w2 = cp[3].z;

    CONST vec2 A = (P1-P2)*w0 + (P2-P0)*w1 + (P0-P1)*w2;
    CONST vec2 B = (P2-2*P1)*w0 + P0*(2*w1-w2);
    CONST vec2 C = P1*w0 - P0*w1;

    passVectorCoefficients(A, B, C);
}
#endif

// Version of computeRationalQuadraticVectorCoefficients specialized for
// w0=1, w1=k where k>0, w2=1 (so no projective transform applied).
void StrokeTessellatorMETHOD computeConicVectorCoefficients(float k)
{
    // Yes, conics require solving a quadratic equation during evaluation...
    assert(k > 0);  // When k<=0, generality of computeRationalQuadraticVectorCoefficients is needed.
    assert(isNotVeryNearlyOne(k));  // Should have been demoted to a non-rational quadratic if k==1.

    CONST vec2 P0 = cp[0].xy;
    CONST vec2 P1 = cp[1].xy;
    CONST vec2 P2 = cp[3].xy;

    // Intersect a normalized rational quadratic Bezier segment P(t) with a polar line.
    // 
    // A polar line is a line through the origin at angle theta:
    //
    //  0 = sin(theta) * Px(t) - cos(theta) * Py(t)
    //
    // Solving for t results in a quadratic polynomial that for a given angle theta can be solved for t.
    
    // The coefficients (A,B,C) of the quadratic polynomial are observed to be "rotations" by theta
    // of the form:
    //
    //   A = sin(theta)*Ax - cos(theta)*Ay
    //   B = sin(theta)*Bx - cos(theta)*By
    //   C = sin(theta)*Cx - cos(theta)*Cy
    // 
    // Seperating out the vectors [Ax,Ay], [Bx,By], [Cx,Cy] allows us to construct quickly
    // coefficients for any theta.
    //
    // Later in the evaluation process, we can construct coefficients A, B, & C from uniformly
    // steps in theta and solve for t at each theta step.

    // Maple development...
    //
    //> Px:=((1-t)^2*P0x + 2*(1-t)*t*P1x + t^2*P2x)/
    //> ((1-t)^2 + 2*(1-t)*t*k + t^2);
    //> Py:=((1-t)^2*P0y + 2*(1-t)*t*P1y + t^2*P2y)/
    //> ((1-t)^2 + 2*(1-t)*t*k + t^2);
    //
    //                         2                          2
    //                  (1 - t)  P0x + 2 (1 - t) t P1x + t  P2x
    //            Px := ---------------------------------------
    //                              2                    2
    //                       (1 - t)  + 2 (1 - t) t k + t
    //
    //
    //                         2                          2
    //                  (1 - t)  P0y + 2 (1 - t) t P1y + t  P2y
    //            Py := ---------------------------------------
    //                              2                    2
    //                       (1 - t)  + 2 (1 - t) t k + t
    //
    //> dPx:=radsimp(diff(Px,t));
    //> dPy:=radsimp(diff(Py,t));
    //
    //               2                          2          2
    //  dPx := -2 (-t  P2x k - 2 P0x t k + P0x t  k - P0x t  + P0x k
    //
    //                                            2        /
    //         + P0x t - P1x + 2 t P1x - t P2x + t  P2x)  /
    //                                                   /
    //
    //                       2              2   2
    //        (-1 + 2 t - 2 t  - 2 t k + 2 t  k)
    //
    //
    //               2                          2                  2
    //  dPy := -2 (-t  P2y k - 2 P0y t k + P0y t  k + P0y t - P0y t
    //
    //                                            2        /
    //         + P0y k - P1y - t P2y + 2 t P1y + t  P2y)  /
    //                                                   /
    //
    //                       2              2   2
    //        (-1 + 2 t - 2 t  - 2 t k + 2 t  k)
    //
    //> eq := sinT*dPx - cosT*dPy = 0;
    //
    //                   2                          2          2
    //  eq := -2 sinT (-t  P2x k - 2 P0x t k + P0x t  k - P0x t  + P0x k
    //
    //                                            2        /
    //         + P0x t - P1x + 2 t P1x - t P2x + t  P2x)  /
    //                                                   /
    //
    //                       2              2   2             2
    //        (-1 + 2 t - 2 t  - 2 t k + 2 t  k)  + 2 cosT (-t  P2y k
    //
    //                            2                  2
    //         - 2 P0y t k + P0y t  k + P0y t - P0y t  + P0y k - P1y
    //
    //                              2        /
    //         - t P2y + 2 t P1y + t  P2y)  /
    //                                     /
    //
    //                       2              2   2
    //        (-1 + 2 t - 2 t  - 2 t k + 2 t  k)  = 0
    //
    //> poly := collect(numer(radsimp(lhs(eq)))/(2),t);
    //
    //  poly := (cosT P0y k - sinT P2x - cosT P2y k - sinT P0x k + sinT P0x
    //
    //                                              2
    //         + sinT P2x k - cosT P0y + cosT P2y) t  + (2 sinT P0x k
    //
    //         + cosT P0y - sinT P0x - 2 sinT P1x + sinT P2x - cosT P2y
    //
    //         + 2 cosT P1y - 2 cosT P0y k) t - sinT P0x k + cosT P0y k
    //
    //         - cosT P1y + sinT P1x
    //
    //> A := collect(coeff(poly,t,2),[sinT,cosT]);
    //> 
    //
    //  A := (-P2x + P0x + P2x k - P0x k) sinT
    //
    //         + (-P2y k + P0y k - P0y + P2y) cosT
    //
    //> B := collect(coeff(poly,t,1),[sinT,cosT]);
    //> 
    //
    //  B := (2 P0x k + P2x - P0x - 2 P1x) sinT
    //
    //         + (2 P1y + P0y - 2 P0y k - P2y) cosT
    //
    //> C := collect(coeff(poly,t,0),[sinT,cosT]);
    //
    //            C := (-P0x k + P1x) sinT + (P0y k - P1y) cosT
    //
    //> Ax := factor(op(1,A)/(sinT));
    //> Ay := factor(op(2,A)/(-cosT));
    //
    //                      Ax := -(k - 1) (P0x - P2x)
    //
    //
    //                     Ay := -(-P2y + P0y) (k - 1)
    //
    //> Bx := collect(op(1,B)/sinT,[P0x]);
    //> By := collect(op(2,B)/(-cosT),[P0y]);
    //
    //                  Bx := (2 k - 1) P0x + P2x - 2 P1x
    //
    //
    //                  By := (2 k - 1) P0y - 2 P1y + P2y
    //
    //> Cx := op(1,C)/(sinT);
    //> Cy := op(2,C)/(-cosT);
    //
    //                          Cx := -P0x k + P1x
    //
    //
    //                          Cy := -P0y k + P1y

    CONST vec2 A = (k-1)*(P2 - P0);           // A = [Ax,Ay] above
    CONST vec2 B = (2*k-1)*P0 - 2*P1 + P2;    // B = [Bx,By] above
    CONST vec2 C = P1 - P0*k;                 // C = [Cx,Cy] above

#if AFFINE
    assert(k > 0);
#else
    if (k < 0) {
        // When k<0, the orientation of the conic reverses and its gradient
        // vectors are reflected through the origin so negate the gradient
        // coefficients.
        passVectorCoefficients(-A, -B, -C);
        return;
    }
#endif

    passVectorCoefficients(A, B, C);
}

void StrokeTessellatorMETHOD computeParabolicVectorCoefficients(CONST vec2 delta01alt, CONST vec2 delta23alt)
{
    // The same Maple development as in computeConicVectorCoefficients specifialized for case of k==1...

    // However if nearly a parabola (k not one but very close to one), the quadratic
    // vector coefficient A might be nearly but not quite a zero vector
    // and lead to unstable quadratic roots.
    CONST vec2 A = vec2(0);
    CONST vec2 B = delta23alt - delta01alt;
    CONST vec2 C = delta01alt;

    passVectorCoefficients(A, B, C);
}

bool StrokeTessellatorMETHOD polarAffineParabolicStroking()
{
    CONST CPtype v4 = v[4];
    CONST CPtype v5 = v[5];

    CONST vec2 delta01alt = v4.xy - cp[0].xy;
    CONST vec2 delta23alt = cp[3].xy - v5.xy;

    if (sqdist(delta01alt) <= 0 || sqdist(delta23alt) <= 0) {
        return true;
    }

    CONST float affine_orientation = determinant(delta01alt, delta23alt);
    CONST float orientation = affine_orientation;

    PER_THREAD row_columns_orientationBitMask_steps[TID] = (orientation < 0) ? 1 : 0;

    // theta[0] = delta01alt.x
    // theta[1] = delta23alt.x
    PER_THREAD /*if (TID==0)*/ theta[TID] = delta01alt.x;
    PER_THREAD if (TID>=1) theta[TID] = delta23alt.x;

    // delta_alt[0] = delta01alt.y
    // delta_alt[1] = delta23alt.y
    // delta_alt[2] = delta23alt.y
    // delta_alt[3] = delta23alt.y
    PER_THREAD /*if (TID==0)*/ delta_alt[TID] = delta01alt.y;
    PER_THREAD if (TID>=1) delta_alt[TID] = delta23alt.y;

    // THREAD PARALLEL, two unique arctangents done at once!
    // theta[0] == theta0
    // theta[1] == theta1
    // theta[2] == theta1
    // theta[3] == theta1
    PER_THREAD theta[TID] = my_atan2(delta_alt[TID], theta[TID]);

    // tbounds[0] = 0
    // tbounds[1] = 1
    // tbounds[2] = 1
    // tbounds[3] = 1
    PER_THREAD tbounds[TID] = float(TID>=1 ? 1 : 0);

    CONST float D01 = angleDelta(theta[0],theta[1],orientation);
    // delta[0] = 0
    // delta[1] = D01
    // delta[2] = D01
    // delta[3] = max(1,D01)
    FOR_THREAD(0) delta[TID] = 0;
    PER_THREAD if (TID >= 1) delta[TID] = D01;
    stashed_D0AB1 = D01;  // XXX fix me! perspective stroking needs this

    computeParabolicVectorCoefficients(delta01alt, delta23alt);

    // XXX deal with angleDelta_stash
    FOR_THREAD(0) angleDelta_stash[TID] = D01;
    FOR_THREAD(1) angleDelta_stash[TID] = 0;
    FOR_THREAD(2) angleDelta_stash[TID] = 0;

    CONST float effective_max_turning_angle_delta = max_turning_angle_delta;

    float D01_steps = ceil(D01/effective_max_turning_angle_delta);

#if SUPPORT_POLAR_VISUALIZE
    if (D01_steps > 0) {
        D01_steps += extra_steps;
    }
#endif

    FOR_THREAD(0) dtheta_dsteps[TID] = D01/D01_steps;
    const float bogus_unreferenced_value = -666;
    FOR_THREAD(1) dtheta_dsteps[TID] = bogus_unreferenced_value;
    FOR_THREAD(2) dtheta_dsteps[TID] = bogus_unreferenced_value;
    FOR_THREAD(3) dtheta_dsteps[TID] = bogus_unreferenced_value;

    // XXX quadratics don't really need half steps, but this code works.
    CONST bool use_half_steps = useHalfSteps();
    if (use_half_steps) {
        CONST float half_effective_max_turning_angle_delta = effective_max_turning_angle_delta*0.5f;
        if (D01 > half_effective_max_turning_angle_delta) {
            D01_steps += 2;
        }
    }

    FOR_THREAD(0) delta_steps[TID] = 0;
    FOR_THREAD(1) delta_steps[TID] = D01_steps;
    FOR_THREAD(2) delta_steps[TID] = D01_steps;
    CONST float total_steps = D01_steps;
    // total_steps of 0 is a line (or point) but lines need at least one step.
    FOR_THREAD(3) delta_steps[TID] = total_steps > 1 ? total_steps : 1;

    return false;  // Don't discard patch.
}

#if !AFFINE
bool StrokeTessellatorMETHOD polarRationalQuadraticStroking()
{
    CONST CPtype v4 = v[4];
    CONST CPtype v5 = v[5];

    CONST vec2 delta01alt = (v4.xy*cp[0].z - cp[0].xy*v4.z)*cp[3].z*v5.z;
    CONST vec2 delta23alt = (cp[3].xy*v5.z - v5.xy*cp[3].z)*cp[0].z*v4.z;

    if (sqdist(delta01alt) <= 0 || sqdist(delta23alt) <= 0) {
        return true;
    }

    CONST float affine_orientation = determinant(delta01alt, delta23alt);
    CONST float orientation = (v[4].z >= 0) ? affine_orientation : -affine_orientation;

    // Rational quadratic (ellipses, hyperbolas) require solving a quadratic equation during evaluation...
    computeRationalQuadraticVectorCoefficients();

    CONST float w0 = cp[0].z;
    CONST float w1 = cp[1].z;
    assert(w1 == cp[2].z);
    CONST float w2 = cp[3].z;

    CONST float a = w0-2*w1+w2;
    CONST float b = 2*(w1-w0);
    CONST float c = w0;

    CONST float discriminant = b*b - 4*a*c;
    if (discriminant > 0) {
        // Solve quadratic equation
        //   0 = A^2*t + B*t + C
        // 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 * (b + (b>=0 ? +sqrt_discr : -sqrt_discr));
        CONST float ta = q/a;
        CONST float tb = c/q;

        // These rotation reversals may be outside curve's [0,1] interval so saturate to [0,1].
        CONST float tA = saturate(safe_min(ta,tb));  // max(0, min(1, min(ta,tb)))
        CONST float tB = saturate(safe_max(ta,tb));  // max(0, min(1, max(ta,tb)))

        // Now inflection parametric values sorted in parametric [0,1] range:
        FOR_THREAD(0) tbounds[TID] = 0;
        FOR_THREAD(1) tbounds[TID] = tA;
        FOR_THREAD(2) tbounds[TID] = tB;
        FOR_THREAD(3) tbounds[TID] = 1;
    } else {
        FOR_THREAD(0) tbounds[TID] = 0;

        // split_t_A & split_t_B indicate where in the curve's parametric space to split the interval for one
        // of two reasons:
        // 1) the curve spans more than 180 degrees of gradient angle change and we need to be able to bound the
        //    solution for the polar line intersecting the curve's gradient to a hemicicle, or
        // 2) the ratiaonl quadartic curve has an undefined gradient where the denominator of its gradient
        //    becomes zero.  We split the interval around these one or two points to avoid an undefined gradient
        //    due to division by zero.
        //
        // Is there a t where the denominator of the gradient components of the rational quadratic curve is zero?
        CONST float gradient_denominator_discriminant = w1*w1 - w0*w2;  // think: B^2 - 4*A*C
        if (gradient_denominator_discriminant >= 0) {
            // So handle reason #2...

            // By solving a quadratic equation, determine the two t values 
            CONST float sqrt_gradient_denominator_discriminant = sqrt(gradient_denominator_discriminant);
            CONST float wtn = (w0-w1-sqrt_gradient_denominator_discriminant)/(w0-2*w1+w2);
            CONST float wtp = (w0-w1+sqrt_gradient_denominator_discriminant)/(w0-2*w1+w2);

            CONST float split_epsilon = 0.00001f;
            CONST float tA = saturate(safe_min(wtn,wtp) - split_epsilon);
            CONST float tB = saturate(safe_max(wtn,wtp) + split_epsilon);

            FOR_THREAD(1) tbounds[TID] = tA;
            FOR_THREAD(2) tbounds[TID] = tB;
        } else {
            // Otherwise reason #1 applies...
            FOR_THREAD(1) tbounds[TID] = 0.5f;
            FOR_THREAD(2) tbounds[TID] = 0.5f;
        }

        FOR_THREAD(3) tbounds[TID] = 1;
    }

    // Compute mid-point of each interval
    // This should be a "stable" point to compute the signed curvature of that interval
    // Intervals of zero or nearly zero length should also have insignifiant curvature
    CONST float t0A = 0.5f*tbounds[1];              // Logically 0.5f*(0 + tA)
    CONST float tAB = 0.5f*(tbounds[1] + tbounds[2]);
    CONST float tB1 = 0.5f*tbounds[2] + 0.5f;  // Logically 0.5f*(tB + 1)

    FOR_THREAD(0) delta[TID] = (a*t0A + b)*t0A + c;
    FOR_THREAD(1) delta[TID] = (a*tAB + b)*tAB + c;
    FOR_THREAD(2) delta[TID] = (a*tB1 + b)*tB1 + c;

    // Advantageous to store the signed curvature sign as a bitmask to convey less data.
    uint orientationBitMask = 0;
    if (delta[0]*orientation < 0) {
        orientationBitMask |= 1;
    }
    if (delta[1]*orientation < 0) {
        orientationBitMask |= 2;
    }
    if (delta[2]*orientation < 0) {
        orientationBitMask |= 4;
    }

    FOR_THREAD(2) row_columns_orientationBitMask_steps[TID] = orientationBitMask;

    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 tA = tbounds[1];
    CONST float tB = tbounds[2];
    // XXX consider stable Bezier evaluation
    vec2 QA = (A*tA + B)*tA + C;
    vec2 QB = (A*tB + B)*tB + C;

    CONST vec2 C0 = cp[0].xy;
    CONST vec2 C3 = cp[3].xy;
    CONST vec2 C1alt = v4.xy;
    CONST vec2 C2alt = v5.xy;
    // Homogeneous difference
    vec2 D01alt = C1alt*cp[0].z - C0*v4.z;
    vec2 D23alt = C3*v5.z - C2alt*cp[3].z;

    if (tA == 0) {
        QA = D01alt;
    }
    if (tB == 0) {
        QB = D01alt;
    }
    if (tA == 1) {
        QA = D23alt;
    }
    if (tB == 1) {
        QB = D23alt;
    }

    CONST vec2 parabola_direction = cp[0].xy - 2*cp[1].xy + cp[3].xy;
    assert(dot(parabola_direction,parabola_direction) > 0);
    if (dot(D01alt,D01alt) == 0) {
        D01alt = parabola_direction;
    }
    if (dot(QA,QA) == 0) {
        QA = parabola_direction;
    }
    if (dot(QB,QB) == 0) {
        QB = parabola_direction;
    }
    if (dot(D23alt,D23alt) == 0) {
        D23alt = parabola_direction;
    }

    // XXX run this atan2 code across threads.
    CONST float theta0 = my_atan2(D01alt.y, D01alt.x);
    CONST float thetaA = my_atan2(QA.y, QA.x);
    CONST float thetaB = my_atan2(QB.y, QB.x);
    CONST float theta1 = my_atan2(D23alt.y, D23alt.x);

    // Degenerate normals grief atan2.
    assert(Cg::isfinite(theta0));
    assert(Cg::isfinite(thetaA));
    assert(Cg::isfinite(thetaB));
    assert(Cg::isfinite(theta1));

    FOR_THREAD(0) theta[TID] = theta0;
    FOR_THREAD(1) theta[TID] = thetaA;
    FOR_THREAD(2) theta[TID] = thetaB;
    FOR_THREAD(3) theta[TID] = theta1;

    CONST float D0A = absoluteAngleDifference(theta[0], theta[1]);
    CONST float DAB = absoluteAngleDifference(theta[1], theta[2]);
    CONST float DB1 = absoluteAngleDifference(theta[2], theta[3]);

    // XXX deal with angleDelta_stash
    FOR_THREAD(0) angleDelta_stash[TID] = D0A;
    FOR_THREAD(1) angleDelta_stash[TID] = DAB;
    FOR_THREAD(2) angleDelta_stash[TID] = DB1;

    float effective_max_turning_angle_delta = max_turning_angle_delta;

    float D0A_steps = ceil(D0A/effective_max_turning_angle_delta);
    float DAB_steps = ceil(DAB/effective_max_turning_angle_delta);
#if 0 // XXX this seemed like a good idea but does not work
    // Does the conic have a middle "negative w" section?
    if (discriminant > 0 && w1 < 0) {
        // Yes, pointless to draw this as it will be discarded so force steps to zero.
        DAB_steps = 0;
    }
#endif
    float DB1_steps = ceil(DB1/effective_max_turning_angle_delta);

    CONST float D0AB = D0A + DAB;
    CONST float D0AB1 = D0A + DAB + DB1;
    stashed_D0AB1 = D0AB1;

    FOR_THREAD(0) delta[TID] = 0;
    FOR_THREAD(1) delta[TID] = D0A;
    FOR_THREAD(2) delta[TID] = D0AB;
    FOR_THREAD(3) delta[TID] = D0AB1;

    // Increment each positive interval's steps by 1 so facet angle is better bounded.
    //
    // The facet angle closely approximates the effective_max_turning_angle_delta
    // but the facet angle can and does fluctuate above and below
    // effective_max_turning_angle_delta.  Adding extra steps, typically 1 step, is
    // sufficient to reduces the situations where the facet angle noticable
    // exceeds effective_max_turning_angle_delta.
    //
#if SUPPORT_POLAR_VISUALIZE
    if (D0A_steps > 0) {
        D0A_steps += extra_steps;
    }
    if (DAB_steps > 0) {
        DAB_steps += extra_steps;
    }
    if (DB1_steps > 0) {
        DB1_steps += extra_steps;
    }
#endif

    // Avoid hitting the max_steps hardware limit by adjusting D0A_steps,
    // DAB_steps, and DB1_steps so their sum (plus an extra steps that
    // might be added for HALF_STEPS), will not exceed the GPU max_steps
    // limit.
#if HALF_STEPS
    float worst_case_step_expansion = 2*3;
#else
    float worst_case_step_expansion = 0;
#endif
    float steps_upper_bound = D0A_steps + DAB_steps + DB1_steps + worst_case_step_expansion;
    // A single iteration might not get the total steps (plus expansion)
    // under max_steps so repeat until it does.
    int loops = 0;
    while (steps_upper_bound > max_steps) {
        CONST float step_scale_back = max_steps/steps_upper_bound;
        D0A_steps = floor(D0A_steps * step_scale_back);
        DAB_steps = floor(DAB_steps * step_scale_back);
        DB1_steps = floor(DB1_steps * step_scale_back);
        steps_upper_bound = D0A_steps + DAB_steps + DB1_steps + worst_case_step_expansion;
        loops++;
    }
    assert(loops <= 2);  // Two iterations has always been sufficient.

    // Compute dtheta_dsteps before adjusting D*_steps for half segments.
    FOR_THREAD(0) dtheta_dsteps[TID] = (D0A_steps > 0) ? D0A/D0A_steps : 0;
    FOR_THREAD(1) dtheta_dsteps[TID] = (DAB_steps > 0) ? DAB/DAB_steps : 0;
    FOR_THREAD(2) dtheta_dsteps[TID] = (DB1_steps > 0) ? DB1/DB1_steps : 0;
    const float bogus_unreferenced_value = -666;
    FOR_THREAD(3) dtheta_dsteps[TID] = bogus_unreferenced_value;

    assert(Cg::isfinite(dtheta_dsteps[0]));
    assert(Cg::isfinite(dtheta_dsteps[1]));
    assert(Cg::isfinite(dtheta_dsteps[2]));
    assert(Cg::isfinite(dtheta_dsteps[3]));

#if HALF_STEPS
    // An exception where the facet angle can significantly (by as much as ~30%)
    // exceed effective_max_turning_angle_delta is near the two inflection points
    // of a cubic curve.
    //
    // The zero curvature at an inflection introduces a significant "unbalance"
    // at steps immediately adjacent to such inflections.
    //
    // The algorithm here compensates for this "overshooting" of the facet angle
    // near inflection points on a cubic by introducing a half-step on either
    // side of each positive interval.
    //
    // Technically introducing a half step just adjacent to the inflection would
    // be sufficient but it becomes more involved/expensive in the evaluation
    // shader to identify such a asymmetric half-steps so we add a symetric half-step
    // on either end of the interval.
    // 
    // Adjust each interval's steps for a initial and final half step.
    //
    CONST bool use_half_steps = useHalfSteps();
    if (use_half_steps) {
        CONST float half_effective_max_turning_angle_delta = effective_max_turning_angle_delta*0.5f;
        if (D0A > half_effective_max_turning_angle_delta) {
            D0A_steps += 2;
        }
        if (DAB > half_effective_max_turning_angle_delta) {
            DAB_steps += 2;
        }
        if (DB1 > half_effective_max_turning_angle_delta) {
            DB1_steps += 2;
        }
    }
#endif

    FOR_THREAD(0) delta_steps[TID] = 0;
    FOR_THREAD(1) delta_steps[TID] = D0A_steps;
    FOR_THREAD(2) delta_steps[TID] = D0A_steps + DAB_steps;
    CONST float total_steps = D0A_steps + DAB_steps + DB1_steps;
    // total_steps of 0 is a line (or point) but lines need at least one step.
    FOR_THREAD(3) delta_steps[TID] = total_steps > 1 ? total_steps : 1;

    return false;  // Don't discard patch.
}
#endif

bool StrokeTessellatorMETHOD polarQuadraticStroking()
{
    if (cp[0].z == 1 && cp[3].z == 1 && v[4].z == 1) {
        // Other control points that should also be unity because they are duplicates.
        assert(cp[1].z == 1);
        assert(v[5].z == 1);
        return polarAffineParabolicStroking();
    } else {
        return polarRationalQuadraticStroking();
    }
}

#endif // __CommonStrokeTessellatorQuadratic_hpp__
