#include "ninox.h"

static void do_HNoiseFilter(struct Image *img);
static void do_HNoiseFilter_window(unsigned char *buffer, int w, int h, int depth, int x1, int y1, int x2, int y2);

// Scan the region given by (x1,y1) - (x2,y2) and return the barycentre (centre of brightness)

static int
_FindCentre_Barycentre(struct Image *img, int x1, int y1, int x2, int y2, int *x_avg, int *y_avg)
	{
	int depth = img->depth;
	int img_width = img->width;
	int x,y,bpp = img->depth/8;
	int count=0;	// count of significant pixels
	int x_total=0, y_total=0;
	int r,g,b,w;
	double l,RealThreshHold;

        if (depth == 24) {
                RealThreshHold = ThreshHold;    // 8 bit data
                }
        else if (depth == 16) {
                RealThreshHold = ThreshHold * 256;              // 16 bit data
                }
        else if (depth == 8)  {
                RealThreshHold = ThreshHold;    // 8 bit data
                }
        else {
                Print("FindCentre: Unsupported depth: %d\n",depth);
                return 0;
                }

        for(y=y1; y<=y2; ++y) {
	   int rowcount=0;

           if (bpp == 3) {
	   	unsigned char  *uptr = (unsigned char  *)img->data + (y * img_width + x1) * bpp;
		for (x=x1; x<=x2; ++x) {
                   b = *(uptr++); g = *(uptr++); r = *(uptr++);
                   l = 0.299 * (double)r + 0.587 * (double)g + 0.114 * (double) b;
                   if (l >= RealThreshHold) {
		      //w = l/RealThreshHold; x_total += x*w; y_total += y*w; count+=w; rowcount++; 
		      x_total += x; y_total += y; count++; rowcount++; 
		      }
                   }
		}
           else if (bpp == 2) {
	   	unsigned short *iptr = (unsigned short *)img->data + y * img_width + x1;
		for (x=x1; x<=x2; ++x,++iptr) {
                   if (*iptr >= RealThreshHold) {
		      //w = *iptr/RealThreshHold; x_total += x*w; y_total += y*w; count+=w; rowcount++; 
		      x_total += x; y_total += y; count++; rowcount++; 
		      }
                   }
		}
           else if (bpp == 1) {
	   	unsigned char  *uptr = (unsigned char  *)img->data + (y * img_width + x1) * bpp;
		for (x=x1; x<=x2; ++x,++uptr) {
                   if (*uptr >= RealThreshHold) {
		      //w = *uptr/RealThreshHold; x_total += x*w; y_total += y*w; count+=w; rowcount++; 
		      x_total += x; y_total += y; count++; rowcount++; 
		      }
                   }
		}

           if (ForceProcess == 0 && (y==y1 || y==y2) && rowcount > 100) {
                Print("Warning: image data found near edge (row %d, %d pixels). Procesing cancelled\n",
                        y,rowcount);
                return 0;
                }
           }

        if (count == 0 && ForceProcess==0) {
           Print("[no image] ");
           return(0);
           }

        if (count < MinPixels && ForceProcess==0) {
           if (!Quiet) Print("[Not enough pixels. Found %d, require %d] ",count,MinPixels);
           return(0);
           }

        *x_avg = x_total / count;
        *y_avg = y_total / count;

        if (!Quiet) Print("centre=(%d,%d), %d px\n",*x_avg,*y_avg,count);

	return(1);
	}

