#include "ninox.h"

// Detect bright or dark pixels and replace them. We assume that
// most noise causes pixels to become brighter, so we err on the side
// of darkness when replacing a bright pixel.

static unsigned int filter_upixel(register unsigned char *p, register int w, int *changed);
static unsigned short filter_ipixel(register unsigned short *p, register int w, int *changed);
static unsigned int filter_cpixel(register unsigned char *p, register int w, unsigned char rgb[3]);

#define BLACK_RGB(d) (*(d)=0,*((d)+1)=0,*((d)+2)=0)
#define WHITE_RGB(d) (*(d)=255,*((d)+1)=255,*((d)+2)=255)

#define SET_RGB(d,s) (*(d)=*(s),*((d)+1)=*((s)+1),*((d)+2)=*((s)+2))
#define ADD_RGB(d,s) (*(d)+=*(s),*((d)+1)+=*((s)+1),*((d)+2)+=*((s)+2))
#define LUM_RGB(p)   (*(p)*0.299 + *((p)+1)*0.587 + *((p)+2)*0.114)

// divide the image into tiles measuring S x S pixels, within each tile replace the
// brightest pixel with the average of it's neighbors.
// Avoid processing too close to the edge of the window so that the neighbor calculation is easier
int
pop_filter(struct Image *img)
	{
	int w = img->width;
	int h = img->height;
	int d = img->depth;
	unsigned short *udata = (unsigned short *)img->data;
	unsigned char *data = (unsigned char *)img->data;
	int v,total,count,r1,r2,tmp;
	int o,oo,x,x1,y1,xx,y,yy,S,npix,max,omax;
	int *map;

	if (! PopFilter) return(1);

	S=PopFilter;
	npix = S * S;

	if (! Quiet) Print("POP filter = %d\n",S);

	// Allocate the x_offset and y_offset maps so that we scan each subregion
	// in a random way
	map = Malloc(sizeof(int) * npix);

	for(y=o=0; y<S; ++y) for(x=0; x<S; ++x,++o) {
	   map[o] = y * w + x;
	   }

	// shuffle the lookup map
	for(o=0; o<npix/2; ++o) {
	   r1 = Random(npix);
	   r2 = Random(npix);

	   tmp=map[r1]; map[r1]=map[r2]; map[r2]=tmp;
	   }

	switch(d) {
	   case 8:
		for(y=1+Random(S); y<h-S; y+=S) {
		   int offset = 1 + Random(S);
		   for(x=offset,o=y*w+offset; x<w-S; x+=S,o+=S) {
			int n,max=-1;
			for(n=0; n<npix; ++n) {
			   int O = o + map[n];
			   v=data[O];
			   if (v>max) {max=v; omax=O; }
			   }

			total =  data[omax-1];
			total += data[omax+1];
			total += data[omax-w-1];
			total += data[omax-w];
			total += data[omax-w+1];
			total += data[omax+w-1];
			total += data[omax+w];
			total += data[omax+w+1];
			data[omax] = total/8;

			// shuffle a bit more
	   		r1 = Random(npix);
	   		r2 = Random(npix);
	   		tmp=map[r1]; map[r1]=map[r2]; map[r2]=tmp;
			}
		   }
		break;
	   case 16:
		for(y=1+Random(S); y<h-S; y+=S) {
		   int offset = 1 + Random(S);
		   for(x=offset,o=y*w+offset; x<w-S; x+=S,o+=S) {
			int n,max=-1;
			for(n=0; n<npix; ++n) {
			  int O = o + map[n];
			  v=udata[O]; if (v>max) {max=v; omax=O;}
			  }

			total =  udata[omax-1];
			total += udata[omax+1];
			total += udata[omax-w-1];
			total += udata[omax-w];
			total += udata[omax-w+1];
			total += udata[omax+w-1];
			total += udata[omax+w];
			total += udata[omax+w+1];

			udata[omax] = total/8;

			// shuffle a bit more
	   		r1 = Random(npix);
	   		r2 = Random(npix);
	   		tmp=map[r1]; map[r1]=map[r2]; map[r2]=tmp;
			}
		   }
		break;
	   default:
		Print("pop filtering not supported on images of depth %d\n",d);
		return 0;
		break;
	   }
		
	free(map);
	return 1;
	}

