#include "ip.h"
#include "main.h"
#include "filter.h"
#include "warp.h"
#include <algorithm>
#include <stdlib.h>

using std::sort;


double nearestInterpolate (Image* img, double x, double y, int c)
{
  int x0 = (int) x;
  int x1 = (int) x + 1;
  int y0 = (int) y;
  int y1 = (int) y + 1;

  return img->getPixel_((x - x0) > 0.5 ? x1 : x0,
                        (y - y0) > 0.5 ? y1 : y0,
                         c);
}


double bilinearInterpolate (Image* img, double x, double y, int c)
{
  if ((x - (int) x) == 0.5 && (y - (int) y) == 0.5)
    return img->getPixel((int) x, (int) y, c);

  int x0 = (x < 0.5)  ? (int) (x - 1.5) : (int) (x - 0.5);
  int x1 = (x < -0.5) ? (int) (x - 0.5) : (int) (x + 0.5);
  int y0 = (y < 0.5)  ? (int) (y - 1.5) : (int) (y - 0.5);
  int y1 = (y < -0.5) ? (int) (y - 0.5) : (int) (y + 0.5);

  x -= 0.5;
  y -= 0.5;

  double v0 = (x - x0) * img->getPixel_(x1,y0,c)
            + (x1 - x) * img->getPixel_(x0,y0,c);
  double v1 = (x - x0) * img->getPixel_(x1,y1,c)
            + (x1 - x) * img->getPixel_(x0,y1,c);

  return ((y - y0) * v1 + (y1 - y) * v0);
}


double gaussianInterpolate (Image* img, double x, double y, int c, int n, 
                            double sigma)
{

if ((x - (int) x == 0) && (y - (int) y == 0))
    return img->getPixel((int) x, (int) y, c);

double weights=0;
double newColor=0;

for (int i = (int) x - n/2; i <= (int) x + 1 + n/2; i++)
    {
      for (int j= (int) y - n/2; j<= (int) y + 1 + n/2; j++) 
	{
	    if (i >=0 && i< img->getWidth() && j>=0 && j< img->getHeight())
		{
		double exponent = -((x-i)*(x-i) + (y-j)*(y-j))/(sigma*sigma);
	     	double factor = pow(M_E, exponent);
	 	weights += factor;	
		newColor += img->getPixel(i,j,c)*factor;
		}
	}
     }
newColor /= weights;
if (newColor <0 || newColor>1) {
	cout << newColor << endl;
	newColor=1;
}
return newColor;


/*
  if ((x - (int) x) == 0.5 && (y - (int) y) == 0.5)
    return img->getPixel((int) x, (int) y, c);

  int k = (n - 1) / 2;
  int x0 = (x < 0.5) ? (int) (x - 1.5) : (int) (x - 0.5);
  int y0 = (y < 0.5) ? (int) (y - 1.5) : (int) (y - 0.5);

  x -= 0.5;
  y -= 0.5;

  double dx = x - x0;
  double dy = y - y0;
  double v  = 0;
  double f  = 0;
  double s  = 0;

  for (int i = -k; i <= k; ++i)
    for (int j = -k; j <= k; ++j)
    {
      f  = pow(M_E, fabs(sqr(dx-i) + sqr(dy-j)) / sqr(sigma) * -1.0);
      s += f;
      v += img->getPixel_(i+x0, j+y0, c) * f;
    }

  return v / s;
*/
}



Image* interpolate (Image* a, Image* b, double alpha)
{
  assert((a->getWidth()    == b->getWidth())    && 
         (a->getHeight()   == b->getHeight())   &&
         (a->getChannels() == b->getChannels()));

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

  for (int j = 0; j < result->getHeight(); ++j)
    for (int i = 0; i < result->getWidth(); ++i)
      for (int c = 0; c < result->getChannels(); ++c)
        result->setPixel(i,j,c,
                         clamp(((1-alpha) * a->getPixel(i,j,c) +
                                alpha * b->getPixel(i,j,c)
                               ),0.0,1.0)
                        );

  return result;
}


