/* Bezier Spline methods from https://ovpwp.wordpress.com/2008/12/17/how-to-draw-a-smooth-curve-through-a-set-of-2d-points-with-bezier-methods/ */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <GL/glut.h>		// GLUT, include glu.h and gl.h

#define DIMENSIONS 2		/* add as many as you need ... */
#define X 0
#define Y 1

#define FIRST 0			/* control points per knot */
#define SECOND 1

#define MAX_POINTS 300

double letter[DIMENSIONS*MAX_POINTS];   // [x][y]   [x][y]   [x][y]   [x][y]  ...
#define getletter(dim, index) letter[((index)*2)+(dim)]
#define setletter(dim, index, val) letter[((index)*2)+(dim)] = (val)

double ControlPoint[2*DIMENSIONS*MAX_POINTS];  // [first x][first y][second x][second y]   [first x][first y][second x][second y]   ...
#define getpt(fs, dim, index) ControlPoint[((index)*4)+((fs)*2)+(dim)]
#define setpt(fs, dim, index, val) ControlPoint[((index)*4)+((fs)*2)+(dim)] = (val)

double xy[DIMENSIONS*MAX_POINTS];
#define getxy(dim, index) xy[(index)*2+(dim)]
#define setxy(dim, index, val) xy[(index)*2+(dim)] = (val)

double rhs[MAX_POINTS];
double tmp[MAX_POINTS];
double cpt[DIMENSIONS];
double b, t1,t2,t3,t4;

int i, n, step, dim, order;

// used in bezier procedure:
double s;
double AB;
double BC;
double CD;
double ABC;
double BCD;
double bezier_delta;
double bezier_result;

/* The procedure below is based on */
/* http://stackoverflow.com/questions/5443653/opengl-coordinates-from-bezier-curves */
/* - it's a 1D Bezier than can be called twice to create a 2D Bezier spline and */
/* 3 times to create a 3D Bezier spline etc.  Each axis is independent. */

void bezier (double A,	/* Start value */
             double B,	/* First control value */
	     double C,	/* Second control value */
	     double D,	/* Ending value */
	     int t, int max_slots)	/* eg 0..255/256 slots */
{				/* Parameter 0 <= t <= 1 */
   s = (max_slots - 1) - t;
   AB = (A * s + B * t) / max_slots;
   BC = (B * s + C * t) / max_slots;
   CD = (C * s + D * t) / max_slots;
   // ABC = (AB*s + CD*t)/max_slots; - this 'bug fix' was wrong...
   ABC = (AB * s + BC * t) / max_slots;	// original code was correct...
   BCD = (BC * s + CD * t) / max_slots;
   bezier_delta = ABC - BCD;
   bezier_result = (ABC * s + BCD * t) / max_slots;
}

void curve (int p1, int p2)
{
   for (step = 0; step < 512; step++) {
     for (dim = X; dim < DIMENSIONS; dim++) {
        t1 = getletter(dim, p1);
        t2 = getpt(FIRST, dim, p1);
        t3 = getpt(SECOND, dim, p1);
        t4 = getletter(dim, p2);
	bezier (t1, t2, t3, t4, step, 512);
        cpt[dim] = bezier_result;
      }
      glVertex2f (cpt[X], cpt[Y]);
      // fprintf(stderr, " (%f,%f)\n", xy[X], xy[Y]);
   }
}