static unsigned int
filter_upixel(register unsigned char *p, register int w, int *changed)
	{
	int min,max;
	int d,val,avg;
	double diff;
	int s1,s2,s3,s4,s5,s6,s7,s8,cp;

	p -= (w+1);
	s1 = *(p++); s2 = *(p++); s3 = *p; p += w;
	s4 = *(p--); cp = *(p--); s5 = *p; p+=w;
	s6 = *(p++); s7 = *(p++); s8 = *p;

	min=s1; max=min;
	if (s2 < min) min=s2; if (s2 > max) max=s2;
	if (s3 < min) min=s3; if (s3 > max) max=s3;
	if (s4 < min) min=s4; if (s4 > max) max=s4;
	if (s5 < min) min=s5; if (s5 > max) max=s5;
	if (s6 < min) min=s6; if (s6 > max) max=s6;
	if (s7 < min) min=s7; if (s7 > max) max=s7;
	if (s8 < min) min=s8; if (s8 > max) max=s8;

	// Calculate modified average of neighbors,
	// discard brightest and darkest pixels pixel
	avg = s1+s2+s3+s4+s5+s6+s7+s8 - max - min;
	avg /= 6;

	// Hot pixel
	if (cp-avg > 50) {
	   *changed=1;
	   return avg;				// light pixel
	   }

	if (avg>0) diff = (double)(cp - avg) / (double)avg;
	else 	   diff = cp - avg;

	if (diff < -InputFilter_ThreshHold) {
		*changed=1;
		return (cp + cp + avg) / 3;		// dark pixel
		}
	if (diff < InputFilter_ThreshHold) {
		*changed=0;
		return cp;			// ok pixel
		}

	*changed=1;
	return avg;				// light pixel
	}

// Data is RGB, 3 bytes per pixel
// "w" is bytes per line in the image
// THIS CODE DOESN'T WORK YET
static unsigned int
filter_cpixel(register unsigned char *p, register int w, unsigned char rgb[3])
	{
	int i,val;
	double diff;
	int cp[3],avg[3];
	double d,a,lum[9],lum_cp;

	p -= (w+3);
	SET_RGB(avg,p); lum[0]=LUM_RGB(p); p+=3;
	ADD_RGB(avg,p); lum[1]=LUM_RGB(p); p+=3;
	ADD_RGB(avg,p); lum[2]=LUM_RGB(p); p += w;
	ADD_RGB(avg,p); lum[3]=LUM_RGB(p); p-=3;

	SET_RGB(cp,p);  lum_cp=LUM_RGB(p); p-=3;

	ADD_RGB(avg,p); lum[4]=LUM_RGB(p); p+=w;
	ADD_RGB(avg,p); lum[5]=LUM_RGB(p); p+=3; 
	ADD_RGB(avg,p); lum[6]=LUM_RGB(p); p+=3;
	ADD_RGB(avg,p); lum[7]=LUM_RGB(p);

	// Calculate average
	for(i=a=0; i<8; ++i) a+=lum[i]; a/=8;

	// How different are we from our neighbors?
	d = lum_cp - a; 

	if (d<0) d=-d;
	diff = d / a;

	if (diff >= InputFilter_ThreshHold) {
		*rgb = avg[0]/8; *(rgb+1)=avg[1]/8; *(rgb+2)=avg[2]/8;
		return(1);
		}

	return(0);
	}

