
// hodograph.cpp - implementation of drawing hodograph of curve for "hairline" stroking algorithm testbed

// Hodograph is from Greek "hodos" and "grapho," meaning path-writing.  It is a XY plot of curve's tangent vector.

#include <assert.h>
#include <stdio.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/dot.hpp>
#include <Cg/normalize.hpp>
#include <Cg/length.hpp>
#include <Cg/inverse.hpp>
#include <Cg/radians.hpp>
#include <Cg/sign.hpp>

#include "polar/StrokeTessellator.hpp"

using namespace Cg;
using std::vector;

#include "cg4cpp_xform.hpp"

#include "gl_framebuffer.hpp"

#include "polar/cpu/CPUStrokeTessellator.hpp"

#include "PathContainer.hpp"

#include "Paint.hpp"

#include "mouse.hpp"
#include "hairline.hpp"
#include "polar/global_polar_settings.hpp"
#include "simple_path.hpp"
#include "keyinput.hpp"
#include "cmdline.hpp"
#include "glcolor.hpp"
#include "polar/eval_bezier.hpp"

#include "request_vsync.h"

GLFramebuffer *hodo_fb;

MouseState ms_hodo;

int hodograph_window = -1;

bool draw_ref_t = false;
float ref_t = 0.5;

bool show_hodograph = true;
bool show_true_hodograph = false;
bool show_second_order_hodograph = false;
bool show_curvature = false;
bool show_perspective_hodograph = true;
bool show_polar_grid_lines = true;
bool show_hodograph_control_points = true;
bool show_hodograph_inflection_points = false;
bool show_gradient_magnitude_change = false;

void initHodographViewMatrices()
{
    ms_hodo.view = float3x3(1,0,0,
        0,1,0,
        0,0,1);
    ms_hodo.inverse_view = inverse(ms_hodo.view);
}

void hodographConfigureProjection()
{
    float3x3 inverse_projection, viewport;

    viewport = ortho(0,ms_hodo.window_width, 0,ms_hodo.window_height);
    float left = -canvas_width, right = canvas_width, top = -canvas_height, bottom = canvas_height;
#if 0
    if (aspect_ratio > 1) {
        top *= aspect_ratio;
        bottom *= aspect_ratio;
    } else {
        left /= aspect_ratio;
        right /= aspect_ratio;
    }
#endif
    glMatrixLoadIdentityEXT(GL_PROJECTION);
    glMatrixOrthoEXT(GL_PROJECTION, left, right, bottom, top,
        -1, 1);
    ms_hodo.projection = ortho(left, right, bottom, top);
    inverse_projection = inverse_ortho(left, right, bottom, top);
    ms_hodo.view_width = right - left;
    ms_hodo.view_height = bottom - top;
    ms_hodo.win2view = mul(inverse_projection, viewport);
    ms_hodo.view2win = inverse(ms_hodo.win2view);
}

void hodographReshape(int w, int h)
{
    glViewport(0,0,w,h);

    ms_hodo.window_width = w;
    ms_hodo.window_height = h;
    ms_hodo.aspect_ratio = ms_hodo.window_height/ms_hodo.window_width;

    establishFramebuffer(hodo_fb, w, h);

    hodographConfigureProjection();
}

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

static void drawHodographGridLines()
{
    hodo_fb->lineWidth(1.0f);

    glBegin(GL_LINES); {
        const float size = 1000;
        sRGB_glColor3fv(dim_green);

        for (float a=10; a<90; a+=10) {
            float r = radians(a);
            float x = cos(r);
            float y = sin(r);

            glVertex4f(-x*size,-y*size,0,0);
            glVertex2f(0,0);
            glVertex4f(+x*size,+y*size,0,0);
            glVertex2f(0,0);

            glVertex4f(-y*size,+x*size,0,0);
            glVertex2f(0,0);
            glVertex4f(+y*size,-x*size,0,0);
            glVertex2f(0,0);
        }

        sRGB_glColor3fv(green);
        glVertex4f(-size,0,0,0);
        glVertex2f(0,0);
        glVertex4f(+size,0,0,0);
        glVertex2f(0,0);

        glVertex4f(0,-size,0,0);
        glVertex2f(0,0);
        glVertex4f(0,+size,0,0);
        glVertex2f(0,0);
    } glEnd();

    hodo_fb->pointSize(1.0);
    sRGB_glColor3fv(yellow);
    glBegin(GL_POINTS); {
        const float size = 100;
        for (float a=10; a<=90; a+=10) {
            float r = radians(a);
            float x = cos(r);
            float y = sin(r);

            glVertex2f(-x*size,-y*size);
            glVertex2f(+x*size,+y*size);

            glVertex2f(-y*size,+x*size);
            glVertex2f(+y*size,-x*size);
        }
    } glEnd();
}

static void drawCubicHodographReferencePoint(const float3 q[3])
{
    const float t = ref_t;
    const float3 p = t*t*(q[0]-2*q[1]+q[2]) + t*(2*q[1]-2*q[0]) + q[0];

    sRGB_glColor3fv(lime_green);
    glBegin(GL_POINTS); {
        glVertex2f(p.x, p.y);
    } glEnd();
}

