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

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

#include "GICBayer.h"

//#define SHOW_STATS

#define MAXITTER 20
#define MINITTER 4
#define MINITTER_FAST 2
#define ITTER_INCR 1
  
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;
}

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

// returns c1
inline short extrapolate(short cmin, short cmax, short c2, 
			 short omin, short omax, short o1, short o2) { 
  if (o1==o2)
    return c2;
  else if (o1<o2)
    return c2 + (c2-cmin)*(o1-o2)/(o2-omin);
  else
    return c2 + (cmax-c2)*(o1-o2)/(omax-o2);
}

// returns c1
inline short interpolate(short c2, short c3, short o1, short o2, short o3) {
  if ((o1>o2)&&(o1>o3)) 
    return MAX(c2,c3);
  else if ((o1<o2)&&(o1<o3)) 
    return MIN(c2,c3);
  else if ((o2-o3)*(c2-c3)<=0)
    return (c2+c3)/2;
  else 
    return (c3*(o1-o2)+c2*(o3-o1)) / (o3-o2);
}


/* ================ GICBayer_Green methods  ================ */

inline void green_better_middle(short* tl, short* id, int y, int w) {
  for (int x=2-(y&1); x<w-1; x+=2)
    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);
}
inline void green_better_edges(short* tl, short* id, int w, int h) {
  // top row
  for (int x=2; x<w; x+=2)
    tl[x] = boundcheck((8*id[x]-tl[x-1]-tl[x+1])/6);
  // bottom row
  for (int x=1; x<w-1; x+=2)
    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);
  // left column
  for (int y=2; y<h; y+=2)
    tl[y*w] = boundcheck((8*id[y*w]-tl[(y-1)*w]-tl[(y+1)*w])/6);
  // right column
  for (int y=1; y<h-1; y+=2)
    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);
}
inline void green_make_up(short* tl, short* id, int y, int w, int h) {
  for (int x=1-(y&1); x<w; x+=2) {
    short gmin=32767;
    short gmax=-32768;
    short cmin=id[y*w+x];
    short cmax=id[y*w+x];
    if (x>0) {
      gmin = MIN(gmin,tl[y*w+x-1]);
      gmax = MAX(gmax,tl[y*w+x-1]);
      if (x>1) {
	cmin = MIN(cmin,id[y*w+x-2]);
	cmax = MAX(cmax,id[y*w+x-2]);
      }
    }
    if (x<w-1) {
      gmin = MIN(gmin,tl[y*w+x+1]);
      gmax = MAX(gmax,tl[y*w+x+1]);
      if (x<w-2) {
	cmin = MIN(cmin,id[y*w+x+2]);
	cmax = MAX(cmax,id[y*w+x+2]);
      }
    }
    if (y>0) {
      gmin = MIN(gmin,tl[y*w+x-w]);
      gmax = MAX(gmax,tl[y*w+x-w]);
      if (y>1) {
	cmin = MIN(cmin,id[(y-2)*w+x]);
	cmax = MAX(cmax,id[(y-2)*w+x]);
      }
    }
    if (y<h-1) {
      gmin = MIN(gmin,tl[y*w+x+w]);
      gmax = MAX(gmax,tl[y*w+x+w]);
      if (y<h-2) {
	cmin = MIN(cmin,id[(y+2)*w+x]);
	cmax = MAX(cmax,id[(y+2)*w+x]);
      }
    }
    
    int count=0;
    int pixel=0;
    if (x>1) {
      pixel += extrapolate(gmin,gmax,tl[y*w+x-1],cmin,cmax,id[y*w+x],id[y*w+x-2]);
      ++count;
    }
    if (x<w-2) {
      pixel += extrapolate(gmin,gmax,tl[y*w+x+1],cmin,cmax,id[y*w+x],id[y*w+x+2]);
      ++count;
    }
    if (y>1) {
      pixel += extrapolate(gmin,gmax,tl[(y-1)*w+x],cmin,cmax,id[y*w+x],id[(y-2)*w+x]);
      ++count;
    }
    if (y<h-2) {
      pixel += extrapolate(gmin,gmax,tl[(y+1)*w+x],cmin,cmax,id[y*w+x],id[(y+2)*w+x]);
      ++count;
    }
    tl[y*w+x]=pixel/count;
  }
}
inline long green_error_top(short* tl, short* id, int w) {
  long error=0;
  for (int x=2; x<w; x+=2) {
    int residule = (tl[x-1]+6*tl[x]+tl[x+1])/8 - id[x];
    error += (residule*residule+8)>>4;
  }
  return error;
}
inline long green_error_middle(short* tl, short* id, int y, int w) {
  long error=0;
  if (y&1) {
    // right column
    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*residule+8)>>4;
  }
  else {
    // left column
    int residule = (6*tl[y*w]+tl[(y-1)*w]+tl[(y+1)*w])/8 - id[y*w];
    error += (residule*residule+8)>>4;
  }
  // middle pixels
  for (int x=2-(y&1); x<w-1; x+=2) {
    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*residule+8)>>4;
  }
  return error;
}
inline long green_error_bottom(short* tl, short* id, int w, int h) {
  long error=0;
  for (int x=1; x<w-1; x+=2) {
    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*residule+8)>>4;
  }
  return error;
}

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

  /* initialize arrays
   */
  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;

    // initialize made-up points
    if (y>1)
      green_make_up(tl,id,y-2,w,h);
  }
  green_make_up(tl,id,h-2,w,h);
  green_make_up(tl,id,h-1,w,h);
  
  /* main loop to calculate points
   */
  int merror=(int)(maxerror*maxerror*128*128*w*h/2)>>4;
  long error;
  int itter=1;
  do {
    ++itter;
    error=0;

    // calculate better answer (edges)
    green_better_edges(tl,id,w,h);

    // do all three steps in one loop for better cache use
    for (int i=0; i<h; ++i) {
      
      // better answer (middle pixels)
      if (i<h-2)
	green_better_middle(tl,id,i+1,w);
      
      // made up points (guessed from red and blue)
      int y=i;
      green_make_up(tl,id,i,w,h);
      
      // calculate error
      if (itter>=MINITTER) {
	if (i==1)
	  // top row
	  error += green_error_top(tl,id,w);
	else if (i>1)
	  // middle rows
	  error += green_error_middle(tl,id,i-1,w);
      }
    }

    // quit if we want it done fast
    if ((maxerr==-1)&&(itter>=MINITTER_FAST))
      break;

    if (itter>=MINITTER) {
      // error (bottom row)
      error+= green_error_bottom(tl,id,w,h);

      // really aweful
      if (error<0)
	continue;

      // solution is good enough => quit
      if (error<merror)
	break;
    }
    
  } while (itter<MAXITTER);

  // calculate better answer one more time for luck :)
  green_better_edges(tl,id,w,h);
  for (int y=1; y<h-1; ++y)
    green_better_middle(tl,id,y,w);

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


