/* Class which describes a single colour channel of an image in a generalized
 * way allowing for easy and accurate scaling to any size.
 * 
 * Written by: Chris Studholme
 * Copyright:  GPL (http://www.fsf.org/copyleft/gpl.html)
 * $Id: GImageComponent.cpp,v 1.8 2003/05/28 02:30:01 cvs Exp $
 */

#include <config.h>
#include <math.h>
#include <stdio.h>

#include "GImageComponent.h"

//#define SHOW_STATS

#define DEFAULT_ERROR 0.1

#define MAXITTER 20
#define MINITTER 3

inline short boundcheck(int x) {
  if (x<-32768) return -32768;
  if (x>32767) return 32767;
  return x;
}

template<class T> inline T MAX(T a, T b) {
  return (a>b ? a : b);
}

template<class T> inline T MAX(T a, T b, T c) {
  if (b>a) a=b;
  if (c>a) a=c;
  return a;
}

template<class T> inline T MAX(T a, T b, T c, T d) {
  if (b>a) a=b;
  if (c>a) a=c;
  if (d>a) a=d;
  return a;
}

template<class T> inline T MAX(T a, T b, T c, T d, T e, T f, T g, T h) {
  if (b>a) a=b;
  if (c>a) a=c;
  if (d>a) a=d;
  if (e>a) a=e;
  if (f>a) a=f;
  if (g>a) a=g;
  if (h>a) a=h;
  return a;
}

template<class T> inline T MIN(T a, T b) {
  return (a<b ? a : b);
}

template<class T> inline T MIN(T a, T b, T c) {
  if (b<a) a=b;
  if (c<a) a=c;
  return a;
}

template<class T> inline T MIN(T a, T b, T c, T d) {
  if (b<a) a=b;
  if (c<a) a=c;
  if (d<a) a=d;
  return a;
}

template<class T> inline T MIN(T a, T b, T c, T d, T e, T f, T g, T h) {
  if (b<a) a=b;
  if (c<a) a=c;
  if (d<a) a=d;
  if (e<a) a=e;
  if (f<a) a=f;
  if (g<a) a=g;
  if (h<a) a=h;
  return a;
}

/* ================ GImageComponent methods  ================ */

GImageComponent::GImageComponent(int width, int height, float error) {
  w=width; h=height;
  tl = new short[w*h];
  
  maxerror = error>0 ? error : DEFAULT_ERROR;
  
  totalerror=0;
}

/* x1<x2, y1<y2,
 * coords are float coords * 256 converted to integer, for example:
 *  (10,12.5) passed as (2560,3200)
 */
int GImageComponent::getPixelI(int x1, int y1, int x2, int y2) {
  /**** DO NOT USE!  Not ready yet. ****/

  // transform coordinates
  x1-=128;  y1-=128;
  x2-=128;  y2-=128;

  int ix1 = x1>>8;
  int ix2 = x2>>8;
  int iy1 = y1>>8;
  int iy2 = y2>>8;
  
  if (ix1<0) ix1=0;
  else if (ix1>=w-1) ix1=w-2;
  if (iy1<0) iy1=0;
  else if (iy1>=h-1) iy1=h-2;
  if (ix2<0) ix2=0;
  else if (ix2>=w-1) ix2=w-2;
  if (iy2<0) iy2=0;
  else if (iy2>=h-1) iy2=h-2;

  int logx=0;
  for (int i=ix2-ix1; i!=0; ++logx,i>>=1) ;   // calculate log2(ix2-ix1)
  int logy=0;
  for (int i=iy2-iy1; i!=0; ++logy,i>>=1) ;   // calculate log2(iy2-iy1)
  logy = logy>logx ? logy-logx : 0;
  
  int result=0;
  int by1 = y1 - (iy1<<8);
  for (int y=iy1; y<=iy2; ++y) {
    int by2 = y2 - (y<<8);
    if ((y<(h-2))&&(by2>256)) by2=256;
    int sy = by2+by1;
    int bx1 = x1 - (ix1<<8);

    int r1=0;
    for (int x=ix1; x<=ix2; ++x) {
      int bx2 = x2 - (x<<8);
      if ((x<(w-2))&&(bx2>256)) bx2=256;
      int sx = bx2+bx1;

      int p1 = (512-sx)*tl[y*w+x]     +sx*tl[y*w+x+1];      // 24 bits
      int p2 = (512-sx)*tl[(y+1)*w+x] +sx*tl[(y+1)*w+x+1];  // 24 bits
      int p3 = (p1*512>>4) + ((8+p2-p1)>>4)*sy;    // net shift right 4, 31 bits

      int dx = bx2-bx1;
      r1 += (((256+p3)>>9)*dx)>>logx;
      bx1=0;
    }

    int dy = by2-by1;
    result += (r1/(x2-x1)*dy)>>logy;     // 31-logx bits
    by1=0;
  }
  return ( result/(y2-y1) )/4;
}

