#include "warp.h"
#include "draw.h"
#include "ip.h"
#include "geometry.h"
#include <stdlib.h>
#include <vector>


/*
 * very bad style - don't try this at home!
 */
extern double nearestInterpolate  (Image* img, double x, double y, int c);
extern double bilinearInterpolate (Image* img, double x, double y, int c);
extern double gaussianInterpolate (Image* img, double x, double y, int c, 
                                   int n, double sigma);


static inline void checkStream (const istream& in)
{
  if (in.fail())
  {
    cerr << "Fatal error: stream failed!" << endl;
    exit(-1);
  }
} 



Image* warp_mosaic (Image* src)
{
  int size;
  cerr << "Tile size (int): ";
  cin  >> size;

  if (size <= 0)
  {
    cerr << "Tile size must be a positive integer!" << endl;
    return NULL;
  }

  checkStream(cin);

  Image* result = new Image(src->getWidth(),
                            src->getHeight(),
                            src->getChannels(),
                            src->getBits()); 

  int xS = ((src->getWidth()  + size - 1) / size) + 1;
  int yS = ((src->getHeight() + size - 1) / size) + 1;
      
  point** mesh = new point*[xS];
  for (int i = 0; i < xS; i++)
      mesh[i] = new point[yS];
      
  for (int i = 0; i < xS; i++)
    for (int j = 0; j < yS; j++)
    {
      point p;
      p.x = round((i * size) + (rand() * size / RAND_MAX / 2)
                  - (size / 4));
      p.y = round((j * size) + (rand() * size / RAND_MAX / 2)
                  - (size / 4));

      p.x = clamp(p.x, 0, src->getWidth()  - 1);
      p.y = clamp(p.y, 0, src->getHeight() - 1);

      mesh[i][j] = p;
    }
  
  for (int i = 0; i < xS; i++)
  {
    mesh[i][0].y    = 0;
    mesh[i][yS-1].y = src->getHeight() - 1;
  }
  for (int j = 0; j < yS; j++)
  {
    mesh[0][j].x    = 0;
    mesh[xS-1][j].x = src->getWidth() - 1;
  }   

  for (int i = 0; i < (xS-1); i++)
    for (int j = 0; j < (yS-1); j++)
    {
      int x = (i * size) + (size / 2);
      int y = (j * size) + (size / 2);

      x = clamp(x, 0, src->getWidth()  - 1);
      y = clamp(y, 0, src->getHeight() - 1);

      Pixel p;
      src->getPixel_(x, y, p);

      drawQuad(result, mesh[i][j], mesh[i][j+1], 
               mesh[i+1][j+1], mesh[i+1][j], p);
    }
                              
  for (int i = 0; i < xS; i++)
    delete[] mesh[i];
  delete[] mesh;

  return result;
}


Image* warp_dots (Image* src)
{
  cerr << "Frequency (double): ";

  double frequency;
  cin  >> frequency;

  if (frequency <= 0 || frequency >= 1)
  {
    cerr << "Frequency must be between 0 and 1" << endl;
    return NULL;
  }

  checkStream(cin);

  Image* result = new Image(src->getWidth(),
                            src->getHeight(),
                            src->getChannels(),
                            src->getBits());

  int    radius = (int) (1.0 / frequency);
  double done   = frequency * src->getWidth() * src->getHeight();
  int x, y;
      
  for (int i = 0; i < done; i++)
  {
    x = rand() * result->getWidth()  / RAND_MAX;
    y = rand() * result->getHeight() / RAND_MAX;

    Pixel p;
    src->getPixel(x,y,p);

    drawAntiAliasedCircle(result, x, y, radius, p);
  }

  return result;
}