/* ================ GICBayer_Red methods  ================ */

inline void red_make_up_horz(short* tl, const short* gtl, int y, int w) {
  // green pixels between two reds horizontally
  // x=0
  tl[y*w] = tl[y*w+1];  // could be better, but it's on the left edge so who cares
  for (int x=2; x<w; x+=2)
    tl[y*w+x] = interpolate(tl[y*w+x-1],tl[y*w+x+1],gtl[y*w+x],gtl[y*w+x-1],gtl[y*w+x+1]);
}
inline void red_make_up_vert(short* tl, const short* gtl, int y, int w, int h) {
  // green pixels between two reds vertically
  if (y<h-1)
    for (int x=1; x<w; x+=2)
      tl[y*w+x] = interpolate(tl[y*w+x-w],tl[y*w+x+w],gtl[y*w+x],gtl[y*w+x-w],gtl[y*w+x+w]);
  else
    for (int x=1; x<w; x+=2)
      tl[y*w+x] = tl[(y-1)*w+x];  // could be better, but it's on the bottom so who cares
}
inline void red_make_up_blue(short* tl, const short* gtl, int y, int w, int h) {
  // blue pixels (note: we only use blue data indirectly through green)
  for (int x=0; x<w; x+=2) {
    short gmin = MIN(gtl[y*w+x],gtl[y*w+x-w+1]);
    short gmax = MAX(gtl[y*w+x],gtl[y*w+x-w+1]);
    short rmin = tl[y*w+x-w+1];
    short rmax = tl[y*w+x-w+1];
    if (x>0) {
      gmin = MIN(gmin,gtl[y*w+x-w-1]);
      gmax = MAX(gmax,gtl[y*w+x-w-1]);
      rmin = MIN(rmin,tl[y*w+x-w-1]);
      rmax = MAX(rmax,tl[y*w+x-w-1]);
    }
    if (y<h-1) {
      gmin = MIN(gmin,gtl[y*w+x+w+1]);
      gmax = MAX(gmax,gtl[y*w+x+w+1]);
      rmin = MIN(rmin,tl[y*w+x+w+1]);
      rmax = MAX(rmax,tl[y*w+x+w+1]);
      if (x>0) {
	gmin = MIN(gmin,gtl[y*w+x+w-1]);
	gmax = MAX(gmax,gtl[y*w+x+w-1]);
	rmin = MIN(rmin,tl[y*w+x+w-1]);
	rmax = MAX(rmax,tl[y*w+x+w-1]);
      }
    }
    
    int count=1;
    int pixel = extrapolate(rmin,rmax,tl[y*w+x-w+1],gmin,gmax,gtl[y*w+x],gtl[y*w+x-w+1]);
    if (x>0) {
      pixel += extrapolate(rmin,rmax,tl[y*w+x-w-1],gmin,gmax,gtl[y*w+x],gtl[y*w+x-w-1]);
      ++count;
    }
    if (y<h-1) {
      pixel += extrapolate(rmin,rmax,tl[y*w+x+w+1],gmin,gmax,gtl[y*w+x],gtl[y*w+x+w+1]);
      ++count;
      if (x>0) {
	pixel += extrapolate(rmin,rmax,tl[y*w+x+w-1],gmin,gmax,gtl[y*w+x],gtl[y*w+x+w-1]);
	++count;
      }
    }
    tl[y*w+x]=pixel/count;
  }
}
inline long red_error_top(short* tl, short* id, int w) {
  long error=0;
  for (int x=1; x<w-1; x+=2) {
    int residule = (tl[x-1]+6*tl[x]+tl[x+1])/8 - id[x/2];
    error += (residule*residule+8)>>4;
  }
  return error;
}
inline long red_error_middle(short* tl, short* id, int y, int w) {
  long error=0;
  id += y*w/4;
  tl += y*w;
  // middle pixels
  for (int x=1; x<w-1; x+=2) {
    int residule = (36*tl[x]+
		    6*tl[-w+x] +6*tl[w+x]+
		    6*tl[x-1]  +6*tl[x+1]+
		    tl[-w+x-1] +tl[-w+x+1]+
		    tl[w+x-1]  +tl[w+x+1])/64 - id[x/2];
    error += (residule*residule+8)>>4;
  }
  
  // right column
  int residule = (6*tl[w-1]+tl[-1]+tl[2*w-1])/8 - id[w/2-1];
  error += (residule*residule+8)>>4;
  return error;
}
inline void red_better_top(short* tl, short* id, int w, int h) {
  // top row
  for (int x=1; x<w-1; x+=2) 
    tl[x] = boundcheck((8*id[x/2]-tl[x-1]-tl[x+1])/6);
}
inline void red_better_middle(short* tl, short* id, int y, int w) {
  id += y*w/4;
  tl += y*w;
  // middle
  for (int x=1; x<w-1; x+=2)
    tl[x] = boundcheck((64*id[x/2] 
			-6*tl[-w+x] -6*tl[w+x]
			-6*tl[x-1]  -6*tl[x+1]
			-tl[-w+x-1] -tl[-w+x+1]
			-tl[w+x-1]  -tl[w+x+1])/36);
  // right column
  tl[w-1] = boundcheck((8*id[w/2-1]-tl[-1]-tl[2*w-1])/6);
}

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

  if ((green.getWidth()!=w)||(green.getHeight()!=h)) {
    fprintf(stderr,"GICBayer_Red: FATAL ERROR! green component has different dimensions\n");
    return;
  }

  const short* gtl = green.getPointArray();
  short* id = new short[w*h/4];

  int itter=0;
  int nitter=2*(maxerr==-1?MINITTER_FAST:MINITTER);
  long error;
  int merror=(maxerr==-1?-1:((int)(maxerror*maxerror*128*128*w*h/4))>>4);

  // main loop
  do {
    // do nitter/2 itterations in batch
    error=0;
    for (int head=0; head<h+nitter; head+=2) {
      for (int i=0; i<nitter; i+=2) {
	int y=head-i;
	
	// initialize arrays (first itteration only)
	if ((itter+i==0)&&(y<h))
	  for (int x=0; x<w; x+=2)
	    tl[y*w+x+1] = id[y*w/4+x/2] = (((short)imagedata[y*w+x+1])-128)<<7;
	
	y-=2;
	if ((y>=0)&&(y<h)) {
	  // calculate made up points (guessed using green as guide)
	  red_make_up_horz(tl,gtl,y,w);
	  red_make_up_vert(tl,gtl,y+1,w,h);
	  red_make_up_blue(tl,gtl,y+1,w,h);

	  // calculate error if necessary
	  if ((i==nitter-2)&&(merror>0)) {
	    if (y==0)
	      error+=red_error_top(tl,id,w);
	    else
	      error+=red_error_middle(tl,id,y,w);
	  }
	  
	  // calculate better answer
	  if (y==0)
	    red_better_top(tl,id,w,h);
	  else
	    red_better_middle(tl,id,y,w);
	}
      }

    }

    itter+=nitter/2;
    nitter=2*ITTER_INCR;

    // quit if we want it done fast
    if (merror<=0)
      break;

    // solution is good enough => quit
    if ((0<=error)&&(error<merror))
      break;

  } while (itter<MAXITTER);

  delete[] id;

  totalerror=sqrt(((double)error)/(128*128/16*w*h/4));