void applyFilter (Image* o, Image* n, Filter* f, int x, int y, int c)
{
  int k   = (f->size() - 1) / 2;
  int lft = x - k;
  int rht = x + k;
  int bot = y - k;
  int top = y + k;

  double v = 0;
  for (int j = bot; j <= top; ++j)
    for (int i = lft; i <= rht; ++i)
      v += (o->getPixel_(i,j,c) * f->get(i-x+k, j-y+k));

  n->setPixel(x,y,c,clamp(v,0.0,1.0));
}


Image* convolve (Image* old, Filter* f)
{
  Image* result = new Image(old->getWidth(),
                            old->getHeight(),
                            old->getChannels(),
                            old->getBits());

  for (int j = 0; j < result->getHeight(); ++j)
    for (int i = 0; i < result->getWidth(); ++i)
      for (int c = 0; c < result->getChannels(); ++c)
        applyFilter(old, result, f, i, j, c);

  return result;
}


/*
 * convolve with a box filter
 */
Image* ip_blur_box (Image* src, int size)
{
  if (size % 2 != 1 || size < 1)
  {
    cerr << "Filter size must be an odd, positive integer!" << endl;
    return NULL;
  }

  Filter* f = new Filter(size);

  for (int j = 0; j < size; ++j)
    for (int i = 0; i < size; ++i)
      f->get(i,j) = 1;
  f->normalize();

  Image* result = convolve(src, f);

  delete f;
  return result;
}


/*
 * convolve with a gaussian filter
 */
Image* ip_blur_gaussian (Image* src, int size, double sigma)
{
  if (size % 2 != 1 || size < 1)
  {
    cerr << "Filter size must be an odd, positive integer!" << endl;
    return NULL;
  }
  if (sigma == 0)
  {
    cerr << "Sigma must be non-zero!" << endl;
    return NULL;
  }

  Filter* f = new Filter(size);

  int k = (size - 1) / 2;
  double d;

  for (int j = 0; j < size; ++j)
    for (int i = 0; i < size; ++i)
    {
      d = sqr(i-k) + sqr(j-k);
      f->get(i,j) = pow(2.0, -d/sqr(sigma));
    }

  f->normalize();

  Image* result = convolve(src, f);

  delete f;
  return result;
}


/*
 * convolve with a triangle filter
 */
Image* ip_blur_triangle (Image* src, int size)
{
  if (size % 2 != 1 || size < 1)
  {
    cerr << "Filter size must be an odd, positive integer!" << endl;
    return NULL;
  }

  Filter* f = new Filter(size);
  int k = (size - 1) / 2;
  int v = 1;

  int* samples = new int[size];
  for (int i = 0; i < size; ++i)
  {
    samples[i] = v;

    if (i < k)
      v++;
    else
      v--;
  }

  for (int j = 0; j < size; ++j)
    for (int i = 0; i < size; ++i)
      f->get(i,j) = samples[i]*samples[j];

  f->normalize();

  Image* result = convolve(src, f);

  delete[] samples;
  delete f;
  return result;
}


/*
 * interpolate with a black image
 */
Image* ip_brighten (Image* src, double alpha)
{
  Image* black = new Image(src->getWidth(),
                           src->getHeight(),
                           src->getChannels());

  Pixel blackpixel = { 0.0, 0.0, 0.0};

  for (int j = 0; j < black->getHeight(); ++j)
    for (int i = 0; i < black->getWidth(); ++i)
      black->setPixel(i,j,blackpixel);

  Image* result = interpolate(black, src, alpha);

  delete black;
  return result;
}


/*
 * interpolate with the average intensity of the src image
 */