// Scan the image from top->bottom, left -> right and stop when we find the first (leftmost) object
static int
_FindCentre_Barycentre_LeftRight(int dir, struct Image *img, int x1, int y1, int x2, int y2, int *x_avg, int *y_avg)
	{
	int depth = img->depth;
	int img_width = img->width;
	int x,y,bpp = img->depth/8;
	int count=0;	// count of significant pixels
	int x_total=0, y_total=0;
	int r,g,b,w, xend,xinc=1;
	double l,RealThreshHold;
	int MinColPixels = 5;
	int state, in_object=0, left_object=0;

	*x_avg = *y_avg = 0;

        if (depth == 24) {
                RealThreshHold = ThreshHold;    // 8 bit data
                }
        else if (depth == 16) {
                RealThreshHold = ThreshHold * 256;              // 16 bit data
                }
        else if (depth == 8)  {
                RealThreshHold = ThreshHold;    // 8 bit data
                }
        else {
                Print("FindCentre: Unsupported depth: %d\n",depth);
                return 0;
		}

	if (dir == CENTRE_LEFT) {
	   // scan left -> right
	   x = x1;
	   xend = x2;
	   xinc = 1;
	   }
	else if (dir == CENTRE_RIGHT) {
	   // scan right -> left
	   x = x2;
	   xend = x1;
	   xinc = -1;
	   }
	else {
	   Print("_FindCentre_Barycentre_LeftRight: Unknown direction %d\n",dir);
	   return 0;
	   }

	state = 0; // 0 = left of object, 1 = on object, 2 = right of object
	for(; x != xend; x+=xinc) {
	   int colcount = 0;

           if (bpp == 3) {
	   	unsigned char  *uptr = (unsigned char  *)img->data + (y1 * img_width + x) * bpp;
		int rowbytes = img_width * 3;

		for (y=y1; y<=y2; ++y) {
                   b = *(uptr); g = *(uptr+1); r = *(uptr+2);
		   uptr += rowbytes;
                   l = 0.299 * (double)r + 0.587 * (double)g + 0.114 * (double) b;
                   if (l >= RealThreshHold) {
		      x_total += x; y_total += y; count++; colcount++; 
		      }
                   }
		}
           else if (bpp == 2) {
	   	unsigned short *iptr = (unsigned short *)img->data + y1 * img_width + x;

		for (y=y1; y<=y2; ++y,iptr += img_width) {
                   if (*iptr >= RealThreshHold) {
		      x_total += x; y_total += y; count++; colcount++; 
		      }
                   }
		}
           else if (bpp == 1) {
	   	unsigned char  *uptr = (unsigned char  *)img->data + y1 * img_width + x;

		for (y=y1; y<=y2; ++y,uptr += img_width) {
                   if (*uptr >= RealThreshHold) {
		      x_total += x; y_total += y; count++; colcount++; 
		      }
                   }
		}

           if (ForceProcess == 0 && (x==x1 || x==x2) && colcount > 50) {
                Print("Warning: image data found near edge (col %d, %d pixels). Procesing cancelled\n",
                        x,colcount);
                return 0;
                }

	   if (state==0) {
		// seeking from left, not yet acquired
	      	if (colcount >= MinColPixels) ++in_object; else in_object=0;
	      	if (in_object > 3) {state=1; left_object=0; }
	 	}
	   else if (state==1) {
		// on object
	      	if (colcount < MinColPixels) ++left_object; else left_object=0;
	      	if (left_object > 3) state=2;
		}
	   else if (state==2) {
		// now right of the object, stop scanning
		break;
		}
	   }

        *x_avg = x_total / count;
        *y_avg = y_total / count;

        if (!Quiet) Print("centre(%s)=(%d,%d), %d px\n",CentreMode==CENTRE_LEFT?"left":"right",*x_avg,*y_avg,count);

	return(1);
	}

static int
_FindCentre(struct Image *img, int x1, int y1, int x2, int y2, int *x_avg, int *y_avg)
	{

	if (CentreMode == CENTRE_BARYCENTRE) 
	   return _FindCentre_Barycentre(img,x1,y1,x2,y2,x_avg,y_avg);

	if (CentreMode == CENTRE_LEFT) 
	   return _FindCentre_Barycentre_LeftRight(CENTRE_LEFT,img,x1,y1,x2,y2,x_avg,y_avg);

	if (CentreMode == CENTRE_RIGHT) 
	   return _FindCentre_Barycentre_LeftRight(CENTRE_RIGHT,img,x1,y1,x2,y2,x_avg,y_avg);

	return 0;
	}


// find the centre of brightness of the whole image
int
FindCentre(struct Image *img, int *x_avg, int *y_avg)
	{
	int x1 = 2;
	int x2 = img->width - 3;
	int y1 = 0;
	int y2 = img->height-1;

	return _FindCentre(img,x1,y1,x2,y2,x_avg,y_avg);
	}

