
// simple_path.cpp - test bed for faster GPU-based algorithms to draw stroked paths

#include <assert.h>
#include <stdio.h>
#include <stdarg.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <GL/glew.h>
#ifdef __APPLE__
# include <GLUT/glut.h>
#else
# include <GL/glut.h>
# ifdef _WIN32
#  include <windows.h>  // for wglGetProcAddress
# else
#  include <GL/glx.h>  // for glXGetProcAddress
# endif
#endif

#include <vector>

#include <Cg/double.hpp>
#include <Cg/vector/xyzw.hpp>
#include <Cg/matrix/1based.hpp>
#include <Cg/matrix/rows.hpp>
#include <Cg/vector.hpp>
#include <Cg/matrix.hpp>
#include <Cg/mul.hpp>
#include <Cg/distance.hpp>
#include <Cg/determinant.hpp>
#include <Cg/length.hpp>
#include <Cg/dot.hpp>
#include <Cg/isfinite.hpp>
#include <Cg/normalize.hpp>
#include <Cg/round.hpp>
#include <Cg/upper2x2.hpp>
#include <Cg/iostream.hpp>
#include <Cg/inverse.hpp>
#include <Cg/degrees.hpp>
#include <Cg/radians.hpp>
#include <Cg/saturate.hpp>
#include <Cg/sqrt.hpp>
#include <Cg/transpose.hpp>
#include <Cg/any.hpp>
#include <Cg/all.hpp>

using namespace Cg;
using std::vector;

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

#include "sRGB_vector.hpp"
#include "countof.h"
#include "request_vsync.h"
#include "showfps.h"
#include "cg4cpp_xform.hpp"
#include "parse_fms_mode.hpp"
#include "gl_framebuffer.hpp"
#include "gl_program.hpp"
#include "gl_debug_callback.h"  // for GLRegisterDebugCallback
#include "cg4cpp_xform.hpp"
#include "gl_texture.hpp"
#include "stats.hpp"
#include "textures.hpp"

#include "sg/path.hpp"

#include "polar/cpu/CPUStrokeTessellator.hpp"
#include "polar/global_polar_settings.hpp"
#include "polar/StrokingBufferState_global.hpp"
#include "polar/StrokerShaderConfig_util.hpp"
#include "polar/polar_shader_state.hpp"

#include "CurveConfig.hpp"

#include "color_space.hpp"

#include "Paint.hpp"

#include "hairline.hpp"
#include "main_window.hpp"
#include "keyinput.hpp"
#include "simple_path.hpp"
#include "hodograph.hpp"
#include "menus.hpp"
#include "decoration.hpp"
#include "cmdline.hpp"
#include "glcolor.hpp"

#ifndef GLAPIENTRYP
# ifdef _WIN32
#  define GLAPIENTRYP __stdcall *
# else
#  define GLAPIENTRYP *
# endif
#endif

typedef StrokeConfigTypes::StrokingSpace StrokingSpace;
typedef StrokeConfigTypes::ShadingMode ShadingMode;
typedef StrokeConfigTypes::CurveType CurveType;

bool curveTypeSelect = true;
char splineType;
vector<float3> inputPoints;

PathContainer simple;

int simple_color_ndx = 0;
float simple_opacity = 0.7;

int cpu_path_count = 0;

static bool draw_annotations = false;  // F9 toggles
static bool draw_normal_extensions = false;  // 'M' toggles

int current_link = 0;

GLSLIncludeProgram stroke_fragment_program[2];

void initSimplePathShader(GLSLIncludeProgram& program, bool aliased)
{
    std::string defines("#version 430 compatibility\n");
    if (aliased) {
        defines += "#define ALIASED 1\n";
    } else {
        defines += "#define ALIASED 0\n";
    }
    program.setCommonDefines(defines);
    program.compileFragmentShaderOnlyFromFiles("shaders/polar/stroke_quad.frag");
}

static bool simple_path_shaders_valid = false;

void initSimplePathShaders()
{
    if (!simple_path_shaders_valid) {
        for (size_t i=0; i<countof(stroke_fragment_program); i++) {
            bool aliased = (i&1);
            initSimplePathShader(stroke_fragment_program[i], aliased);
        }
        simple_path_shaders_valid = true;
    }
}

void invalidateSimplePathShaders()
{
    simple_path_shaders_valid = false;
}

void advanceCurrentLink()
{
    current_link++;
    printf("current_link = %d\n", current_link);
    postGraphingWindowsRedisplay();
}

void reverseCurrentLink()
{
    current_link--;
    printf("current_link = %d\n", current_link);
    postGraphingWindowsRedisplay();
}

float3 getSimpleColor()
{
    switch (simple_color_ndx) {
    case 0:
        return float3(0.8f, 0.1f, 0.3f);
    case 1:
        return float3(1.0f, 1.0f, 1.0f);
    case 2:
        return float3(0.0f, 0.0f, 0.0f);
    case 3:
        return float3(0.7f, 0.7f, 0.3f);
    default:
        assert(!"unknown color");
        return float3(0);
    }
}

void setSimpleColor()
{
    float3 color = getSimpleColor();
    simple.setColor(color);
}

void
initPoints(vector<float3>& inputPoints)
{
    inputPoints.clear();

    int num_coords = config_list[curveType].num_coords;
    splineType = config_list[curveType].spline_type;
    const float* coord = config_list[curveType].coord;
    if (splineType == 'G') {
        // Index into SVG grammar string
        int svg_string_ndx = num_coords;
        const char* svg_string = svg_strings[svg_string_ndx];
        simple.load(1, &svg_string);
    } else if (splineType == 'C') {
        for (int i = 0; i < num_coords; i+=2 ){
            const float3 cp = float3(coord[i+0], coord[i+1], 1);
            inputPoints.push_back(cp);
        }
        simple.loadCubicSpline(inputPoints);
    } else if (splineType == 'O') {
        for (int i = 0; i < num_coords; i+=2 ){
            const float3 cp = float3(coord[i+0], coord[i+1], 1);
            inputPoints.push_back(cp);
        }
        simple.loadCatmullRomLoop(inputPoints);
    } else if (splineType == 'S') {
        for (int i = 0; i < num_coords; i+=2 ){
            const float3 cp = float3(coord[i+0], coord[i+1], 1);
            inputPoints.push_back(cp);
        }
        simple.loadCatmullRomSpline(inputPoints);
    } else if (splineType == 'Q') {
        for (int i = 0; i < num_coords; i+=2 ){
            const float3 cp = float3(coord[i+0], coord[i+1], 1);
            inputPoints.push_back(cp);
        }
        simple.loadQuadraticSpline(inputPoints);
    } else if (splineType == 'L') {
        for (int i = 0; i < num_coords; i+=2 ){
            const float3 cp = float3(coord[i+0], coord[i+1], 1);
            inputPoints.push_back(cp);
        }
        simple.loadLineSpline(inputPoints);
    } else if (splineType == 'W') {
        const float3 cp = float3(coord[0], coord[1], 1);
        inputPoints.push_back(cp);
        for (int i = 2; i < num_coords; i+=5 ){
            float w = coord[i+4];
            w *= w_scale;
            const float3 cp0 = float3(coord[i+0], coord[i+1], w);
            const float3 cp1 = float3(coord[i+2], coord[i+3], 1);
            inputPoints.push_back(cp0);
            inputPoints.push_back(cp1);
        }
        simple.loadConicSpline(inputPoints);
    } else {
        assert(!"unknown splineType");
    }
    setSimpleColor();
    simple.setOpacity(simple_opacity);
    curveTypeSelect = false;
}

void dumpInputPoints()
{
    int scalars = 0;
    for (size_t i=0; i<inputPoints.size(); i++) {
        float3 cp = inputPoints[i];
        float w = cp[2];
        scalars += 2;
        if (w != 1) {
            scalars += 1;
        }
    }
    printf(" { \"Snapshot\", '%c', %d,{ ", splineType, scalars);
    for (size_t i=0; i<inputPoints.size(); i++) {
        float3 cp = inputPoints[i];
        float x = cp.x;
        float y = cp.y;
        float w = cp[2];
        printf("%.9g,%.9g, ", x, y);
        if (w != 1) {
            assert(splineType == 'W');
            i++;
            float3 cp2 = inputPoints[i];
            float x2 = cp2.x;
            float y2 = cp2.y;
            assert(cp2[2] == 1);
            printf("%.9g,%.9g, ", x2, y2);
            printf("%.9g, ", w);
        }
    }
    printf(" }},\n");
}

void glVertex2f(const float2& v)
{
    glVertex2f(v.x, v.y);
}

void glVertex2fw(const float3& v)
{
    float w = v.z;
    glVertex4f(v.x, v.y, 0, w);
}

void glVertex2f_mul(const float2& v)
{
    float3x3 mvp = mul(getMainProjectionMatrix(), modelview);
    float3 p = mul(mvp, float3(v,1));
    glVertex4f(p.x, p.y, 0, p.z);
}

