
// CommonStrokeTessellatorCubic.hpp - common GPU and CPU polar stroking miter handling

#ifndef __CommonStrokeTessellatorCubic_hpp__
#define __CommonStrokeTessellatorCubic_hpp__

#include "CommonStrokeTessellatorPrivate.hpp"

#include "CommonStrokeTessellatorLinkType.hpp"

// Must be a method in order to access accessing stable_control_quadratic_solution.
vec2 StrokeTessellatorMETHOD findInflectionPoints(vec2 delta01, vec2 delta12, INOUTPARAM(bool) is_loop)
{
    // Convert control points to cubic monomial polynomial coefficients
#if AFFINE
#if 1  // use double precision
    VECTOR Ai = VECTOR(cp[3]) - VECTOR(cp[0]) - 3*VECTOR(delta12),
           Ci = VECTOR(delta01),
           Bi = VECTOR(delta12) - VECTOR(delta01);
#else  // use single-precision (can't enable due to bug 1800876)
    VECTOR Ai = VECTOR(cp[3] - cp[0] - 3*delta12),
           Ci = VECTOR(delta01),
           Bi = VECTOR(delta12 - delta01);
#endif
#else
#if 1  // use double precision
    VECTOR Ai = VECTOR(cp[1].z)*cp[2].z*cp[3].xy*cp[0].z - VECTOR(cp[1].z)*cp[2].z*cp[0].xy*cp[3].z - 3*VECTOR(delta12),
           Ci = VECTOR(delta01),
           Bi = VECTOR(delta12) - VECTOR(delta01);
#else  // use single-precision (can't enable due to bug 1800876)
    VECTOR Ai = VECTOR(cp[1].z*cp[2].z*cp[3].xy*cp[0].z - cp[1].z*cp[2].z*cp[0].xy*cp[3].z - 3*delta12),
           Ci = VECTOR(delta01),
           Bi = VECTOR(delta12 - delta01);
#endif
#endif  // AFFINE

    // "Inflection points of a cubic Bezier"
    // http://www.caffeineowl.com/graphics/2d/vectorial/cubic-inflexion.html
    // 0 = cross(P'(t),P"(t)) = Px*Pyy-Py*Pxx
    //
    // Developed in Maple:
    //
    // with(linalg):
    // > unprotect(D);
    // > A:=vector([Ax,Ay,0]):
    // > B:=vector([Bx,By,0]):
    // > C:=vector([Cx,Cy,0]):
    // > D:=vector([Dx,Dy,0]):
    // > P:=A*t^3+B*t^2+C*t+D;
    // P := A*t^3+B*t^2+C*t+D
    // > Pt:=(diff(P,t));Ptt:=(diff(P,t,t));
    // Pt := 3*A*t^2+2*B*t+C
    // Ptt := 6*A*t+2*B
    // > v:=crossprod(Pt,Ptt);
    // v := vector([0, 0, (3*t^2*Ax+2*t*Bx+Cx)*(6*t*Ay+2*By)-(3*t^2*Ay+2*t*By+Cy)*(6*t*Ax+2*Bx)])
    // > collect(v[3]/2,t);  % dividing by 2 just scales coefficients without changing polynomial's roots
    // (3*Bx*Ay-3*Ax*By)*t^2+(3*Cx*Ay-3*Cy*Ax)*t+Cx*By-Cy*Bx
    //
    // So the inflection points are at the solutions (roots) of
    //
    //   0 = (3*Bx*Ay-3*Ax*By)*t^2+(3*Cx*Ay-3*Cy*Ax)*t+Cx*By-Cy*Bx
    //     = a*t^2+b*t+c
    //
    // where a, b, & c are:
    //   a = 3*(Bx*Ay-Ax*By)
    //   b = 3*(Cx*Ay-Cy*Ax)
    //   c = Cx*By-Cy*Bx
    // These are determinants!  Scaled by 3 for a & b.

    // v1[0] = Ci
    // v1[1] = Bi
    // v1[2] = Ai
    // v1[3] = Ci
    PER_THREAD v1[TID] = Ci;  // Threads 0 and 3
    FOR_THREAD(1) v1[TID] = VECTOR(Bi);
    FOR_THREAD(2) v1[TID] = Ai;

#if ThreadSequential
    SCALAR a = determinant(Bi,Ai);
    SCALAR b = determinant(Ai,Ci);  // actually negative b coefficient
    SCALAR c = determinant(Ci,Bi);
#else  // ThreadParallel
    // v2[0] = c
    // v2[1] = a/3
    // v2[2] = b/3
    // v2[3] = -b/3
    PER_THREAD_COUNT(3) v2[TID] = determinant(v1[TID],v1[TID+1]);
    FOR_THREAD(3) v2[TID] = -v2[2];
#endif

    // Is the quadratic really a degenerate line?
    // If (Bx*Ay-Ax*By) is approximately 0...
    if (isAlmostZero(float(v2[1]))) {
        // Is the line really a degenerate point?
        // If (Cx*Ay-Cy*Ax) is approximately 0...
        if (v2[2] == 0) {
            // Is the situation super-degenerate due to the segment being a line with
            // duplicated control ponints for the start and end?
            if (v2[0] == 0) {
                // Yes, return zero valued solutions that will clamp to [0,1] range.
                return vec2(0);
            } else {
                // Look for either X or Y component reversing its derivative.
                //
                // Take second derivative of X & Y (colinear) equations...
                //
                // > eq:=Ax*(1-t)^3 + 3*Bx*(1-t)^2*t + 3*Cx*(1-t)*t^2 + Dx*t^3;
                //                  3               2                   2       3
                //  eq := Ax (1 - t)  + 3 Bx (1 - t)  t + 3 Cx (1 - t) t  + Dx t
                // > diff(diff(eq,t),t);
                //  6 Ax (1 - t) + 6 Bx t - 12 Bx (1 - t) - 12 Cx t + 6 Cx (1 - t) + 6 Dx t
                // > collect(%/6,t);
                // (-Ax + 3 Bx - 3 Cx + Dx) t + Ax + Cx - 2 Bx
#if AFFINE
#if ThreadSequential
                VECTOR AA = 9*Bi - Ai - 9*Ci + VECTOR(cp[0]);
                VECTOR BB = Ai + 3*Ci - 6*Bi;
                return vec2(-BB/AA);
#else  // ThreadParallel
                if (TID < 2) v2[TID] = -(v1[2][TID] + 3*v1[0][TID] - 6*v1[1][TID]) /
                                        (9*v1[1][TID] - v1[2][TID] - 9*v1[0][TID] + cp[0][TID]);
                return vec2(v2[0], v2[1]);
#endif
#else
#if ThreadSequential
                VECTOR AA = 9*Bi - Ai - 9*Ci + VECTOR(cp[0].xy)*cp[1].z*cp[2].z*cp[3].z;
                VECTOR BB = Ai + 3*Ci - 6*Bi;
                return vec2(-BB/AA);
#else  // ThreadParallel
                PER_THREAD if (TID < 2) v2[TID] = -(v1[2][TID] + 3*v1[0][TID] - 6*v1[1][TID]) /
                                                 (9*v1[1][TID] - v1[2][TID] - 9*v1[0][TID] + cp[1].z*cp[2].z*cp[0][TID]*cp[3].z);
                return vec2(v2[0], v2[1]);
#endif
#endif // AFFINE
            }
        } else {
            CONST float t = float(v2[0]/v2[2]);
            return vec2(t);
        }
    } else {
        // Exact cusp has a discriminant of zero; needs to be treated as zero solutions.
        CONST SCALAR discriminant = v2[2]*v2[2] - 4.0f*v2[1]*v2[0];
        if (discriminant > 0) {
#if 1  // use double precision
#if ThreadSequential
            CONST SCALAR sqrt_discr = sqrt(discriminant);
            CONST float ta = float((v2[2] - sqrt_discr)/(2.0f*v2[1]));
            CONST float tb = float((v2[2] + sqrt_discr)/(2.0f*v2[1]));
            return vec2(ta,tb);
#else  // ThreadParallel
            PER_THREAD {
                SCALAR sqrt_discr = sqrt(discriminant);
                if (TID > 0) sqrt_discr = -sqrt_discr;
                if (TID < 2) v2[TID] = float((v2[2] + sqrt_discr)/(2.0f*v2[1]));
            }
            return vec2(v2[0], v2[1]);
#endif
#else  // use single-precision
            CONST float sqrt_discr = sqrt(float(discriminant));
#if ThreadSequential
            float ta = (float(v2[2]) - sqrt_discr)/(2.0f*float(v2[1]));
            float tb = (float(v2[2]) + sqrt_discr)/(2.0f*float(v2[1]));
            return vec2(ta,tb);
#else  // ThreadParallel
            if (TID > 0) sqrt_discr = -sqrt_discr;
            if (TID < 2) v2[TID] = float(float(v2[2] + sqrt_discr)/(2.0f*float(v2[1])));
            return vec2(v2[0], v2[1]);
#endif
#endif
        } else {
            is_loop = true;
            return vec2(1,1);
        }
    }
}