Image* warp_ripple (Image* src)
{
  int mag, period, horiz, iMethod, filtersize, sigma;
  double phase;

  cerr << "Magnitude (int): ";
  cin >> mag;

  if (mag <= 0)
  {
    cerr << "Magnitude must be a positive integer!" << endl;
    return NULL;
  }

  cerr << "Period (int): ";
  cin >> period;

  if (period <= 0)
  {
    cerr << "Period must be a positive integer!" << endl;
    return NULL;
  }

  cerr << "Phase shift (double): ";
  cin >> phase;

  cerr << "0 for horizontal, 1 for vertical: ";
  cin >> horiz;

  if (horiz != 0 && horiz != 1)
  {
    cerr << "Must be either 0 or 1!" << endl;
    return NULL;
  }

  cerr << "Interpolation mode (int) : ";
  cin  >> iMethod;

  if (iMethod != I_NEAREST && iMethod != I_BILINEAR && 
      iMethod != I_GAUSSIAN)
  {
    cerr << "Mode must be one of 0 (linear), 1 (bilinear), or "
         << "2 (gaussian)!" << endl;
    return NULL;
  }

  if (iMethod == I_GAUSSIAN)
  {
    cerr << "Filter size (int) : ";
    cin  >> filtersize;

    if ((filtersize % 2) != 1 || filtersize < 1)
    {
      cerr << "Filter size must be an odd, positive integer!" << endl;
      return NULL;
    }

    cerr << "Sigma (double) : ";
    cin  >> sigma;

    if (sigma == 0)
    {
      cerr << "Sigma must be non-zero!" << endl;
      return NULL;
    }                    
  }

  checkStream(cin);

  Image* result = new Image(src->getWidth(),
                            src->getHeight(),
                            src->getChannels(),
                            src->getBits());

  for (int i = 0; i < result->getWidth(); i++)
    for (int j = 0; j < result->getHeight(); j++)
    {
      double x, y;
      if (horiz == 0)
      {
        x = i + mag * sin(j * M_PI_2 / period + phase);
        y = j;
      }
      else
      {
        x = i;
        y = j + mag * sin(i * M_PI_2 / period + phase);
      }

      for (int c = 0; c < result->getChannels(); c++)
      {
        double v = 0;

        if ((floor(x) == x) && (floor(y) == y))
          v = src->getPixel_(x,y,c);
        else if (iMethod == I_NEAREST)
          v = nearestInterpolate(src, x, y, c);
        else if (iMethod == I_BILINEAR)
          v = bilinearInterpolate(src, x, y, c);
        else if (iMethod == I_GAUSSIAN)
          v = gaussianInterpolate(src, x, y, c, filtersize, sigma);

        result->setPixel(i, j, c, clamp(v,0,1));
      }
    }

  return result;
}


Image* warp_emboss (Image* src)
{
  Image* grey = ip_grey(src);

  Image* result = new Image(src->getWidth(),
                            src->getHeight(),
                            1);

  int hoff = -2, voff = 2;
  double v;
  for (int i = 0; i < grey->getWidth(); ++i)
    for (int j = 0; j < grey->getHeight(); ++j)
    {
      v = grey->getPixel(i, j, 0) - grey->getPixel_(i+hoff, j+voff, 0) + 0.5;
      result->setPixel(i, j, 0, clamp(v,0,1));
    }

  delete grey;

  return result;
}


Image* warp_lines (Image* src)
{
  int maxthickness;

  cerr << "Maximum thickness (int): ";
  cin  >> maxthickness;

  if (maxthickness < 1)
  {
    cerr << "Maximum thickness must be a positive integer!" << endl;
    return NULL;
  }

  checkStream(cin);

  Image* grey = ip_grey(src);

  Image* result = new Image(src->getWidth(),
                            src->getHeight(),
                            1, 8);

  double v, line;
  int mid = (maxthickness - 1) / 2;

  for (int i = 0; i < grey->getWidth(); ++i)
    for (int j = 0; 
         j < grey->getHeight() + maxthickness;
         j += maxthickness + 1)
    {
      /*
       * get average luminance for sample column.  count pixels outside 
       * the image as white, so they don't increase the width of the
       * lines artificially.
       */
      v = 0;
      for (int k = 0; k < maxthickness + 1; ++k)
      {
        if (k + j >= grey->getHeight())
          v += 1.0;
        else
          v += grey->getPixel_(i,k+j,0);
      }
      v /= (maxthickness + 1);

      // determine half the line's thickness
      line  = (1.0 - v) * maxthickness;
      line  = clamp(line, 0, maxthickness);
      line /= 2.0;

      /*
       * draw black for pixels entirely in the line, white for pixels
       * entirely out, and grey for partial pixels
       */
      for (int k = 0; k < maxthickness; ++k)
      {
        v = fabs((double) (mid - k)) + 0.5 - line;
        v = clamp(v, 0, 1);

        result->setPixel_(i, j+k, 0, v);
      }

      // add in the white margin between rows
      result->setPixel_(i, j+maxthickness, 0, 1.0);
    }

  delete grey;

  return result;
}


