/* Copyright (C) 2007 Michael Goffioul
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License as published by
** the Free Software Foundation; either version 2 of the License, or
** (at your option) any later version.
**
** This program is distributed in the hope that it will be useful,
** but WITHOUT ANY WARRANTY; without even the implied warranty of
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
** GNU General Public License for more details.
**
** You should have received a copy of the GNU General Public License
** along with this program; if not, write to the Free Software
** Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301
*/

#include "gl-graphics-renderer.h"

#include <GL/gl.h>

enum {
  AXE_ANY_DIR   = 0,
  AXE_DEPTH_DIR = 1,
  AXE_HORZ_DIR  = 2,
  AXE_VERT_DIR  = 3
};

inline Matrix xform_matrix (void)
{
  Matrix m (4, 4, 0.0);
  for (int i = 0; i < 4; i++)
    m(i,i) = 1;
  return m;
}

inline ColumnVector xform_vector (void)
{
  ColumnVector v (4, 0.0);
  v(3) = 1;
  return v;
}

inline ColumnVector xform_vector (double x, double y, double z)
{
  ColumnVector v (4, 1.0);
  v(0) = x; v(1) = y; v(2) = z;
  return v;
}

inline ColumnVector transform (const Matrix& m, double x, double y, double z)
{
  return (m * xform_vector (x, y, z));
}

inline Matrix xform_scale (double x, double y, double z)
{
  Matrix m (4, 4, 0.0);
  m(0,0) = x; m(1,1) = y; m(2,2) = z; m(3,3) = 1;
  return m;
}

inline Matrix xform_translate (double x, double y, double z)
{
  Matrix m = xform_matrix ();
  m(0,3) = x; m(1,3) = y; m(2,3) = z; m(3,3) = 1;
  return m;
}

inline void scale (Matrix& m, double x, double y, double z)
{
  m = m * xform_scale (x, y, z);
}

inline void translate (Matrix& m, double x, double y, double z)
{
  m = m * xform_translate (x, y, z);
}

inline void xform (ColumnVector& v, const Matrix& m)
{
  v = m*v;
}

inline void scale (ColumnVector& v, double x, double y, double z)
{
  v(0) *= x;
  v(1) *= y;
  v(2) *= z;
}

inline void translate (ColumnVector& v, double x, double y, double z)
{
  v(0) += x;
  v(1) += y;
  v(2) += z;
}

inline void normalize (ColumnVector& v)
{
  double fact = 1.0/sqrt(v(0)*v(0)+v(1)*v(1)+v(2)*v(2));
  scale (v, fact, fact, fact);
}

inline double dot (const ColumnVector& v1, const ColumnVector& v2)
{
  return (v1(0)*v2(0)+v1(1)*v2(1)+v1(2)*v2(2));
}

inline double norm (const ColumnVector& v)
{
  return sqrt (dot (v, v));
}

inline ColumnVector cross (const ColumnVector& v1, const ColumnVector& v2)
{
  ColumnVector r = xform_vector ();
  r(0) = v1(1)*v2(2)-v1(2)*v2(1);
  r(1) = v1(2)*v2(0)-v1(0)*v2(2);
  r(2) = v1(0)*v2(1)-v1(1)*v2(0);
  return r;
}

inline Matrix unit_cube (void)
{
  static double data[32] = {
    0,0,0,1,
    1,0,0,1,
    0,1,0,1,
    0,0,1,1,
    1,1,0,1,
    1,0,1,1,
    0,1,1,1,
    1,1,1,1};
  Matrix m (4, 8);
  memcpy (m.fortran_vec (), data, sizeof(double)*32);
  return m;
}