void glVertex2fw_mul(const float3& v)
{
    float3x3 mvp = mul(getMainProjectionMatrix(), modelview);
    float3 p = mul(mvp, v);
    float w = p.z;
    glVertex4f(p.x, p.y, 0, w);
}

#if SUPPORT_POLAR_VISUALIZE
static void drawUniformQuadStrip(const StrokingResults& results)
{
    const vector<float3> &Ppos = results.Ppos;
    const vector<float3> &Pneg = results.Pneg;

    assert(Ppos.size() == Pneg.size());
    glBegin(GL_QUAD_STRIP); {
        for (size_t i = 0; i<results.Ppos.size(); i++){
            glVertex2fw(Ppos[i]);
            glVertex2fw(Pneg[i]);
        }
    } glEnd();
}

static float2 ortho(float2 v)
{
    return float2(v.y, -v.x);
}

static void evalUniformCubicStrokingCPU(StrokingResults& results, const float3 inputPoints[6])
{
    // Generate cubic Bezier
    float2 deltaIn[3];
    deltaIn[0] = inputPoints[1].xy - inputPoints[0].xy;
    deltaIn[1] = inputPoints[2].xy - inputPoints[1].xy;
    deltaIn[2] = inputPoints[3].xy - inputPoints[2].xy;

    vector<float3>& Pmid = results.Pmid;
    vector<float2>& normal = results.normal;
    vector<float3>& Ppos = results.Ppos;
    vector<float3>& Pneg = results.Pneg;

    float u = 0;
    const unsigned int uniform_segments = stroke_tessellator_settings.UniformSegments();
    float delta_u = 1/float(uniform_segments);
    for (unsigned int i=0; i<=uniform_segments; i++){	//parameter control
        // Hard-coded for 4 points
        Pmid.push_back((1-u)*(1-u)*(1-u)*inputPoints[0] + 
            3*u*(1-u)*(1-u)*inputPoints[1] +
            3*u*u*(1-u)*inputPoints[2] +
            u*u*u*inputPoints[3]);
        const float2 gradient_direction = 3*(1-u)*(1-u)*deltaIn[0] + 3*2*(1-u)*u* deltaIn[1] + 3*u*u*deltaIn[2];
        const float2 norm_grad = normalize(gradient_direction);
        const float2 this_normal = ortho(norm_grad);
        normal.push_back(this_normal);
        const float stroke_radius = current_paint.StrokeRadius();
        const float w = Pmid[i].z;
        Ppos.push_back(float3(Pmid[i].xy + w*stroke_radius*this_normal,w));
        Pneg.push_back(float3(Pmid[i].xy - w*stroke_radius*this_normal,w));
        u += delta_u;
    }
}

static void evalUniformQuadraticStrokingCPU(StrokingResults& results, const float3 inputPoints[6])
{
    float3 cp[3];
    cp[0] = inputPoints[0];
    cp[1] = inputPoints[1];
    cp[2] = inputPoints[3];

    vector<float3>& Pmid = results.Pmid;
    vector<float2>& normal = results.normal;
    vector<float3>& Ppos = results.Ppos;
    vector<float3>& Pneg = results.Pneg;

    float u = 0;
    const unsigned int uniform_segments = stroke_tessellator_settings.UniformSegments();
    float delta_u = 1/float(uniform_segments);
    for (unsigned int i=0; i<=uniform_segments; i++){	//parameter control
        // Hard-coded for 4 points
        const float3 p = evalRationalBezierQuadratic(cp, u);
        Pmid.push_back(p);
        const float3 dp3 = evalRationalBezierQuadraticGradient(cp, u);
        const float2 dp = dp3.xy/dp3.z;
        const float2 gradient_direction = dp;
        float2 norm_grad = normalize(gradient_direction);
        float2 this_normal = ortho(norm_grad);
        normal.push_back(this_normal);
        const float stroke_radius = current_paint.StrokeRadius();
        const float w = Pmid[i].z;
        Ppos.push_back(float3(Pmid[i].xy + w*stroke_radius*this_normal,w));
        Pneg.push_back(float3(Pmid[i].xy - w*stroke_radius*this_normal,w));
        u += delta_u;
    }
}

// Cull if returns true
static bool evalUniformStrokingCPU(StrokingResults& results, const float3 inputPoints[6], const GLuint ndx[6])
{
    results.reset();

    bool is_quadratic = (ndx[2]==0);

    if (is_quadratic) {
        evalUniformQuadraticStrokingCPU(results, inputPoints);
    } else {
        evalUniformCubicStrokingCPU(results, inputPoints);
    }
    return false;
}

static void drawUniformStrokingGeneratingLines(const StrokingResults& results)
{
    sRGB_glColor3fv(white);
    fb->disableCoverageModulation();
    const vector<float3>& Pmid = results.Pmid;
    glBegin(GL_LINE_STRIP); {
        for (size_t i=0; i<Pmid.size();i++){
            glVertex2fw(Pmid[i]);
        }
    } glEnd();
    fb->enableCoverageModulation();
}

static void drawUniformStrokingGeneratingNormals(const StrokingResults& results)
{
    sRGB_glColor3fv(whitishBlue);
    fb->disableCoverageModulation();
    const vector<float3>& Ppos = results.Ppos;
    const vector<float3>& Pneg = results.Pneg;
    assert(Ppos.size() == Pneg.size());
    glBegin(GL_LINES); {
        for (size_t i = 0; i<Ppos.size(); i++){
            glVertex2fw(Ppos[i]);
            glVertex2fw(Pneg[i]);
        }
    } glEnd();
    fb->enableCoverageModulation();
}

static void drawUniformStrokingGeneratingPoints(const StrokingResults& results)
{
    sRGB_glColor3fv(red);
    fb->pointSize(5);
    const vector<float3>& Pmid = results.Pmid;
    glBegin(GL_POINTS); {
        for (size_t i = 0; i < Pmid.size(); i++){
            glVertex2fw(Pmid[i]);
        }
    } glEnd();
}
#endif

void sRGB_glColor4f(float3 color, float opacity)
{
    sRGB_glColor4f(color.r*opacity,color.g*opacity,color.b*opacity,opacity);
}

static void sendShadingModeOpacity(ShadingMode mode, float opacity)
{
    if (skip_paint_update) {
        return;
    }
    switch (mode) {
    case StrokeConfigTypes::CONSTANT_COLOR:
        {
            float3 color = getSimpleColor();
            sRGB_glColor4f(color, opacity);
        }
        break;
    case StrokeConfigTypes::LINEAR_GRADIENT:
    case StrokeConfigTypes::CENTER_RADIAL_GRADIENT:
    case StrokeConfigTypes::FOCAL_RADIAL_GRADIENT:
    case StrokeConfigTypes::IMAGE_GRADIENT:
    case StrokeConfigTypes::ARC_LENGTH_TEXTURE_2D:
    case StrokeConfigTypes::ARC_LENGTH_TEXTURE_2D_EARLY_FRAGMENT_TEST:
        sRGB_glColor4f(opacity,opacity,opacity,opacity);
        break;
    default:
        assert(!"unknown ShadingMode");
    }
}

#if SUPPORT_POLAR_VISUALIZE
static void coverUniformStrokingCPU(const StrokingResults& results)
{
    // Note: uniform stroking lacks arc length texture support.
    sendShadingModeOpacity(StrokeConfigTypes::CONSTANT_COLOR, simple_opacity);
    drawUniformQuadStrip(results);
}

static void annotateUniformStrokingCPU(const StrokingResults& results)
{
    // Draw the curve:
    if (decoration & (DRAW_GENERATING_LINES|DRAW_GENERATING_NORMALS)) {
        fb->lineWidth(1.0f);
    }
    if (decoration & DRAW_GENERATING_LINES) {
        drawUniformStrokingGeneratingLines(results);
    }

    if (decoration & DRAW_GENERATING_NORMALS) {
        drawUniformStrokingGeneratingNormals(results);
    }

    if (decoration & DRAW_GENERATING_POINTS) {
        drawUniformStrokingGeneratingPoints(results);
    }
}

static void stencilUniformStrokingCPU(const StrokingResults& results)
{
    float opacity = 0.6;
    sRGB_glColor4f(0.7*opacity,0.7*opacity,0.3*opacity,opacity);  // yellow
    drawUniformQuadStrip(results);
}
#endif

static void GoodOKBadColor(float angle)
{
    const float max_turning_angle_delta = stroke_tessellator_settings.MaxTurningAngleDelta();
    if (angle > 1.1f * max_turning_angle_delta) {
        sRGB_glColor3fv(red);  // BAD, under-tessellated
    } else if (angle > max_turning_angle_delta) {
        sRGB_glColor3fv(dark_yellow);  // SLIGHTLY BAD, slightly under-tessellated
    } else if (angle < 0.9f * max_turning_angle_delta) {
        sRGB_glColor3fv(yellow);  // OK, but over-tessellated
    } else {
        sRGB_glColor3fv(green);  // GOOD, not too little or too much tessellation
    }
}