Image* warp_stamp (Image* src)
{
  int whichmask, count, filtersize;
  double alpha, sigma;

  cerr << "Which mask (0 = SPAM, 1 = Gaussian): ";
  cin  >> whichmask;

  if (whichmask != 0 && whichmask != 1)
  {
    cerr << "Invalid mask selection!" << endl;
    return NULL;
  }

  if (whichmask == 1)
  {
    cerr << "Filter size (int): ";
    cin  >> filtersize;

    if (filtersize < 1 || (filtersize % 2) != 1)
    {
      cerr << "Filter size must be odd positive integer!" << endl;
      return NULL;
    }

    cerr << "Sigma (double): ";
    cin  >> sigma;

    if (sigma == 0)
    {
      cerr << "Sigma must be non-zero!" << endl;
      return NULL;
    }
  }

  cerr << "Number of stamps (int): ";
  cin  >> count;

  if (count < 1)
  {
    cerr << "Number of stamps must be a positive integer!" << endl;
    return NULL;
  }

  cerr << "Mask alpha (double): ";
  cin  >> alpha;

  if (alpha < 0 || alpha > 1)
  {
    cerr << "Alpha must be between 0 and 1!" << endl;
    return NULL;
  }

  checkStream(cin);

  int maskWidth;
  int maskHeight;
  double* mask;
    
  if (whichmask == 0)
  {
    // SPAM mask

    maskWidth  = 19;
    maskHeight = 7;

    double spam[133] = {1,1,1,0,0,1,0,0,0,0,1,0,0,1,0,1,0,0,1,
                        0,1,1,0,0,1,0,0,0,0,1,1,1,1,0,1,0,0,1,
                        0,0,1,1,0,1,0,0,0,0,1,1,1,1,0,1,0,0,1,
                        0,0,1,1,0,1,1,1,0,0,1,0,0,1,0,1,1,1,1,
                        0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,1,1,
                        1,1,0,0,0,1,0,1,1,0,0,1,1,0,0,1,0,0,1,
                        0,1,1,1,0,1,1,1,0,0,0,1,1,0,0,0,0,0,0};

    mask = new double[133];
    for (int i = 0; i < 133; ++i)
      mask[i] = spam[i];
  }
  else
  {
    // Noise weighted by gaussian mask
  
    maskWidth  = filtersize;
    maskHeight = filtersize;
          
    mask = new double[sqr(filtersize)];

    int mid = (filtersize - 1) / 2;
    double dist, randNum;

    for (int i = -mid; i <= mid; ++i)
      for (int j = -mid; j <= mid; ++j)
      {
        randNum = ((double) rand()) / RAND_MAX;
        dist = sqr(i) + sqr(j);

        mask[(i + mid) + (j + mid) * maskWidth] = 
          pow(2.0, -dist / sqr(sigma)) * randNum;
      }
  }

  Image* result = new Image(*src);

  int midW = (maskWidth  - 1) / 2;
  int midH = (maskHeight - 1) / 2;

  for (int n = 0; n < count; ++n)
  {
    int x = rand() * src->getWidth()  / RAND_MAX;
    int y = rand() * src->getHeight() / RAND_MAX;

    for (int j = -midH; j <= midH; ++j)
      for (int i = -midW; i <= midW; ++i)
      {
        double weight = alpha * mask[(midH - j) * maskWidth + (midW + i)];

        int curX = clamp(x + i, 0, src->getWidth()  - 1);
        int curY = clamp(y + j, 0, src->getHeight() - 1);

        for (int c = 0; c < src->getChannels(); ++c)
        {
          double oldPixel = src->getPixel_(x, y, c);
          double curPixel = result->getPixel_(curX, curY, c);
          double value    = oldPixel * weight + curPixel * (1.0 - weight);

          result->setPixel_(curX, curY, c, clamp(value, 0, 1));
        }
      }
  }


  delete[] mask;

  return result;
}