void
gl_graphics_renderer::draw (const figure::properties& props)
{
  const double *color = props.get_color ().rgb ();
  Matrix children = props.get_children ();

  if (! error_state)
    {
      glClearColor (color[0], color[1], color[2], 1);
      glClear (GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    }

  for (int i = 0; i < children.length (); i++)
    graphics_renderer::draw (gh_manager::get_object (children(i)));
}

void
gl_graphics_renderer::draw (const axes::properties& props)
{
  setup_xform (props);

  Matrix xlim = sx.scale (props.get_xlim ().matrix_value ());
  Matrix ylim = sy.scale (props.get_ylim ().matrix_value ());
  Matrix zlim = sz.scale (props.get_zlim ().matrix_value ());
  double xmin = xlim(0), xmax = xlim(1);
  double ymin = ylim(0), ymax = ylim(1);
  double zmin = zlim(0), zmax = zlim(1);
  
  double xd = (props.get_xdir ().string_value () == "normal" ? 1 : -1);
  double yd = (props.get_ydir ().string_value () == "normal" ? 1 : -1);
  double zd = (props.get_zdir ().string_value () == "normal" ? 1 : -1);

  ColumnVector p1, p2, xv (3), yv (3), zv (3);
  int xstate, ystate, zstate;

  p1 = transform (x_render, xmin, (ymin+ymax)/2, (zmin+zmax)/2);
  p2 = transform (x_render, xmax, (ymin+ymax)/2, (zmin+zmax)/2);
  xv(0) = xround (p2(0)-p1(0)); xv(1) = xround (p2(1)-p1(1)); xv(2) = (p2(2)-p1(2));
  if (xv(0) == 0 && xv(1) == 0)
    xstate = AXE_DEPTH_DIR;
  else if (xv(2) == 0)
    {
      if (xv(0) == 0)
        xstate = AXE_VERT_DIR;
      else if (xv(1) == 0)
        xstate = AXE_HORZ_DIR;
    }
  double xPlane;
  if (xv(2) == 0)
    if (xv(1) == 0)
      xPlane = (xv(0) > 0 ? xmax : xmin);
    else
      xPlane = (xv(1) < 0 ? xmax : xmin);
  else
    xPlane = (xv(2) < 0 ? xmin : xmax);
  double xPlaneN = (xPlane == xmin ? xmax : xmin);
  double fx = (xmax-xmin)/sqrt(xv(0)*xv(0)+xv(1)*xv(1));

  p1 = transform (x_render, (xmin+xmax)/2, ymin, (zmin+zmax)/2);
  p2 = transform (x_render, (xmin+xmax)/2, ymax, (zmin+zmax)/2);
  yv(0) = xround (p2(0)-p1(0)); yv(1) = xround (p2(1)-p1(1)); yv(2) = (p2(2)-p1(2));
  if (yv(0) == 0 && yv(1) == 0)
    ystate = AXE_DEPTH_DIR;
  else if (yv(2) == 0)
    {
      if (yv(0) == 0)
        ystate = AXE_VERT_DIR;
      else if (yv(1) == 0)
        ystate = AXE_HORZ_DIR;
    }
  double yPlane;
  if (yv(2) == 0)
    if (yv(1) == 0)
      yPlane = (yv(0) > 0 ? ymax : ymin);
    else
      yPlane = (yv(1) < 0 ? ymax : ymin);
  else
    yPlane = (yv(2) < 0 ? ymin : ymax);
  double yPlaneN = (yPlane == ymin ? ymax : ymin);
  double fy = (ymax-ymin)/sqrt(yv(0)*yv(0)+yv(1)*yv(1));

  p1 = transform(x_render, (xmin+xmax)/2, (ymin+ymax)/2, zmin);
  p2 = transform(x_render, (xmin+xmax)/2, (ymin+ymax)/2, zmax);
  zv(0) = xround(p2(0)-p1(0)); zv(1) = xround (p2(1)-p1(1)); zv(2) = (p2(2)-p1(2));
  if (zv(0) == 0 && zv(1) == 0)
    zstate = AXE_DEPTH_DIR;
  else if (zv(2) == 0)
  {
    if (zv(0) == 0)
      zstate = AXE_VERT_DIR;
    else if (zv(1) == 0)
      zstate = AXE_HORZ_DIR;
  }
  double zPlane;
  if (zv(2) == 0)
    if (zv(1) == 0)
      zPlane = (zv(0) > 0 ? zmin : zmax);
    else
      zPlane = (zv(1) < 0 ? zmin : zmax);
  else
    zPlane = (zv(2) < 0 ? zmin : zmax);
  double zPlaneN = (zPlane == zmin ? zmax : zmin);
  double fz = (zmax-zmin)/sqrt(zv(0)*zv(0)+zv(1)*zv(1));

  bool mode2d = (((xstate > AXE_DEPTH_DIR ? 1 : 0) +
        (ystate > AXE_DEPTH_DIR ? 1 : 0) +
        (zstate > AXE_DEPTH_DIR ? 1 : 0)) == 2);
  // TODO: tickdirmode
  if (true /*TickDirMode.is("auto")*/)
  {
    //autoMode++;
    //TickDir.set(mode2d ? "in" : "out", true);
    //autoMode--;
  }

  double xticklen = 7, yticklen = 7, zticklen = 7;

  // TODO: use tickdir property
  double tickdir = (mode2d ? -1 : 1) /*(TickDir.is("in") ? -1 : 1)*/;
  double xtickoffset = (mode2d && tickdir < 0 ? 0 : xticklen) + 5;
  double ytickoffset = (mode2d && tickdir < 0 ? 0 : yticklen) + 5;
  double ztickoffset = (mode2d && tickdir < 0 ? 0 : zticklen) + 5;

  bool xySym = (xd*yd*(xPlane-xPlaneN)*(yPlane-yPlaneN) > 0);
  bool x2Dtop = false;
  bool y2Dright = false;
  double zpTick = zPlane;

  /* 2D mode */
  if (xstate == AXE_HORZ_DIR && ystate == AXE_VERT_DIR)
  {
    if (props.get_xaxislocation ().string_value () == "top")
    {
      double tmp = yPlane;
      yPlane = yPlaneN;
      yPlaneN = tmp;
      x2Dtop = true;
    }
    if (props.get_yaxislocation ().string_value () == "right")
    {
      double tmp = xPlane;
      xPlane = xPlaneN;
      xPlaneN = tmp;
      y2Dright = true;
    }
    if (props.get_layer ().current_value () == "top")
      zpTick = zPlaneN;
  }

  color_property axe_color = props.get_color ();
  bool visible = (props.get_visible ().string_value () == "on");
  bool box = (props.get_box ().string_value () == "on");

  // Axes planes

  if (! axe_color.is_radio () && visible)
    {
      set_color (axe_color);
      set_polygon_offset (true, 2.5);

      glBegin (GL_QUADS);

      // X plane
      glVertex3d (xPlane, ymin, zmin);
      glVertex3d (xPlane, ymax, zmin);
      glVertex3d (xPlane, ymax, zmax);
      glVertex3d (xPlane, ymin, zmax);

      // Y plane
      glVertex3d (xmin, yPlane, zmin);
      glVertex3d (xmax, yPlane, zmin);
      glVertex3d (xmax, yPlane, zmax);
      glVertex3d (xmin, yPlane, zmax);

      // Z plane
      glVertex3d (xmin, ymin, zPlane);
      glVertex3d (xmax, ymin, zPlane);
      glVertex3d (xmax, ymax, zPlane);
      glVertex3d (xmin, ymax, zPlane);

      glEnd ();

      set_polygon_offset (false);
    }

  // Axes box

  set_linestyle ("-", true);
  // TODO: use linewidth property
  set_linewidth (0.5f);

  if (visible)
    {
      glBegin (GL_LINES);

      // X box
      set_color (props.get_xcolor ());
      glVertex3d (xPlaneN, yPlaneN, zPlane);
      glVertex3d (xPlane, yPlaneN, zPlane);
      if (box)
        {
          glVertex3d (xPlaneN, yPlane, zPlane);
          glVertex3d (xPlane, yPlane, zPlane);
          glVertex3d (xPlaneN, yPlane, zPlaneN);
          glVertex3d (xPlane, yPlane, zPlaneN);
          glVertex3d (xPlaneN, yPlaneN, zPlaneN);
          glVertex3d (xPlane, yPlaneN, zPlaneN);
        }

      // Y box
      set_color (props.get_ycolor ());
      glVertex3d (xPlaneN, yPlaneN, zPlane);
      glVertex3d (xPlaneN, yPlane, zPlane);
      if (box)
        {
          glVertex3d (xPlane, yPlaneN, zPlane);
          glVertex3d (xPlane, yPlane, zPlane);
          glVertex3d (xPlane, yPlaneN, zPlaneN);
          glVertex3d (xPlane, yPlane, zPlaneN);
          glVertex3d (xPlaneN, yPlaneN, zPlaneN);
          glVertex3d (xPlaneN, yPlane, zPlaneN);
        }

      // Z box
      set_color (props.get_zcolor ());
      if (xySym)
        {
          glVertex3d (xPlaneN, yPlane, zPlaneN);
          glVertex3d (xPlaneN, yPlane, zPlane);
        }
      else
        {
          glVertex3d (xPlane, yPlaneN, zPlaneN);
          glVertex3d (xPlane, yPlaneN, zPlane);
        }
      if (box)
        {
          glVertex3d (xPlane, yPlane, zPlaneN);
          glVertex3d (xPlane, yPlane, zPlane);
          if (xySym)
            {
              glVertex3d (xPlane, yPlaneN, zPlaneN);
              glVertex3d (xPlane, yPlaneN, zPlane);
            }
          else
            {
              glVertex3d (xPlaneN, yPlane, zPlaneN);
              glVertex3d (xPlaneN, yPlane, zPlane);
            }
          glVertex3d (xPlaneN, yPlaneN, zPlaneN);
          glVertex3d (xPlaneN, yPlaneN, zPlane);
        }

      glEnd ();
    }
}

void
gl_graphics_renderer::set_viewport (int w, int h)
{
  graphics_renderer::set_viewport (w, h);
  glViewport (0, 0, w, h);
}

void
gl_graphics_renderer::setup_xform (const axes::properties& props)
{
  sx = props.get_xscale ();
  sy = props.get_yscale ();
  sz = props.get_zscale ();

  double xd = (props.get_xdir ().string_value () == "normal" ? 1 : -1);
  double yd = (props.get_ydir ().string_value () == "normal" ? 1 : -1);
  double zd = (props.get_zdir ().string_value () == "normal" ? 1 : -1);

  Matrix xlim = sx.scale (props.get_xlim ().matrix_value ());
  Matrix ylim = sy.scale (props.get_ylim ().matrix_value ());
  Matrix zlim = sz.scale (props.get_zlim ().matrix_value ());

  double xo = xlim(xd > 0 ? 0 : 1);
  double yo = ylim(yd > 0 ? 0 : 1);
  double zo = zlim(zd > 0 ? 0 : 1);

  Matrix pb (1, 3, 1.0);
  // TODO: should come from axes properties: plotboxaspectratio
  pb(0) = pb(1) = pb(2) = 1;

  // TODO: should come from axes properties: camerapositionmode,
  // cameratargetmode, cameraupvectormode, cameraviewanglemode
  bool autocam = true;
  // TODO: should come from axes properties: plotboxaspectratiomode
  bool dowarp = (autocam && (props.get_dataaspectratiomode ().string_value () == "auto"));

  ColumnVector c_eye (4, 0.0), c_center (4, 0.0), c_upv (4, 0.0);
  c_eye(3) = c_center(3) = c_upv(3) = 1;

  // TODO: cameratargetmode
  if (true)
    {
      c_center(0) = (xlim(0)+xlim(1))/2;
      c_center(1) = (ylim(0)+ylim(1))/2;
      c_center(2) = (zlim(0)+zlim(1))/2;
      // TODO: set cameratarget
    }
  else
    {
      // TODO: cameratarget
    }

  // TODO: camerapositionmode
  if (true)
    {
      Matrix view = props.get_view ().matrix_value ();
      double az = view(0), el = view(1);
      double d = 5*sqrt(pb(0)*pb(0)+pb(1)*pb(1)+pb(2)*pb(2));

      if (el == 90 || el == -90)
        c_eye(2) = d*signum(el);
      else
        {
          az *= M_PI/180.0;
          el *= M_PI/180.0;
          c_eye(0) = d*cos(el)*sin(az);
          c_eye(1) = -d*cos(el)*cos(az);
          c_eye(2) = d*sin(el);
        }
      c_eye(0) = c_eye(0)*(xlim(1)-xlim(0))/(xd*pb(0))+c_center(0);
      c_eye(1) = c_eye(1)*(ylim(1)-ylim(0))/(yd*pb(1))+c_center(1);
      c_eye(2) = c_eye(2)*(zlim(1)-zlim(0))/(zd*pb(2))+c_center(2);
      // TODO: set cameraposition
    }
  else
    {
      // TODO: cameraposition
    }

  // TODO: cameraupvectormode
  if (true)
    {
      Matrix view = props.get_view ().matrix_value ();
      double az = view(0), el = view(1);

      if (el == 90 || el == -90)
        {
          c_upv(0) = -sin(az*M_PI/180.0)*(xlim(1)-xlim(0))/pb(0);
          c_upv(1) = cos(az*M_PI/180.0)*(ylim(1)-ylim(0))/pb(1);
        }
      else
        c_upv(2) = 1;
      // TODO: set cameraupvector
    }
  else
    {
      // TODO: cameraupvector
    }

  // TODO: these matrices should be stored in the axes properties
  x_view = xform_matrix ();
  x_projection = xform_matrix ();
  x_viewport = xform_matrix ();
  x_normrender = xform_matrix ();
  x_render = xform_matrix ();

  Matrix x_pre = xform_matrix ();
  Matrix x_mat1, x_mat2;

  scale (x_pre, pb(0), pb(1), pb(2));
  translate (x_pre, -0.5, -0.5, -0.5);
  scale (x_pre, xd/(xlim(1)-xlim(0)), yd/(ylim(1)-ylim(0)), zd/(zlim(1)-zlim(0)));
  translate (x_pre, -xo, -yo, -zo);

  xform (c_eye, x_pre);
  xform (c_center, x_pre);
  scale (c_upv, pb(0)/(xlim(1)-xlim(0)), pb(1)/(ylim(1)-ylim(0)), pb(2)/(zlim(1)-zlim(0)));
  translate (c_center, -c_eye(0), -c_eye(1), -c_eye(2));

  ColumnVector F (c_center), f (F), UP (c_upv);
  normalize (f);
  normalize (UP);

  if (abs (dot (f, UP)) > 1e-15)
    {
      double fa = 1/sqrt(1-f(2)*f(2));
      scale (UP, fa, fa, fa);
    }

  ColumnVector s = cross (f, UP);
  ColumnVector u = cross (s, f);

  scale (x_view, 1, 1, -1);
  Matrix l = xform_matrix ();
  l(0,0) = s(0); l(0,1) = s(1); l(0,2) = s(2);
  l(1,0) = u(0); l(1,1) = u(1); l(1,2) = u(2);
  l(2,0) = -f(0); l(2,1) = -f(1); l(2,2) = -f(2);
  x_view = x_view * l;
  translate (x_view, -c_eye(0), -c_eye(1), -c_eye(2));
  scale (x_view, pb(0), pb(1), pb(2));
  translate (x_view, -0.5, -0.5, -0.5);

  Matrix x_cube = x_view * unit_cube ();
  ColumnVector cmin = x_cube.row_min (), cmax = x_cube.row_max ();
  double xM = cmax(0)-cmin(0);
  double yM = cmax(1)-cmin(1);

  Matrix bb = get_boundingbox (props);
  double v_angle;

  // TODO: cameraviewanglemode
  if (true)
    {
      double af;

      if (dowarp)
        af = 1.0 / (xM > yM ? xM : yM);
      else
        {
          if ((bb(2)/bb(3)) > (xM/yM))
            af = 1.0 / yM;
          else
            af = 1.0 / xM;
        }
      v_angle = 2 * (180.0 / M_PI) * atan (1 / (2 * af * norm (F)));
      // TODO: set cameraviewangle
    }

  double pf = 1 / (2 * tan ((v_angle / 2) * M_PI / 180.0) * norm (F));
  scale (x_projection, pf, pf, 1);

  if (dowarp)
    {
      xM *= pf;
      yM *= pf;
      translate (x_viewport, bb(0)+bb(2)/2, get_height()-(bb(1)+bb(3)/2)+1, 0);
      scale (x_viewport, bb(2)/xM, -bb(3)/yM, 1);
    }
  else
    {
      double pix = 1;
      if (autocam)
        {
          if ((bb(2)/bb(3)) > (xM/yM))
            pix = bb(3);
          else
            pix = bb(2);
        }
      else
        pix = (bb(2) < bb(3) ? bb(2) : bb(3));
      translate (x_viewport, bb(0)+bb(2)/2, get_height()-(bb(1)+bb(3)/2)+1, 0);
      scale (x_viewport, pix, -pix, 1);
    }

  x_normrender = x_viewport * x_projection * x_view;

  x_cube = x_normrender * unit_cube ();
  cmin = x_cube.row_min ();
  cmax = x_cube.row_max ();
  x_zmin = cmin(2);
  x_zmax = cmax(2);

  x_render = x_normrender;
  scale (x_render, xd/(xlim(1)-xlim(0)), yd/(ylim(1)-ylim(0)), zd/(zlim(1)-zlim(0)));
  translate (x_render, -xo, -yo, -zo);
  std::cout << "x_render" << std::endl << x_render << std::endl;

  x_mat1 = x_view;
  scale (x_mat1, xd/(xlim(1)-xlim(0)), yd/(ylim(1)-ylim(0)), zd/(zlim(1)-zlim(0)));
  translate (x_mat1, -xo, -yo, -zo);
  x_mat2 = x_viewport * x_projection;

  // setup OpenGL transformation

  double xZ1, xZ2;

  xZ1 = x_zmin-(x_zmax-x_zmin)/2;
  xZ2 = x_zmax+(x_zmax-x_zmin)/2;

  glMatrixMode (GL_MODELVIEW);
  glLoadIdentity ();
  glScaled(1, 1, -1);
  glMultMatrixd (x_mat1.data ());
  glMatrixMode (GL_PROJECTION);
  glLoadIdentity ();
  glOrtho (0, get_width (), get_height (), 0, xZ1, xZ2);
  glMultMatrixd (x_mat2.data ());
  glMatrixMode (GL_MODELVIEW);

  glClear (GL_DEPTH_BUFFER_BIT);
}

void
gl_graphics_renderer::set_color (const color_property& p)
{
  const double *rgb = p.rgb ();
  if (! error_state)
    glColor3dv (rgb);
}

void
gl_graphics_renderer::set_polygon_offset (bool enable, double offset)
{
  if (enable)
    {
      glPolygonOffset (offset, offset);
      glEnable (GL_POLYGON_OFFSET_FILL);
    }
  else
    glDisable (GL_POLYGON_OFFSET_FILL);
}

void
gl_graphics_renderer::set_linewidth (float w)
{
  glLineWidth (w);
}

void
gl_graphics_renderer::set_linestyle (const std::string& s, bool use_stipple)
{
  bool solid = false;

  if (s == "-")
    {
      glLineStipple (1, (unsigned short)0xFFFF);
      solid = true;
    }
  else if (s == ":")
    glLineStipple (1, (unsigned short)0x8888);
  else if (s == "--")
    glLineStipple (1, (unsigned short)0x0FFF);
  else if (s == "-.")
    glLineStipple (1, (unsigned short)0x020F);
  else
    glLineStipple (1, (unsigned short)0x0000);

  if (solid && ! use_stipple)
    glDisable (GL_LINE_STIPPLE);
  else
    glEnable (GL_LINE_STIPPLE);
}