static unsigned short
filter_ipixel(register unsigned short *p, register int w, int *changed)
	{
	int min,max;
	int val,avg;
	double diff;
	int s1,s2,s3,s4,s5,s6,s7,s8,cp;

	p -= (w+1);
	s1 = *(p++); s2 = *(p++); s3 = *p; p += w;
	s4 = *(p--); cp = *(p--); s5 = *p; p+=w;
	s6 = *(p++); s7 = *(p++); s8 = *p;

	min=s1; max=min;
	if (s2 < min) min=s2; if (s2 > max) max=s2;
	if (s3 < min) min=s3; if (s3 > max) max=s3;
	if (s4 < min) min=s4; if (s4 > max) max=s4;
	if (s5 < min) min=s5; if (s5 > max) max=s5;
	if (s6 < min) min=s6; if (s6 > max) max=s6;
	if (s7 < min) min=s7; if (s7 > max) max=s7;
	if (s8 < min) min=s8; if (s8 > max) max=s8;

	// Calculate modified average of neighbors,
	// discard brightest and darkest pixels
	avg = s1+s2+s3+s4+s5+s6+s7+s8 - max - min;
	avg /= 6;

	// Noisy pixel, significantly different from the average of its neighbours
	// is replaced by the average of the neighbours.
	if (cp-avg > 50 * 256) {
	   *changed=1;
	   return avg;
	   }


	if (avg==0) avg=1;
	diff = (double)(cp - avg) / (double)avg;

	if (diff < -InputFilter_ThreshHold) {
		// pixel significantly darker than neighbors, replace with average of
		// this pixel and neighbours. Weights toward darker pixels being ok, noise
		// tends to only add to pixels so darker pixels are probably more "correct".
		// This will change the contrast of the image slightly.
		*changed=1;
		return (cp + avg) / 2;		// dark pixel
		}
	if (diff > InputFilter_ThreshHold) {
		// Pixel is significantly lighter than neighbours, replace
		// with the average of the neighbours.
		*changed=1;
		return avg;
		}

	// -InputFilter_ThreshHold <= 0 >= InputFilter_ThreshHold is ok
	// to leave unchanged
	*changed=0;
	return cp;
	}