STATIC bool closeToQuadratic(vec2 delta01, vec2 delta23, vec2 cp0, vec2 cp3)
{
    CONST vec2 D1 = 1.5f*delta01;
    CONST vec2 D2 = 1.5f*delta23;

    CONST vec2 P1 = cp0+D1;
    CONST vec2 P2 = cp3-D2;

    // Euclidean distance/length (L2 norm) involve computing square roots.
    // Instead use L1 norm, which is cheaper to calculate. It should produce
    // the same values within a factor of sqrt(2), which is close enough for
    // the ratio we compute.
    CONST vec2 absPdiff = abs(P1-P2);
    CONST float mag = absPdiff.x + absPdiff.y;

    // Compute sum of L1 norms (Manhattan distance) of D1 and D2.
    CONST float edges = abs(D1.x) + abs(D1.y) + abs(D2.x) + abs(D2.y);

    // 4096 is good enough for the 1984intro.svg to get most of its cubic curves reclassified.
    // Using <= makes sure mag & edges both exactly zero looks like a quadratic
    CONST bool closeEnough = (mag * 4096 <= edges);
    return closeEnough;
}

float StrokeTessellatorMETHOD angleDifference(float a, float b, float c, bool is_loop)
{
    if (is_loop) {
        return angleDelta(a, b, c);
    } else {
        return absoluteAngleDifference(a, b);
    }
}