#ifdef SHOW_STATS
  fprintf(stderr,"red: error is %ld after %d itterations\n",error,itter);
#endif
}


/* ================ GICBayer_Blue methods  ================ */

inline void blue_make_up_horz(short* tl, const short* gtl, int y, int w) {
  // green pixels between two blues horizontally
  for (int x=1; x<w-1; x+=2)
    tl[y*w+x] = interpolate(tl[y*w+x-1],tl[y*w+x+1],gtl[y*w+x],gtl[y*w+x-1],gtl[y*w+x+1]);
  // x=w-1
  tl[y*w+w-1] = tl[y*w+w-2];  // could be better, but it's on the right edge so who cares
}
inline void blue_make_up_vert(short* tl, const short* gtl, int y, int w, int h) {
  // green pixels between two blues vertically
  for (int x=0; x<w; x+=2) {
    if (y==0)
      tl[y*w+x] = tl[(y+1)*w+x];  // could be better, but it's on the top so who cares
    else
      tl[y*w+x] = interpolate(tl[y*w+x-w],tl[y*w+x+w],gtl[y*w+x],gtl[y*w+x-w],gtl[y*w+x+w]);
  }
}
inline void blue_make_up_red(short* tl, const short* gtl, int y, int w) {
  // red pixels (note: we only use red data indirectly through green)
  for (int x=1; x<w; x+=2) {
    short gmin = MIN(gtl[y*w+x],gtl[y*w+x+w-1]);
    short gmax = MAX(gtl[y*w+x],gtl[y*w+x+w-1]);
    short bmin = tl[y*w+x+w-1];
    short bmax = tl[y*w+x+w-1];
    if (x<w-1) {
      gmin = MIN(gmin,gtl[y*w+x+w+1]);
      gmax = MAX(gmax,gtl[y*w+x+w+1]);
      bmin = MIN(bmin,tl[y*w+x+w+1]);
      bmax = MAX(bmax,tl[y*w+x+w+1]);
    }
    if (y>0) {
      gmin = MIN(gmin,gtl[y*w+x-w-1]);
      gmax = MAX(gmax,gtl[y*w+x-w-1]);
      bmin = MIN(bmin,tl[y*w+x-w-1]);
      bmax = MAX(bmax,tl[y*w+x-w-1]);
      if (x<w-1) {
	gmin = MIN(gmin,gtl[y*w+x-w+1]);
	gmax = MAX(gmax,gtl[y*w+x-w+1]);
	bmin = MIN(bmin,tl[y*w+x-w+1]);
	bmax = MAX(bmax,tl[y*w+x-w+1]);
      }
    }
    
    int count=1;
    int pixel = extrapolate(bmin,bmax,tl[y*w+x+w-1],gmin,gmax,gtl[y*w+x],gtl[y*w+x+w-1]);
    if (x<w-1) {
      pixel += extrapolate(bmin,bmax,tl[y*w+x+w+1],gmin,gmax,gtl[y*w+x],gtl[y*w+x+w+1]);
      ++count;
    }
    if (y>0) {
      pixel += extrapolate(bmin,bmax,tl[y*w+x-w-1],gmin,gmax,gtl[y*w+x],gtl[y*w+x-w-1]);
      ++count;
      if (x<w-1) {
	pixel += extrapolate(bmin,bmax,tl[y*w+x-w+1],gmin,gmax,gtl[y*w+x],gtl[y*w+x-w+1]);
	++count;
      }
    }
    tl[y*w+x]=pixel/count;
  }
}
inline void blue_better_middle(short* tl, short* id, int y, int w) {
  id += (y/2)*w/2;
  tl += y*w;
  // left column
  tl[0] = boundcheck((8*id[0]-tl[-w]-tl[w])/6);
  // middle
  for (int x=2; x<w; x+=2)
    tl[x] = boundcheck((64*id[x/2] 
			-6*tl[-w+x] -6*tl[w+x]
			-6*tl[x-1]  -6*tl[x+1]
			-tl[-w+x-1] -tl[-w+x+1]
			-tl[w+x-1]  -tl[w+x+1])/36);
}
inline void blue_better_bottom(short* tl, short* id, int w, int h) {
  id += (h-2)*w/4;
  tl += (h-1)*w;
  // bottom row
  for (int x=2; x<w; x+=2) 
    tl[x] = boundcheck((8*id[x/2]-tl[x-1]-tl[x+1])/6);
}
inline long blue_error_middle(short* tl, short* id, int y, int w) {
  long error=0;
  id += (y/2)*w/2;
  tl += y*w;
  // left column
  int residule = (6*tl[0]+tl[-w]+tl[w])/8 - id[0];
  error += (residule*residule+8)>>4;
  
  // middle pixels
  for (int x=2; x<w; x+=2) {
    int residule = (36*tl[x]+
		    6*tl[-w+x] +6*tl[w+x]+
		    6*tl[x-1]  +6*tl[x+1]+
		    tl[-w+x-1] +tl[-w+x+1]+
		    tl[w+x-1]  +tl[w+x+1])/64 - id[x/2];
    error += (residule*residule+8)>>4;
  }
  return error;
}
inline long blue_error_bottom(short* tl, short* id, int w, int h) {
  long error=0;
  id += (h-2)*w/4;
  tl += (h-1)*w;
  // bottom row
  for (int x=2; x<w; x+=2) {
    int residule = (tl[x-1]+6*tl[x]+tl[x+1])/8 - id[x/2];
    error += (residule*residule+8)>>4;
  }
  return error;
}

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

  if ((green.getWidth()!=w)||(green.getHeight()!=h)) {
    fprintf(stderr,"GICBayer_Red: FATAL ERROR! green component has different dimensions\n");
    return;
  }

  const short* gtl = green.getPointArray();
  short* id = new short[w*h/4];

  int itter=0;
  int nitter=2*(maxerr==-1?MINITTER_FAST:MINITTER);
  long error;
  int merror=(maxerr==-1?-1:((int)(maxerror*maxerror*128*128*w*h/4))>>4);

  // main loop
  do {
    // do nitter/2 itterations in batch
    error=0;
    for (int head=0; head<h+nitter; head+=2) {
      for (int i=0; i<nitter; i+=2) {
	int y=head-i;
	
	// initialize arrays (first itteration only)
	if ((itter+i==0)&&(y<h))
	  for (int x=0; x<w; x+=2)
	    tl[y*w+w+x] = id[y*w/4+x/2] = (((short)imagedata[y*w+w+x])-128)<<7;
	
	// calculate made up points (guessed using green as guide)
	if ((y>=0)&&(y<h)) {
	  blue_make_up_horz(tl,gtl,y+1,w);
	  blue_make_up_vert(tl,gtl,y,w,h);
	  blue_make_up_red(tl,gtl,y,w);
	}

	y-=2;
	if (y>=0) {
	  // calculate error if necessary
	  if ((i==nitter-2)&&(merror>0)) {
	    if (y<h-2)
	      error += blue_error_middle(tl,id,y+1,w);
	    else if (y==h-2)
	      error += blue_error_bottom(tl,id,w,h);
	  }
	  
	  // calculate better answer
	  if (y<h-2)
	    blue_better_middle(tl,id,y+1,w);
	  else if (y==h-2)
	    blue_better_bottom(tl,id,w,h);
	}
      }

    }

    itter+=nitter/2;
    nitter=2*ITTER_INCR;

    // quit if we want it done fast
    if (merror<=0)
      break;

    // solution is good enough => quit
    if ((0<=error)&&(error<merror))
      break;

  } while (itter<MAXITTER);

  delete[] id;

  totalerror=sqrt(((double)error)/(128*128/16*w*h/4));
#ifdef SHOW_STATS
  fprintf(stderr,"blue: error is %ld after %d itterations\n",error,itter);
#endif
}