// Filter the source image, limit to rectangle given by img->cutout
// Change the source image.
double
input_filter(struct Image *img)
	{
	int width = img->width;
	int height = img->height;
	int depth = img->depth;
	int i,x,y,o,os,od,filter_count=0,threshhold,changed;
	int filter_maybe=0;
	unsigned short *isrc,*idst;
	unsigned char *usrc,*udst;
	int w,h;
	unsigned char *tmpbuf = NULL;
	static struct Image *itmp = NULL;
	int x1,y1,x2,y2;

	x1 = img->cutout.x1; x2 = img->cutout.x2;
	y1 = img->cutout.y1; y2 = img->cutout.y2;

	// clip the bounds
	if (x1<0) x1=0; if (x1>=width) x1=width-1;
	if (x2<0) x2=0; if (x2>=width) x2=width-1;
	if (y1<0) y1=0; if (y1>=height-1) y1=height-1;
	if (y2<0) y2=0; if (y2>=height-1) y2=height-1;

	w = x2 - x1 + 1;
	h = y2 - y1 + 1;

// printf("(%d,%d) - (%d,%d) w=%d h=%d\n",x1,y1,x2,y2,w,h); fflush(stdout);

	tmpbuf = (unsigned char *)malloc(w * h * depth/8);
	if (tmpbuf==NULL) {
 	   Print("input_filter: out of memory\n");
	   return 0;
	   }

	switch(depth) {
	   case 8:
			usrc = (unsigned char *) img->data;
			udst = (unsigned char *) tmpbuf;

			// set the threshhold so we only consider
			// pixels greater than 5% in brightness
			threshhold = 256 / 20;

			if (DisplayFrames && InputFilter_Display) {
			   itmp = ConvertImage(img,IMG_FIT,16);
			   itmp->cutout.x1 = x1; itmp->cutout.y1 = y1;
			   itmp->cutout.x2 = x2; itmp->cutout.y2 = y2;
			   }
			   
			// copy row 0 with no filtering
			o = y1 * width + x1;
			memcpy(udst,usrc+o,w);

			// filter rows 1 .. height-2
			for(y=y1+1; y<y2; ++y) {
			  os = y * width + x1 + 1;  // offset into source
			  od = (y-y1) * w;  // offset into destination

			  udst[od] = usrc[os]; ++od; ++os;
			  for(x=x1+1; x<x2; ++x,++os,++od)
				if (usrc[os] < threshhold || usrc[os-1] < threshhold || usrc[os+1]<threshhold) udst[od]=usrc[os];
				else {
				   udst[od] = filter_upixel(usrc+os,width,&changed);
				   filter_count += changed;
				   filter_maybe++;
				   if (changed && DisplayFrames && InputFilter_Display) {
					*((unsigned short *)itmp->data + os) = 65535;
					}
				   }
			  udst[od] = usrc[os];
			  }
			
			if (DisplayFrames && InputFilter_Display) {
			   ShowImage(itmp);
			   DestroyImage(itmp);
			   }

			// copy last row with no filtering
			os = y * width + x1;  // offset into source
			od = (y-y1) * w;  // offset into destination
			memcpy(udst+od,usrc+os,w);

			// copy it back
			os=0;
			od = y1 * width + x1;
			for(y=0; y<h; ++y,os+=w,od+=width) {
			   memcpy(usrc+od,udst+os,w);
			   }

			break;

	   case 16:
			isrc = (unsigned short *) img->data;
			idst = (unsigned short *) tmpbuf;

			if (DisplayFrames && InputFilter_Display) {
			   itmp = ConvertImage(img,IMG_FIT,16);
			   itmp->cutout.x1 = x1; itmp->cutout.y1 = y1;
			   itmp->cutout.x2 = x2; itmp->cutout.y2 = y2;
			   }
			   
			// set the threshhold so we only consider
			// pixels greater than 5% in brightness
			threshhold = 65535 / 20;

			// copy row 0 with no filtering
			o = y1 * width + x1;
			memcpy(idst,isrc+o,w * 2);

			// filter rows 1 .. height-2
			for(y=y1+1; y<y2; ++y) {
			  os = y * width + x1 + 1;  // offset into source
			  od = (y-y1) * w;  // offset into destination

			  idst[od++] = isrc[os++]; 
			  for(x=x1+1; x<x2; ++x,++os,++od)
				if (isrc[os] < threshhold || isrc[os-1] < threshhold || isrc[os+1]<threshhold) idst[od]=isrc[os];
				else {
				   idst[od] = filter_ipixel(isrc+os,width,&changed);
				   filter_count += changed;
				   filter_maybe++;
				   if (changed && DisplayFrames && InputFilter_Display) {
					*((unsigned short *)itmp->data + os) = 65535;
					}
				   }
			  idst[od] = isrc[os];
			  }
			
			if (DisplayFrames && InputFilter_Display) {
			   ShowImage(itmp);
			   DestroyImage(itmp);
			   }

			os = y * width + x1;  // offset into source
			od = (y-y1) * w;  // offset into destination
			// copy last row with no filtering
			memcpy(idst+od,isrc+os,w * 2);

			// copy it back
			os=0;
			od = y1 * width + x1;
			for(y=0; y<h; ++y,os+=w,od+=width) {
			   memcpy(isrc+od,idst+os,w*2);
			   }

			break;

	   case 24:
			usrc = (unsigned char *) img->data;
			udst = (unsigned char *) tmpbuf;

			// set the threshhold so we only consider
			// pixels greater than 5% in brightness
			threshhold = 256 / 20;

			// copy row 0 with no filtering
			o = (y1 * width + x1) * 3;
			memcpy(udst,usrc+o,w * 3);

			// filter rows 1 .. height-2
			for(y=y1+1; y<y2; ++y) {
			  os = (y * width + x1) * 3;  // offset into source
			  od = (y-y1) * w * 3;  // offset into destination

			  SET_RGB(udst+od,usrc+os); od+=3; os+=3;
			  for(x=x1+1; x<x2; ++x,os+=3,od+=3)
				if (LUM_RGB(usrc+os) < threshhold) SET_RGB(udst+od,usrc+os);
				else {
				   unsigned char p[3];
				   ++filter_maybe;
				   if (filter_cpixel(usrc+os,width*3,p)) {
				      SET_RGB(udst+od,p); ++filter_count;
				      }
				   else {
				      SET_RGB(udst+od,usrc+os);
				      }
				   }
			  SET_RGB(udst+od,usrc+os);
			  }
			
			// copy last row with no filtering
			os = (y * width + x1)*3;  // offset into source
			od = (y-y1) * w * 3;  // offset into destination
			memcpy(udst+od,usrc+os,w * 3);

			// copy it back
			os=0;
			od = (y1 * width + x1) * 3;
			for(y=0; y<h; ++y,os+=w*3,od+=width*3) {
			   memcpy(usrc+od,udst+os,w*3);
			   }

			break;
		}

	free(tmpbuf);
	if (! filter_maybe) return 0.0;
	return (double)filter_count / (double)filter_maybe;
	}

// Generate the straight 3x3 average, but weight the centre pixel
static unsigned short
average_small_ipixel(register unsigned short *p, register int w)
	{
	int s1,s2,s3,s4,s5,s6,s7,s8,cp;

	p -= (w+1);
	s1 = *(p++); s2 = *(p++); s3 = *p; p += w;
	s4 = *(p--); cp = *(p--); s5 = *p; p+=w;
	s6 = *(p++); s7 = *(p++); s8 = *p;

	return (s1+s2+s3+s4+s5+s6+s7+s8+cp+cp) / 10;
	}

