// lagrange - Hip - 2009-12-13
// provlacci Lagrangeov interpolacijski polinom kroz zadane toccke

#include <GL/glut.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

// maksimalni broj toccaka
#define NTMAX 30

// polumer toccke
#define TPOLUMJER 0.15

GLUquadricObj *disk;

double ktocka[NTMAX][3]; // tabela s tocckama
int ntocaka = 0; // broj toccaka
int akt = -1; // "aktivna" toccka

// Nevilleov algoritam (Press et al: Numerical Recipes)
double lagrange(double ktocka[][3], int n, double x) {
  double den, dif, dift, ho, hp, w, c[NTMAX + 1], d[NTMAX + 1], y, dy;
  int i, m, ns;

  ns = 1;
  dif = fabs(x - ktocka[0][0]);
  for(i = 1; i <= n; i++) {
    dift = fabs(x - ktocka[i - 1][0]);
    if(dift < dif) { ns = i; dif = dift; }
    c[i] = ktocka[i - 1][1];
    d[i] = ktocka[i - 1][1];
  }
  y = ktocka[ns - 1][1]; ns--;
  for(m = 1; m <= n - 1; m++) {
    for(i = 1; i <= n - m; i++) {
      ho = ktocka[i - 1][0] - x;
      hp = ktocka[i + m - 1][0] - x;
      w = c[i + 1] - d[i];
      den = ho - hp;
      if(den == 0.0) { printf("Division by zero!\n"); exit(0); }
      den = w / den;
      d[i] = hp * den;
      c[i] = ho * den;
    }
    if(2 * ns < n - m) dy = c[ns + 1];
    else { dy = d[ns]; ns--; }
    y += dy;
  }
  return y;
} // lagrange

double udaljenost(double t1[], double t2[]) {
  return sqrt(pow(t1[0] - t2[0], 2.0) +
    pow(t1[1] - t2[1], 2.0) + pow(t1[2] - t2[2], 2.0));
} // udaljenost

void crtajTocku(double t[], double r) {
  glPushMatrix();
    glTranslated(t[0], t[1], t[2]);
    gluDisk(disk, 0.0, r, 16, 2);
  glPopMatrix();
} // crtajTocku

void iscrtaj(void) {
  double x, y;
  int i, j;

  glClear(GL_COLOR_BUFFER_BIT);

  if(ntocaka > 1) {
    // crtkane linije koje spajaju toccke
    glLineStipple(1, 0x00FF);
    glEnable(GL_LINE_STIPPLE);
    glColor3f(1.0, 1.0, 0.0);
    glBegin(GL_LINE_STRIP);
      for(i = 0; i < ntocaka; i++) glVertex3dv(&ktocka[i][0]);
    glEnd();
    glDisable(GL_LINE_STIPPLE);
    
    // Lagrangeov interpolacijski polinom
    glColor3f(0.0, 1.0, 1.0); // cijan
    glBegin(GL_LINE_STRIP);
      for(x = ktocka[0][0]; x < ktocka[ntocaka - 1][0]; x += 0.1) {
         y = lagrange(ktocka, ntocaka, x); // printf("%lf %lf\n", x, y);
         glVertex3d(x, y, 0.0);
      }
      glVertex3d(ktocka[ntocaka - 1][0], ktocka[ntocaka - 1][1], 0.0);
    glEnd();      
  }

  // crtanje toccaka
  for(i = 0, j = 0; i < ntocaka; i++) { 
    if(i == akt) glColor3f(0.0, 0.0, 1.0);
    else glColor3f(1.0, 0.0, 0.0); 
    crtajTocku(&ktocka[i][0], TPOLUMJER);
  }

  glutSwapBuffers();
} // iscrtaj

void skaliraj(int w, int h) {
  float a = 5.0, b;

  glViewport(0, 0, w, h);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  if(w > h) {
    b = a * (float)w / (float)h;
    glOrtho(-b, b, -a, a, -a, a);
  } else {
    b = a * (float)h / (float)w;
    glOrtho(-a, a, -b, b, -a, a);
  }

  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
} // skaliraj

void mish(int tipka, int stanje, int x, int yy) {
  static double zadnjiWX = -100.0;
  GLint viewport[4];
  GLdouble mvmatrix[16], projmatrix[16];
  int y, i, j;
  double wx, wy, wz, w[3], d;

  if(tipka == GLUT_LEFT_BUTTON) {
    if(stanje == GLUT_DOWN) {
      glGetIntegerv(GL_VIEWPORT, viewport);
      glGetDoublev(GL_MODELVIEW_MATRIX, mvmatrix);
      glGetDoublev(GL_PROJECTION_MATRIX, projmatrix);

      // viewport[2] - shirina prozora; viewport[3] - visina prozora
      y = viewport[3] - yy - 1; // korekcija ishodishta
      gluUnProject((GLdouble) x, (GLdouble) y, 0.0, 
        mvmatrix, projmatrix, viewport, &wx, &wy, &wz);
      w[0] = wx; w[1] = wy; w[2] = 0.0;

      akt = -1;
      for(i = 0; i < ntocaka; i++) {
        d = udaljenost(ktocka[i], w); // printf("%d udaljenost: %lf\n", i, d);
        if(d <= TPOLUMJER) akt = i;
      }
      // dodavanje nove toccke
      if((akt == -1) && (wx > zadnjiWX)) {
        ktocka[ntocaka][0] = wx;
        ktocka[ntocaka][1] = wy;
        ktocka[ntocaka][2] = 0.0;
        ntocaka++;
        zadnjiWX = wx;
      }
      glutPostRedisplay();
    }

    if(stanje == GLUT_UP) {
      glGetIntegerv(GL_VIEWPORT, viewport);
      glGetDoublev(GL_MODELVIEW_MATRIX, mvmatrix);
      glGetDoublev(GL_PROJECTION_MATRIX, projmatrix);

      // viewport[2] - shirina prozora; viewport[3] - visina prozora
      y = viewport[3] - yy - 1; // korekcija ishodishta
      gluUnProject ((GLdouble) x, (GLdouble) y, 0.0, 
        mvmatrix, projmatrix, viewport, &wx, &wy, &wz);

      //printf ("(%4d, %4d) -> (%f, %f)\n", x, y, wx, wy);
      if(akt >= 0) {
        ktocka[akt][1] = wy; // dozvoljena je samo promjena y
        akt = -1;
      }
      glutPostRedisplay();
    }
  }

  if((tipka == GLUT_MIDDLE_BUTTON) && (stanje == GLUT_UP)) exit(0);
} // mish

int main(int argc, char** argv) {
  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);
  glutInitWindowSize(500, 500);
  glutInitWindowPosition(0, 0);
  glutCreateWindow(argv[0]);

  glClearColor(0.0, 0.0, 0.0, 0.0);
  disk = gluNewQuadric();

  glutDisplayFunc(iscrtaj);
  glutReshapeFunc(skaliraj);
  glutMouseFunc(mish);

  glutMainLoop();
  return 0;
} // main