bool StrokeTessellatorMETHOD polarCubicStroking(int link_type)
{
    // Having v5 & v4 is a hack since v[] is (logically) const.
    CPtype v4 = v[4];
    CPtype v5 = v[5];

    if (isJoinOrCap(link_type)) {
        if (isCap(link_type)) {
            v5 = v4;
        }
    }

    CONST vec2 C0 = cp[0].xy;
    CONST vec2 C3 = cp[3].xy;
    CONST vec2 C1alt = v4.xy;
    CONST vec2 C2alt = v5.xy;

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

    if (sqdist(D01alt) <= 0 || sqdist(D23alt) <= 0) {
        return true;  // Cull this patch!
    }

    // Control point differences.
    // FYI: the cubic's gradient control points are 3 times the control point differences.
    // But for our purposes, this uniform scaling by 3 is irrelevant to our computations.
    vec2 dC[3];
    for (uint i=0; i<3; i++) {
        CONST uint curr = i;
        CONST uint next = i+1;
#if AFFINE
        dC[curr] = cp[next].xy - cp[curr].xy;
#else
        CONST uint other1 = (i+2) % 4;
        CONST uint other2 = (i+3) % 4;
        dC[curr] = (cp[next].xy*cp[curr].z - cp[curr].xy*cp[next].z);
        dC[curr] *= cp[other1].z*cp[other2].z;
#endif
    }

    // hodograph.cpp accesses these through GradCp method for graphing.
    for (int i=0; i<3; i++) {
        grad_cp[i] = dC[i];
    }

    CONST vec2 D01 = dC[0];
    CONST vec2 D12 = dC[1];
    CONST vec2 D23 = dC[2];

#if AFFINE
    CONST bool is_masquerading_quadratic = closeToQuadratic(D01, D23,
        cp[0],
        cp[3]);
#else
    // Control point differences.
    // FYI: the cubic's gradient control points are 3 times the control point differences.
    // But for our purposes, this uniform scaling by 3 is irrelevant to our computations.
    CONST bool is_masquerading_quadratic = closeToQuadratic(D01, D23,
        cp[1].z*cp[2].z*cp[0].xy*cp[3].z,
        cp[1].z*cp[2].z*cp[3].xy*cp[0].z);
#endif

    // Compute monomial coefficients for (one-third) cubic Bezier gradient
    //
    // Maple math:
    //
    // > Tt:=collect((1-t)^2*D01 + 2*(1-t)*t*D12 + t^2*D23,t);
    //                               2
    //    Tt := (D01 - 2 D12 + D23) t  + (-2 D01 + 2 D12) t + D01
    //
    // > A:=coeff(Tt,t,2);
    //                    A := D01 - 2 D12 + D23
    // > B:=coeff(Tt,t,1);
    //                    B := -2 D01 + 2 D12
    // > C:=coeff(Tt,t,0);
    //                    C := D01
    //
    CONST vec2 B = 2*(D12 - D01);
    CONST vec2 C = D01;
    vec2 A;
    bool is_loop = false;
    if (is_masquerading_quadratic) {
        // Trick: force A to be zero when a quadratic is maseuerading as a cubic.
        // This ensures N2 quadratic coefficient is zero during evaluation, forcing a single linear solution. 
        A = vec2(0);
        CONST float theta0 = my_atan2(D01alt.y, D01alt.x);
        CONST float theta1 = my_atan2(D23alt.y, D23alt.x);

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

        FOR_THREAD(0) tbounds[TID] = 0;
        FOR_THREAD(1) tbounds[TID] = 1;
        FOR_THREAD(2) tbounds[TID] = 1;
        FOR_THREAD(3) tbounds[TID] = 1;
    } else {
        // Non-quadratic monomial coefficient for (one-third) cubic Bezier gradient 
        A = D01 - 2 * D12 + D23;

        // Solve for (two or zero) inflection points of cubic Bezier segment...
        CONST vec2 t = findInflectionPoints(D01, D12, is_loop);
        float ta = t[0];
        float tb = t[1];

        CONST float tA = saturate(min(ta,tb));  // max(0, min(1, min(ta,tb)))
        CONST float tB = saturate(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;

        // Stable Bezier evaluation of (saturated) inflection point positions.
        vec2 QA = (1-tA)*(1-tA)*D01 + 2*(1-tA)*tA*D12 + tA*tA*D23;
        vec2 QB = (1-tB)*(1-tB)*D01 + 2*(1-tB)*tB*D12 + tB*tB*D23;

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

        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);

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

    // A_C[0] = A.x
    // A_C[1] = C.x
    // A_C[2] = A.y
    // A_C[3] = C.y
#if PEDANTIC
    FOR_THREAD(0) A_C[TID] = A.x;
    FOR_THREAD(1) A_C[TID] = C.x;
    FOR_THREAD(2) A_C[TID] = A.y;
    FOR_THREAD(3) A_C[TID] = C.y;
#else
    PER_THREAD if ((TID & 1) == 0) A_C[TID] = A[TID>>1];
    PER_THREAD if ((TID & 1) != 0) A_C[TID] = C[TID>>1];
#endif

    // 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

    // Assignments below to delta[TID] could be writing when predicated by TID < 3,
    // but easier to skip predication as the current value of delta[3] is not used later.

    // Joins and masquerading cubics won't have gradient curvature reversals.
    if (isJoin(link_type) || is_masquerading_quadratic) {
        // Cross product of delta01alt & delta23alt is an oriented (signed) area.
        // Compute as a determinant
        CONST float join_curvature = determinant(D01alt, D23alt);

        PER_THREAD /*if (TID < 3)*/ delta[TID] = join_curvature;
    } else {
        // 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)

        // https://stackoverflow.com/questions/16140281/how-to-determine-curvature-of-a-cubic-bezier-path-at-an-end-point
        CONST vec2 ca = D01;
        CONST vec2 cb = D12 - ca;
        CONST vec2 cc = D23 - 2*cb - ca;

        CONST float bc = determinant(cb, cc);
        CONST float ac = determinant(ca, cc);
        CONST float ab = determinant(ca, cb);

        // r = ab + ac*t + bc*t*t;
        // Evaluate in Horner form for numerical accuracy
        FOR_THREAD(0) delta[TID] = (bc*t0A + ac)*t0A + ab;
        FOR_THREAD(1) delta[TID] = (bc*tAB + ac)*tAB + ab;
        FOR_THREAD(2) delta[TID] = (bc*tB1 + ac)*tB1 + ab;
    }

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

    FOR_THREAD(2) row_columns_orientationBitMask_steps[TID] = orientationBitMask;

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

    // 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;
    if (isJoinOrCap(link_type)) {
#if 0
        // Round joins shouldn't have too many steps.
        // XXX this should be adapted to screen space
        effective_max_turning_angle_delta *= 4;
#endif
    }

    float D0A_steps = ceil(D0A/effective_max_turning_angle_delta);
    float DAB_steps = ceil(DAB/effective_max_turning_angle_delta);
    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;

    // Specific to supporting triangular joins: keep triangular join from degenerating
    // to a bevel.
    if (D0A_steps > 0) {
        if (isJoin(link_type) && max_join_steps == 2) {
            D0A_steps = 3;
            assert(DAB_steps == 0);
            assert(DB1_steps == 0);
        }
    }

    // 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) {
        assert(link_type == LT_CUBIC_SEGMENT);
        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 // __CommonStrokeTessellatorCubic_hpp__