Image* ip_contrast (Image* src, double alpha)
{
  Image* grey = new Image(src->getWidth(),
                          src->getHeight(),
                          src->getChannels());

  double avg[4] = { 0.0, 0.0, 0.0, 0.0 };

  for (int j = 0; j < src->getHeight(); ++j)
    for (int i = 0; i < src->getWidth(); ++i)
      for (int c = 0; c < src->getChannels(); ++c)
        avg[c] += src->getPixel(i,j,c);

  for (int c = 0; c < src->getChannels(); ++c)
    avg[c] /= (src->getWidth() * src->getHeight());

  Pixel greypixel;
  switch (grey->getChannels())
  {

    case 3:
      greypixel.r = avg[0] * 0.3 + avg[1] * 0.59 + avg[2] * 0.11;
      greypixel.g = avg[0] * 0.3 + avg[1] * 0.59 + avg[2] * 0.11;
      greypixel.b = avg[0] * 0.3 + avg[1] * 0.59 + avg[2] * 0.11;
      break;

    case 1:
      greypixel.r = avg[0];
      break;

    default:
      break;
  }

  for (int j = 0; j < grey->getHeight(); ++j)
    for (int i = 0; i < grey->getWidth(); ++i)
      grey->setPixel(i,j,greypixel);

  Image* result = interpolate(grey, src, alpha);

  delete grey;
  return result;
}


/*
 * use a mask image for a per-pixel alpha value to perform
 * interpolation with a second image
 */
Image* ip_composite (Image* src, const char* imageName, 
                     const char* maskName)
{
  Image* secondImage = new Image(imageName);
  if (secondImage->bad())
  {
    cerr << "Couldn't read " << imageName << "!" << endl;
    delete secondImage;
    return NULL;
  }

  if (secondImage->getWidth()    != src->getWidth() ||
      secondImage->getHeight()   != src->getHeight() ||
      secondImage->getChannels() != src->getChannels())
  {
    cerr << "Images do not match in size or number of channels!" << endl;
    delete secondImage;
    return NULL;
  }

  Image* maskImage   = new Image(maskName);
  if (maskImage->bad())
  {
    cerr << "Couldn't read " << maskName << "!" << endl;
    delete secondImage;
    delete maskImage;
    return NULL;
  }

  if (maskImage->getWidth()  != src->getWidth() ||
      maskImage->getHeight() != src->getHeight())
  {
    cerr << "Mask image is the wrong size!" << endl;
    delete secondImage;
    delete maskImage;
    return NULL;
  }

  if (maskImage->getChannels() != 1)
  {
    cerr << "Warning: only using red channel of mask image" << endl;
  }

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

  double alpha, v;
  for (int j = 0; j < src->getHeight(); ++j)
    for (int i = 0; i < src->getWidth(); ++i)
    {
      alpha = maskImage->getPixel(i, j, 0);

      for (int c = 0; c < src->getChannels(); ++c)
      {
        v = alpha * src->getPixel(i, j, c)
          + (1 - alpha) * secondImage->getPixel(i, j, c);
        result->setPixel(i, j, c, v);
      }
    }

  delete secondImage;
  delete maskImage;
  return result;
}


/*
 * cut away all but a subset of the image
 */
Image* ip_crop (Image* src, int x0, int y0, int x1, int y1)
{
  if (x0 < 0 || y0 < 0 || x1 < 0 || y1 < 0 ||
      x0 >= src->getWidth()  ||
      y0 >= src->getHeight() ||
      x1 >= src->getWidth()  ||
      y1 >= src->getHeight())
  {
    cerr << "The upper-left and lower-right corners must be within the"
         << " image bounds!" << endl;
    return NULL;
  }
  if (x0 >= x1 || y0 >= y1)
  {
    cerr << "The first point must be above and to the left of the"
         << " second point!" << endl;
    return NULL;
  }

  if ((x1 - x0 + 1) == src->getWidth() &&
      (y1 - y0 + 1) == src->getHeight())
    return NULL;

  Image* subImage = new Image(x1 - x0 + 1, 
                              y1 - y0 + 1,
                              src->getChannels(),
                              src->getBits());

  Pixel pixel;
  for (int j = y0; j <= y1; ++j)
    for (int i = x0; i <= x1; ++i)
      subImage->setPixel(i-x0, j-y0, src->getPixel(i, j, pixel));

  return subImage;
}