static void drawDoublePoint(float3 p,
                            const GLfloat inner_color[3], float inner_size,
                            const GLfloat outer_color[3], float outer_size)
{
    // First, draw outer point
    sRGB_glColor3fv(outer_color);
    hodo_fb->pointSize(outer_size);
    glBegin(GL_POINTS); {
        glVertex2fw(p);
    } glEnd();
    // Second, next inner point
    sRGB_glColor3fv(inner_color);
    hodo_fb->pointSize(inner_size);
    glBegin(GL_POINTS); {
        glVertex2fw(p);
    } glEnd();
}

static void drawPolarHodographReferencePoint(const CPUStrokeTessellator& stroker)
{
    const float2 A = stroker.A();
    const float2 B = stroker.B();
    const float2 C = stroker.C();

    const float t = ref_t;

    if (show_hodograph) {
        const float3 p = float3(A*t*t + B*t + C,1);
        drawDoublePoint(p, black, 3, lime_green, 7);
    }

    if (show_true_hodograph) {
        if (stroker.LinkTypeIsQuadratic()) {
            float3 cps[3];
            stroker.GetQuadraticControlPoints(cps);

            const float3 p = evalRationalBezierQuadraticGradient(cps, t);
            drawDoublePoint(p, black, 3, lime_green, 7);
        } else {
            float3 cps[4];
            stroker.GetCubicControlPoints(cps);

            sRGB_glColor3fv(whitish_blue);
            const float3 p = evalRationalBezierCubicGradient(cps, t);
            drawDoublePoint(p, black, 3, lime_green, 7);
        }
    }

    if (show_second_order_hodograph) {
        const float3 p2 = float3(2*A*t + B,1);
        drawDoublePoint(p2, red, 3, blue, 7);
    }
}

#if SUPPORT_POLAR_VISUALIZE
static void drawUniformCubicHodographCurve(const float3 q[3])
{
    sRGB_glColor3fv(whitish_blue);
    hodo_fb->lineWidth(1.0f);
    glBegin(GL_LINE_STRIP); {
        const unsigned int uniform_segments = stroke_tessellator_settings.UniformSegments();
        assert(uniform_segments >= 2);
        const float dt = 1/float(uniform_segments);
        float t = 0;
        for (unsigned int i = 0; i<=uniform_segments; i++) {
            // Hodograph from quadratic bezier control points:
            const float3 p = t*t*(q[0]-2*q[1]+q[2]) + t*(2*q[1]-2*q[0]) + q[0];
            glVertex2f(p.x, p.y);
            t += dt ;
        }
    } glEnd();
}

static void drawUniformCubicHodographControlPoints(const float3 q[3])
{
    hodo_fb->pointSize(3.0f);
    sRGB_glColor3fv(whitish_blue);
    glBegin(GL_POINTS); {
        const unsigned int uniform_segments = stroke_tessellator_settings.UniformSegments();
        assert(uniform_segments >= 2);
        const float dt = 1/float(uniform_segments);
        float t = 0;
        for (unsigned int i = 0; i<=uniform_segments; i++){
            // Hodograph from quadratic bezier control points:
            const float3 p = t*t*(q[0]-2*q[1]+q[2]) + t*(2*q[1]-2*q[0]) + q[0];
            glVertex2f(p.x, p.y);
            t += dt ;
        }
    } glEnd();

    hodo_fb->pointSize(5.0f);
    glBegin(GL_POINTS); {
        sRGB_glColor3fv(green);
        glVertex2f(q[0].x, q[0].y);
        sRGB_glColor3fv(yellow);
        glVertex2f(q[1].x, q[1].y);
        sRGB_glColor3fv(red);
        glVertex2f(q[2].x, q[2].y);
    } glEnd();
}

static void drawUniformCubicHodograph(const StrokingResults& results)
{
    assert(inputPoints.size() == 4);

    // Assert is an integral cubic.
    assert(inputPoints[0].z == 1);
    assert(inputPoints[1].z == 1);
    assert(inputPoints[2].z == 1);
    assert(inputPoints[3].z == 1);

    // Control point differences of the cubic bezier:
    const float3 q[3] = { inputPoints[1] - inputPoints[0], 
                          inputPoints[2] - inputPoints[1], 
                          inputPoints[3] - inputPoints[2] };

    drawUniformCubicHodographCurve(q);

    drawUniformCubicHodographControlPoints(q);

    if (draw_ref_t) {
        drawCubicHodographReferencePoint(q);
    }
}
#endif

static float2 evalQuadraticHodograph(const float3 q[2], float t)
{
    const float3 p = (1-t)*q[0] + t*q[1];
    return p.xy;
}

#if SUPPORT_POLAR_VISUALIZE
static void drawUniformQuadraticHodographCurve(const float3 q[2])
{
    sRGB_glColor3fv(whitish_blue);
    hodo_fb->lineWidth(1.0f);
    glBegin(GL_LINE_STRIP); {
        const unsigned int uniform_segments = stroke_tessellator_settings.UniformSegments();
        assert(uniform_segments >= 2);
        const float dt = 1/float(uniform_segments);
        float t = 0;
        for (unsigned int i = 0; i<=uniform_segments; i++){
            // Hodograph from linear bezier control points:
            const float2 p = evalQuadraticHodograph(q, t);
            glVertex2f(p.x, p.y);
            t += dt ;
        }
    } glEnd();
}