int
CalculateCutout(struct Image *img)
	{
	int x_avg,y_avg;
	int x1,y1,x2,y2;
	int err;

	// Locate the centre of the object in the frame
	if (img->sub_x && img->sub_y) 
	     err = _FindCentre(img,img->sub_x-CutX/2, img->sub_y - CutY/2, img->sub_x + CutX/2, img->sub_y + CutX/2, &x_avg, &y_avg);
	else err =  FindCentre(img,&x_avg,&y_avg);

	if (! err) return(0);

	// calculate the cutout from the source image (x1,y1) - (x2,y2).
	x1 = x_avg - CutX/2;
	x2 = x1 + CutX - 1;
	y1 = y_avg - CutY/2;
	y2 = y1 + CutY - 1;

//Print("Cutout (%d,%d) - (%d,%d)\n",x1,y1,x2,y2);

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

	return(1);
	}

// Cutout the rectangle (x1,y1)-(x2,y2) and move it to the centre of the destination.
// Assume no resizing, ie both *in and *out have width and height as given
int
CutOut(struct Image *img)
	{
	int width = img->width;
	int height = img->height;
	int bpp = img->depth/8;
	int xc=width/2;
	int yc=height/2;
	int y,i,cutx,cuty,yo,rowbytes;
	int X1,X2,Y1,Y2;  // output 
	int x1,y1,x2,y2;
	unsigned char *src,*dst,*ptr;
	unsigned short *uptr;

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

	src = img->data;
	dst = ZeroMalloc(width * height * bpp);
	if (!dst) {
          Print("Out of Memory!\n");
          return 0;
          }

	if (ZBlack) switch(bpp) {
	   case 16:
		uptr = (unsigned short *)dst;
		for(i=0; i<width*height; ++i) uptr[i]=ZBlack<<8;
		break;
	   }

        // Switch to the destination
        img->data = dst;

	// calculate new cutout, after possible clipping
	cutx = x2 - x1 + 1;
	cuty = y2 - y1 + 1;

	// Calculate where our cutout will end up so we can clip it
	// and adjust the src boundaries if required
	X1 = xc - cutx/2; X2 = X1 + cutx - 1;
	Y1 = yc - cuty/2; Y2 = Y1 + cuty - 1;

	// Clip the src rectangle
	if (x1<0) { X1 -= x1; x1=0; }
	if (y1<0) { Y1 -= y1; y1=0; }
	if (x2>=width)  {X2 -= (x2-width+1);  x2=width-1;}
	if (y2>=height) {Y2 -= (y2-height+1); y2=height-1;}

	// Clip the dst rectangle
	if (X1<0) {x1 -= X1; X1=0;}
	if (Y1<0) {y1 -= Y1; Y1=0;}
	if (X2>=width)  {x2 -= (X2-width+1);  X2=width-1;}
	if (Y2>=height) {y2 -= (Y2-height+1); Y2=height-1;}

	// Now (x1,y1)-(x2,y2) is guaranteed to fit in the destination
	// recalculate the cutout size in case the bounds have changed
	cutx = x2 - x1 + 1;
	cuty = y2 - y1 + 1;

	rowbytes = cutx * bpp;
	dst = img->data + (Y1 * width + X1) * bpp;
	ptr = src  + (y1 * width + x1) * bpp;

	for(y = y1; y<=y2; ++y) {
	   if (ChangeGamma || ChangeGain) gammacpy(dst,ptr,cutx,bpp,Gamma);
	   else memcpy(dst,ptr,rowbytes);
	   dst += width * bpp;
	   ptr += width * bpp;
	   }

	img->cutout.x1 = X1; img->cutout.x2 = X2;
	img->cutout.y1 = Y1; img->cutout.y2 = Y2;
	free(src);
	return(1);
	}

int
ApplyGamma(struct Image *img, double g)
	{
	int y;
	int width = img->width;
	int height = img->height;
	int bpp = img->depth/8;
	int rowbytes = width * bpp;
	int o;
	unsigned char *src = img->data;
	unsigned char *buf = Malloc(width * height * bpp);

	for(o=y=0; y<height; ++y, o+= rowbytes)
	   gammacpy(buf+o,src+o,width,bpp,Gamma);

	free(img->data);
	img->data = buf;
	return 1;
	}