static void drawPolarQuadStripAnnotationsOffsetPosition(const float3x3& xform, const vector<float3> &Ppos, size_t i)
{
    float3 a = Ppos[i-1];
    float3 b = Ppos[i];
    float3 c = Ppos[i+1];

    float3 a2 = mul(xform, a);
    float3 b2 = mul(xform, b);
    float3 c2 = mul(xform, c);

    float2 aw = project(a2);
    float2 bw = project(b2);
    float2 cw = project(c2);

    float2 dif1 = bw-aw;
    float2 dif2 = cw-bw;

    float2 t1 = normalize(dif1);
    float2 t2 = normalize(dif2);

    float cos_angle = dot(t1,t2);
    if (cos_angle <= 0) {
        // Cusp point!
        sRGB_glColor3fv(purple);
    } else {
        float abs_cos_angle = fabs(cos_angle);
        float angle = acos(abs_cos_angle);

        GoodOKBadColor(angle);
    }
    glVertex2fw(Ppos[i]);
}

static void drawPolarQuadStripAnnotationsMinMaxPoints(const StrokingResults& results, bool black_and_white)
{
    if (results.Ppos_max_ndx != ~0u) {
        if (black_and_white) {
            sRGB_glColor3fv(gray);
        } else {
            GoodOKBadColor(results.Ppos_max_angle);
        }
        glVertex2fw(results.Ppos[results.Ppos_max_ndx]);
    }
    if (results.Ppos_min_ndx != ~0u) {
        if (black_and_white) {
            sRGB_glColor3fv(black);
        } else {
            GoodOKBadColor(results.Ppos_min_angle);
        }
        glVertex2fw(results.Ppos[results.Ppos_min_ndx]);
    }
    if (results.Pneg_max_ndx != ~0u) {
        if (black_and_white) {
            sRGB_glColor3fv(gray);
        } else {
            GoodOKBadColor(results.Pneg_max_angle);
        }
        glVertex2fw(results.Pneg[results.Pneg_max_ndx]);
    }
    if (results.Pneg_min_ndx != ~0u) {
        if (black_and_white) {
            sRGB_glColor3fv(black);
        } else {
            GoodOKBadColor(results.Pneg_min_angle);
        }
        glVertex2fw(results.Pneg[results.Pneg_min_ndx]);
    }
}

static void drawPolarQuadStripAnnotations(const StrokingResults& results)
{
    const vector<float3> &Ppos = results.Ppos;
    const vector<float3> &Pneg = results.Pneg;

    fb->pointSize(5.0);
    glBegin(GL_POINTS); {
        const float3x3 xform = scale(getMainWindowWidth(), getMainWindowHeight());

        for (size_t i = 1; i<Ppos.size()-1; i++){
            drawPolarQuadStripAnnotationsOffsetPosition(xform, Ppos, i);
            drawPolarQuadStripAnnotationsOffsetPosition(xform, Pneg, i);
        }
    } glEnd();

    fb->pointSize(9.0);
    glBegin(GL_POINTS); {
        drawPolarQuadStripAnnotationsMinMaxPoints(results, false);
    } glEnd();

    fb->pointSize(3.0);
    glBegin(GL_POINTS); {
        drawPolarQuadStripAnnotationsMinMaxPoints(results, true);
    } glEnd();
}

static void drawNormalExtensions(const StrokingResults& results)
{
    glBegin(GL_LINES); {
        const vector<float3> &Pmid = results.Pmid;
        const vector<float3> &Pradius = results.Pradius;
        const vector<float3> &Pinfinite = results.Pinfinite;

        for (size_t i = 0; i<Pmid.size(); i++){
            if (i&1) {
                sRGB_glColor3fv(gray);
            } else {
                sRGB_glColor3fv(white);
            }
            glVertex2fw(Pmid[i]);
            glVertex2fw(Pinfinite[i]);
        }
        sRGB_glColor3fv(purple);
        for (size_t i = 0; i<Pmid.size(); i++){
            glVertex2fw(Pmid[i]);
            glVertex2fw(Pradius[i]);
        }
    } glEnd();
}

static void configShadingMode(StrokeConfigTypes::ShadingMode mode)
{
    switch (mode) {
    case StrokeConfigTypes::CONSTANT_COLOR:
        break;
    case StrokeConfigTypes::LINEAR_GRADIENT:
    case StrokeConfigTypes::CENTER_RADIAL_GRADIENT:
    case StrokeConfigTypes::FOCAL_RADIAL_GRADIENT:
        glEnable(GL_TEXTURE_1D);
        break;
    case StrokeConfigTypes::IMAGE_GRADIENT:
    case StrokeConfigTypes::ARC_LENGTH_TEXTURE_2D:
        glEnable(GL_TEXTURE_2D);
        break;
    case StrokeConfigTypes::ARC_LENGTH_TEXTURE_2D_EARLY_FRAGMENT_TEST:
        glEnable(GL_TEXTURE_2D);
        break;
    case StrokeConfigTypes::NUM_SHADING_MODES:
    default:
        assert(!"unknown ShadingMode");
    }
}

static void deconfigShadingMode(ShadingMode mode)
{
    switch (mode) {
    case StrokeConfigTypes::CONSTANT_COLOR:
        break;
    case StrokeConfigTypes::LINEAR_GRADIENT:
    case StrokeConfigTypes::CENTER_RADIAL_GRADIENT:
    case StrokeConfigTypes::FOCAL_RADIAL_GRADIENT:
        glDisable(GL_TEXTURE_1D);
        break;
    case StrokeConfigTypes::IMAGE_GRADIENT:
    case StrokeConfigTypes::ARC_LENGTH_TEXTURE_2D:
        glDisable(GL_TEXTURE_2D);
        break;
    case StrokeConfigTypes::ARC_LENGTH_TEXTURE_2D_EARLY_FRAGMENT_TEST:
        glDisable(GL_TEXTURE_2D);
        break;
    case StrokeConfigTypes::NUM_SHADING_MODES:
    default:
        assert(!"unknown ShadingMode");
    }
}

// Debugging mode, color each strip from red (start) to green (end).
int color_by_segment_length = 0;
// Debugging mode, to draw only a range of quads in the strip
int splice_strip = 0;
size_t splice_start;
size_t splice_stop;

int do_decoration = 0;

static float determinant(float2 a, float2 b)
{
    return a.x*b.y - a.y*b.x;
}

static float2 homodiff(const float3& a, const float3& b)
{
    return (a.xy*b.z - b.xy*a.z)/(a.z*b.z);
}

static float3 homoadd(const float3& a, const float2& b)
{
    return float3(a.xy + b.xy*a.z, a.z);
}

static float homodistance(const float3& a, const float3& b)
{
    float2 diff = homodiff(a, b);
    return length(diff);
}

static float sqdist(const float2& a, const float2& b)
{
    const float2 d = a-b;
    return dot(d,d);
}

static float biggerDistance(const float3& Pa, const float2& Na,
                            const float3& Pb, const float2& Nb,
                            const float stroke_width_fraction)
{
#if 0
    assert(dot(Na, Nb) > 0);
#endif

    const float2 Pa_neg = (Pa.xy - stroke_width_fraction*Na*Pa.z)/Pa.z;
    const float2 Pb_neg = (Pb.xy - stroke_width_fraction*Nb*Pb.z)/Pb.z;

    const float2 Pa_pos = (Pa.xy + stroke_width_fraction*Na*Pa.z)/Pa.z;
    const float2 Pb_pos = (Pb.xy + stroke_width_fraction*Nb*Pb.z)/Pb.z;

    const float sq_dist1 = sqdist(Pa_neg, Pb_neg);
    const float sq_dist2 = sqdist(Pa_pos, Pb_pos);

    if (sq_dist1 > sq_dist2) {
        return sqrt(sq_dist1);
    } else {
        return sqrt(sq_dist2); 
    }
}

static float biggerDistance(const float3& Pa_neg, const float3& Pb_neg,
                            const float3& Pa_pos, const float3& Pb_pos,
                            const float stroke_width_fraction)
{
    const float sq_dist1 = sqdist(project(Pa_neg), project(Pb_neg));
    const float sq_dist2 = sqdist(project(Pa_pos), project(Pb_pos));

    if (sq_dist1 > sq_dist2) {
        return sqrt(sq_dist1);
    } else {
        return sqrt(sq_dist2); 
    }
}


static float segmentArcLength(const StrokingResults& results, size_t i)
{
    const vector<float3>& P = results.p;
    const vector<float2>& normal = results.normal;

    const float2 prior_normal = normal[i-1];
    const float2 this_normal = normal[i];
    const float3 prior_P = P[i-1];
    const float3 this_P = P[i];
    if (prior_P.z < 0 || this_P.z < 0) {
        return 0;
    }
    const float stroke_width_fraction = current_paint.ArcLengthExtension() * current_paint.StrokeRadius();

    float wrong = biggerDistance(prior_P, prior_normal, this_P, this_normal, stroke_width_fraction);
    return wrong;

    const vector<float3>& Ppos_obj = results.Ppos_obj;
    const vector<float3>& Pneg_obj = results.Pneg_obj;

    return biggerDistance(Pneg_obj[i-1], Pneg_obj[i],
                          Ppos_obj[i-1], Ppos_obj[i], stroke_width_fraction);
}