/*
 * convolve with an edge detection kernel
 */
Image* ip_edge_detect (Image* src)
{
  Filter* f = new Filter(3);

  f->get(0,0) = -1;
  f->get(0,1) = -1;
  f->get(0,2) = -1;
  f->get(1,0) = -1;
  f->get(1,1) =  8;
  f->get(1,2) = -1;
  f->get(2,0) = -1;
  f->get(2,1) = -1;
  f->get(2,2) = -1;

  Image* result = convolve(src, f);

  delete f;
  return result;
}


/*
 * extract a channel from the image
 */
Image* ip_extract (Image* src, int channel)
{
  if (channel < 0 || channel >= src->getChannels())
  {
    cerr << "Invalid channel number!" << endl;
    return NULL;
  }
  if (src->getChannels()==1)
  {
	  cerr << "This is a one channel image!" << endl;
	  return NULL;
  }
  Image* monochrome = new Image(src->getWidth(),
                                src->getHeight(),
                                3,
                                src->getBits());


  Pixel newval;
  int factor[3]={0,0,0};
  factor[channel]=1;
  


  for (int j = 0; j < src->getHeight(); ++j)
	  for (int i = 0; i < src->getWidth(); ++i) {
		  newval = src->getPixel(i,j);
		  newval.r *= factor[0];
		  newval.g *= factor[1];
		  newval.b *= factor[2];
		  monochrome->setPixel(i,j,newval);
	  }

  return monochrome;
}


/*
 * create a new one channel image with the psychosomatic intensities
 * of the source image
 */
Image* ip_grey (Image* src)
{
  if (src->getChannels() == 1)
    return NULL;

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

  double v;
  for (int j = 0; j < src->getHeight(); ++j)
    for (int i = 0; i < src->getWidth(); ++i)
    {
      v  = src->getPixel(i,j,RED)   * 0.2126;
      v += src->getPixel(i,j,GREEN) * 0.7152;
      v += src->getPixel(i,j,BLUE)  * 0.0722;

      greyImage->setPixel(i,j,0,v);
    }

  return greyImage;
}


/*
 * subtract the image from a white image
 */
Image* ip_invert (Image* src)
{
  Image* white = new Image(src->getWidth(),
                           src->getHeight(),
                           src->getChannels());

  Pixel whitepixel = { 0.5, 0.5, 0.5 };

  for (int j = 0; j < white->getHeight(); ++j)
    for (int i = 0; i < white->getWidth(); ++i)
      white->setPixel(i,j,whitepixel);

  Image* result = interpolate(white, src, -1.0);

  delete white;
  return result;
}


void getVals (Image* img, int x, int y, int size, double* vals)
{
  int k = (size - 1) / 2;

  for (int j = 0; j < size; ++j)
    for (int i = 0; i < size; ++i)
      for (int c = 0; c < img->getChannels(); ++c)
        vals[sqr(size)*c + size*j + i] = img->getPixel_(x-k+i, y-k+j, c);
}


/*
 * replace each pixel with the median of the values surrounding it
 */
Image* ip_median (Image* src, int size)
{
  if (size % 2 != 1 || size < 1)
  {
    cerr << "Filter size must be an odd, positive integer!" << endl;
    return NULL;
  }

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

  double* vals = new double[sqr(size) * src->getChannels()];
  int     k    = (sqr(size) - 1) / 2;

  for (int j = 0; j < src->getHeight(); ++j)
    for (int i = 0; i < src->getWidth(); ++i)
    {
      getVals(src, i, j, size, vals);

      for (int c = 0; c < src->getChannels(); ++c)
      {
        //sort(vals +  c*sqr(size), vals + (c+1)*sqr(size));
		  cerr << "Warning: median filter is not correct." << endl;
        result->setPixel(i, j, c, vals[k + c*sqr(size)]);
      }
    }

  delete[] vals;
  return result;
}


