
// CPUStrokeTessellatorEvalMain.cpp - CPU-centric "main" code for evaluation shader implementation of polar stroking geometry

// Stuff here will not appear in the GPU-centric version

#include <stdio.h>
#include <math.h>

#include "CPUStrokeTessellator.hpp"

#include <Cg/dot.hpp>
#include <Cg/min.hpp>
#include <Cg/mul.hpp>
#include <Cg/normalize.hpp>
#include <Cg/distance.hpp>
#include <Cg/determinant.hpp>
#include <Cg/isfinite.hpp>
#include <Cg/max.hpp>
#include <Cg/isnan.hpp>
#include <Cg/lerp.hpp>
#include <Cg/sign.hpp>
#include <Cg/upper2x2.hpp>
#include <Cg/inverse.hpp>
#include <Cg/iostream.hpp>

#include "countof.h"

#include "../common/CommonStrokeTessellatorPrivate.hpp"

// Make using iostream easier.
using std::cout;
using std::endl;

#include "../common/CommonStrokeTessellatorMath.hpp"

#include "../common/CommonStrokeTessellatorEvalMain.hpp"

bool CPUStrokeTessellator::transformIsAffine() const
{
    if (modelview_projection[3][0] == 0 &&
        modelview_projection[3][1] == 0 &&
        modelview_projection[3][3] == 1) {
        return true;
    } else {
        return false;
    }
}

float3 CPUStrokeTessellator::adjustEdge(const float3 Pedge,  // clip space
                                        const float3 Pmid,   // clip space
                                        const float3 p,      // object space version of Pmid
                                        const float2 effective_offset,
                                        float stroke_radius) const
{
    const float Pmid_w = Pmid.z;
    assert(Pmid_w > 0);  // Guaranteed before adjustEdge is called.
    const float Pedge_w = Pedge.z;
    if (Pedge_w > 0) {
        const float2 clip_to_window_scale = viewport_scale.xy;
        const float3 winPmid = float3(Pmid.xy*clip_to_window_scale, Pmid_w);
        const float3 winPedge = float3(Pedge.xy*clip_to_window_scale, Pedge_w);

        // Is the distance from generating curve to offset edges less than
        // window-space stroke radius threshold?
        const float d = distance(project(winPmid), project(winPedge));
        if (d < min_window_stroke_radius) {
            // Yes, needs to compute adjusted object-space edge position...
            const float scale = min_window_stroke_radius / d;
            const float dilated_stroke_radius = stroke_radius * scale;
            const float p_w = p.z;
            const float3 dilated_Pedge = float3(p.xy + dilated_stroke_radius*effective_offset, p_w);
            const float3 new_Pedge = mul(modelview_projection, dilated_Pedge);
            return new_Pedge;
        }
    }
    // No adjustment needed; return unchanged Pedge.
    return Pedge;
}

void StrokeTessellatorMETHOD stashOffset(uint i, vec2 offset)
{
    results->polar_offset.push_back(offset);
    assert(results->polar_offset[i].x == offset.x);
    assert(results->polar_offset[i].y == offset.y);
}

void StrokeTessellatorMETHOD stashNormal(uint i, vec2 normal)
{
    results->normal.push_back(normal);
    assert(results->normal[i].x == normal.x);
    assert(results->normal[i].y == normal.y);
}

void StrokeTessellatorMETHOD stashTheta(uint i, float theta)
{
    results->theta.push_back(theta);
}

void StrokeTessellatorMETHOD stashParametricT(uint i, float t)
{
    results->t.push_back(t);
}