Image* warp_swirl (Image* src)
{
  int x0, y0, iMethod, filtersize;
  double innertheta, constant, linear, quadratic, sigma;

  cerr << "Current rect : 0 0 " << src->getWidth() - 1 << " "
       << src->getHeight() - 1 << endl;
  cerr << "Rotate around (int int) : ";
  cin  >> x0 >> y0;

  if (x0 < 0 || y0 < 0 || 
      x0 >= src->getWidth() ||
      y0 >= src->getHeight())
  {  
    cerr << "Center of rotation must be within the image!" << endl;
    return NULL; 
  }

  cerr << "Inner theta (degrees, double): ";
  cin  >> innertheta;
  innertheta = deg2rad(innertheta);

  cerr << "Const, linear, quad dropoffs (double x 3): ";
  cin  >> constant >> linear >> quadratic;

  cerr << "Interpolation mode (int) : ";
  cin  >> iMethod;
          
  if (iMethod != I_NEAREST && iMethod != I_BILINEAR &&
      iMethod != I_GAUSSIAN)
  {     
    cerr << "Mode must be one of 0 (linear), 1 (bilinear), or "
         << "2 (gaussian)!" << endl;
    return NULL;
  }      
  
  if (iMethod == I_GAUSSIAN)
  {    
    cerr << "Filter size (int) : ";
    cin  >> filtersize;
 
    if ((filtersize % 2) != 1 || filtersize < 1)
    {  
      cerr << "Filter size must be an odd, positive integer!" << endl;
      return NULL;
    }
      
    cerr << "Sigma (double) : ";
    cin  >> sigma;

    if (sigma == 0)
    {
      cerr << "Sigma must be non-zero!" << endl;
      return NULL;
    }
  }    

  checkStream(cin);

  Image* result = new Image(src->getWidth(),
                            src->getHeight(),
                            src->getChannels(),
                            src->getBits());

  double x, y, theta, dist;

  for (int j = 0; j < src->getHeight(); ++j)
    for (int i = 0; i < src->getWidth(); ++i)
    {
      dist  = sqrt((double) (sqr(i - x0) + sqr(j - y0)));
      theta = innertheta / (constant + linear * dist + quadratic * sqr(dist));
      x = (i - x0) * cos(theta) - (j - y0) * sin(theta) + x0;
      y = (i - x0) * sin(theta) + (j - y0) * cos(theta) + y0;

      for (int c = 0; c < src->getChannels(); ++c)
      {
        switch (iMethod)
        {
          case I_NEAREST:
            result->setPixel(i, j, c,
                nearestInterpolate(src, x, y, c));
            break;
     
          case I_BILINEAR:
            result->setPixel(i, j, c,
                bilinearInterpolate(src, x, y, c));
            break;
     
          case I_GAUSSIAN:
            result->setPixel(i, j, c,
                gaussianInterpolate(src, x, y, c, filtersize, sigma));
            break;           
        }
      }
    }

  return result;
}


Image* warp_spherical (Image* src)
{
  int x0, y0, iMethod, filtersize;
  double sigma;

  cerr << "Current rect : 0 0 " << src->getWidth() - 1 << " "
       << src->getHeight() - 1 << endl;
  cerr << "Sphere center (int int) : ";
  cin  >> x0 >> y0;

  if (x0 < 0 || y0 < 0 || 
      x0 >= src->getWidth() ||
      y0 >= src->getHeight())
  {  
    cerr << "Center of sphere must be within the image!" << endl;
    return NULL; 
  }

  cerr << "Interpolation mode (int) : ";
  cin  >> iMethod;
            
  if (iMethod != I_NEAREST && iMethod != I_BILINEAR &&
      iMethod != I_GAUSSIAN)
  {
    cerr << "Mode must be one of 0 (linear), 1 (bilinear), or "
         << "2 (gaussian)!" << endl;
    return NULL;
  }
     
  if (iMethod == I_GAUSSIAN)
  {
    cerr << "Filter size (int) : ";
    cin  >> filtersize;
     
    if ((filtersize % 2) != 1 || filtersize < 1)
    {
      cerr << "Filter size must be an odd, positive integer!" << endl;
      return NULL;
    }
                            
    cerr << "Sigma (double) : ";
    cin  >> sigma;
    
    if (sigma == 0)
    {       
      cerr << "Sigma must be non-zero!" << endl;
      return NULL;
    }
  }

  checkStream(cin);

  Image* result = new Image(src->getWidth(),
                            src->getHeight(),
                            src->getChannels(),
                            src->getBits());

  double xs, ys, x, y, d, r, theta;

  for (int j = 0; j < result->getHeight(); ++j)
    for (int i = 0; i < result->getWidth(); ++i)
    {
      xs = (i - x0) * 2.0 / result->getWidth();
      ys = (j - y0) * 2.0 / result->getHeight();
      d = sqrt(sqr(xs) + sqr(ys));

      if (d > 1)
      {
        for (int c = 0; c < result->getChannels(); ++c)
          result->setPixel(i, j, c, 0.0);
      }
      else
      {
        r = asin(d) / M_PI_2;
        theta = atan2(ys, xs);

        x = r * cos(theta) * result->getWidth()  / 2.0 + x0;
        y = r * sin(theta) * result->getHeight() / 2.0 + y0;

        for (int c = 0; c < result->getChannels(); ++c)
        {
          switch (iMethod)
          {
            case I_NEAREST:   
              result->setPixel(i, j, c,
                  nearestInterpolate(src, x, y, c));
              break;
    
            case I_BILINEAR:
              result->setPixel(i, j, c,
                  bilinearInterpolate(src, x, y, c));
              break;
    
            case I_GAUSSIAN:
              result->setPixel(i, j, c,
                  gaussianInterpolate(src, x, y, c, filtersize, sigma));
              break;
          }
        }
      }
    }

  return result;
}