/*
 * interpolate with random noise
 */
Image* ip_noisify (Image* src, double alpha)
{
  Image* result = new Image(src->getWidth(),
                            src->getHeight(),
                            src->getChannels(),
                            src->getBits());

  double v;

  for (int j = 0; j < src->getHeight(); ++j)
    for (int i = 0; i < src->getWidth(); ++i)
      for (int c = 0; c < src->getChannels(); ++c)
      {
        v = (1 - alpha) * src->getPixel(i,j,c)
          + alpha * drand48();
        result->setPixel(i,j,c, clamp(v,0,1));
      }

  return result;
}


/*
 * round each pixel to the nearest value in the new number of bits
 */
Image* ip_quantize_simple (Image* src, int bitsPerChannel)
{
  if (bitsPerChannel < 1 || bitsPerChannel > 8)
  {
    cerr << "Bits per channel must be between 1 and 8!" << endl;
    return NULL;
  }

  if (bitsPerChannel == src->getBits())
    return NULL;

  Image* quantized = new Image(src->getWidth(),
                               src->getHeight(),
                               src->getChannels(),
                               bitsPerChannel);

  Pixel pixel;
  for (int j = 0; j < src->getHeight(); ++j)
    for (int i = 0; i < src->getWidth(); ++i)
      quantized->setPixel(i,j, src->getPixel(i,j,pixel));

  return quantized;
}


static double Bayer2[4][4] = { { 16,  8, 14,  6 },
                               {  4, 12,  2,  10 },
                               { 13,  5, 15,  7 },
                               {  1,  9,  3, 11 } };


/*
 * dither each pixel to the nearest value in the new number of bits
 * using a static 4x4 matrix
 */
Image* ip_quantize_ordered (Image* src, int bitsPerChannel)
{
  if (bitsPerChannel < 1 || bitsPerChannel > 8)
  {
    cerr << "Bits per channel must be between 1 and 8!" << endl;
    return NULL;
  }

  if (bitsPerChannel == src->getBits())
    return NULL;

  Image* quantized = new Image(src->getWidth(),
                               src->getHeight(),
                               src->getChannels(),
                               bitsPerChannel);


  double levels = pow(2, bitsPerChannel);
  double v;
  for (int j = 0; j < src->getHeight(); ++j)
    for (int i = 0; i < src->getWidth(); ++i)
      for (int c = 0; c < src->getChannels(); ++c)
      {
        v = src->getPixel(i,j,c)
          + (double) ((2*Bayer2[i%4][j%4]-1))/((levels-1.0)*32.0)
          -0.5/(levels-1) ;
/*
	if (i==50)
		cout << "source=" << src->getPixel(i,j,c) << " noise=" 
          		(double) ((2*Bayer2[i%4][j%4]-1))/((levels-1.0)*32.0)
			<< " output=" << v << endl;*/
        quantized->setPixel(i,j,c, clamp(v, 0.0, 1.0));
      }

/*  old
  double scale = 1.0 / pow(2, bitsPerChannel + 1);
  double v;
  for (int j = 0; j < src->getHeight(); ++j)
    for (int i = 0; i < src->getWidth(); ++i)
      for (int c = 0; c < src->getChannels(); ++c)
      {
        v = src->getPixel(i,j,c)
          + (scale * Bayer2[i%4][j%4] / 15.0)
          - (scale * 0.5);
        quantized->setPixel(i,j,c, clamp(v, 0.0, 1.0));
      }
*/

  return quantized;
}


/*
 * dither each pixel to the nearest value in the new number of bits
 * using error diffusion
 */