static void drawUniformQuadraticHodographControlPoints(const float3 q[2])
{
    hodo_fb->pointSize(3.0f);
    sRGB_glColor3fv(whitish_blue);
    glBegin(GL_POINTS); {
        const unsigned int uniform_segments = stroke_tessellator_settings.UniformSegments();
        assert(uniform_segments >= 2);
        const float dt = 1/float(uniform_segments);
        float t = 0;
        for (unsigned int i = 0; i<=uniform_segments; i++){
            // Hodograph from linear bezier control points:
            const float2 p = evalQuadraticHodograph(q, t);
            glVertex2f(p.x, p.y);
            t += dt ;
        }
    } glEnd();

    hodo_fb->pointSize(5.0f);
    glBegin(GL_POINTS); {
        sRGB_glColor3fv(green);
        glVertex2f(q[0].x, q[0].y);
        sRGB_glColor3fv(red);
        glVertex2f(q[1].x, q[1].y);
    } glEnd();
}
#endif

static void drawQuadraticHodographReferencePoint(const float3 q[3])
{
    const float t = ref_t;
    const float2 p = evalQuadraticHodograph(q, t);

    sRGB_glColor3fv(lime_green);
    glBegin(GL_POINTS); {
        glVertex2f(p);
    } glEnd();
}

#if SUPPORT_POLAR_VISUALIZE
static void drawUniformQuadraticHodograph(const StrokingResults& results)
{
    assert(inputPoints.size() == 3);

    // Control point differences of the cubic bezier:
    const float3 q[2] = { inputPoints[1] - inputPoints[0], 
                          inputPoints[2] - inputPoints[1], };

    drawUniformQuadraticHodographCurve(q);

    drawUniformQuadraticHodographControlPoints(q);

    if (draw_ref_t) {
        drawQuadraticHodographReferencePoint(q);
    }
}

static void drawUniformHodograph(const StrokingResults& results)
{
    if (inputPoints.size() == 3) {
        drawUniformQuadraticHodograph(results);
    } else {
        drawUniformCubicHodograph(results);
    }
}
#endif

static void drawPolarStrokingHodographSamples(const StrokingResults& results)
{
    sRGB_glColor3fv(whitish_blue);
    const vector<float2> &gradient = results.gradient;
    glBegin(GL_LINE_STRIP); {
        for (size_t i = 0; i<gradient.size(); i++){
            glVertex2f(gradient[i]);
        }
    } glEnd();

    hodo_fb->pointSize(3.0f);
    sRGB_glColor3fv(whitish_green);
    glBegin(GL_POINTS); {
        for (size_t i = 1; i<gradient.size()-1; i++){
            glVertex2f(gradient[i]);
        }
    } glEnd();
}

static void drawPolarStrokingSecondOrderHodographSamples(const StrokingResults& results)
{
    sRGB_glColor3fv(whitish_red);
    const vector<float2> &gradient2 = results.gradient2;
    glBegin(GL_LINE_STRIP); {
        for (size_t i = 0; i<gradient2.size(); i++){
            glVertex2f(gradient2[i]);
        }
    } glEnd();

    hodo_fb->pointSize(3.0f);
    sRGB_glColor3fv(whitish_purple);
    glBegin(GL_POINTS); {
        for (size_t i = 1; i<gradient2.size()-1; i++){
            glVertex2f(gradient2[i]);
        }
    } glEnd();
}

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

static void drawPolarStrokingCurvatureSamples(const StrokingResults& results)
{
    const vector<float> &curvature = results.curvature;
    const vector<float2> &normal = results.normal;
    assert(curvature.size() == normal.size());
    glBegin(GL_LINE_STRIP); {
        for (size_t i = 0; i<curvature.size(); i++){
            float k = curvature[i];
            if (k > 0) {
                sRGB_glColor3fv(whitish_green);
            } else if (k < 0) {
                sRGB_glColor3fv(whitish_red);
            } else {
                sRGB_glColor3fv(whitish_yellow);
            }
            const float2 gradient_direction = ortho(normal[i]);
            float2 curvature_point = k * gradient_direction;
            curvature_point *= 500.0f;
            glVertex2f(curvature_point);
        }
    } glEnd();

    hodo_fb->pointSize(3.0f);
    sRGB_glColor3fv(whitish_purple);
    glBegin(GL_POINTS); {
        for (size_t i = 0; i<curvature.size(); i++){
            const float k = curvature[i];
            if (k > 0) {
                sRGB_glColor3fv(green);
            } else if (k < 0) {
                sRGB_glColor3fv(red);
            } else {
                sRGB_glColor3fv(yellow);
            }
            const float2 gradient_direction = ortho(normal[i]);
            float2 curvature_point = k * gradient_direction;
            curvature_point *= 500.0f;
            glVertex2f(curvature_point);
        }
    } glEnd();
}

// Compute 2D squared distance
static float sqdist(float2 v)
{ 
    return dot(v,v);
}

// Compute 2D squared distance
static float sqdist(float3 v3)
{ 
    float2 v = v3.xy/v3.z;
    return dot(v,v);
}