const float too_flat_threshold = 0.001;  // asin(0.001) = 0.0573 degrees

static void drawDecorations(const StrokingResults& results)
{
    glDisable(GL_STENCIL_TEST);
    glPointSize(10);
    const vector<float3>& Ppos_obj = results.Ppos_obj;
    const vector<float3>& Pneg_obj = results.Pneg_obj;
    const vector<float3>& Pmid_obj = results.p;

    const vector<float2>& normal = results.normal;

    for (size_t i = 1; i<Ppos_obj.size(); i++){
        const size_t this_i = i;
        const size_t prev_i = i-1;

        const float3 Pa = Pmid_obj[prev_i];
        const float3 Pb = Pmid_obj[this_i];

        const float2 Na = normal[prev_i];
        const float2 Nb = normal[this_i];

        float Pa_pos_scale = 1;
        float Pb_pos_scale = 1;
        float Pa_neg_scale = 1;
        float Pb_neg_scale = 1;

        const float lambda_denom = determinant(Nb, Na);
        // Are normals sufficiently distinct from each other?
        if (fabs(lambda_denom) >= too_flat_threshold) {
            const float lambda_numer = determinant(Nb, homodiff(Pb,Pa));
            const float lambda  = lambda_numer/lambda_denom;
            // Pi = intersection point for rays from Pa and Pb going in directions Na and Nb respectively.
            const float3 Pi = homoadd(Pa, lambda*Na);

            if (i&1) {
                glColor3f(1,0,0);
            } else {
                glColor3f(0,1,0);
            }
            glBegin(GL_LINES); {
                glVertex2fw_mul(Pneg_obj[prev_i]);
                glVertex2fw_mul(Pi);

                glVertex2fw_mul(Ppos_obj[this_i]);
                glVertex2fw_mul(Pi);
            } glEnd();

            glBegin(GL_POINTS); {
                glColor3f(1,1,1);
                glVertex2fw_mul(Pi);
            } glEnd();
        }
    }
}

static int do_hack = 0;

static void drawPolarQuadStrip(const StrokingResults& results, PathDrawState* state)
{
    const float texture_scale = g_stroker_context.getArcLengthTexGenScale();
    const float texture_offset = g_stroker_context.getArcLengthOffset() * texture_scale;

    const float scaled_arc_length_offset = texture_offset;
    const vector<float3>& Ppos = results.Ppos;
    const vector<float3>& Pneg = results.Pneg;
    const vector<float3>& Pmid = results.Pmid;
    const vector<float3>& P = results.p;
    const vector<float3>& Ppos_obj = results.Ppos_obj;
    const vector<float3>& Pneg_obj = results.Pneg_obj;
    const vector<float2>& frieze_height = results.frieze_height;

    const vector<float2>& normal = results.polar_offset;

    assert(Ppos.size() == Pneg.size());
    assert(Ppos_obj.size() == Ppos.size());
    assert(Pneg_obj.size() == Pneg.size());
    const float initial_accumulated_path_length = state ? state->accumulated_path_length : 0;
    float accumulated_path_length = initial_accumulated_path_length;
    const float scale = getCurrentTextureScale();
    const float stroke_width = current_paint.StrokeWidth();
    if (state) {
        bindCurrentTexture();
        ShadingMode mode = getShadingMode();
        configShadingMode(mode);
    }
    float strip_accumulated_path_length = 0;
    float total_strip_accumulated_path_length = 0;
    if (color_by_segment_length) {
        for (size_t i = 1; i<Ppos.size(); i++){
            total_strip_accumulated_path_length += segmentArcLength(results, i);
        }
    }

    strip_accumulated_path_length = 0;
    bool do_quads = g_stroker_context.usesArcLengthExtension();
    if (do_quads) {
        initSimplePathShaders();

        const vector<float3>& Ppos_obj = results.Ppos_obj;
        const vector<float3>& Pneg_obj = results.Pneg_obj;
        const vector<float3>& Pmid_obj = results.p;

        const vector<float2>& normal = results.normal;

        if (fb->CoverageSamples() == 1) {
            stroke_fragment_program[1].use();
        } else {
            stroke_fragment_program[0].use();
        }

        glBegin(GL_QUADS); {
            for (size_t i = 1; i<Ppos.size(); i++){
                const size_t this_i = i;
                const size_t prev_i = i-1;

                const float3 Pa = Pmid_obj[prev_i];
                const float3 Pb = Pmid_obj[this_i];

                const float2 Na = normal[prev_i];
                const float2 Nb = normal[this_i];

                float Pa_pos_scale = 1;
                float Pb_pos_scale = 1;
                float Pa_neg_scale = 1;
                float Pb_neg_scale = 1;

                const float lambda_denom = determinant(Nb, Na);
                // Are normals sufficiently distinct from each other?
                if (fabs(lambda_denom) >= too_flat_threshold) {
                    const float lambda_numer = determinant(Nb, homodiff(Pb,Pa));
                    const float lambda  = lambda_numer/lambda_denom;
                    // Pi = intersection point for rays from Pa and Pb going in directions Na and Nb respectively.
                    const float3 Pi = homoadd(Pa, lambda*Na);

                    const float2 v = homodiff(Pb, Pa);
                    const float2 w = homodiff(Pi, Pa);
                    const float winding = determinant(w,v);
                    const float sgn = winding >= 0 ? +1 : -1;
                    const float c1 = dot(w,v);
                    const float c2 = dot(v,v);
                    const float b = (c2 != 0) ? c1/c2 : 0;

                    const float stroke_radius = current_paint.StrokeRadius();

                    // Pc is point closed to Pi on the line formed by Pa and Pb.
                    const float3 Pc = homoadd(Pa, b*v);

                    const float CI_dist = homodistance(Pi, Pc);
                    const float CA_dist = homodistance(Pa, Pc);
                    const float CB_dist = homodistance(Pb, Pc);

                    const float AI_dist = homodistance(Pi, Pa);
                    const float BI_dist = homodistance(Pi, Pb);

                    // Division-by-zero is OK here because biggest will be NaN and fall into "else" case below.
                    const float A_extension = stroke_radius*CI_dist/AI_dist;
                    const float B_extension = stroke_radius*CI_dist/BI_dist;

                    Pa_pos_scale = CI_dist - sgn*A_extension;
                    Pb_pos_scale = CI_dist - sgn*B_extension;
                    Pa_neg_scale = CI_dist + sgn*A_extension;
                    Pb_neg_scale = CI_dist + sgn*B_extension;

                    // Find biggest magnitude of all four.
                    const float biggest = max(max(fabs(Pa_pos_scale),fabs(Pa_neg_scale)),
                                              max(fabs(Pb_pos_scale),fabs(Pb_neg_scale)));
                    if (fabs(biggest) > 0) {
                        // Everything OK.
                        const float one_over_biggest = 1/biggest;
                        Pa_pos_scale *= one_over_biggest;
                        Pa_neg_scale *= one_over_biggest;
                        Pb_pos_scale *= one_over_biggest;
                        Pb_neg_scale *= one_over_biggest;
                    } else {
                        // Could get here because biggest is either zero or NaN!
                        Pa_pos_scale = 1;
                        Pa_neg_scale = -1;
                        Pb_pos_scale = 1;
                        Pb_neg_scale = -1;
                    }
                }

                const float last_accumulated_path_length = accumulated_path_length;
                accumulated_path_length += segmentArcLength(results, i);

                float negA = frieze_height[prev_i][0];
                float posA = frieze_height[prev_i][1];
                float negB = frieze_height[this_i][0];
                float posB = frieze_height[this_i][1];

                if (negA == 0.5f || posA == 0.5f) {
                    Pa_pos_scale = 1;
                    Pb_pos_scale = 1;
                    Pa_neg_scale = 1;
                    Pb_neg_scale = 1;
                }

                const float s0 = last_accumulated_path_length*scale/stroke_width + scaled_arc_length_offset;
                glTexCoord4f(s0*Pa_pos_scale, Pa_pos_scale, posA, -1*Pa_pos_scale);
                glVertex2fw(Ppos[i-1]);
                glTexCoord4f(s0*Pa_neg_scale, Pa_neg_scale, negA, -1*Pa_neg_scale);
                glVertex2fw(Pneg[i-1]);

                const float s1 = accumulated_path_length*scale/stroke_width + scaled_arc_length_offset;
                glTexCoord4f(s1*Pb_neg_scale, Pb_neg_scale, negB, +1*Pb_neg_scale);
                glVertex2fw(Pneg[i]);
                glTexCoord4f(s1*Pb_pos_scale, Pb_pos_scale, posB, +1*Pb_pos_scale);
                glVertex2fw(Ppos[i]);

            }
        } glEnd();

        glUseProgram(0);
    } else {
        if (do_hack) {
            printf(">> start\n");
        }
        glBegin(GL_QUAD_STRIP); {
            for (size_t i = 0; i<Ppos.size(); i++){
                if (i > 0) {
                    const float segment_len = segmentArcLength(results, i);
                    accumulated_path_length += segment_len;
                    if (color_by_segment_length) {
                        strip_accumulated_path_length += segment_len;
                    }
                    if (do_hack) {
                        printf("%d/%d: len = %g (%g), [%g,%g]\n", int(i), int(Ppos.size()), accumulated_path_length, segment_len, frieze_height[i][0], frieze_height[i][1]);
                    }
                } else {
                    if (do_hack) {
                        printf("%d/%d: len = %g (%g), [%g,%g]\n", 0, int(Ppos.size()), accumulated_path_length, 0.0f, frieze_height[i][0], frieze_height[i][1]);
                    }
                }

                float neg = frieze_height[i][0];
                float pos = frieze_height[i][1];

                if (splice_strip) {
                    if (i < splice_start || i > splice_stop) {
                        continue;
                    }
                }

                if (color_by_segment_length) {
                    const float t = strip_accumulated_path_length/total_strip_accumulated_path_length;
                    const float3 red = float3(1,0,0);
                    const float3 green = float3(0,1,0);
#define lerp(a,b,t) ((a) + (t)*((b)-(a)))
                    const float3 color = lerp(red, green, t);
                    glColor3fv(&color[0]);
                }
                const float s = accumulated_path_length*scale/stroke_width + scaled_arc_length_offset;
                glTexCoord2f(s, pos);
                glVertex2fw(Ppos[i]);
                glTexCoord2f(s, neg);
                glVertex2fw(Pneg[i]);
            }
        } glEnd();
        if (do_hack) {
            printf("<< end\n");
        }
    }

    ShadingMode mode = getShadingMode();
    deconfigShadingMode(mode);
    if (state) {
        state->accumulated_path_length = accumulated_path_length;
    }

    if (do_decoration) {
        drawDecorations(results);
    }
}

