$ g++ -o spline spline.cpp -lGL -lGLU -lglut

#include <GL/gl.h>
#include <GL/glu.h>
#include <GL/glut.h>
#include <stdlib.h>

// ---- spline stuff ----------------------------------------------------------
//
// inspired by Paul Bourke (http://astronomy.swin.edu.au/~pbourke)
//

typedef struct
{
    float x;
    float y;
    float z;
} point_t;

void spline_knots(int * u, int n, int t)
{
    int m = n+t+1;
    for (int j = 0; j < m; ++j)
    {
        if (j < t) u[j] = 0;
        else u[j] = (j <= n) ? (j - t + 1) : (n - t + 2);
    }
}

double spline_blend(int k, int t, int * u, double v)
{
    if (t == 1) return ((u[k] <= v) && (v < u[k+1])) ? 1.0 : 0.0;
    if ((u[k+t-1] == u[k]) && (u[k+t] == u[k+1]))
        return 0.0;
    else if (u[k+t-1] == u[k])
        return (u[k+t] - v) / (u[k+t] - u[k+1]) * spline_blend(k+1, t-1, u, v);
    else if (u[k+t] == u[k+1])
        return (v - u[k]) / (u[k+t-1] - u[k]) * spline_blend(k , t-1, u, v);
    else return (v - u[k]) / (u[k+t-1] - u[k]) * spline_blend(k, t-1, u, v)
        + (u[k+t] - v) / (u[k+t] - u[k+1]) * spline_blend(k+1, t-1, u, v);
}

void spline_point(int * u, int n, int t, double v, point_t * ctrl, point_t * out)
{
    out->x = out->y = out->z = 0.0;
    for (int k = 0; k <= n; ++k)
    {
        double b = spline_blend(k, t, u, v);
        out->x += ctrl[k].x * b;
        out->y += ctrl[k].y * b;
        out->z += ctrl[k].z * b;
    }
}

void spline_curve(point_t * in, int n, int * knot, int t, point_t * out, int res)
{
    double interval = 0.0;
    double increment = (n-t+2)/(double)(res-1);
    for (int i = 0; i < res; ++i)
    {
        spline_point(knot, n, t, interval, in, &(out[i]));
        interval += increment;
    }
    out[res-1] = in[n];
}

point_t pi[] =
{
    {  0.0,  0.0, 0.0 },
    {  1.0,  2.0, 0.0 },
    {  2.0,  1.0, 0.0 },
    {  4.0,  4.0, 2.0 },
    {  4.0,  0.0, 1.0 },
    {  2.0,  0.5, 0.0 },
    { -0.5, -1.0, 0.0 },
    {  0.0,  0.0, 0.0 },
};

#define N (sizeof(pi)/sizeof(point_t)-1)

#define T 3
int knot[N+T+1];

#define RESOLUTION 200
point_t po[RESOLUTION];

// ---- opengl stuff ----------------------------------------------------------

bool grid = true;
bool line_between_points = true;
bool points = true;

void draw_grid(void)
{
    glColor4f(0.2, 0.2, 0.2, 1.0);
    glBegin(GL_LINES);
    for (int x = 0; x < 21; ++x)
    {
        glVertex3f(-5.0 + 0.5 * x, -5.0, 0.0);
        glVertex3f(-5.0 + 0.5 * x,  5.0, 0.0);
    }
    for (int y = 0; y < 21; ++y)
    {
        glVertex3f(-5.0, -5.0 + 0.5 * y, 0.0);
        glVertex3f( 5.0, -5.0 + 0.5 * y, 0.0);
    }
    glEnd();
}

void draw_line_between_points(void)
{
    glColor4f(0.5, 0.5, 0.5, 1.0);
    glBegin(GL_LINE_STRIP);
    for (int i = 0; i <= N; ++i)
        glVertex3fv(reinterpret_cast<const GLfloat *>(&pi[i]));
    glEnd();
}

void draw_points(void)
{
    glColor4f(1.0, 1.0, 1.0, 1.0);
    glPointSize(3.0);
    glBegin(GL_POINTS);
    for (int i = 0; i <= N; ++i)
        glVertex3fv(reinterpret_cast<const GLfloat *>(&pi[i]));
    glEnd();
}

void draw_spline(void)
{
    glColor4f(1.0, 0.5, 0.5, 1.0);
    glBegin(GL_LINE_STRIP);
    for (int i = 0; i < RESOLUTION; ++i)
        glVertex3fv(reinterpret_cast<const GLfloat *>(&po[i]));
    glEnd();
}

void display(void)
{
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    glLoadIdentity();
    glTranslatef(0.0, 0.0, -15.0);
    glRotatef(-75.0, 1.0, 0.0, 0.0);

    if (grid) draw_grid();
    if (line_between_points) draw_line_between_points();
    if (points) draw_points();
    draw_spline();

    glFlush();
    glutSwapBuffers();
}

void reshape(int w, int h)
{
    glViewport(0, 0, w, h);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluPerspective(45.0, (double)w / (double)h, 0.1, 100.0);
    glMatrixMode(GL_MODELVIEW);
}

void keyboard(unsigned char key, int x, int y)
{
    switch (key)
    {
        case 'g':
        case 'G':
            grid = !grid;
            glutPostRedisplay();
            break;

        case 'l':
        case 'L':
            line_between_points = !line_between_points;
            glutPostRedisplay();
            break;

        case 'p':
        case 'P':
            points = !points;
            glutPostRedisplay();
            break;

        case 27: exit(0);
    }
}

int main(int argc, char ** argv)
{
    spline_knots(knot, N, T);
    spline_curve(pi, N, knot, T, po, RESOLUTION);

    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH | GLUT_ALPHA);
    glutInitWindowPosition(0, 0);
    glutInitWindowSize(600, 500);
    glutCreateWindow("Spline 0");
    glutDisplayFunc(display);
    glutReshapeFunc(reshape);
    glutKeyboardFunc(keyboard);
    glClearColor(0.0, 0.0, 0.0, 0.0);
    glEnable(GL_DEPTH_TEST);
    glutMainLoop();
    return 0;
}