static void drawPolarStrokingCurvatureReferencePoint(const CPUStrokeTessellator& stroker)
{
    if (draw_ref_t) {
        float curvature;
        float3 gradient;

        if (stroker.LinkTypeIsQuadratic()) {
            float3 cp[3];
            stroker.GetQuadraticControlPoints(cp);
            gradient = evalRationalBezierQuadraticGradient(cp, ref_t);
            float3 gradient2 = evalRationalBezierQuadraticSecondGradient(cp, ref_t);

            const float curvature_numer = gradient.x*gradient2.y - gradient.y*gradient2.x;
            const float curvature_denom = pow(sqdist(gradient), 1.5f);
            curvature = curvature_numer/curvature_denom;
        } else {
            float3 cp[4];
            stroker.GetCubicControlPoints(cp);
            gradient = evalRationalBezierCubicGradient(cp, ref_t);
            float3 gradient2 = evalRationalBezierCubicSecondGradient(cp, ref_t);
            const float curvature_numer = gradient.x*gradient2.y - gradient.y*gradient2.x;
            const float curvature_denom = pow(sqdist(gradient), 1.5f);
            curvature = curvature_numer/curvature_denom;
        }
        hodo_fb->pointSize(5.0f);

        const float k = curvature;
        if (k > 0) {
            sRGB_glColor3fv(green);  // positive curvature has green points
        } else if (k < 0) {
            sRGB_glColor3fv(red);    // negative curvature has red points
        } else {
            sRGB_glColor3fv(yellow); // zero (or NaN) curvature has yellow points
        }
        const float2 gradient_direction = normalize(gradient.xy/gradient.z);
        float2 curvature_point = k * gradient_direction;
        curvature_point *= 500.0f;
        glBegin(GL_POINTS); {
            glVertex2f(curvature_point);
        } glEnd();
    }
}

static void drawPolarStrokingGradientMagnitudeChange(const CPUStrokeTessellator& stroker, const StrokingResults& results)
{
    if (stroker.LinkTypeIsQuadratic()) {
        float3 cps[3];
        float3 icps[3];
        stroker.GetQuadraticControlPoints(cps);
        stroker.GetQuadraticControlPoints(icps);

        icps[0].z = 1;
        icps[1].x /= icps[1].z;
        icps[1].y /= icps[1].z;
        icps[1].z /= icps[1].z;
        icps[2].z = 1;

        const float w[3] = { cps[0].z, cps[1].z, cps[2].z };

        sRGB_glColor3fv(orange);
        const vector<float> &tarray = results.t;
        const vector<float2> &normal = results.normal;
        glBegin(GL_LINE_STRIP); {
            for (size_t i = 0; i<tarray.size(); i++){
                float t = tarray[i];
                float3 rational_gradient3 = evalRationalBezierQuadraticGradient(cps, t);
                float2 rational_gradient = rational_gradient3.xy/rational_gradient3.z;
                float3 integral_gradient3 = evalRationalBezierQuadraticGradient(icps, t);
                float2 integral_gradient = integral_gradient3.xy/integral_gradient3.z;
                const float d = length(rational_gradient)/length(integral_gradient);
                const float denom = pow(w[0]+(-2*w[0]+2*w[1]+(w[0]-2*w[1]+w[2])*t)*t,2.0f);
                const float2 gradient_direction = ortho(normal[i]);
                glVertex2f(gradient_direction * 100 * d);
            }
        } glEnd();

        hodo_fb->pointSize(3.0f);
        sRGB_glColor3fv(orange_red);
        glBegin(GL_POINTS); {
            for (size_t i = 1; i<tarray.size()-1; i++){
                float t = tarray[i];
                float3 rational_gradient3 = evalRationalBezierQuadraticGradient(cps, t);
                float2 rational_gradient = rational_gradient3.xy/rational_gradient3.z;
                float3 integral_gradient3 = evalRationalBezierQuadraticGradient(icps, t);
                float2 integral_gradient = integral_gradient3.xy/integral_gradient3.z;
                const float d = length(rational_gradient)/length(integral_gradient);
                const float denom = pow(w[0]+(-2*w[0]+2*w[1]+(w[0]-2*w[1]+w[2])*t)*t,2.0f);
                const float2 gradient_direction = ortho(normal[i]);
                glVertex2f(gradient_direction * 100 * d);
            }
        } glEnd();
    } else {
        // XXX do something different for non-quadratics?
    }
}