static void stencilPolarQuadStrip(const StrokingResults& results)
{
    const vector<float3> &Ppos = results.Ppos;
    const vector<float3> &Pneg = results.Pneg;
    const vector<float3> &P = results.p;

    assert(Ppos.size() == Pneg.size());
    glBegin(GL_QUAD_STRIP); {
        for (size_t i = 0; i<Ppos.size(); i++){
            glVertex2fw(Ppos[i]);
            glVertex2fw(Pneg[i]);
        }
    } glEnd();
}

static void loadStrokerWithPaintState(CPUStrokeTessellator& stroker, const Paint& paint)
{
    stroker.setStrokeWidth(paint.StrokeWidth());
    stroker.setStrokeRadiusFraction(paint.ArcLengthExtension() * paint.StrokeWidth() * 0.5f);
    stroker.setMaxJoinSteps(StrokeTessellator::JoinEnumToMaxJoinSteps(paint.StrokeTessellatorJoinStyle()));
    stroker.setEffectiveMiterLimit(paint.EffectiveMiterLimit());
    stroker.setCapMask(paint.StrokeTessellatorCapMask());
}

bool prepStroker(CPUStrokeTessellator &stroker, const Paint& paint, const float3 link_control_points[], const GLuint ndx[6], unsigned int primitive_id)
{
    stroker.setPrimitiveID(primitive_id);
    loadStrokerWithPaintState(stroker, paint);
    if (g_stroker_context.isWindowSpace()) {
        stroker.setStrokeSpace(CPUStrokeTessellator::WINDOW_SPACE);
    } else {
        stroker.setStrokeSpace(CPUStrokeTessellator::EYE_SPACE);
    }
    // Establish transforms.
    {
        float3x3 mvp = mul(getMainProjectionMatrix(), modelview);
        float2 vs = float2(getMainWindowWidth(), -getMainWindowHeight());
#if 0
        float3x3 identity = float3x3(1,0,0, 0,1,0, 0,0,1);
        // XXX grr, weird that projection matrix is really mvp.
        stroker.setTransforms(identity, mvp, mvp, vs);
#else
        stroker.setTransforms(modelview, getMainProjectionMatrix(), mvp, vs);
#endif
        if (verbose) {
            cout << "modelview = " << modelview << endl;
            cout << "projection = " << getMainProjectionMatrix() << endl;
            cout << "mvp = " << mvp << endl;
            cout << "vs = " << vs << endl;
        }
        if (stroker.invalidTransforms()) {
            return true;  // Cull everything when transforms contain NaN values.
        }
    }
    // The Setup Math: appropriate for a tessellation control shader
    stroker.setupLink(link_control_points, ndx);
    bool culled = stroker.isCulled();
    return culled;
}

bool needsPerspectivePolarStroking(const CPUStrokeTessellator &stroker)
{
    return stroker.LinkTypeIsCurved();
}

static float ratio_exponent = 1.25;

static const float Klinear = (1+0.05);
static const float Kquadratic = 0.025;

static const float KKlinear = (1+0.005);
static const float KKquadratic = 0.0025;

// Cull if returns true
static bool evalPolarStrokingCPU(CPUStrokeTessellator &stroker,
                                 StrokingResults& results,
                                 const float3 link_control_points[],
                                 const GLuint ndx[6],
                                 unsigned int primitive_id)
{
    results.reset();

    bool cull = prepStroker(stroker, current_paint, link_control_points, ndx, primitive_id);
    if (cull) {
        // Skip evaluation.
    } else {
        // The Evaluation Math: appropriate for a tessellation evaluation shader
        if (verbose) {
            stroker.print(stdout);
        }
        stroker.evalRibs(results);
    }
    return cull;
}

static void annotatePolarStrokingCPU(const CPUStrokeTessellator &stroker, const StrokingResults& results)
{
    const bool is_curved = stroker.LinkTypeIsCurved();
    if (draw_annotations && is_curved) {
        drawPolarQuadStripAnnotations(results);
    }
    if (draw_normal_extensions) {  // 'M' toggles
        drawNormalExtensions(results);
    }

    if (decoration & (DRAW_GENERATING_LINES|DRAW_GENERATING_NORMALS)) {
        fb->lineWidth(1.0f);
    }
    if (decoration & DRAW_GENERATING_LINES) {
        sRGB_glColor3fv(white);
        fb->disableCoverageModulation();
        const std::vector<float3>& Pmid = results.Pmid;
        glBegin(GL_LINE_STRIP); {
            for (size_t i = 0; i < Pmid.size(); i++){
                glVertex2fw(Pmid[i]);
            }
        } glEnd();
        fb->enableCoverageModulation();
    }

    if (decoration & DRAW_GENERATING_NORMALS) {
        sRGB_glColor3fv(whitishBlue);
        fb->disableCoverageModulation();
        const std::vector<float3>& Ppos = results.Ppos;
        const std::vector<float3>& Pneg = results.Pneg;
        assert(Ppos.size() == Pneg.size());
        glBegin(GL_LINES); {
            for (size_t i = 0; i<Ppos.size(); i++){
                glVertex2fw(Ppos[i]);
                glVertex2fw(Pneg[i]);
            }
        } glEnd();
        fb->enableCoverageModulation();
    }
    if (decoration & DRAW_GENERATING_POINTS) {
        sRGB_glColor3fv(red);
        fb->pointSize(5);
        const std::vector<float3>& Pmid = results.Pmid;
        glBegin(GL_POINTS); {
            for (size_t i = 0; i < Pmid.size(); i++){
                glVertex2fw(Pmid[i]);
            }
        } glEnd();
    }

    // Is a cubic?
    if (stroker.LinkType() == -3) {  // XXX
        if (show_hodograph_inflection_points) {
            const float3x3 mvp = stroker.ModelviewProjectionMatrix();
            for (int i=0; i<4; i++) {
                float t = stroker.Tbound(i);
                if (t > 0.0f && t < 1.0f) {
                    const float3 p = (1-t)*(1-t)*(1-t)*stroker.Cp(0) +
                        3*(1-t)*(1-t)*t*stroker.Cp(1) +
                        3*(1-t)*t*t*stroker.Cp(2) +
                        t*t*t*stroker.Cp(3);
                    const float3 clip_p = mul(mvp, p);
                    // Draw black point with red outline.
                    fb->pointSize(5);
                    sRGB_glColor3fv(blue);
                    glBegin(GL_POINTS); {
                        glVertex2fw(clip_p);
                    } glEnd();        
                    fb->pointSize(3);
                    sRGB_glColor3fv(black);
                    glBegin(GL_POINTS); {
                        glVertex2fw(clip_p);
                    } glEnd();  
                }
            }
        }
        if (draw_ref_t) {
            // Draw the reference point.
            const float3x3 mvp = stroker.ModelviewProjectionMatrix();

            const float t = ref_t;
            const float3 p = (1-t)*(1-t)*(1-t)*stroker.Cp(0) +
                3*(1-t)*(1-t)*t*stroker.Cp(1) +
                3*(1-t)*t*t*stroker.Cp(2) +
                t*t*t*stroker.Cp(3);
            const float3 clip_p = mul(mvp, p);
            // Draw black point with red outline.
            fb->pointSize(5);
            sRGB_glColor3fv(red);
            glBegin(GL_POINTS); {
                glVertex2fw(clip_p);
            } glEnd();        
            fb->pointSize(3);
            sRGB_glColor3fv(black);
            glBegin(GL_POINTS); {
                glVertex2fw(clip_p);
            } glEnd();     
        }
    }

    // Is a quadratic?
    if (stroker.LinkType() == -2) {  // XXX
        // Quadratics don't have inflection points
        // XXX conics sort of do
        if (draw_ref_t) {
            // Draw the reference point.
            const float3x3 mvp = stroker.ModelviewProjectionMatrix();

            const float t = ref_t;
            const float3 p = (1-t)*(1-t)*stroker.Cp(0) +
                2*(1-t)*t*stroker.Cp(1) +
                t*t*stroker.Cp(3);
            const float3 clip_p = mul(mvp, p);
            // Draw black point with red outline.
            fb->pointSize(5);
            sRGB_glColor3fv(red);
            glBegin(GL_POINTS); {
                glVertex2fw(clip_p);
            } glEnd();        
            fb->pointSize(3);
            sRGB_glColor3fv(black);
            glBegin(GL_POINTS); {
                glVertex2fw(clip_p);
            } glEnd();     
        }
    }
}