// If we rescale then reset width and height
int process(struct Image *img)
	{
	int x1,x2,y1,y2;
	int i;

	if (feof(stdin)) {
	   Print("EOF, quitting\n");
	   exit(1);
	   }

	// Locate the object and create the cutout rectangle from (x1,y1) to (x2,y2)
	if (DoCutout && !CalculateCutout(img))
	      return(0);

	// Apply the Horizontal Noise filter
	if (HNoiseFilter) {
	   if (!Quiet) Print("Hnoise: ");
	   do_HNoiseFilter(img);
	   if (!Quiet) Print("\n");
	   }

	// Apply input filter if required
	if (InputFilter && !Quiet) Print("IFilter: ");
        for(i=0; i<InputFilter; ++i) {
           double pc = input_filter(img);
	   if (!Quiet) {
		if (i) Print("+");
		Print("%-2.0f%%",pc*100.0);
		}
	   }
	if (InputFilter && !Quiet) Print("\n");

	if (AMean)
	   AlphaTrimmedMean(img,1,2);

	if (StreamFilter)
	   if (! do_StreamFilter(img))
		return(0);

	// Cut out the source image and copy into destination buffer
	if (DoCutout) CutOut(img);
	else if (ChangeGamma || ChangeGain)
	   ApplyGamma(img,Gamma);

	// If we have specified a -subregion then we now cut this out and use it
	// instead of the whole centered image
	if (EnableSubRegion) {
	   if (SR_X1 >= img->width || SR_X2 >= img->width || SR_Y1 >= img->height ||
		SR_Y2 >= img->height) {
		Print("Error: Subregion (%d,%d,%d,%d) outside image dimensions of %d x %d\n",
			SR_X1,SR_Y1,SR_X2,SR_Y2,img->width, img->height);
		exit(1);
		}

	   // Now set the cutout region and recreate the image
	   img->cutout.x1 = SR_X1;
	   img->cutout.x2 = SR_X2;
	   img->cutout.y1 = SR_Y1;
	   img->cutout.y2 = SR_Y2;

	   CutX = newWidth  = SR_X2-SR_X1+1;
	   CutY = newHeight = SR_Y2-SR_Y1+1;

	   CutOut(img);
	   }

	if (do_Rotate) {
	   if (! RotateImage(img,RotateAngle)) {
		Print("processing cancelled\n");
		return 0;
		}
	   Print("Rotated by %lf\n",RotateAngle);
	   }

	// If we have to upscale the data then do it from 
	// bottom -> top so we can just copy the existing 
	// image data in place.
	if (UpScale > 1) {
	   upscale_image(img,UpScale);

	  if (UpScale_Smoothing && UpScale_Smoothing_When == SMOOTHING_BEFORE)
	     smooth_image(img,UpScale);
	  }

	if (DownScale > 1) {
	   downscale_image(img,DownScale);
	   }

	if (Do3x3Smooth) smooth_image(img,1);

        if (newWidth > 0)  img->dst_width = (newWidth * UpScale) / DownScale;
        if (newHeight > 0) img->dst_height = (newHeight * UpScale) / DownScale;

        // Width must be a multiple of 4 bytes
        img->dst_width -= img->dst_width & 3;

        if (img->dst_width < 0) img->dst_width = img->width;
        if (img->dst_height < 0) img->dst_height = img->height;

        if (img->dst_width > img->width) img->dst_width = img->width;
        if (img->dst_height > img->height) img->dst_height = img->height;

	if (Morphing) {
	   if (!morph_image(img,Morphing)) {
		Print("Morphing failed on %s\n",img->src_fname);
		return 0;
		}
	   }

        // Possibly merge data from a MergeFile
        if (MergeFile) {
           merge_data(MergeFile,img);
           }

	return(1);
	}