// p points to top-left corner of a large pixel that is sXs in size.
// buf is a buffer allocated by our caller that we can use for storing our
// intermediate results
// rowlen is the number of pixels in each row.
// s is the scale factor.

static void
smooth_large_ipixel(unsigned short *dst,unsigned short *src, int rowlen, int s)
	{
	int n,i,first=1;
	unsigned short *src_top,*src_bottom,*src_right;
	unsigned short *dst_bottom,*dst_right;
	unsigned short pixbuf[100];

	src_top = src;

	while(s>1) {
	   src_bottom = src_top + rowlen * (s-1);
	   src_right = src_top + s - 1;

	   dst_bottom = dst + rowlen * (s-1);
	   dst_right = dst + s - 1;

	   // horizontal pass (top and bottom)
	   for(i=n=0; i<s; ++i) {
		pixbuf[n++] = average_small_ipixel(src_top+i,rowlen);
		pixbuf[n++] = average_small_ipixel(src_bottom+i,rowlen);
		}

	   // vertical pass (left and right)
	   for(i=1; i<s-1; ++i) {
		pixbuf[n++] = average_small_ipixel(src_top+rowlen*i,rowlen);
		pixbuf[n++] = average_small_ipixel(src_right + rowlen*i,rowlen);
		}

	   // Copy the pixels to the destination
	   for(i=n=0; i<s; ++i) {
		*(dst+i) = pixbuf[n++];
		*(dst_bottom+i) = pixbuf[n++];
		}
	   for(i=1; i<s-1; ++i) {
		*(dst+rowlen*i) = pixbuf[n++];
		*(dst_right+rowlen*i) = pixbuf[n++];
		}

	   // Walk inward toward centre of large pixel (down and right) and repeat
	   s-=2; dst += rowlen+1;

	   // subsequent passes use dst as the source so we get a real blending
	   src_top = dst;
	   }
	}

// Input is an image that has just been upscaled by a factor of "scale".
// Each original pixel has been replaced by a scaleXscale square block of pixels
// Now we work out how to smooth this block to reduce edge artifacts
// Input width and height refer to the upscaled image, not the original image

int
smooth_image(struct Image *img, int scale)
	{
	int width = img->width;
	int height = img->height;
	int depth = img->depth;

	if (scale<2) return(0);

	switch(depth) {
	   case 8:
		Print("smoothing not implemented for depth 8\n");
		return 0;
		break;
	   case 16:
		return smooth_image_16(img->data,width,height,scale);
		break;
	   case 24:
		Print("smoothing not implemented for depth 24\n");
		return 0;
		break;
	   default:
		Print("smoothing not implemented for unknown depth %d\n",depth);
		return 0;
		break;
	   }
	return(0);
	}