/* x1<x2, y1<y2 */
int GImageComponent::getPixel(float x1, float y1, float x2, float y2) {
  TransformCoordinates(x1,y1);
  TransformCoordinates(x2,y2);
  
  int ix1 = (int)floor(x1);
  int ix2 = (int)floor(x2);
  int iy1 = (int)floor(y1);
  int iy2 = (int)floor(y2);
  
  if (ix1<0) ix1=0;
  else if (ix1>=w-1) ix1=w-2;
  if (iy1<0) iy1=0;
  else if (iy1>=h-1) iy1=h-2;
  if (ix2<0) ix2=0;
  else if (ix2>=w-1) ix2=w-2;
  if (iy2<0) iy2=0;
  else if (iy2>=h-1) iy2=h-2;
  
  double result=0;
  double by1 = y1-iy1;
  for (int y=iy1; y<=iy2; ++y) {
    double by2 = y2-y;
    if ((y<(h-2))&&(by2>1)) by2=1;
    double sy = by2+by1;
    double dy = by2-by1;
    double bx1 = x1-ix1;

    for (int x=ix1; x<=ix2; ++x) {
      double bx2 = x2-x;
      if ((x<(w-2))&&(bx2>1)) bx2=1;
      double sx = bx2+bx1;
      double dx = bx2-bx1;

      result += ( (sx-2)*(sy-2)*tl[y*w+x]
		  -sx*(sy-2)*tl[y*w+x+1]
		  -(sx-2)*sy*tl[(y+1)*w+x]
		  +sx*sy*tl[(y+1)*w+x+1] 
		  )*dx*dy;
      bx1=0;
    }
    by1=0;
  }
  return ((int)(result/(x2-x1)/(y2-y1)))/4;
}

/* y1<y2
 */
void GImageComponent::getScanLineHorz(short* pixels, float x1, float y1,
				      float x2, float y2, int npixels) {
  float delta=(x2-x1)/npixels;
  float start=x1;
  if (delta>0)
    for (int i=0; i<npixels; ++i) {
      float end=start+delta;
      pixels[i] = getPixel(start,y1,end,y2);
      start=end;
    }
  else
    for (int i=0; i<npixels; ++i) {
      float end=start+delta;
      pixels[i] = getPixel(end,y1,start,y2);
      start=end;
    }

}

/* x1<x2
 */
void GImageComponent::getScanLineVert(short* pixels, float x1, float y1,
				      float x2, float y2, int npixels) {
  float delta=(y2-y1)/npixels;
  float start=y1;
  if (delta>0)
    for (int i=0; i<npixels; ++i) {
      float end=start+delta;
      pixels[i] = getPixel(x1,start,x2,end);
      start=end;
    }
  else
    for (int i=0; i<npixels; ++i) {
      float end=start+delta;
      pixels[i] = getPixel(x1,end,x2,start);
      start=end;
    }

}