int
gammacpy(unsigned char *dst, unsigned char *src, int npix, int bpp, double gamma)
	{
	static double last_gamma = 1.0;
	static int last_bpp = -1;
	static unsigned char *table = NULL;
	static unsigned short *itable;
	unsigned short *idst = (unsigned short *) dst;
	unsigned short *isrc = (unsigned short *) src;

	// Do we have to init the lookup table?
	if (last_bpp != bpp || last_gamma != gamma || table == NULL) {
	   if (table) {
		Print("\nfree table\n");
		free(table);
		}

	   Print("\nCreating gamma table %f for %dbpp\n",gamma,bpp);

	   switch(bpp) {
		case 1:
		case 3:
		table = (unsigned char *) malloc(256); if (! table) {
		   Print("adjust_gamma: Out of memory\n");
		   exit(1);
		   }
		create_gamma_table_8bpp(table,gamma);
		break;

		case 2:
		table = (unsigned char *) malloc(65536 * sizeof(short)); if (! table) {
		   Print("adjust_gamma: Out of memory\n");
		   exit(1);
		   }
		itable = (unsigned short *) table;
		create_gamma_table_16bpp(itable,gamma);
		break;

		default:
		Print("adjust_gamma: unknown bpp %d\n",bpp);
		exit(1);
		}
	   }

	switch (bpp) {
	   	case 1:
		while(npix--) *(dst++) = table[*(src++)];
		break;

		case 2:
		while(npix--) *(idst++) = itable[*(isrc++)];
		break;

		case 3:
		while(npix--) {
		   *(dst++) = table[*(src++)];
		   *(dst++) = table[*(src++)];
		   *(dst++) = table[*(src++)];
		   }
		break;
		}

	last_gamma = gamma;
	last_bpp = bpp;
	return(1);
	}

int
create_gamma_table_8bpp(unsigned char *table, double gamma)
	{
	int i;

	for(i=0; i<256; ++i) {
	   double v = (double)i * Gain; if (v>255) v=255;
	   v = pow(v / 255.0, 1.0/gamma) * 255.0;
	   if (v>255) v=255;
	   table[i] = (int)v;
	   }

	return(1);
	}

int
create_gamma_table_16bpp(unsigned short *table, double gamma)
	{
	int i;

	for(i=0; i<65536; ++i) {
	   double v = (double)i * Gain; if (v>65535) v=65535;
	   v = pow(v / 65535.0, 1.0/gamma) * 65535.0;
	   if (v>65535) v=65535;
	   table[i] = (int)v;
	   }

	return(1);
	}

static void
do_HNoiseFilter(struct Image *img)
	{
	int window=50;
	int x1 = img->cutout.x1;
	int x2 = img->cutout.x2;
	int y1 = img->cutout.y1;
	int y2 = img->cutout.y2;

	// Process the image in slices to focus on transient noise
	while(x1 < x2) {
	   do_HNoiseFilter_window(img->data, img->width, img->height, img->depth, x1, y1, x1+window-1, y2);
	   x1 += window;
	   }
	}
	   
// Force (x1,y1) and (x2,y2) to be a rectangle within (0,0) - (w-1,h-1)
int Clip(int w, int h, int *x1, int *y1, int *x2, int *y2)
	{
	if (*x2 < *x1 || *y2 < *y1) {
	   Print("\n\nClipBouondaries: (%d,%d) - (%d,%d) invalid region\n",*x1,*y1,*x2,*y2);
	   exit(1);
	   }

	if (*x1<0) *x1=0; if (*x1>=w) *x1=w-1;
	if (*x2<0) *x2=0; if (*x2>=w) *x2=w-1;

	if (*y1<0) *y1=0; if (*y1>=h) *y1=h-1;
	if (*y2<0) *y2=0; if (*y2>=h) *y2=h-1;

	return(1);
	}