static void drawPolarStrokingTrueHodograph(const CPUStrokeTessellator& stroker, const StrokingResults& results)
{
    if (stroker.LinkTypeIsQuadratic()) {
        float3 cps[3];
        stroker.GetQuadraticControlPoints(cps);

        sRGB_glColor3fv(whitish_blue);
        const vector<float> &t = results.t;
        glBegin(GL_LINE_STRIP); {
            for (size_t i = 0; i<t.size(); i++){
                const float3 gradient = evalRationalBezierQuadraticGradient(cps, t[i]);
                glVertex2fw(gradient);
            }
        } glEnd();

        hodo_fb->pointSize(3.0f);
        sRGB_glColor3fv(whitish_green);
        glBegin(GL_POINTS); {
            for (size_t i = 1; i<t.size()-1; i++){
                const float3 gradient = evalRationalBezierQuadraticGradient(cps, t[i]);
                glVertex2fw(gradient);
            }
        } glEnd();
    } else {
        // Cubic case
        float3 cps[4];
        stroker.GetCubicControlPoints(cps);

        sRGB_glColor3fv(whitish_blue);
        const vector<float> &t = results.t;
        glBegin(GL_LINE_STRIP); {
            for (size_t i = 0; i<t.size(); i++){
                const float3 gradient = evalRationalBezierCubicGradient(cps, t[i]);
                glVertex2fw(gradient);
            }
        } glEnd();

        hodo_fb->pointSize(3.0f);
        sRGB_glColor3fv(whitish_green);
        glBegin(GL_POINTS); {
            for (size_t i = 1; i<t.size()-1; i++){
                const float3 gradient = evalRationalBezierCubicGradient(cps, t[i]);
                glVertex2fw(gradient);
            }
        } glEnd();
    }
}

static void drawPolarStrokingHodographWindowSpaceSamples(const CPUStrokeTessellator& stroker)
{
    const float scale = 1/3.0f;
    const float3x3 M = stroker.ObjectToWindowSpaceMatrix();

    if (stroker.LinkTypeIsQuadratic()) {
        // Quadratic Bezier or conic case:
        float3 cp_orig[3];
        stroker.GetQuadraticControlPoints(cp_orig);

        float3 cp[3];
        for (size_t i=0; i<countof(cp); i++) {
            cp[i] = mul(M, cp_orig[i]);
        }

        sRGB_glColor3fv(cyan);

        const int steps = stroker.Steps();
        hodo_fb->lineWidth(1.0f);
        glBegin(GL_LINE_STRIP); {
            for (int i=0; i<=steps; i++) {
                const float t = i/float(steps);
                const float3 p = evalRationalBezierQuadraticGradient(cp, t) * scale;
                glVertex2fw(p);
            }
        } glEnd();

        hodo_fb->pointSize(3.0f);
        glBegin(GL_POINTS); {
            for (int i=0; i<=steps; i++) {
                const float t = i/float(steps);
                const float3 p = evalRationalBezierQuadraticGradient(cp, t) * scale;
                glVertex2fw(p);
            }
        } glEnd();

        if (draw_ref_t) {
            const float t = ref_t;
            const float3 p = evalRationalBezierQuadraticGradient(cp, t) * scale;

            drawDoublePoint(p, black, 3, lime_green, 7);
        }
    } else {
        // Cubic case:
        float3 cp[4];
        for (int i=0; i<icountof(cp); i++) {
            cp[i] = mul(M, stroker.Cp(i));
        }

        double x[4], y[4], w[4];

        for (int i=0; i<4; i++) {
            x[i] = cp[i].x;
            y[i] = cp[i].y;
            w[i] = cp[i].z;
        }

        sRGB_glColor3fv(cyan);

        hodo_fb->lineWidth(1.0f);
        glBegin(GL_LINE_STRIP); {
            const int steps = 100;
            for (int i=0; i<=steps; i++) {
                const float t = i/float(steps);
                const float3 p = evalRationalBezierCubicGradient(x, y, w, t) * scale;
                glVertex2fw(p);
            }
        } glEnd();

        hodo_fb->pointSize(3.0f);
        glBegin(GL_POINTS); {
            const int steps = 100;
            for (int i=0; i<=steps; i++) {
                const float t = i/float(steps);
                const float3 p = evalRationalBezierCubicGradient(x, y, w, t) * scale;
                glVertex2fw(p);
            }
        } glEnd();

        if (draw_ref_t) {
            const float t = ref_t;
            const float3 p = evalRationalBezierCubicGradient(x, y, w, t) * scale;

            drawDoublePoint(p, black, 3, lime_green, 7);
        }
    }
}

static void drawRawPolarStrokingHodographInflectionPoints(const CPUStrokeTessellator& stroker)
{
    // Quadratic gradient control points
    const float2& v0 = stroker.GradCp(0);
    const float2& v1 = stroker.GradCp(1);
    const float2& v2 = stroker.GradCp(2);

    glBegin(GL_POINTS); {
        for (int i=0; i<4; i++) {
            float t = stroker.Tbound(i);
            if (t > 0.0f && t < 1.0f) {
                const float2 v = (1-t)*(1-t)*v0 + 2*(1-t)*t*v1 + t*t*v2;
                glVertex2f(v);
            }
        }
    } glEnd();
}

static void drawPolarStrokingHodographInflectionPoints(const CPUStrokeTessellator& stroker, const StrokingResults& results)
{
    if (stroker.LinkTypeIsQuadratic()) {
        return;
    }

    // First, draw inflection points LARGE.
    hodo_fb->pointSize(5.0f);
    sRGB_glColor3fv(white);
    drawRawPolarStrokingHodographInflectionPoints(stroker);

    // Second, overdraw inflection points small for halo effect.
    hodo_fb->pointSize(3.0f);
    sRGB_glColor3fv(blue);
    drawRawPolarStrokingHodographInflectionPoints(stroker);

    // Quadratic gradient control points
    const float2& v0 = stroker.GradCp(0);
    const float2& v1 = stroker.GradCp(1);
    const float2& v2 = stroker.GradCp(2);

    // Draw polar lines for each inflection point.
    for (int i=0; i<4; i++) {
        float t = stroker.Tbound(i);
        if (t <= 0) {
            sRGB_glColor3fv(light_green);
        } else if (t >= 1) {
            sRGB_glColor3fv(light_red);
        } else {
            sRGB_glColor3fv(light_cyan);
        }
        glBegin(GL_LINE_STRIP); {
            const float2 v = (1-t)*(1-t)*v0 + 2*(1-t)*t*v1 + t*t*v2;
            glVertex4f(v.x, v.y, 0, 0);
            glVertex4f(0,0, 0, 1);
            glVertex4f(-v.x, -v.y, 0, 0);
        } glEnd();
    }
}