int GImageComponent::ReduceGrain(bool extra) {

  int noadjust=0;
  for (int y=1; y<h-1; ++y)
    for (int x=1; x<w-1; ++x) {

      // check if central pixel is too bright
      int max = (extra ? 
		 MAX(tl[y*w+x-1],tl[y*w+x+1],tl[(y-1)*w+x],tl[(y+1)*w+x]) :
		 MAX(tl[y*w+x-1],tl[y*w+x+1],tl[(y-1)*w+x],
		     tl[(y+1)*w+x],tl[(y-1)*w+x-1],tl[(y-1)*w+x+1],
		     tl[(y+1)*w+x-1],tl[(y+1)*w+x+1]));
      if (tl[y*w+x]>max) {
	tl[y*w+x]=max;
	continue;
      }

      // check if central pixel is too dim
      int min = (extra ? 
		 MIN(tl[y*w+x-1],tl[y*w+x+1],tl[(y-1)*w+x],tl[(y+1)*w+x]) :
		 MIN(tl[y*w+x-1],tl[y*w+x+1],tl[(y-1)*w+x],
		     tl[(y+1)*w+x],tl[(y-1)*w+x-1],tl[(y-1)*w+x+1],
		     tl[(y+1)*w+x-1],tl[(y+1)*w+x+1]));
      if (tl[y*w+x]<min) {
	tl[y*w+x]=min;
	continue;
      }

      // keep track of pixels that were not changed
      ++noadjust;
    }

  return w*h-noadjust;
}



/* ================ GImageComponent0 methods  ================ */