static void
do_HNoiseFilter_window(unsigned char *buffer, int w, int h, int depth, int x1, int y1, int x2, int y2)
	{
	int i,nrows,y,r,count;
	int pixels_per_row;
	int row,diff;
	unsigned short *ubuffer = (unsigned short *)buffer;
	double *row_avg,*diff_avg;
	int minval = 20;
	int uminval = minval*256;
	int threshold = 3;
	unsigned int *temp_buffer;
	unsigned int *modified_row;
	int useupper = 1;
	int uselower = 1;

	// start 1 row down since we take diff with previous row
	if (y1==0) y1=1;

	Clip(w,h,&x1,&y1,&x2,&y2);

//printf("Clip %d %d %d %d\n",x1,y1,x2,y2);

	pixels_per_row = x2-x1+1;
	nrows = y2-y1+1;
	row_avg = (double *)malloc(sizeof(double) * nrows);
	diff_avg = (double *)malloc(sizeof(double) * nrows);
	if (row_avg == NULL || diff_avg == NULL) {
	   Print("Out of memory\n");
	   exit(1);
	   }

	modified_row = (unsigned int *)malloc(sizeof(int) * nrows);
	temp_buffer = (unsigned int *)malloc(sizeof(int) * pixels_per_row * nrows);
	if (temp_buffer == NULL) {
	   Print("Out of memory\n");
	   exit(1);
	   }

	switch(depth) {
	   case 8:
		for(r=0,y=y1; y<=y2; ++y,++r) {
		   int o = y*w + x1;
		   row=diff=count=0;
		   for(i=0; i<pixels_per_row; ++i,++o)
			if (buffer[o] >= minval) {
			   int v = buffer[o];
			   row += v; diff += v;
			   diff -= (int)buffer[o-w];
			   count++;
			   }
		   if (count<10) row_avg[r]=-1;
		   else {
			row_avg[r] = (double)row / (double)count;
			diff_avg[r] = (double)diff / (double)count;
			diff_avg[r] /= row_avg[r]; diff_avg[r]*=100;
//printf("row %d: avg=%lf  diff=%lf\n",r,row_avg[r],diff_avg[r]);
			}
		   }
		break;

	   case 16:
		for(r=0,y=y1; y<=y2; ++y,++r) {
		   int o = y*w + x1;
		   row=diff=count=0;
		   for(i=0; i<pixels_per_row; ++i,++o)
			if (ubuffer[o] >= uminval) {
			   int v = (int)ubuffer[o]>>8;
			   row += v; diff += v;
			   diff -= (int)ubuffer[o-w]>>8;
			   count++;
			   }
		   if (count<10) row_avg[r]=-1;
		   else {
			row_avg[r] = (double)row / (double)count;
			diff_avg[r] = (double)diff / (double)count;
			diff_avg[r] /= row_avg[r]; diff_avg[r]*=100;
			}
		   }
		break;

	   default:
		Print("Bit depth not handled\n");
		exit(1);
		}

	// Any row that is significantly brighter than its previous and following rows is
	// suspicious
	count=0;
	for(i=0; i<nrows; ++i) modified_row[i]=0;

	for(r=1; r<nrows-1; ++r) 
	   if (row_avg[r-1]>0 && row_avg[r]>0 && row_avg[r+1]>0) {
	   	int d1 = diff_avg[r] - diff_avg[r-1];
	   	int d2 = diff_avg[r] - diff_avg[r+1];

	   	if ((uselower && d1<=-threshold && d2<=-threshold) 
		   || (useupper && d1>=threshold && d2>=threshold)) {
		  double br = ((row_avg[r]/row_avg[r-1]) + (row_avg[r]/row_avg[r+1]))/2;
		  int o_dst = r * pixels_per_row;
		  int o_src = (y1+r) * w + x1;

		  if (depth==8) for(i=0; i<pixels_per_row; ++i,++o_src,++o_dst)
		   	   temp_buffer[o_dst] = (
				(buffer[o_src] / br) * 4 + 
				buffer[o_src-w] + 
			   	buffer[o_src+w]) / 6;
		  else if (depth==16) for(i=0; i<pixels_per_row; ++i,++o_src,++o_dst)
		   	   temp_buffer[o_dst] = (
				(ubuffer[o_src] / br) * 4 + 
				ubuffer[o_src-w] + 
			   	ubuffer[o_src+w]) / 6;
		   else { Print("do_HNoiseFilterWindow: Depth %d not handled\n",depth); exit(1);}

		  modified_row[r] = 1;
		  ++count;
		  }
		}

	if (!Quiet) Print("%d ",count);

	// Copy back the modified rows
	for(r=1; r<nrows-1 && count>0; ++r) if (modified_row[r]) {
	   int o_src = r * pixels_per_row;
	   int o_dst = (y1+r) * w + x1;
	   switch(depth) {
		case 8:
	   		for(i=0; i<pixels_per_row; ++i,++o_src,++o_dst)
				buffer[o_dst] = temp_buffer[o_src];
			break;
		case 16:
	   		for(i=0; i<pixels_per_row; ++i,++o_src,++o_dst)
				ubuffer[o_dst] = temp_buffer[o_src];
			break;
		}
	   --count;
	   }

	free(temp_buffer);
	free(row_avg);
	free(diff_avg);
	free(modified_row);
	}