Image* ip_quantize_fs (Image* src, int bitsPerChannel)
{
  if (bitsPerChannel < 1 || bitsPerChannel > 8)
  {
    cerr << "Bits per channel must be between 1 and 8!" << endl;
    return NULL;
  }

  if (bitsPerChannel >= src->getBits())
    return NULL;

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

  double*** error = new double**[src->getWidth() + 1];
  for (int i = 0; i < src->getWidth() + 1; ++i)
  {
    error[i] = new double*[src->getHeight()];

    for (int j = 0; j < src->getHeight(); ++j)
    {
      error[i][j] = new double[src->getChannels()];

      for (int c = 0; c < src->getChannels(); ++c)
        error[i][j][c] = 0;
    }
  }

  double v;
  double e;
  for (int j = 0; j < src->getHeight(); ++j)
    for (int i = 0; i < src->getWidth(); ++i)
      for (int c = 0; c < src->getChannels(); ++c)
      {
        v = src->getPixel(i,j,c) + error[i][j][c];
        result->setPixel(i,j,c, clamp(v,0.0,1.0));

        e = v - result->getPixel(i,j,c);

        error[i+1][j][c] += (e * 7.0 / 16.0);
        if (j+1 < src->getHeight())
        {
          if (i-1 >= 0)
            error[i-1][j+1][c] += (e * 3.0 / 16.0);

          error[i][j+1][c]   += (e * 5.0 / 16.0);
          error[i+1][j+1][c] += (e * 1.0 / 16.0);
        }
      }

  for (int i = 0; i < src->getWidth() + 1; ++i)
  {
    for (int j = 0; j < src->getHeight(); ++j)
      delete[] error[i][j];
    delete[] error[i];
  }
  delete[] error;

  return result;
}


/*
 * rotate image using one of three sampling techniques
 */
Image* ip_rotate (Image* src, double theta, int x, int y, int mode, 
                  int size, double sigma)
{
  if (x < 0 || y < 0 || x >= src->getWidth() ||
      y >= src->getHeight())
  {
    cerr << "Center of rotation must be within the image!" << endl;
    return NULL;
  }

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

  if (mode == I_GAUSSIAN)
  {
    if (size % 2 != 0 || size < 1)
    {
      cerr << "Filter size must be an even, positive integer!" << endl;
      return NULL;
    }
    if (sigma == 0)
    {
      cerr << "Sigma must be non-zero!" << endl;
      return NULL;
    }
  }

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

  double x0, y0;
  double cosTheta = cos(deg2rad(theta));
  double sinTheta = sin(deg2rad(theta));

  for (int j = 0; j < src->getHeight(); ++j)
    for (int i = 0; i < src->getWidth(); ++i)
    {
      x0 = (i - x) * cosTheta - (j - y) * sinTheta + x;
      y0 = (i - x) * sinTheta + (j - y) * cosTheta + y;

      for (int c = 0; c < src->getChannels(); ++c)
      {
        switch (mode)
        {
          case I_NEAREST:
            rotated->setPixel(i, j, c,
                nearestInterpolate(src, x0, y0, c));
            break;

          case I_BILINEAR:
            rotated->setPixel(i, j, c,
                bilinearInterpolate(src, x0, y0, c));
            break;

          case I_GAUSSIAN:
            rotated->setPixel(i, j, c,
                gaussianInterpolate(src, x0, y0, c, size, sigma));
            break;
        }
      }
    }

  return rotated;
}


/*
 * interpolate with the greyscale version of the image
 */