GImageComponent0::GImageComponent0(const unsigned char* imagedata, 
				   int width, int height, 
				   float maxerr)
  : GImageComponent(width,height,maxerr) {

  int itter;
  long error;

  //int merror = (int)(maxerror*width*height*128);
  int werror = (int)(maxerror*width*128);
  int herror = (int)(maxerror*height*128);

#ifdef SHOW_STATS
  fprintf(stderr," error: werror is %ld, herror is %ld\n",werror,herror);
#endif  

  /* initialize array
   */
  short* id = new short[w*h];
  for (int y=0; y<h; ++y)
    for (int x=0; x<w; ++x)
      tl[y*w+x] = id[y*w+x] = (((short)imagedata[y*w+x])-128)<<7;
    
  /* top row 
   */
  itter=0;
  error=werror+1;
  do {
    ++itter;
    
    // calculate better answer
    for (int x=1; x<w-1; ++x) 
      tl[x] = boundcheck((8*id[x] - tl[x-1] - tl[x+1])/6);

    if (itter<MINITTER)
      continue;
    
    // calculate error
    error=0;
    for (int x=1; x<w-1; ++x) {
      int residule = (tl[x-1]+6*tl[x]+tl[x+1])/8 - id[x];
      error += residule<0?-residule:residule;
    }
  } while ((error>werror)&&(itter<MAXITTER));
  
  totalerror+=error;
#ifdef SHOW_STATS
  fprintf(stderr,"   top: error is %ld after %d itterations\n",error,itter);
#endif  

  /* bottom row
   */
  itter=0;
  error=werror+1;
  do {
    ++itter;
    
    // calculate better answer
    for (int x=1; x<w-1; ++x) 
      tl[(h-1)*w+x] = boundcheck((8*id[(h-1)*w+x] 
				   -tl[(h-1)*w+x-1] 
				   -tl[(h-1)*w+x+1])/6);
    
    if (itter<MINITTER)
      continue;
    
    // calculate error
    error=0;
    for (int x=1; x<w-1; ++x) {
      int residule = (tl[(h-1)*w+x-1]+6*tl[(h-1)*w+x]+tl[(h-1)*w+x+1])/8
	- id[(h-1)*w+x];
      error += residule<0?-residule:residule;
    }
  } while ((error>werror)&&(itter<MAXITTER));
  
  totalerror+=error;
#ifdef SHOW_STATS
  fprintf(stderr,"bottom: error is %ld after %d itterations\n",error,itter);
#endif
  
  /* left column 
   */
  itter=0;
  error=herror+1;
  do {
    ++itter;
    
    // calculate better answer
    for (int y=1; y<h-1; ++y) 
      tl[y*w] = boundcheck((8*id[y*w] 
			     - tl[(y-1)*w] - tl[(y+1)*w])/6);
    
    if (itter<MINITTER)
      continue;
    
    // calculate error
    error=0;
    for (int y=1; y<h-1; ++y) {
      int residule = (6*tl[y*w]+tl[(y-1)*w]+tl[(y+1)*w])/8
	- id[y*w];
      error += residule<0?-residule:residule;
    }
  } while ((error>herror)&&(itter<MAXITTER));
  
  totalerror+=error;
#ifdef SHOW_STATS
  fprintf(stderr,"left  : error is %ld after %d itterations\n",error,itter);
#endif
  
  /* right column 
   */
  itter=0;
  error=herror+1;
  do {
    ++itter;
    
    // calculate better answer
    for (int y=1; y<h-1; ++y) 
      tl[y*w+w-1] = boundcheck((8*id[y*w+w-1] 
				 - tl[(y-1)*w+w-1] - tl[(y+1)*w+w-1])/6);
    
    if (itter<MINITTER)
      continue;
    
    // calculate error
    error=0;
    for (int y=1; y<h-1; ++y) {
      int residule = (6*tl[y*w+w-1]+tl[(y-1)*w+w-1]+tl[(y+1)*w+w-1])/8
	- id[y*w+w-1];
      error += residule<0?-residule:residule;
    }
  } while ((error>herror)&&(itter<MAXITTER));
  
  totalerror+=error;
#ifdef SHOW_STATS
  fprintf(stderr,"right : error is %ld after %d itterations\n",error,itter);
#endif
  
  /* the middle pixels
   */
  itter=0;
  int mitter=2;   // maximum itterations
  for (int head=1, tail=1; tail<h-1;) {
    ++itter;

    // calculate better answer (head to tail inclusive)
    for (int y=head; y>=tail; --y)
      for (int x=1; x<w-1; ++x)
	tl[y*w+x] = boundcheck((64*id[y*w+x] 
				-6*tl[(y-1)*w+x] -6*tl[(y+1)*w+x]
				-6*tl[y*w+x-1]   -6*tl[y*w+x+1]
				-tl[(y-1)*w+x-1] -tl[(y-1)*w+x+1]
				-tl[(y+1)*w+x-1] -tl[(y+1)*w+x+1])/36);

    // check answer at tail (fast method)
    if ((maxerr<0)&&(itter>=2)) {
      ++tail;
      --itter;
      totalerror+=werror;
    }

    // check answer at tail (accurate method)
    else while ((itter>=MINITTER)&&(tail<h-1)) {
      error=0;
      int y=tail;
      for (int x=1; x<w-1; ++x) {
	int residule = (36*tl[y*w+x]+
			6*tl[(y-1)*w+x]+6*tl[(y+1)*w+x]+
			6*tl[y*w+x-1]  +6*tl[y*w+x+1]+
			tl[(y-1)*w+x-1]+tl[(y-1)*w+x+1]+
			tl[(y+1)*w+x-1]+tl[(y+1)*w+x+1])/64 -id[y*w+x];
	error += residule<0?-residule:residule;
      }
      
      // not good enough, calculate tail again
      if ((error>werror)&&(itter<MAXITTER))
	break;
      
      // advance tail and check next line
      ++tail;
      totalerror+=error;
      if (itter>mitter)
	mitter=itter;
      --itter;
    }
    
    // advance head
    if (head<h-2)
      ++head;
  }

  delete[] id;
  
#ifdef SHOW_STATS
  fprintf(stderr,"middle: error is %ld after %d itterations\n",error,mitter);
#endif
  
  totalerror=totalerror/(w*h*128);
}