static void drawPolarStrokingHodographControlPoints(const CPUStrokeTessellator& stroker)
{
    if (stroker.LinkTypeIsQuadratic()) {
        float3 cp[3];
        stroker.GetQuadraticControlPoints(cp);
        // Linear gradient control points
        const float3 v0 = float3(cp[1].z*(cp[1].xy*cp[0].z - cp[0].xy*cp[1].z), cp[0].z*cp[1].z);
        const float3 v1 = float3(cp[1].z*(cp[2].xy*cp[1].z - cp[1].xy*cp[2].z), cp[1].z*cp[2].z);

        // Draw quadratic gradient control points.
        hodo_fb->pointSize(5.0f);
        glBegin(GL_POINTS); {
            sRGB_glColor3fv(green);
            glVertex2fw(v0);

            sRGB_glColor3fv(red);
            glVertex2fw(v1);
        } glEnd();
    } else {

        // Quadratic gradient control points
        const float2 v0 = stroker.GradCp(0);
        const float2 v1 = stroker.GradCp(1);
        const float2 v2 = stroker.GradCp(2);

        // Draw quadratic gradient control points.
        hodo_fb->pointSize(5.0f);
        glBegin(GL_POINTS); {
            sRGB_glColor3fv(green);
            glVertex2f(v0);

            sRGB_glColor3fv(yellow);
            glVertex2f(v1);

            sRGB_glColor3fv(red);
            glVertex2f(v2);
        } glEnd();
    }
}

static void GoodOKBadColor(float dtheta_dTheta)
{
    if (dtheta_dTheta > 1.05) {
        sRGB_glColor3fv(red);  // BAD, under-tessellated
    } else if (dtheta_dTheta < 0.95) {
        sRGB_glColor3fv(yellow);  // OK, but over-tessellated
    } else {
        sRGB_glColor3fv(green);  // GOOD, not too little or too much tessellation
    }
}

static void drawPolarStrokingHodograph(const CPUStrokeTessellator& stroker, const StrokingResults& results)
{
    if (show_hodograph) {
        drawPolarStrokingHodographSamples(results);
    }

    if (show_true_hodograph) {
        drawPolarStrokingTrueHodograph(stroker, results);
    }

    if (show_second_order_hodograph) {
        drawPolarStrokingSecondOrderHodographSamples(results);
    }

    if (show_curvature) {
        drawPolarStrokingCurvatureSamples(results);
        drawPolarStrokingCurvatureReferencePoint(stroker);
    }

    if (show_hodograph_inflection_points) {
        drawPolarStrokingHodographInflectionPoints(stroker, results);
    }

    if (show_hodograph_control_points) {
        drawPolarStrokingHodographControlPoints(stroker);
    }

    if (show_gradient_magnitude_change) {
        drawPolarStrokingGradientMagnitudeChange(stroker, results);
    }

    if (draw_ref_t) {
        drawPolarHodographReferencePoint(stroker);
    }

    if (show_perspective_hodograph) {
        drawPolarStrokingHodographWindowSpaceSamples(stroker);
    }
}

static void drawHodograph()
{
    if (show_polar_grid_lines) {
        drawHodographGridLines();
    }

    // Is there at least one link?
    if (simple.indexed_path_state.numLinks() > 0) {
        // Yes, draw hodograph for first link.
        const StrokeLink& link = simple.indexed_path_state.ClampedLink(current_link);

        // Collected "patch" worth of 6 indexed control points.
        float3 link_control_points[6];
        patchControlPoints(link_control_points, link);

        CPUStrokeTessellator stroker(stroke_tessellator_settings);
        stroker.setTrivialCulling(false);
        const unsigned int primitive_id = 0;
        bool cull = prepStroker(stroker, current_paint, link_control_points, link.ndx, primitive_id);
        if (cull) {
            printf("CPU culled hodograph\n");
            return;
        }

        CPUStrokeTessellator stroker32 = stroker;
        stroker32.setSteps(31);
        StrokingResults results32;
        stroker32.evalRibs(results32);

        StrokingResults results;
        stroker.evalRibs(results);

#if SUPPORT_POLAR_VISUALIZE
        // Draw hodograph (technically gradient magnitude is off by one-third).
        if (stroke_tessellator_settings.isUniformStroking()){
            // Uniform parametric evaluation:
            if (show_hodograph) {
                drawUniformHodograph(results);
            }
        } else
#endif
        {
            drawPolarStrokingHodograph(stroker, results);
        }
    }
}