int
smooth_image_16(unsigned char *data,int width,int height,int scale)
	{
	int x,y,buffer_size,large_width,large_height;
	int cutx = CutX,cuty = CutY;
	static unsigned char *buffer = NULL;
	static int old_buffer_size = -1;
	unsigned short *isrc,*ibuf;

	if (! DoCutout) {
	   cutx = width;
	   cuty = height;
	   }

	buffer_size = width * height * 2;

	if (buffer == NULL || old_buffer_size != buffer_size) {
	   if (buffer) free(buffer);
	   buffer = (unsigned char *)malloc(buffer_size);
	   if (buffer==NULL) {
		Print("Out of memory allocating buffer in smooth_image\n");
		exit(1);
		}
	   old_buffer_size = buffer_size;
	   }


	// Image dimensions in "large pixels" 
	large_height = height / scale;
	large_width = width / scale;

	// For ODD scaling values we use a simple 3x3 mean function
	// working clockwise around each large pixel, from edge to centre.
	// Leave the centre pixel alone.
	if (scale & 1) {
	   int o,x1,x2,y1,y2;
	   int v_inc = width * (scale-1);

	   // Start at top left corner of large pixel at (1,1)

	   y1 = (large_height - cuty)/2 - 1; if (y1<0) y1=0;
	   y2 = y1 + cuty+2; if (y2>large_height) y2=large_height;
	   x1 = (large_width - cutx)/2 -1; if (x1<0) x1=0;
	   x2 = x1 + cutx+2; if (x2>large_width) x2=large_width;

	   y1 *= scale; y2 *= scale;
	   x1 *= scale; x2 *= scale;

	   // Copy the src image into the dst to help in our
	   // smooth_large_ipixel routine
	   for(y=y1; y<y2; ++y) {
	     o = y * width;
	     isrc = (unsigned short *)data + o;
	     ibuf = (unsigned short *)buffer + o;
	     for(x=x1; x<x2; ++x) ibuf[x] = isrc[x];
	     }

	   y1 = (large_height - cuty)/2; if (y1<1) y1=1;
	   y2 = y1 + cuty; if (y2>large_height-2) y2=large_height-2;
	   x1 = (large_width - cutx)/2; if (x1<1) x1=1;
	   x2 = x1 + cutx; if (x2>large_width-2) x2=large_width-2;

	   for(y = y1; y <= y2; ++y) {
	     o = width * scale * y;
	     isrc = (unsigned short *)data + o;
	     ibuf = (unsigned short *)buffer + o;
	     for(x=x1; x<=x2; ++x)
		smooth_large_ipixel(isrc+x*scale,ibuf+x*scale,width,scale);
	     }

	   }
	else if (scale==2 || scale == 4) {
	  int o;
	  isrc = (unsigned short *)data;
	  ibuf = (unsigned short *)buffer;

	  memcpy(ibuf,isrc,width*height*2);
	   for(y=1; y<height-1; ++y) 
		for(x=1,o=y*width+1; x<width-1; ++x, ++o)
		   isrc[o] = average_small_ipixel(ibuf+o,width);
	   }
	else {
	   Print("Scale factor %d not supported\n",scale);
	   exit(1);
	   }

	return(0);
	}

// Normal Average
#define AVERAGE(a,b) (((a)+(b))>>1)

// weighted average, favours the dark side by 25%
#define DARK_AVERAGE(a,b) ((a)>(b)?((a)*3+(b)*5)>>3:((b)*3+(a)*5)>>3)

// iN-place scale existing widthxheight image by factor n.
// Work from bottom->top right->left duplicating pixel data.
// width & height refer to the src image dimensions.
int
upscale_image(struct Image *img, int n)
	{
	unsigned char *s,*src = img->data;
	int width = img->width;
	int height = img->height;
	int bpp = img->depth/8;
	int i,j,k,l,x,y;
	unsigned char *dst;
	unsigned short *isrc,*idst;
	int rowbytes = width * bpp;
	int imagebytes = rowbytes * height;

	// short-cut, nothing to do
	if (n<=1) return 1;

	// Allocate new image buffer for upscaled data
	dst = malloc(imagebytes * n * n);
	if (dst==NULL) {
	   Print("Out of Memory\n");
	   return 0;
	   }

	img->data = dst;

	switch(bpp) {
	   case 1:
		s = src;
		for(y=0; y<height; ++y,s+=rowbytes)
	   	   // Duplicate each input row n times, reset src position
	   	   // each time to end of current input row
	   	   for(j=0; j<n; ++j,s-=rowbytes)
	      		// copy a row with expansion
	      		for(x=0; x<rowbytes; ++x,++s)
			   for(i=0; i<n; ++i) *(dst++) = *s;
		break;

	   case 2:
		isrc = (unsigned short *) src;
		idst = (unsigned short *) dst;

		for(y=0; y<height; ++y, isrc+=width) {
	   	   // Duplicate each input row n times, reset src position
	   	   // each time to end of current input row

	   	   for(j=0; j<n; ++j,isrc-=width)
			// copy 1 pixel in row
	      		for(x=0; x<width; ++x,++isrc)
			   for(i=0; i<n; ++i) *(idst++) = *isrc;
		   }

		break;
	   case 3:
		s = src;
		for(y=0; y<height; ++y,s+=rowbytes)
	   	   // Duplicate each input row n times, reset src position
	   	   // each time to end of current input row
	   	   for(j=0; j<n; ++j,s-=rowbytes)
	      		// copy a row with expansion
	      		for(x=0; x<rowbytes; x+=bpp,s+=bpp)
			   for(i=0; i<n; ++i,dst+=bpp) {
				*dst = *s; *(dst+1) = *(s+1); *(dst+2) = *(s+2); 
				}
		break;
	   }

	free(src);
	img->width  *= n; img->height *= n;
	img->dst_width  *= n; img->dst_height *= n;

        img->cutout.x1 *= n; img->cutout.x2 *= n;
        img->cutout.y1 *= n; img->cutout.y2 *= n;

	return(0);
	}