Image* ip_saturate (Image* src, double alpha)
{
  if (src->getChannels() == 1)
  {
    cerr << "Can't adjust saturation on a greyscale image!" << endl;
    return NULL;
  }

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

  double v;
  Pixel  p;
  for (int j = 0; j < src->getHeight(); ++j) 
    for (int i = 0; i < src->getWidth(); ++i)
    {                      
      v  = src->getPixel(i,j,RED)   * 0.3;
      v += src->getPixel(i,j,GREEN) * 0.59;
      v += src->getPixel(i,j,BLUE)  * 0.11;

      p.r = v;
      p.g = v;
      p.b = v;

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

  Image* result = interpolate(greyImage, src, alpha);

  delete greyImage;
  return result;
}


/*
 * scale image using one of three sampling techniques
 */
Image* ip_scale (Image* src, double xFac, double yFac, int mode, 
                 int size, double sigma)
{
  if (xFac == 0 || yFac == 0)
  {
    cerr << "Scale factors must be non-zero!" << endl;
    return NULL;
  }

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

  if (mode == I_GAUSSIAN)
  {
    if (size % 2 != 1 || size < 1)
    {
      cerr << "Filter size must be an even, positive integer!" << endl;
      return NULL;
    }
    if (sigma == 0)
    {
      cerr << "Sigma must be non-zero!" << endl;
      return NULL;
    }
  }

  int newWidth  = (int) (src->getWidth()  * xFac);
  int newHeight = (int) (src->getHeight() * yFac);

  Image* scaled = new Image(newWidth,
                            newHeight,
                            src->getChannels(),
                            src->getBits());

  double x, y;
  for (int j = 0; j < newHeight; ++j)
  {
    y = j / yFac;

    for (int i = 0; i < newWidth; ++i)
    {
      x = i / xFac;

      for (int c = 0; c < src->getChannels(); ++c)
      {
        switch (mode)
        {
          case I_NEAREST:
            scaled->setPixel(i, j, c,
                nearestInterpolate(src, x, y, c));
            break;

          case I_BILINEAR:
            scaled->setPixel(i, j, c,
                bilinearInterpolate(src, x, y, c));
            break;

          case I_GAUSSIAN:
            scaled->setPixel(i, j, c,
                gaussianInterpolate(src, x, y, c, size, sigma));
            break;
        }
      }
    }
  }

  return scaled;
}


/*
 * create a one bit per channel image, full color where the intensity is
 * above the threshold, no color where it isn't
 */
Image* ip_threshold (Image* src, double cutoff)
{
  if (cutoff < 0.0 || cutoff > 1.0)
  {
    cerr << "Cutoff must be between 0.0 and 1.0!" << endl;
    return NULL;
  }
  

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

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

      switch (src->getChannels())
      {
        case 3:
		  if (pixel.b>cutoff)
			  pixel.b=1;
		  else 
			  pixel.b=0;
		  if (pixel.g>cutoff)
			  pixel.g=1;
		  else 
			  pixel.g=0;

        case 1:
          if (pixel.r > cutoff)
			  pixel.r=1;
		  else
			  pixel.r=0;
          thresh->setPixel(i, j, pixel);
          break;

        default:
          break;
      }
    }

  return thresh;
}


/*
 * perform do what you like
 */
Image* ip_warp (Image* src)
{
  cerr << "Select a warp" 		<< endl 
       					<< endl
       << " 1)  Mosaic"			<< endl
       << " 2)  Dots"			<< endl
       << " 3)  Ripple"			<< endl
       << " 4)  Emboss"			<< endl
       << " 5)  Lines"			<< endl
       << " 6)  Stamp"			<< endl
       << " 7)  Swirl"			<< endl
       << " 8)  Spherical"		<< endl
       << " 9)  Rotate Hue"		<< endl
       << " 10) Voronoi (unfinished)"	<< endl
       					<< endl
       << "> " << flush;

  int selection;
  cin >> selection;

  switch (selection)
  {
  case 1:
    return warp_mosaic(src);

  case 2:
    return warp_dots(src);

  case 3:
    return warp_ripple(src);

  case 4:
    return warp_emboss(src);

  case 5:
    return warp_lines(src);

  case 6:
    return warp_stamp(src);

  case 7:
    return warp_swirl(src);

  case 8:
    return warp_spherical(src);

  case 9:
    return warp_rotate_hue(src);

  case 10:
    return warp_voronoi(src);

  default:
    cerr << "Invalid selection!" << endl;
    return NULL;
  }

  return NULL;
}