void hodographDisplay()
{
    updateColorSpace(hodo_fb);
    hodo_fb->bind(); {
        glClear(GL_COLOR_BUFFER_BIT | GL_STENCIL_BUFFER_BIT);
        MatrixLoadToGL(ms_hodo.view);
        drawHodograph();
    } hodo_fb->copy();
    glBindFramebuffer(GL_FRAMEBUFFER, 0);
    glutSwapBuffers();
}

void
hodographMouse(int button, int state, int mouse_space_x, int mouse_space_y)
{
    mouse_space_y = ms_hodo.window_height - mouse_space_y;
    const int modifiers = glutGetModifiers() & (GLUT_ACTIVE_CTRL|GLUT_ACTIVE_SHIFT);

    if (button == GLUT_LEFT_BUTTON) {
        if (state == GLUT_DOWN) { 
            float2 win = float2(mouse_space_x, mouse_space_y);
            {
                ms_hodo.anchor = project(mul(ms_hodo.win2view, float3(win,1)));

                ms_hodo.rotate_x = mouse_space_x;
                ms_hodo.scale_y = mouse_space_y;
                if (!(glutGetModifiers() & GLUT_ACTIVE_CTRL)) {
                    ms_hodo.scaling = 1;
                } else {
                    ms_hodo.scaling = 0;
                }
                if (!(glutGetModifiers() & GLUT_ACTIVE_SHIFT)) {
                    ms_hodo.zooming = 1;
                } else {
                    ms_hodo.zooming = 0;
                }
            }
        } else {
            assert(state == GLUT_UP);
            ms_hodo.zooming = 0;
            ms_hodo.scaling = 0;
        }
    }
    if (button == GLUT_MIDDLE_BUTTON) {
        if (state == GLUT_DOWN) {
            if (modifiers & GLUT_ACTIVE_CTRL) {
                ms_hodo.slide_scale = 1.0/8.0;
            } else {
                ms_hodo.slide_scale = 1.0;
            }
            ms_hodo.slide_x = mouse_space_x;
            ms_hodo.slide_y = mouse_space_y;
            if (modifiers & GLUT_ACTIVE_SHIFT) {
                // If Shiftheld down, confine translation to Y direction
            } else {
                ms_hodo.x_sliding = true;
            }
            ms_hodo.y_sliding = true;
        } else {
            assert(state == GLUT_UP);
            ms_hodo.x_sliding = false;
            ms_hodo.y_sliding = false;
        }
        if (verbose) {
            printf("sliding%s%s: %f,%f\n",
                ms_hodo.x_sliding ? "X" : "",
                ms_hodo.y_sliding ? "Y" : "",
                ms_hodo.slide_x, ms_hodo.slide_y);
        }
    }
}

void hodographMotion(int mouse_space_x, int mouse_space_y)
{
    mouse_space_y = ms_hodo.window_height - mouse_space_y;
    if (ms_hodo.zooming || ms_hodo.scaling) {
        float3x3 t, r, s, m;
        float angle = 0;
        float zoom = 1;
        if (ms_hodo.scaling) {
            angle = 0.3 * (ms_hodo.rotate_x - mouse_space_x) * canvas_width/ms_hodo.window_width;
        }
        if (ms_hodo.zooming) {
            zoom = pow(1.003f, (ms_hodo.scale_y - mouse_space_y) * canvas_height/ms_hodo.window_height);
        }

        t = translate(ms_hodo.anchor);
        r = rotate(angle);
        s = scale(zoom, zoom);

        r = mul(r, s);
        m = mul(t, r);
        t = translate(-ms_hodo.anchor);
        m = mul(m, t);
        ms_hodo.view = mul(m, ms_hodo.view);
        ms_hodo.inverse_view = inverse(ms_hodo.view);

        validateModelview();

        ms_hodo.rotate_x = mouse_space_x;
        ms_hodo.scale_y = mouse_space_y;
        postAllWindowsRedisplay();
    }
    if (ms_hodo.x_sliding || ms_hodo.y_sliding) {
        float x_offset = 0,
              y_offset = 0;

        float3x3 m;

        if (ms_hodo.x_sliding) {
            x_offset = (mouse_space_x - ms_hodo.slide_x) * ms_hodo.view_width/ms_hodo.window_width;
            x_offset *= ms_hodo.slide_scale;
        }
        if (ms_hodo.y_sliding) {
            y_offset = (ms_hodo.slide_y - mouse_space_y) * ms_hodo.view_height/ms_hodo.window_height;
            y_offset *= ms_hodo.slide_scale;
        }
        m = translate(x_offset, y_offset);
        ms_hodo.view = mul(m, ms_hodo.view);
        ms_hodo.inverse_view = inverse(ms_hodo.view);

        validateModelview();

        ms_hodo.slide_y = mouse_space_y;
        ms_hodo.slide_x = mouse_space_x;
        ms_hodo.scale_y = mouse_space_y;
        postAllWindowsRedisplay();
    }
}

void postHodographWindowRedisplay()
{
    if (hodograph_window >= 0) {
        glutPostWindowRedisplay(hodograph_window);
    }
}

bool isHodographFramebuffer(const GLFramebuffer *fb)
{
    return fb == hodo_fb;
}