static void stencilPolarStrokingCPU(const CPUStrokeTessellator &stroker, const StrokingResults& results, /*unused*/PathDrawState& state)
{
    stencilPolarQuadStrip(results);

    // Stencil testing should still be enabled from "stencil" step.
    assert(glIsEnabled(GL_STENCIL_TEST));
}

static void coverPolarStrokingCPU(const CPUStrokeTessellator &stroker, const StrokingResults& results, PathDrawState& state)
{
    sendShadingModeOpacity(getShadingMode(), simple_opacity);
    drawPolarQuadStrip(results, &state);
}

static void drawPatchSet(GLSLProgram* program)
{
    // Toggle trivial culling with 'C' to avoid these being culled.
    for (int i=0; i<simple.numPatches(); i++) {
        const int x = 100, y = 500;
        const int width = 250, height = 250;
        const int space = 260;
        const int delta = space*i;
        program->set4f("corners",
            x + delta, y,
            x + delta + width, y+height);
        simple.drawOnePatch(program, i);
    }
}

static void drawSimplePathWithNVpr()
{
    // Use NV_path_rendering to draw paths.
    glEnable(GL_STENCIL_TEST);
    const StrokeConfigTypes::CompositeMode compositing = g_stroker_context.getCompositeMode();
    if (compositing >= StrokeConfigTypes::COMPOSITE_BLEND) {
        // Only blended compositing mode for NV_path_rendering is "blend once and only once"
        glBlendFunc(GL_ONE, GL_ONE_MINUS_SRC_ALPHA);  // "over" blending
        glEnable(GL_BLEND);
        float opacity = 0.6;
        sRGB_glColor4f(0.7*opacity,0.7*opacity,0.3*opacity,opacity);  // yellow
    } else {
        sRGB_glColor3fv(red);
    }
    glStencilFunc(GL_NOTEQUAL, 0, 0x1F);
    glStencilOp(GL_KEEP, GL_KEEP, GL_ZERO);

    simple.drawStrokeNVPR();

    if (compositing >= StrokeConfigTypes::COMPOSITE_BLEND) {
        glDisable(GL_BLEND);
    }
    glDisable(GL_STENCIL_TEST);
}

struct LineStripInfo {
    size_t min_ndx;
    size_t max_ndx;

    float min_abs_cos;
    float max_abs_cos;

    float min_angle;
    float max_angle;

    float avg_angle;

    vector<float> angle;
};

static void processDirectionLineStrip(const std::vector<float2> &normal, bool do_normalize, LineStripInfo &result)
{
    double min_abs_cos = 1;
    double max_abs_cos = 0;
    double angle_sum = 0;
    int min_ndx = -1;
    int max_ndx = -1;

    for (size_t i=1; i<normal.size(); i++) {
        const float2 g1 = normal[i-1];
        const float2 g2 = normal[i];

        const float2 a = do_normalize ? normalize(g1) : g1;
        const float2 b = do_normalize ? normalize(g2) : g2;

        const float v = dot(a,b);
        if (v <= 0) {
            // Encountered cusp.
            
            // Ignore this gradient angle because "inside" the elbow of the stroke.
        } else {
            const float abs_v = fabs(v);

            if (abs_v < min_abs_cos) {
                min_ndx = int(i);
                min_abs_cos = abs_v;
            }
            if (abs_v > max_abs_cos) {
                max_ndx = int(i);
                max_abs_cos = abs_v;
            }

            const float angle = degrees(acos(abs_v));
            angle_sum += angle;
            result.angle.push_back(angle);
        }
    }

    result.min_ndx = min_ndx;
    result.min_abs_cos = min_abs_cos;
    result.min_angle = acos(min_abs_cos);

    result.max_ndx = max_ndx;
    result.max_abs_cos = max_abs_cos;
    result.max_angle = acos(max_abs_cos);

    result.avg_angle = float(angle_sum / (normal.size()-1));
}

void processPositionLineStrip(const std::vector<float3> &position, const float3x3 xform, LineStripInfo &result)
{
    bool loud = false;  // XXX debugging, remove me at some point
    double min_abs_cos = 1;
    double max_abs_cos = 0;
    double angle_sum = 0;
    int min_ndx = -1;
    int max_ndx = -1;

    assert(position.size() >= 2);

    for (size_t i=1; i<position.size()-1; i++) {
        const float3 ap = position[i-1];
        const float3 bp = position[i];
        const float3 cp = position[i+1];

        const float2 a1 = project(mul(xform, ap));
        const float2 b1 = project(mul(xform, bp));
        const float2 c1 = project(mul(xform, cp));

        const float2 dif1 = b1-a1;
        const float2 dif2 = c1-b1;

        float a_sq = dot(dif1,dif1);
        float b_sq = dot(dif2,dif2);

        if (a_sq == 0 || b_sq == 0) {
            // Join or cap
        } else {
            const float2 a = normalize(dif1);
            const float2 b = normalize(dif2);
            const float c = determinant(float2x2(b, a));

            const float s = asin(c);
            const float ds = degrees(s);
            const float ads = fabs(ds);

            const float v = dot(a,b);
            if (v <= 0) {
                // Encountered cusp.

                // Ignore this gradient angle because "inside" the elbow of the stroke.
                if (loud) {
                    printf("%d: ignore (%g)\n", int(i), double(degrees(acos(v))));
                }
            } else {
                const float abs_v = fabs(v);
                if (abs_v < min_abs_cos) {
                    min_ndx = int(i);
                    min_abs_cos = abs_v;
                }
                if (abs_v > max_abs_cos) {
                    max_ndx = int(i);
                    max_abs_cos = abs_v;
                }
                const float angle = degrees(acos(abs_v));
                angle_sum += angle;
                result.angle.push_back(angle);
                if (loud) {
                    printf("%d, %g\n", int(i), ads);
                }
            }
        }
    }

    result.min_ndx = min_ndx;
    result.min_abs_cos = min_abs_cos;
    result.min_angle = acos(min_abs_cos);

    result.max_ndx = max_ndx;
    result.max_abs_cos = max_abs_cos;
    result.max_angle = acos(max_abs_cos);

    result.avg_angle = float(angle_sum / (position.size()-2));
}

static void reportStrokeResultInstance(const char* info_name, const LineStripInfo &info)
{
    float c = info.min_abs_cos;
    float r = info.min_angle;
    float d = degrees(r);
    printf("  %s.min_abs_cos = %g, acos(min_abs_cos) = %g, degrees = %g\n", info_name, c, r, d);
    c = info.max_abs_cos;
    r = info.max_angle;
    d = degrees(r);
    printf("    %s.max_abs_cos = %g, acos(max_abs_cos) = %g, degrees = %g\n", info_name, c, r, d);
    printf("    %s.angle_avg = %g\n", info_name, info.avg_angle);
}