// In-place downscale the src image, write it back into
int
downscale_image(struct Image *img, int n)
	{
	unsigned char *data = img->data;
	int width = img->width;
	int height = img->height;
	int bpp = img->depth/8;
	int i,j,o,x,y;
	unsigned char *src,*dst;
	unsigned short *isrc,*idst;
	int rowbytes = width * bpp;
	int nw = width/n;
	int nh = height/n;
	unsigned int *rowbuffer_i = NULL;
	unsigned int *rowbuffer_r = NULL;
	unsigned int *rowbuffer_b = NULL;
	unsigned int *rowbuffer_g = NULL;
	int total,total_r,total_g,total_b;
	int sample_size = n * n;

	src = data;   // first pixel in input image

	// We are shrinking the image, so work left->right and top->bottom
	// so we don't overwrite needed data

	switch(bpp) {
	   case 1:
		// Monochrome 8bpp

		rowbuffer_i = (unsigned int *)ZeroMalloc(nw * sizeof(unsigned int));
		dst = (unsigned char *)data;

		for(y=Random(n/2); y<height-n; y+=n) {
		   for(i=y; i<y+n; ++i) {
		      x = Random(n/2);	// start each row in a different place to remove artifacts
		      src = (unsigned char *)data + i*width +x;

		      for(o=0; x<width-n; x+=n,++o)
			for(j=total=0; j<n; ++j) rowbuffer_i[o] += *(src++);
		      }

		   // Output averages and reset the row buffer
		   for(x=0; x<nw; ++x) {
			*(dst++) = rowbuffer_i[x] / sample_size;
			rowbuffer_i[x]=0;
			}
		   }

		free(rowbuffer_i);
		break;

	   case 3:
		// RGB 24bpp
		rowbuffer_r = (unsigned int *)ZeroMalloc(nw * sizeof(unsigned int));
		rowbuffer_g = (unsigned int *)ZeroMalloc(nw * sizeof(unsigned int));
		rowbuffer_b = (unsigned int *)ZeroMalloc(nw * sizeof(unsigned int));
		dst = (unsigned char *)data;

		for(y=Random(n/2); y<height-n; y+=n) {
		   for(i=y; i<y+n; ++i) {
		      x = Random(n/2);	// start each row in a different place to remove artifacts
		      src = (unsigned char *)data + (i*width+x) * 3;

		      for(o=0; x<width-n; x+=n, ++o)
			for(j=0; j<n; ++j) {
			   rowbuffer_r[o] += *(src++);
			   rowbuffer_g[o] += *(src++);
			   rowbuffer_b[o] += *(src++);
			   }
		      }

		   // Output averages
		   for(x=0; x<nw; ++x) {
			*(dst++) = rowbuffer_r[x] / sample_size; rowbuffer_r[x]=0;
			*(dst++) = rowbuffer_g[x] / sample_size; rowbuffer_g[x]=0;
			*(dst++) = rowbuffer_b[x] / sample_size; rowbuffer_b[x]=0;
			}
		   }

		free(rowbuffer_r);
		free(rowbuffer_g);
		free(rowbuffer_b);
		break;

	   case 2:
		// Monochrome 16bpp

		rowbuffer_i = (unsigned int *)ZeroMalloc(nw * sizeof(unsigned int));
		idst = (unsigned short *)data;

		for(y=Random(n/2); y<height-n; y+=n) {
		   for(i=y; i<y+n; ++i) {
		      x = Random(n/2);	// start each row in a different place to remove artifacts
		      isrc = (unsigned short *)data + i*width +x;

		      for(o=0; x<width-n; x+=n,++o)
			for(j=total=0; j<n; ++j) rowbuffer_i[o] += *(isrc++);
		      }

		   // Output averages and reset the row buffer
		   for(x=0; x<nw; ++x) {
			*(idst++) = rowbuffer_i[x] / sample_size;
			rowbuffer_i[x]=0;
			}
		   }

		free(rowbuffer_i);
		break;
	      }

	img->width = nw; img->height = nh;
	img->dst_width = nw; img->dst_height = nh;

        img->cutout.x1 /= n; img->cutout.x2 /= n;
        img->cutout.y1 /= n; img->cutout.y2 /= n;

	return(1);
	}