/*  Get open-ended Bezier Spline Control Points. */
static void DrawCurve (int Length)
{
   n = Length - 1;
   if (n < 1) return;

   /* Calculate Bezier control points */
   for (dim = X; dim < DIMENSIONS; dim++) {
      /* Set right hand side X or Y values */
      for (i = 1; i < n - 1; ++i) {
        t1 = getletter(dim, i);
	t2 = getletter(dim, i + 1);
        rhs[i] = 4 * t1 + 2 * t2;
      }
      t1 = getletter(dim, 0);
      t2 = getletter(dim, 1); 
      rhs[0] = t1 + 2 * t2;
      t1 = getletter(dim, n - 1);
      rhs[n - 1] = 3 * t1;
      /*  Solves a tridiagonal system for one of coordinates of Bezier control points. */
      b = 2.0;
      setxy(dim, 0, rhs[0] / b);
      for (i = 1; i < n; i++) {	/* Decomposition and forward substitution. */
         tmp[i] = 1 / b; 
         if (i < n - 1) t1 = 4.0; else t1 = 2.0; 
         b = t1 - tmp[i];
         t1 = getxy(dim, i - 1);
         setxy(dim, i, (rhs[i] - t1) / b);
      }
      for (i = 1; i < n; i++) {
	 t1 = getxy(dim, n - i - 1);
	 t2 = getxy(dim, n - i);	/* Backsubstitution. */
	 setxy(dim, n - i - 1, t1 - tmp[n - i] * t2);	/* Backsubstitution. */
      }
   }

   for (i = 0; i <= n; ++i) {
      for (dim = X; dim < DIMENSIONS; dim++) {
	t1 = getxy(dim, i);
	setpt(FIRST, dim, i, t1);
	if (i < n - 1) {
	   t1 = getletter(dim, i + 1);
	   t2 = getxy(dim, i + 1);
	   t3 = 2 * t1 - t2;
	} else {
	   t1 = getletter(dim, n);
	   t2 = getxy(dim, n - 1);
	   t3 = (t1 + t2) / 2;
	}
	setpt(SECOND, dim, i, t3);
      }
      if (i > 0) curve (i-1, i);
   }
}

void bitmap_output (int x, int y, char *string, void *font)
{
   glRasterPos2f (x, y); while (*string) glutBitmapCharacter (font, *string++);
}

void display ()
{
   static int odd = 0;

   odd = (odd + 1) & 15;
   if (odd == 0) {
      glClearColor (0.0f, 0.0f, 0.0f, 1.0f);	// Set background color to black and opaque
      glClear (GL_COLOR_BUFFER_BIT);	// Clear the color buffer (background)
      glLineWidth (6.0);
      glColor3f (1.0, 1.0, 1.0);
      glBegin (GL_LINE_STRIP);
      
      /* Trivial test */
      setletter(X, 0, 10 * 4);
      setletter(Y, 0, 0 * 4 + 20);
      setletter(X, 1, 10 * 4);
      setletter(Y, 1, 50 * 4 + 20);
      setletter(X, 2, 28 * 4);
      setletter(Y, 2, 50 * 4 + 20);
      setletter(X, 3, 30 * 4);
      setletter(Y, 3, 0 * 4 + 20);
      setletter(X, 4, 10 * 4);
      setletter(Y, 4, 0 * 4 + 20);
      setletter(X, 5, 10 * 4);
      setletter(Y, 5, 50 * 4 + 20);

      DrawCurve (4);

      glEnd ();
      bitmap_output (40, 275, "Letter 'n'", GLUT_BITMAP_9_BY_15);
      glFlush ();		// Render now
      glutSwapBuffers ();
   }
}

void idle (void)
{
   glutPostRedisplay ();
}

void reshape (int width, int height)
{
   glViewport (0, 0, 480, 360);
   glMatrixMode (GL_PROJECTION);
   glLoadIdentity ();
   glOrtho (0, 480, 0, 360, -1, 1);
   glMatrixMode (GL_MODELVIEW);
   glLoadIdentity ();
   glEnable (GL_BLEND);
   glEnable (GL_LINE_SMOOTH);
   glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
}

/* Main function: GLUT runs as a console application starting at main()  */
int main (int argc, char **argv)
{
   glutInitDisplayMode (GLUT_DOUBLE | GLUT_RGBA);
   glutInitWindowPosition (50, 50);
   glutInitWindowSize (480, 360);

   glutInit (&argc, argv);	// Initialize GLUT

   glutCreateWindow ("Bezier path through several points");	// Create a window with the given title
   glutInitWindowSize (480, 360);	// Set the window's initial width & height
   glViewport (0, 0, 480, 360);
   glMatrixMode (GL_PROJECTION);
   glLoadIdentity ();
   glOrtho (0, 480, 0, 360, -1, 1);
   glMatrixMode (GL_MODELVIEW);
   glLoadIdentity ();
   glEnable (GL_BLEND);
   glEnable (GL_LINE_SMOOTH);
   glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
   glutInitWindowPosition (50, 50);	// Position the window's initial top-left corner
   glutDisplayFunc (display);	// Register display callback handler for window re-paint
   glutReshapeFunc (reshape);
   glutIdleFunc (idle);
   glutMainLoop ();		// Enter the event-processing loop

   exit (0);
   return (0);
}