// Params are (min,max) (assumed 8 bit, we will scale correctly for 16 bit).
// each pixel is adjusted as follows:
// - subtract min, clamp to zero
// - multiply by 255/(max-min), clamp to MAXPIXEL
int
LevelsAdjust(struct Image *img, int min, int max)
	{
	unsigned char *data = (unsigned char *)img->data;
	unsigned short *udata = (unsigned short *)img->data;
	int o,depth = img->depth;
	int npix = img->width * img->height;
	double scale;
	double val;

	scale = 255.0 / ((double)max - (double)min);

	switch(depth) {
		case 8:
		for(o=0; o<npix; ++o) {
		   val = data[o]; val -= min; if (val<0) val=0;
		   val *= scale; if (val>255) val=255;
		   data[o] = val;
		   }
		break;
		case 16:
		for(o=0; o<npix; ++o) {
		   val = udata[o]; val -= (min<<8); if (val<0) val=0;
		   val *= scale; if (val>65535) val=65535;
		   udata[o] = val;
		   }
		break;
		default:
		printf("Levels not implemented for depth %d\n",depth);
		break;
	  	}

	return 1;
	}

#define A_MINMAX(n) \
	if ((n)>maxx) \
	   { max=maxx;maxx=(n); } \
	else if ((n)>max) \
	   max=(n); \
	if ((n) < minn) \
	   { min=minn;minn=(n); } \
	else if ((n)<min) min=(n)

int
AlphaTrimmedMean(struct Image *img, int radius, int cut)
	{
	int width = img->width;
	int height = img->height;
	int depth = img->depth;
	int bufsize = width * height * depth/8;
	unsigned char *dst,*src = img->data;
	unsigned short *idst,*isrc = (unsigned short *)img->data;
	static unsigned char *Buf = NULL;
	static int BufSize = 0;
	int x,y,o;
	unsigned int minn,min,max,maxx; // least 2 and max 2 values

	if (bufsize != BufSize) {
	   if (Buf != NULL) free(Buf);
	   Buf = Malloc(bufsize);
	   BufSize = bufsize;
	   }

	dst = Buf;
	idst = (unsigned short *)Buf;
	src = img->data;
	isrc = (unsigned short *)img->data;

	switch(depth) {
	   case 8:
		for(y=1; y<height-1; ++y) for(x=1,o=y*width+1; x<width-1; ++x) {
		   unsigned int n,v = src[o];
		   min=minn=255; max=maxx=0;

		   n=src[o-1];       v+=n; A_MINMAX(n);
		   n=src[o+1];       v+=n; A_MINMAX(n);
		   n=src[o-width-1]; v+=n; A_MINMAX(n);
		   n=src[o-width];   v+=n; A_MINMAX(n);
		   n=src[o-width+1]; v+=n; A_MINMAX(n);
		   n=src[o+width-1]; v+=n; A_MINMAX(n);
		   n=src[o+width];   v+=n; A_MINMAX(n);
		   n=src[o+width+1]; v+=n; A_MINMAX(n);

		   v -= min+minn+max+maxx;
		   dst[o++] = v/5;
		   }
		break;

	   case 16:
		for(y=1; y<height-1; ++y) for(x=1,o=y*width+1; x<width-1; ++x) {
		   unsigned int n,v = isrc[o];
		   min=minn=65535; max=maxx=0;

		   n=isrc[o-1];       v+=n; A_MINMAX(n);
		   n=isrc[o+1];       v+=n; A_MINMAX(n);
		   n=isrc[o-width-1]; v+=n; A_MINMAX(n);
		   n=isrc[o-width];   v+=n; A_MINMAX(n);
		   n=isrc[o-width+1]; v+=n; A_MINMAX(n);
		   n=isrc[o+width-1]; v+=n; A_MINMAX(n);
		   n=isrc[o+width];   v+=n; A_MINMAX(n);
		   n=isrc[o+width+1]; v+=n; A_MINMAX(n);

		   v -= minn+maxx;
		   //v -= minn+min+max+maxx;
		   idst[o++] = v/7;
		   }
		break;

	   default:
		Print("Alpha Trimmed Mean not supported for depth %d\n",depth);
		return 0;
		break;
	   }

	memcpy(img->data,Buf,BufSize);
	return 1;
	}