static void analyzeStrokeResult(const CPUStrokeTessellator& stroker, StrokingResults& results)
{
    assert(results.Ppos.size() >= 2);
    assert(results.Pneg.size() >= 2);

    LineStripInfo result_polar_offset;
    processDirectionLineStrip(results.polar_offset, false, result_polar_offset);

    LineStripInfo result_normal;
    processDirectionLineStrip(results.normal, true, result_normal);

    const float3x3 xform = scale(stroker.ViewportScale());

    LineStripInfo result_Ppos;
    processPositionLineStrip(results.Ppos, xform, result_Ppos);

    results.Ppos_min_ndx = result_Ppos.min_ndx;
    results.Ppos_max_ndx = result_Ppos.max_ndx;

    results.Ppos_min_angle = result_Ppos.min_angle;
    results.Ppos_max_angle = result_Ppos.max_angle;

    LineStripInfo result_Pneg;
    processPositionLineStrip(results.Pneg, xform, result_Pneg);

    results.Pneg_min_ndx = result_Pneg.min_ndx;
    results.Pneg_max_ndx = result_Pneg.max_ndx;

    results.Pneg_min_angle = result_Pneg.min_angle;
    results.Pneg_max_angle = result_Pneg.max_angle;

    if (force_analyze || verbose) {
        if (result_Ppos.min_abs_cos > result_Ppos.max_abs_cos) {
            printf("  JOIN/CAP count = %d\n", int(results.polar_offset.size()));
        } else {
            printf("A = %g, B = %g\n", stroker.ObjectToWindowSpaceMatrix()[2][0], stroker.ObjectToWindowSpaceMatrix()[2][0]);
            printf("count = %d\n", int(results.polar_offset.size()));
            reportStrokeResultInstance("polar_offset", result_polar_offset);
            reportStrokeResultInstance("Ppos", result_Ppos);
            reportStrokeResultInstance("Pneg", result_Pneg);
            reportStrokeResultInstance("normal", result_normal);
        }
        force_analyze = false;
    }
}

float dThetaProj_dTheta(double A, double B, double X, double Y, double x, double y, double cos_theta, double sin_theta)
{
    const double Z = A*A + B*B;
    const double numer = 1+A*X-A*y+B*Y+((B*y-A*x)*sin_theta+(B*x+A*y)*cos_theta)*cos_theta;
    const double denom = 1+(2*B+Z*Y)*Y+((2*A*sin_theta-Z*Y)*Y+(2*Z*sin_theta*Y+2*B*sin_theta)*X+(-2*B*Y+(2*A+Z*X)*X)*cos_theta)*cos_theta;
    const double result = numer/denom;
    return float(result);
}

void patchControlPoints(float3 link_control_points[6], const StrokeLink& link)
{
    for (int j=0; j<6; j++) {
        const GLuint ndx = link.ndx[j];
        const float x = simple.indexed_path_state.Coord(ndx+0);
        const float y = simple.indexed_path_state.Coord(ndx+1);
        link_control_points[j] = float3(x,y,1);
    }
    // Is this patch encoding a quadratic Bezier segment AND is the
    // extrapolating control point index odd to indicate
    // a rational extrapolating control point for a conic?
    if ((link.ndx[2] == 0) && (link.ndx[1] & 1)) {
        // Yes, fetch the 3rd component and use it as a homogeneous coordinate.
        const GLuint ndx = link.ndx[1];
        const float w = simple.indexed_path_state.Coord(ndx+2);
        link_control_points[1].z = w;
        // Also fix up the extra directional control points.
        assert(link.ndx[4] == ndx);
        assert(link.ndx[5] == ndx);
        link_control_points[4].z = w;
        link_control_points[5].z = w;
    } else {
        // Otherwise caps and joins may also reference a rational
        // extrapolating quadratic Bezier control point.
        // Check the last two indices of the link of oddness...
        for (size_t i=4; i<countof(link.ndx); i++) {
            // Is either of the last two indexes in the link odd?
            const GLuint ndx = link.ndx[i];
            if (ndx & 1) {
                // Should be a cap or join link.
                assert(link.ndx[1] == 0);
                // Yes, fetch the 3rd component and use it as a homogeneous coordinate.
                const float w = simple.indexed_path_state.Coord(ndx+2);
                link_control_points[i].z = w;
            }
        }
    }
}

void coverSimplePathLinkWithCPU(StrokerContext& context, const StrokeLink& link, PathDrawState& state, LinkRenderStats& stats, unsigned int primitive_id)
{
    // Collected "patch" worth of 6 indexed control points.
    float3 link_control_points[6];
    patchControlPoints(link_control_points, link);

    StrokingResults results;
#if SUPPORT_POLAR_VISUALIZE
    if (stroke_tessellator_settings.isUniformStroking()) {
        if (link.ndx[1] == 0) {
            // Skip joins and caps.
        } else {
            stats.incrementLinksRendered();
            const bool cull = evalUniformStrokingCPU(results, link_control_points, link.ndx);
            if (cull) {
                // Draw nothing if culled.
                stats.incrementLinksTriviallyCulled();
            } else {
                coverUniformStrokingCPU(results);
            }
        }
    } else
#endif
    {
        stats.incrementLinksRendered();
        CPUStrokeTessellator stroker(stroke_tessellator_settings, context.getShadingMode(), state);
        stroker.setPassNumber(1);  // Cover pass is pass one.
        const bool cull = evalPolarStrokingCPU(stroker, results, link_control_points, link.ndx, primitive_id);
        if (cull) {
            // Draw nothing if culled.
            stats.incrementLinksTriviallyCulled();
        } else {
            analyzeStrokeResult(stroker, results);
            // Draw into clip space; assume modelview and projection matrices are identity.
            coverPolarStrokingCPU(stroker, results, state);
        }
    }
    const int quads_generated = results.numQuads();
    stats.addQuadsGenerated(quads_generated);
}

void stencilSimplePathLinkWithCPU(StrokerContext& context, const StrokeLink& link, PathDrawState& state, LinkRenderStats& stats, unsigned int primitive_id)
{
    // Collected "patch" worth of 6 indexed control points.
    float3 link_control_points[6];
    patchControlPoints(link_control_points, link);

    StrokingResults results;
#if SUPPORT_POLAR_VISUALIZE
    if (stroke_tessellator_settings.isUniformStroking()) {
        if (link.ndx[1] == 0) {
            // Skip joins and caps.
        } else {
            stats.incrementLinksRendered();
            const bool cull = evalUniformStrokingCPU(results, link_control_points, link.ndx);
            if (cull) {
                // Draw nothing if culled.
                stats.incrementLinksTriviallyCulled();
            } else {
                stencilUniformStrokingCPU(results);
            }
        }
    } else 
#endif
    {
        stats.incrementLinksRendered();
        CPUStrokeTessellator stroker(stroke_tessellator_settings, context.getShadingMode(), state);
        stroker.setPassNumber(0);  // Stencil pass is pass zero.
        stroker.setPathOffset(cpu_path_count);
        const bool cull = evalPolarStrokingCPU(stroker, results, link_control_points, link.ndx, primitive_id);
        if (cull) {
            // Draw nothing if culled.
            stats.incrementLinksTriviallyCulled();
        } else {
            analyzeStrokeResult(stroker, results);
            // Draw into clip space; assume modelview and projection matrices are identity.
            stencilPolarStrokingCPU(stroker, results, state);
        }
    }
    const int quads_generated = results.numQuads();
    stats.addQuadsGenerated(quads_generated);
}

void annotateSimplePathLinkWithCPU(StrokerContext& context, const StrokeLink& link, /*unused*/PathDrawState& state, LinkRenderStats& stats, unsigned int primitive_id)
{
    // Collected "patch" worth of 6 indexed control points.
    float3 link_control_points[6];
    patchControlPoints(link_control_points, link);

    StrokingResults results;
#if SUPPORT_POLAR_VISUALIZE
    if (stroke_tessellator_settings.isUniformStroking()) {
        if (link.ndx[1] == 0) {
            // Skip joins and caps.
        } else {
            stats.incrementLinksRendered();
            bool cull = evalUniformStrokingCPU(results, link_control_points, link.ndx);
            if (cull) {
                // Draw nothing if culled.
                stats.incrementLinksTriviallyCulled();
            } else {
                annotateUniformStrokingCPU(results);
            }
        }
    } else
#endif
    {
        stats.incrementLinksRendered();
        CPUStrokeTessellator stroker(stroke_tessellator_settings, context.getShadingMode(), state);
        stroker.setPassNumber(-1);  // Treat annotation pass as pass zero.
        const bool cull = evalPolarStrokingCPU(stroker, results, link_control_points, link.ndx, primitive_id);
        if (cull) {
            // Draw nothing if culled.
            stats.incrementLinksTriviallyCulled();
        } else {
            assert(results.Ppos.size() >= 2);
            assert(results.Pneg.size() >= 2);
            analyzeStrokeResult(stroker, results);
            // Draw into clip space; assume modelview and projection matrices are identity.
            annotatePolarStrokingCPU(stroker, results);
        }

    }
    const int quads_generated = results.numQuads();
    stats.addQuadsGenerated(quads_generated);
}

// processPathContainerWithCPU iterates a function over all links in a PathContainer

typedef void (*PathLinkFunc)(StrokerContext& context, const StrokeLink& link, PathDrawState& state, LinkRenderStats& stats, unsigned int primitive_id);

void processPathContainerWithCPU(StrokerContext& context, const PathContainer& path, PathLinkFunc func, PathDrawState& state, LinkRenderStats& stats)
{
    // CPU path:
    const size_t num_links = path.indexed_path_state.numLinks();
    for (size_t i=0; i<num_links; i++) {
        const StrokeLink& link = path.indexed_path_state.Link(i);
        const unsigned int primitive_id = static_cast<unsigned int>(i);
        func(context, link, state, stats, primitive_id);
    }
    cpu_path_count += static_cast<unsigned int>(path.indexed_path_state.numLinks());
}