static const Point4d LUM(0.30, 0.59, 0.11, 1);


/*
 * http://www.sgi.com/grafica/matrix/
 * http://www.sgi.com/software/opengl/advanced98/notes/node187.html
 */

Image* warp_rotate_hue (Image* src)
{
  double theta;

  cerr << "Theta (degrees, double): ";
  cin  >> theta;

  checkStream(cin);

  Image* result = new Image(src->getWidth(),
                            src->getHeight(),
                            src->getChannels(),
                            src->getBits());

  Matrixd Rx = Matrixd::IdentityMatrix();
  Rx(1,1) =  cos(-M_PI / 4.0);
  Rx(1,2) =  sin(-M_PI / 4.0);
  Rx(2,1) = -sin(-M_PI / 4.0);
  Rx(2,2) =  cos(-M_PI / 4.0);

  Matrixd Ry = Matrixd::IdentityMatrix();
  Ry(0,0) =  cos(M_PI / 5.0);
  Ry(0,2) = -sin(M_PI / 5.0);
  Ry(2,0) =  sin(M_PI / 5.0);
  Ry(2,2) =  cos(M_PI / 5.0);

  Matrixd Rxy = Rx * Ry;

  Point4d lum = Rxy * LUM;

  Matrixd S = Matrixd::IdentityMatrix();
//  S(0,2) = -lum[0] / lum[2];
//  S(1,2) = -lum[1] / lum[2];
  S(2,2) = 1e6;

  theta = deg2rad(theta);

  Matrixd Rz = Matrixd::IdentityMatrix();
  Rz(0,0) =  cos(theta);
  Rz(0,1) =  sin(theta);
  Rz(1,0) = -sin(theta);
  Rz(1,1) =  cos(theta);

  Matrixd m = Rxy * S * Rz * S.invert() * Rxy.invert();

  // use this rotation matrix to transform each color in the image
  Pixel p;
  Point4d q;

  for (int j = 0; j < result->getHeight(); ++j)
    for (int i = 0; i < result->getWidth(); ++i)
    {
      src->getPixel(i,j,p);

      q[0] = p.r; 
      q[1] = p.g; 
      q[2] = p.b;
      q[3] = 1;

      q = m * q;

      p.r = clamp(q[0], 0, 1); 
      p.g = clamp(q[1], 0, 1); 
      p.b = clamp(q[2], 0, 1); 

      result->setPixel(i,j,p);
    }

  return result;
}


Image* warp_voronoi (Image* src)
{
  Image* result = new Image(src->getWidth(),
                            src->getHeight(),
                            src->getChannels(),
                            src->getBits());

  glEnable(GL_DEPTH_TEST);
  glDisable(GL_LIGHTING);
  glShadeModel(GL_FLAT);

  glClear(GL_DEPTH_BUFFER_BIT | GL_COLOR_BUFFER_BIT);

  glMatrixMode(GL_PROJECTION);
  glPushMatrix();
    glLoadIdentity();
    glOrtho(0, src->getWidth(), 0, src->getHeight(), -10, 10);

  glMatrixMode(GL_MODELVIEW);
  glPushMatrix();
    glLoadIdentity();

  VECTOR(Vector2i) points;
  points.reserve(100);

  for (unsigned int i = 0; i < points.capacity(); ++i)
  {
    Vector2i p(drand48() * (src->getWidth()  - 1),
               drand48() * (src->getHeight() - 1));
    points.push_back(p);
  }

  Pixel c;

  for (unsigned int i = 0; i < points.size(); ++i)
  {
    src->getPixel_(points[i][0], src->getHeight() - points[i][1] - 1, c);
    glColor3f(c.r, c.g, c.b);

    glPushMatrix();
      glTranslatef(points[i][0], points[i][1], 0);
      glutSolidCone(50, 10, 20, 1);
    glPopMatrix();
  }

  result->glReadPixelsWrapper();

  glDisable(GL_DEPTH_TEST);

  glMatrixMode(GL_PROJECTION);
  glPopMatrix();

  glMatrixMode(GL_MODELVIEW);
  glPopMatrix();

  return result;
}