enum HodographMenuEntry {
    TOGGLE_HODOGRAPH,
    TOGGLE_TRUE_HODOGRAPH,
    TOGGLE_SECOND_ORDER_HODOGRAPH,
    TOGGLE_CURVATURE,
    TOGGLE_PERSPECTIVE_HODOGRAPH,
    TOGGLE_POLAR_GRID_LINES,
    TOGGLE_HODOGRAPH_INFLECTION_POINTS,
    TOGGLE_HODOGRAPH_CONTROL_POINTS,
    TOGGLE_GRADIENT_MAGNITUDE_CHANGE,
    HIDE_HODOGRAPH_WINDOW,
};

static void hideHodographWindow()
{
    assert(hodograph_window >= 0);
    glutSetWindow(hodograph_window);
    delete hodo_fb;
    hodo_fb = NULL;
    glutDestroyWindow(hodograph_window);
    hodograph_window = -1;
}

void hodographMenu(int item)
{
    switch (item) {
    case TOGGLE_HODOGRAPH:
        show_hodograph = !show_hodograph;
        break;
    case TOGGLE_TRUE_HODOGRAPH:
        show_true_hodograph = !show_true_hodograph;
        break;
    case TOGGLE_SECOND_ORDER_HODOGRAPH:
        show_second_order_hodograph = !show_second_order_hodograph;
        break;
    case TOGGLE_CURVATURE:
        show_curvature = !show_curvature;
        break;
    case TOGGLE_PERSPECTIVE_HODOGRAPH:
        show_perspective_hodograph = !show_perspective_hodograph;
        break;
    case TOGGLE_POLAR_GRID_LINES:
        show_polar_grid_lines = !show_polar_grid_lines;
        break;
    case HIDE_HODOGRAPH_WINDOW:
        hideHodographWindow();
        return;
    case TOGGLE_HODOGRAPH_CONTROL_POINTS:
        show_hodograph_control_points = !show_hodograph_control_points;
        break;
    case TOGGLE_HODOGRAPH_INFLECTION_POINTS:
        show_hodograph_inflection_points = !show_hodograph_inflection_points;
        postMainWindowRedisplay();  // Inflection points also drawn on main window.
        break;
    case TOGGLE_GRADIENT_MAGNITUDE_CHANGE:
        show_gradient_magnitude_change = !show_gradient_magnitude_change;
        break;
    default:
        assert(!"bad menu item");
    }
    postHodographWindowRedisplay();
}

static void makeHodographMenu()
{
    int menu_id = glutCreateMenu(hodographMenu);
    glutAddMenuEntry("Toggle hodograph (white)", TOGGLE_HODOGRAPH);
    glutAddMenuEntry("Toggle true hodograph (white)", TOGGLE_TRUE_HODOGRAPH);
    glutAddMenuEntry("Toggle second order hodograph (light purple)", TOGGLE_SECOND_ORDER_HODOGRAPH);
    glutAddMenuEntry("Toggle curvature (red)", TOGGLE_CURVATURE);
    glutAddMenuEntry("Toggle perspective hodograph (cyan)", TOGGLE_PERSPECTIVE_HODOGRAPH);
    glutAddMenuEntry("Toggle hodograph control points (green-yellow-red dots)", TOGGLE_HODOGRAPH_CONTROL_POINTS);
    glutAddMenuEntry("Toggle hodograph inflection points (blue line/dot)", TOGGLE_HODOGRAPH_INFLECTION_POINTS);
    glutAddMenuEntry("Toggle polar grid lines", TOGGLE_POLAR_GRID_LINES);
    glutAddMenuEntry("Toggle gradient magnitude change (orange, quadratic only)", TOGGLE_GRADIENT_MAGNITUDE_CHANGE);
    glutAddMenuEntry("Hide window", HIDE_HODOGRAPH_WINDOW);
    glutAttachMenu(GLUT_RIGHT_BUTTON);
}

static void hodographKeyboard(unsigned char c, int x, int y)
{
    switch (c) {
    case '+':
    case '=':
        advanceCurrentLink();
        return;
    case '-':
    case '_':
        reverseCurrentLink();
        return;
    default:
        keyboard(c, x, y);
    }
}

void toggleHodographWindow()
{
    if (hodograph_window < 0){
        // Next window should appear to the right
        setMainWindow();
        glutInitWindowSize(canvas_width, canvas_height);
        glutInitWindowPosition(glutGet(GLUT_WINDOW_X)+glutGet(GLUT_WINDOW_WIDTH)+20, glutGet(GLUT_WINDOW_Y));
        hodograph_window = glutCreateWindow("Hodograph for the Cubic Bezier.");
        glutDisplayFunc(hodographDisplay);
        glutReshapeFunc(hodographReshape);
        glutKeyboardFunc(hodographKeyboard);
        glutSpecialFunc(special);
        glutMouseFunc(hodographMouse);
        glutMotionFunc(hodographMotion);
        makeHodographMenu();
        requestSynchornizedSwapBuffers(frame_sync);
        glClearColor(0,0,0,0);
    } else{
        hideHodographWindow();
    }
}

bool isHodoFramebuffer(const GLFramebuffer *this_fb)
{
    return this_fb == hodo_fb;
}