void stencilSimplePathWithCPU(StrokerContext& context, PathDrawState& state, LinkRenderStats& stats)
{
    // FIRST "stencil" step
    glColorMask(0,0,0,0);  // disable color writes
    // Write 1 to stencil buffer
    glEnable(GL_STENCIL_TEST);
    const StrokeConfigTypes::CompositeMode compositing = context.getCompositeMode();
    if (compositing >= StrokeConfigTypes::COMPOSITE_BLEND_ONCE_LAST) {
        glStencilOp(GL_KEEP, GL_KEEP, GL_INCR);
    } else {
        glStencilOp(GL_KEEP, GL_KEEP, GL_REPLACE);
    }
    glStencilFunc(GL_ALWAYS, 1, ~0);
    glStencilMask(~0);

    processPathContainerWithCPU(context, simple, stencilSimplePathLinkWithCPU, state, stats);

    assert(glIsEnabled(GL_STENCIL_TEST));
}

void coverSimplePathWithCPU(StrokerContext& context, PathDrawState& state, LinkRenderStats& stats)
{
    // Stencil testing should still be enabled from "stencil" step.
    assert(glIsEnabled(GL_STENCIL_TEST));

    // SECOND "cover" step
    // If stencil equals reference, draw color and reset stencil
    if (view_conservative_rasterization) {
        glStencilFunc(GL_ALWAYS, 1, ~0);
        glStencilOp(GL_KEEP, GL_KEEP, GL_ZERO);
    } else {
        const StrokeConfigTypes::CompositeMode compositing = g_stroker_context.getCompositeMode();
        if (compositing >= StrokeConfigTypes::COMPOSITE_BLEND_ONCE_LAST) {
            glStencilFunc(GL_EQUAL, 1, ~0);
            glStencilOp(GL_DECR /*saturating!*/, GL_KEEP, GL_ZERO);
        } else {
            glStencilFunc(GL_EQUAL, 1, ~0);
            glStencilOp(GL_KEEP, GL_KEEP, GL_ZERO);
        }
    }
    glColorMask(1,1,1,1);  // re-enable color writes

    if (fb->isMixedSamples()) {
        glEnable(GL_CONSERVATIVE_RASTERIZATION_NV);
    }

    processPathContainerWithCPU(context, simple, coverSimplePathLinkWithCPU, state, stats);

    if (fb->isMixedSamples()) {
        glDisable(GL_CONSERVATIVE_RASTERIZATION_NV);
    }

    // Stop stencil testing
    glDisable(GL_STENCIL_TEST);
}

void annotateSimplePathWithCPU(StrokerContext& context, PathDrawState& state, LinkRenderStats& stats)
{
    processPathContainerWithCPU(context, simple, annotateSimplePathLinkWithCPU, state, stats);
}

void drawSimpleStencilThenCoverPathWithCPU(StrokerContext& context, LinkRenderStats& stats)
{
    PathDrawState state;

    const StrokeConfigTypes::CompositeMode compositing = context.getCompositeMode();
    if (compositing >= StrokeConfigTypes::COMPOSITE_BLEND_ONCE) {
        glBlendFunc(GL_ONE, GL_ONE_MINUS_SRC_ALPHA);  // "over" blending
        glEnable(GL_BLEND);

        LinkRenderStats stencil_stats;
        LinkRenderStats annotate_stats;

        stencilSimplePathWithCPU(context, state, stencil_stats);
        coverSimplePathWithCPU(context, state, stats);
        annotateSimplePathWithCPU(context, state, annotate_stats);

        // Make sure stats from different passes are consistent.
        assert(stencil_stats == stats);
        assert(annotate_stats == stats);

        if (dump_buffers && StrokeConfigTypes::shadingModeNeedsArcLength(getShadingMode())) {
            state.dumpBuffers(simple.indexed_path_state.numLinks());
        }
        glDisable(GL_BLEND);
    } else if (compositing == StrokeConfigTypes::COMPOSITE_BLEND) {
        glBlendFunc(GL_ONE, GL_ONE_MINUS_SRC_ALPHA);  // "over" blending
        glEnable(GL_BLEND);
        processPathContainerWithCPU(context, simple, coverSimplePathLinkWithCPU, state, stats);
        glDisable(GL_BLEND);
    } else {
        assert(compositing == StrokeConfigTypes::COMPOSITE_OFF);
        glDisable(GL_BLEND);
        processPathContainerWithCPU(context, simple, coverSimplePathLinkWithCPU, state, stats);
    }
}

static void drawSimplePathWithCPU(LinkRenderStats& stats)
{
#if SUPPORT_POLAR_VISUALIZE
    if (stroke_tessellator_settings.isUniformStroking()) {
        drawSimpleStencilThenCoverPathWithCPU(g_stroker_context, stats);
    } else
#endif
    {
        g_stroker_context.setArcLengthOffset(arc_length_offset);
        g_stroker_context.setArcLengthTexGenScale(getCurrentTextureScale());

        // Draw into clip space; set modelview and projection matrices to identity.
        glMatrixPushEXT(GL_MODELVIEW); {
            glMatrixPushEXT(GL_PROJECTION); {
                glMatrixLoadIdentityEXT(GL_PROJECTION);
                glMatrixLoadIdentityEXT(GL_MODELVIEW);

                drawSimpleStencilThenCoverPathWithCPU(g_stroker_context, stats);

            } glMatrixPopEXT(GL_PROJECTION);
        } glMatrixPopEXT(GL_MODELVIEW);
    }
}

static void drawSimpleControlPoints()
{
    fb->pointSize(10);
    int numOfVertices = int(inputPoints.size());
    float green_scale = 1/float(numOfVertices-1);
    for (size_t i=0; i<inputPoints.size(); i++){
        sRGB_glColor3f(1-green_scale*i,1,0);    // Yellow for first; green for last
        glBegin(GL_POINTS); {
            glVertex2f(inputPoints[i].x, inputPoints[i].y);
        } glEnd();
    }
}

static void drawSimpleSplinePoints()
{
    fb->pointSize(8);
    const std::vector<float2> cp = simple.CatmullRomExtraPoints();
    glBegin(GL_POINTS); {
        for (size_t i=0; i<cp.size(); i++) {
            if (i&1) {
                sRGB_glColor3f(0.5f,0.5f,1);    // white-blue
            } else {
                sRGB_glColor3f(1,0.5f,0.5f);    // white-red
            }
            glVertex2f(cp[i]);
        }
    } glEnd();
}

static void drawSimpleHullLines()
{
    if (isBackgroundWhite()) {
        sRGB_glColor3f(0.4,0.2,0.4);  // magenta
    } else {
        sRGB_glColor3f(1,0.5,1);  // magenta
    }
    fb->disableCoverageModulation();
    fb->lineWidth(1.0f);
    glBegin(GL_LINE_STRIP); {
        for (size_t i=0; i<inputPoints.size();i++){
            glVertex2f(inputPoints[i].x, inputPoints[i].y);
        }
    } glEnd();
    fb->blendConfiguration();
    fb->enableCoverageModulation();
}

void drawSimplePath(LinkRenderStats& stats)
{
    if (curveTypeSelect){
        initPoints(inputPoints);
    }

    if (use_NV_path_rendering) {
        drawSimplePathWithNVpr();
    } else {
        if (outline) {
            fb->lineWidth(1.0f);
            glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
        } else {
            glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
        }

        drawSimplePathWithCPU(stats);
        if (!outline) {
            glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
        }
    }

    if (decoration & DRAW_HULL_LINES) {
        drawSimpleHullLines();
    }
    if (decoration & DRAW_CONTROL_POINTS) {
        drawSimpleControlPoints();
    }
    if (decoration & DRAW_SPLINE_POINTS) {
        drawSimpleSplinePoints();
    }
}

void updateSpline()
{
    if (splineType == 'G') {
        printf("no updates of SVG grammar specifies paths allowed");
    } else if (splineType == 'C') {
        simple.updateCubicSpline(inputPoints);
    } else if (splineType == 'O') {
        simple.updateCatmullRomLoop(inputPoints);
    } else if (splineType == 'S') {
        simple.updateCatmullRomSpline(inputPoints);
    } else if (splineType == 'Q') {
        simple.updateQuadraticSpline(inputPoints);
    } else if (splineType == 'L') {
        simple.updateLineSpline(inputPoints);
    } else if (splineType == 'W') {
        simple.updateConicSpline(inputPoints);
    } else {
        assert(!"unknown splineType");
    }
}

void updatePoint(int ndx, const float2 p)
{
    inputPoints[ndx].xy = p;
    updateSpline();
}

void toggleAnnotations()
{
    draw_annotations = !draw_annotations;
    printf("draw_annotations = %d\n", draw_annotations);
    postMainWindowRedisplay();
}

void toggleNormalExtensions()
{
    draw_normal_extensions = !draw_normal_extensions;
    printf("draw_normal_extensions = %d\n", draw_normal_extensions);
    postMainWindowRedisplay();
}