void CPUStrokeTessellator::evalRib(uint i)
{
    // Expect evalRib to update these members (output variables 
    // in the evaluation shader!):
    //  * pervertex_radius
    //  * pervertex_offset
    //  * pervertex_p

    const RibIndexResult rib(i);
    evalRib(rib);

    float2 effective_offset = pervertex_offset;
    const float3 Pmid_obj = pervertex_p;

    const float Pmid_obj_w = Pmid_obj.z;  // Z is really W!
    effective_offset *= Pmid_obj_w;  // Scale effective_offset by p's homogeneous coordinate.
    const float3 Pneg_obj = float3(Pmid_obj.xy + pervertex_radius[0]*effective_offset, Pmid_obj_w);
    const float3 Ppos_obj = float3(Pmid_obj.xy + pervertex_radius[1]*effective_offset, Pmid_obj_w);

    // Save Ppos and Pneg (in object space) prior to transforming.
    results->p.push_back(Pmid_obj);
    results->Ppos_obj.push_back(Ppos_obj);
    results->Pneg_obj.push_back(Pneg_obj);
    // Is either pervertex_radius.x OR pervertex_radius.y zero?
    if (pervertex_radius.x*pervertex_radius.y == 0) {
        // This is a join or cap.
        results->frieze_height.push_back(0.5f - sign(pervertex_radius)*0.5f);
    } else {
        results->frieze_height.push_back(sign(pervertex_radius)*0.5f + 0.5f);
        // Could this be a round join?  If so, flip coordinates.
        if (pervertex_radius.x > pervertex_radius.y) {
            results->frieze_height.back() = 1 - results->frieze_height.back();
        }
    }

    float3 Pneg, Pmid, Ppos;
    if (stroke_space == WINDOW_SPACE) {
        Pmid = Pmid_obj * float3(viewport_scale.zw,1);
        Ppos = Ppos_obj * float3(viewport_scale.zw,1);
        Pneg = Pneg_obj * float3(viewport_scale.zw,1);
    } else {
        Pmid = mul(modelview_projection, Pmid_obj);
        Ppos = mul(modelview_projection, Ppos_obj);
        Pneg = mul(modelview_projection, Pneg_obj);
        // Homogenous coordinate "w" stored in z component of float3.
        const float Pmid_w = Pmid.z;
        if ((min_window_stroke_radius > 0) && (pervertex_radius[0] != 0) && (Pmid_w > 0)) {
            Pneg = adjustEdge(Pneg, Pmid, Pmid_obj, effective_offset, pervertex_radius[0]);
            Ppos = adjustEdge(Ppos, Pmid, Pmid_obj, effective_offset, pervertex_radius[1]);
        }
    }

    // Clip space positions.
    results->Pmid.push_back(Pmid);
    results->Ppos.push_back(Ppos);
    results->Pneg.push_back(Pneg);

    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 t = results->t[i];
    if (Cg::isnan(t)) {
        // A value of NaN for t indicates theta has no solution.  This can
        // happen for a hyperbola in perspective.
    } else {
        assert(t >= 0);
        assert(t <= 1);
    }

    // Gradient coefficients A, B, and C work for quadratic and cubic Bezier segments.
    const float2 gradient = float2(t*(t*A + B) + C);  // Horner quadratic evaluation
    results->gradient.push_back(gradient);
    const float2 gradient2 = float2(2*t*A + B);  // Second order gradient

    const float curvature_numer = gradient.x*gradient2.y - gradient.y*gradient2.x;
    const float curvature_denom = pow(sqdist(gradient), 1.5f);
    const float curvature = curvature_numer/curvature_denom;

    float3 Pradius_obj;
    float3 Pinfinite_obj;
    if (curvature < 0) {
        Pradius_obj = float3(Pmid_obj.xy * (-curvature) + effective_offset*Pmid_obj_w, Pmid_obj_w*(-curvature));
        Pinfinite_obj = float3(effective_offset*Pmid_obj_w, 0);
    } else {
        Pradius_obj = float3(Pmid_obj.xy * (curvature) - effective_offset*Pmid_obj_w, Pmid_obj_w*(curvature));
        Pinfinite_obj = float3(-effective_offset*Pmid_obj_w, 0);
    }
    const float3 Pradius = mul(modelview_projection, Pradius_obj);;
    const float3 Pinfinite = mul(modelview_projection, Pinfinite_obj);;

    results->gradient2.push_back(gradient2);
    results->curvature.push_back(curvature);
    results->Pradius.push_back(Pradius);
    results->Pinfinite.push_back(Pinfinite);
}

void CPUStrokeTessellator::verifyIncreasingT(StrokingResults& results) const
{
    for (size_t i=1; i<results.t.size(); i++) {
        // Are clip space positions in front of camera?
        if (results.Pmid[i].z > 0 && results.Pmid[i-1].z > 0) {
            // Yes, both positions have positive clip-space W values.

            if (transformIsAffine()) {
                // When the transform is affine, offset positions should have
                // positive clip-space W values too.
                assert(results.Ppos[i-1].z > 0);
                assert(results.Pneg[i-1].z > 0);
                assert(results.Ppos[i].z > 0);
                assert(results.Pneg[i].z > 0);
            }

            if (results.t[i-1] > results.t[i]) {
                const float diff = results.t[i-1] - results.t[i];
                // Small errors here mean the curve "backs up" slightly; if diff is sufficiently small,
                // such "back ups" will not be noticable.
                // Seems acceptable to be off by 1 part in 100,000.
                static float diff_threshold = 0.00001f;
                if (diff >= diff_threshold) {
                    printf("non-increasing t exceeded diff_threshold(%g)! [%d]=%g,[%d]=%g, diff=%g\n",
                        diff_threshold,
                        int(i-1), results.t[i-1],
                        int(i), results.t[i],
                        diff);
                    float percent = 100*float(i)/results.t.size();
                    printf("  array size = %u (percent through = %.3g%%)\n",
                        (unsigned int)results.t.size(), percent);
                }
            }
        } else {
#if 0 // XXX I don't think these are really valid under perspective
            if (results.Pmid[i-1].z > 0) {
                if (isfinite(results.Ppos[i-1].z)) {
                    assert(results.Ppos[i-1].z > 0);
                }
                if (isfinite(results.Pneg[i-1].z)) {
                    assert(results.Pneg[i-1].z > 0);
                }
            } else {
                if (isfinite(results.Pmid[i-1].z)) {
                    assert(results.Pmid[i-1].z <= 0);
                }
                if (isfinite(results.Ppos[i-1].z)) {
                    assert(results.Ppos[i-1].z <= 0);
                }
                if (isfinite(results.Pneg[i-1].z)) {
                    assert(results.Pneg[i-1].z <= 0);
                }
            }
            if (results.Pmid[i].z > 0) {
                if (isfinite(results.Ppos[i].z)) {
                    assert(results.Ppos[i].z > 0);
                }
                if (isfinite(results.Pneg[i].z)) {
                    assert(results.Pneg[i].z > 0);
                }
            } else {
                if (isfinite(results.Pmid[i].z)) {
                    assert(results.Pmid[i].z <= 0);
                }
                if (isfinite(results.Ppos[i].z)) {
                    assert(results.Ppos[i].z <= 0);
                }
                if (isfinite(results.Pneg[i].z)) {
                    assert(results.Pneg[i].z <= 0);
                }
            }
#endif
        }
    }
}

void CPUStrokeTessellator::evalRibs(StrokingResults& results)
{
    setStrokingResults(&results);

    const uint steps = Steps();
    for (uint rib_index=0; rib_index<=steps; rib_index++) {
        evalRib(rib_index);
    }

    setStrokingResults(NULL);

    verifyIncreasingT(results);
}

