#include "ninox.h"
#include <ctype.h>
#include <string.h>

/*
 * Generate a quality estimate for the given region by computing gradients
 * on downsampled images
 */

#define MAX_FILES 100000

typedef struct
	{
	int active;
	int bucket;  // 0..6
	char *fname;
	char *newname;
	double val;
	} Sharp;

static double AnalyseImage16(unsigned short *buf, int x_samples,int y_samples);

static double Histo_Quality(struct Image *img);
static double PixDiff(unsigned char *p1, unsigned char *p2);
inline static double PixDiffLC1(unsigned char *p1, unsigned char *p2, int depth, double l, double c);
static int Compare(const void *, const void *);
static unsigned short * _smooth_image_16(unsigned short *buf, int width, int height);
static void _Trim(int before, int after);

double cutoff = 0.5;
static Sharp *flist = NULL;
static int flist_size = 0;
static int flist_idx = 0;

static int Sorted = 0;

#define MAX_QREGIONS 8
struct QRegion QR[MAX_QREGIONS];
static int nQR = 0;

int
QAddRegion(int x1, int y1, int x2, int y2, int h_low, int h_high)
	{
	if (nQR == MAX_QREGIONS) {
	   printf("Maximum of %d Regions reached\n",MAX_QREGIONS);
	   exit(1);
	   }

	printf("Defined QRegion %d (%d,%d) - (%d,%d) : [ %d - %d]\n",nQR,
		x1,y1,x2,y2,h_low,h_high);
	fflush(stdout);

	QR[nQR].r.x1 = x1;
	QR[nQR].r.y1 = y1;
	QR[nQR].r.x2 = x2;
	QR[nQR].r.y2 = y2;
	QR[nQR].histo_low = h_low;
	QR[nQR].histo_high = h_high;
	nQR++;
	}

// hist[256]
// convert all formats to 8 bit for histogram

static int
MakeHistogram(struct Image *img, struct Region r, int *hist)
	{
	int i,o,x1,x2,y1,y2,x,y,npix;
	unsigned char *data;
	unsigned short *udata;

	for(i=0; i<256; ++i) hist[i]=0;
	npix=0;

	x1 = r.x1; y1 = r.y1;
	x2 = r.x2; y2 = r.y2;

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

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

		for(y=y1; y<=y2; ++y) {
		   o = y * img->width + x1;
		   for(x=x1; x<=x2; ++x,++o) { hist[data[o]]++; npix++;}
		   }

		break;
	   case 16:
		udata = (unsigned short *)img->data;

		for(y=y1; y<=y2; ++y) {
		   o = y * img->width + x1;
		   for(x=x1; x<=x2; ++x,++o) {
			hist[udata[o]>>8]++; npix++;
			//udata[o] += 30<<8; // AEW
			}
		   }

		break;
	   case 24:
		data = (unsigned char *)img->data;
		for(y=y1; y<=y2; ++y) {
		   o = (y * img->width + x1) * 3;
		   for(x=x1; x<=x2; ++x,o+=3) {
		   	int b = data[o];
		   	int g = data[o+1];
		   	int r = data[o+2];

		   	double v = r * 0.299 + g * 0.547 + b * 0.114;
		   	hist[(int)v]++;
			npix++;
			}
		   }
		break;
	   }

//printf("%d\n",npix); fflush(stdout);
	return npix;
	}

static double
Histo_Quality_Region(struct Image *img, struct QRegion R)
	{
	int i,npix,total;
	int hist[256];
	FILE *z;
	static int first=1;

	if (first) {unlink(QHISTO_FILE); first=0;}

	if (! (z = OpenLogfile(QHISTO_FILE,"a"))) {
	   printf("Cannot open %s\n",QHISTO_FILE);
	   exit(1);
	   }

	fprintf(z,"%s (%d,%d)-(%d,%d) ",img->src_fname,R.r.x1,R.r.y1,R.r.x2,R.r.y2);

	npix = MakeHistogram(img,R.r,hist);
	for(i=R.histo_low,total=0; i<=R.histo_high; ++i) {
	   fprintf(z,"%-3d:%-4d ",i,hist[i]);
	   total += hist[i];
	   }

	fprintf(z," [%d]\n",total);
	fclose(z);
	return (double)total;
	}

// Calculate the histogram quality for the given image
static double
Histo_Quality(struct Image *img)
	{
	int i,npix,thresh = ThreshHold;
	int p,min,max,total;

	if (QHISTO_MIN<0) QHISTO_MIN = ThreshHold;
	if (QHISTO_MAX<0) QHISTO_MAX = ThreshHold + 20;
	if (QHISTO_MAX > 255) QHISTO_MAX = 255;

	// If no regions defined then create a default region that covers the whole image
	if (nQR==0)
	   QAddRegion(0,0,img->width-1,img->height-1,QHISTO_MIN,QHISTO_MAX);

	for(i=total=0; i<nQR; ++i)
	   total += Histo_Quality_Region(img,QR[i]);

	return (double) total;
	}

// Set a threshhold value v (8 bit) and add together N-v for all 
// pixels of brightness N where N > v.
// This is good for small objects like moons of Jupiter
static double
BrThresh_Quality(struct Image *img)
	{
	unsigned char *data = (unsigned char *)img->data;
	unsigned short *udata = (unsigned short *)img->data;
	int thresh = ThreshHold;
	int x,y,x1,y1,x2,y2,o;
	int npix=0;
	double total=0;

	// Ignore pixels close to the edge
	y1=2; y2 = img->height-2;
	x1=2; x2 = img->width-2;

	switch(img->depth) {
	   case 8:
		for(y=y1; y<=y2; ++y) {
		   o = y * img->width + x1;
		   for(x=x1; x<=x2; ++x,++o) {
			if (data[o]>thresh) {
			   total += (data[o]-thresh) * (data[o]-thresh);
			   ++npix;
			   }
			}
		   }

		break;
	   case 16:
		thresh <<= 8;

		for(y=y1; y<=y2; ++y) {
		   o = y * img->width + x1;
		   for(x=x1; x<=x2; ++x,++o) {
			if (udata[o] >= thresh) {
			   total += (udata[o]-thresh)*(udata[o]-thresh);
			   ++npix;
			   }
			}
		   }

		total /= 256;
		break;
	   case 24:
		for(y=y1; y<=y2; ++y) {
		   o = (y * img->width + x1) * 3;
		   for(x=x1; x<=x2; ++x,o+=3) {
		   	int b = data[o];
		   	int g = data[o+1];
		   	int r = data[o+2];

		   	int v = r * 0.299 + g * 0.547 + b * 0.114;
			if (v >= thresh) {
		   	   total += (v-thresh) * (v-thresh);
			   ++npix;
			   }
			}
		   }
		break;
		}

	return total;
	}

static int
Compare(const void * a, const void * b)
	{
	Sharp *A = (Sharp *)a;
	Sharp *B = (Sharp *)b;

	if (B->bucket != A->bucket)
	   return (B->bucket>A->bucket)?1:-1;

	if (B->val > A->val) return (1);
	if (B->val < A->val) return (-1);
	return 0;
	}
	
void
QualitySort()
	{
	// Sort the list
	qsort(flist,flist_idx,sizeof(Sharp),Compare);
	Sorted = 1;
	}


void
WriteQFile(char *fname)
	{
	FILE *out = OpenLogfile(fname,"w");
	int i;

	if (! Sorted) QualitySort();

	if (out == NULL) {
	   Print("WriteQFile: cannot open '%s' for writing\n",fname);
	   return;
	   }

	// Sort the list
	qsort(flist,flist_idx,sizeof(Sharp),Compare);

	for(i=0; i<flist_idx; ++i) {
	   if (flist[i].newname) fprintf(out,"%s ",flist[i].newname);
	   fprintf(out,"%s %lf\n",flist[i].fname,flist[i].val);
	   }
	fclose(out);
	}

// Take the sorted list of files and renumber them
// arg is either QRENUMBER_FIRST or QRENUMBER_LAST, selecting which
// numeric portion of the filename to change
void
QRenumber(int which)
	{
	int b,i,j,ndigits;
	char *ptr=NULL,*newname=NULL,*p2=NULL;
	int bk,bcount[10];
	int qcount = 0;

	if (QTrim_Before || QTrim_After) _Trim(QTrim_Before,QTrim_After);

	if (! Sorted) QualitySort();

	// init all buckets to starting value
	for(i=0; i<9; ++i) bcount[i]=0;

	if (!Quiet) Print("\tCalculating new names...\n");
	for(i=0; i<flist_idx; ++i) {

	   // locate the start of the filename portion
	   ptr = flist[i].fname + strlen(flist[i].fname) - 1;
	   while(ptr != flist[i].fname && *(ptr-1) != '/') --ptr;

	   // ptr points to the start of the section that is to be replaced
	   if (which == QRENUMBER_FIRST) {
	      // locate first numeric section in the filename
	      while(*(ptr) && !isdigit(*ptr)) ++ptr;
	      }
	   else if (which == QRENUMBER_LAST) {
	      ptr = flist[i].fname + strlen(flist[i].fname) - 1;
	      while(ptr != flist[i].fname && *(ptr-1) != '/' && !isdigit(*ptr)) --ptr;
	      while (ptr != flist[i].fname && isdigit(*ptr) && isdigit(*(ptr-1))) --ptr;
	      }
	   else {
  	      Print("Invalid argument to -renumber : Aborted\n");
	      exit(1);
	      }

	   // p2 points to the first character past this numeric section
	   if (!isdigit(*ptr)) {
		Print("Filename '%s' has no numeric portion, renumbering cancelled\n",flist[i].fname);
		exit(1);
		}

	   p2 = ptr; while(isdigit(*p2)) ++p2;

	   if (p2 <= ptr) {
		Print("Filename '%s' does not start with digits. Cancelled\n",flist[i].fname);
		return;
		}

	   j = ptr - flist[i].fname; // number of bytes to copy from original name
	   bk = flist[i].bucket;  // category, used for correct numbering
	   b = bcount[bk]++;
	   
	   if (flist_idx <= 100) ndigits=2;
	   else if (flist_idx<=1000) ndigits=3;
	   else if (flist_idx<=10000) ndigits=4;
	   else if (flist_idx<=100000) ndigits=5;
	   else ndigits=6;  // This is not likely to happen

	   newname = (char *)malloc(strlen(flist[i].fname)+15); // make space for longer filename
	   strncpy(newname,flist[i].fname,j);
	   sprintf(newname+j,"q%0*d%s%s",ndigits,b,flist[i].active?"":"_trim",p2);

	   flist[i].newname = newname;
	   }

	// No errors, so now proceed with the renumbering

	// If set this tells us how many good frames we care about.
	// After we rename this many files we can start deleting the rest
	if (QualityMaxCount>0) qcount = QualityMaxCount;
	else qcount = flist_idx;

	if (!qcount) {
	   Print("oops, qcount is 0 is QRenumber!\n");
	   exit(1);
	   }

	if (!Quiet) Print("\tRenaming %d files...\n",qcount);
	for(i=0; i<flist_idx; ++i) {
	   char *old_name = flist[i].fname;
	   char *new_name = flist[i].newname;

	   if (! access(newname,0)) {
	   	if (AllowWrite(new_name)) unlink(new_name);
		else {
		   Print("qrenumber error: will not overwrite existing file '%s'\n",newname);
		   exit(1);
		   }
		}

	   if (i >= qcount) {
		// Past the end of the list we care about
		Print("Deleting [%s]\n",old_name);
		unlink(old_name);
		continue;
		}

	   if (rename(old_name,new_name)) {
		Print("Rename '%s' -> '%s' failed with error: ",old_name,new_name);
		perror(NULL);
		exit(1);
		}
	   }

	if (!Quiet) Print("\tdone.\n");
	return;
	}

/* "trim" N entries from the start and end of every set of buckets
 * by setting the active flag to 0.
 *
 * Note: This function processes the frames in the same order that they were added
 * by QAssociate().
 */

static void
_Trim(int before, int after)
	{
	int i;
	int bchange_idx[10];
	int cur_bucket,b;
	FILE *z;

	// umm, if no associated frames, do nothing.
	if (! flist_idx) return;

	/* build a list of bucket change points, init entries to
	 * -1 to signal "no bucket"
  	 */
	for(i=0; i<10; ++i) bchange_idx[i] = -1;

	cur_bucket = flist[0].bucket;
	for(i=1; i< flist_idx; ++i) {
	   b = flist[i].bucket;
	   if (b != cur_bucket) {
		bchange_idx[b] = i;
		cur_bucket = b;
		}
	   }

	/* Now we have entries in bchange_idx[] that tell us where the
	 * bucket changes occur, disable entries either side of each
	 * change point
	 */

	z = fopen(".ninox/trim.txt","w");

	for(i=0; i<10; ++i) if (bchange_idx[i]>=0) {
	   int p = bchange_idx[i];
	   int start = p - before;
	   int end = p + after;
	   int j;

	   if (start<0) start=0;
	   if (end >= flist_idx) end=flist_idx-1;

printf("Trimming (%d,%d) entries around bucket change %d (%d - %d)\n",before,after,i,start,end);
if (z) fprintf(z,"Trimming (%d,%d) entries around bucket change %d (%d - %d)\n",before,after,i,start,end);
	   for(j=start; j<=end; ++j) {
		flist[j].active=0;
		flist[j].val = -1;
printf("\t%s\n",flist[j].fname);
if (z) fprintf(z,"\t%s\n",flist[j].fname);
		}
	   }
	   
	if (z) fclose(z);

	return;
	}

void
QAssociate(char *fname, double quality)
	{
	char pfx[256],sfx[32];
	char *ptr;
	int n,bucket;

	if (flist_idx == MAX_FILES) {
	   Print("QAssociate: Max file count of %d reached\n",MAX_FILES);
	   exit(1);
	   }

	if (flist_idx == flist_size) {
	   // time to allocate some more entries
	   flist_size += 1000;
	   if (flist == NULL) flist = (Sharp *) ZeroMalloc(sizeof(Sharp) * flist_size);
	   else flist = (Sharp *) realloc(flist,sizeof(Sharp) * flist_size);

	   if (flist == NULL) {
		printf("QAssociate: Cannot allocate memory for %d entries\n",flist_size);
		exit(1);
		}
	   }

	flist[flist_idx].fname = strdup(fname);
	flist[flist_idx].val = quality;
	flist[flist_idx].newname = NULL;

	// Look for filename like NNN-D.fit, D is the bucket number
	ptr = fname + strlen(fname);
	while(ptr != fname && *(ptr-1) != '-') --ptr;

	if (sscanf(ptr,"%1d.%s",&n,sfx) == 2)
	   bucket = n;
	else 
	   bucket=0;
	
	flist[flist_idx].bucket = bucket;
	flist[flist_idx].active = 1;
	flist_idx++;
	}

/*
 * Generate a series of downsampled images and add up their quality
 * 
 * We handle FITS monochrome files only, either 16bpp or 8bpp. Process the region
 * represented by (x1,y1) to (x2,y2)
 */

// How many bright pixels to we average to get the real maximum value
#define MAXP 6

double
QualityEstimate(struct Image *img)
	{
	unsigned char *buffer = img->data;
	int width = img->width;
	int height = img->height;
	int depth = img->depth;
	int x1,y1,x2,y2;
	int subsample,region_w,region_h;
	int i,j,bpp,n,x,y,nbytes,max,maxp[MAXP],x_inc;
	int x_samples,y_samples,y_last;
	unsigned short *ubuffer,*buf,val;
	unsigned short *new=NULL;
	double mult,q,dval=0.0;
	int threshhold = ThreshHold;

	if (img->histo_warning) {
	   Print("Overbright warning, quality set to 0\n");
	   return 0;
	   }

	if (QualityFunction == HISTO) 
	   return Histo_Quality(img);

	if (QualityFunction == BR_THRESH) 
	   return BrThresh_Quality(img);

	if (depth != 16 && depth != 8 && depth != 24) {
	   Print("QualityEstimate: Depth must be 8,16 or 24, not %d\n",depth);
	   return -1;
	   }

	if (depth == 16) threshhold <<= 8;

	// dimensions of the region we want to analyse
	region_w = (img->cutout.x2 - img->cutout.x1 + 1);
	region_h = (img->cutout.y2 - img->cutout.y1 + 1);

	// We can assume the image has already been centered
	x1=(width-region_w)/2; y1=(height-region_h)/2;
	x2=x1 + region_w; y2 = y1 + region_h;

	Clip(width,height,&x1,&y1,&x2,&y2);

	bpp = depth/8;

	// Allocate the intermediate buffer. Will be 16bpp greyscale
	buf = (unsigned short *)malloc(region_w * region_h * 2);
	if (buf == NULL) {
	  Print("QualityEstimate: Out of memory");
	  exit(1);
	  }

	/* Write this image (debugging) */
	if (QWriteIntermediateFiles) {
	      unsigned short *udata = (unsigned short *)img->data;
	      char filename[1024];
	      FILE *out;

	      sprintf(filename,"sample_orig.ppm");
	      out = OpenLogfile(filename,"wb");
	      if (out == NULL) {
		Print("Cannot write subsampled image\n");
		exit(1);
		}

	      fprintf(out,"P5\n%d %d\n255\n",width,height);
	      if (img->depth == 16) for(i=0; i<width*height; ++i) putc(udata[i]>>8,out);
	      else if (img->depth == 8) for(i=0; i<width*height; ++i) putc(img->data[i],out);

	      fclose(out);
	      }

	subsample = QSUBSAMPLE_MIN;
	while(subsample <= QSUBSAMPLE_MAX) {
	   unsigned char *ptr;

	   /*
	    * Number of h & v pixels in subimage
	    */
	   x_samples = region_w / subsample;
	   y_samples = region_h / subsample;

	   if (x_samples < 2 || y_samples < 2) break;

	   y_last = y1 + (y_samples-1) * subsample; // second last row of subsampled output
	   x_inc = subsample * bpp;

	   for(i=0; i<MAXP; ++i) maxp[i]=0;

	   // First row - ignore histo-stretch
	   y=y1; n=0;
	   ptr = buffer + (y*width + x1) * bpp;
	   for(x=0; x < x_samples; ++x, ptr += x_inc) {
		  buf[n++] = SubSample(ptr,width,bpp,subsample,subsample);
		  }

	   // Rows 1 .. y_last-1 additional histo-stretch code
	   for(y+=subsample; y < y_last; y += subsample) {
	      ptr = buffer + (y*width + x1) * bpp;
	      for(x=0; x < x_samples; ++x, ptr += x_inc) {
		  register int v = SubSample(ptr,width,bpp,subsample,subsample);

          	  if (v > maxp[2] && v < 65530) {
			int slot;

              		if (v>maxp[0]) slot=0; else if (v>maxp[1]) slot=1; else slot=2;
                	for(j=MAXP-1; j>slot; --j) maxp[j] = maxp[j-1]; maxp[j]=v;
              		}

		  buf[n++] = v;
		  }
	     }

	   // Last row, ignore histo-stretch
	   ptr = buffer + (y*width + x1) * bpp;
	   for(x=0; x < x_samples; ++x, ptr += x_inc) {
		  buf[n++] = SubSample(ptr,width,bpp,subsample,subsample);
		  }

          // Average the bottom half brightest pixels to get the real max
           // to reduce noise effects
           j = MAXP/2;
           for(i=j,max=0; i<MAXP; ++i) max += maxp[i]; max /= (MAXP-j);

	   // Stretch histogram
	   if (max>0) {
	      mult = (double)60000 / (double)max;
	      for(i=0; i<n; ++i) {
		unsigned int v = buf[i]; v *= mult;
		if (v > 65535) v = 65535;
		buf[i] = v;
		}
	      }

	   // 3x3 smoothing
	   new = _smooth_image_16(buf,x_samples,y_samples);

	   /* Write this image (debugging) */
	   if (QWriteIntermediateFiles) {
	      char filename[1024];
	      FILE *out;

	      sprintf(filename,"sample_%d.ppm",subsample);
	      out = OpenLogfile(filename,"wb");
	      if (out == NULL) {
		Print("Cannot write subsampled image %d\n",subsample);
		exit(1);
		}

	      fprintf(out,"P5\n%d %d\n255\n",x_samples,y_samples);
	      for(i=0; i<n; ++i) putc(new[i]>>8,out);

	      fclose(out);
	      Print("[%d:%s] ",subsample,filename);
	      }

	   q = AnalyseImage16(new,x_samples,y_samples);
	   if (!Quiet) Print("%d:%2.2lf ",subsample,q);
	   dval += q;

	   if (new != NULL) {free(new); new=NULL;}

	   do { subsample += QSUBSAMPLE_INC; }
	   while (width/subsample == x_samples && height/subsample == y_samples);
	   }

	free(buf);
	return dval;
	}

static unsigned short *
_smooth_image_16(unsigned short *buf, int width, int height)
	{
	unsigned short *new = ZeroMalloc(width * height *2);
	int x,y;

	for(y=1; y<height-1; ++y) {
	   int o = y * width + 1;
	   for(x=1; x<width-1; ++x,++o) {
		unsigned int v = (unsigned int)buf[o];
		v += (unsigned int)buf[o-width-1];
		v += (unsigned int)buf[o-width];
		v += (unsigned int)buf[o-width+1];

		v += (unsigned int)buf[o-1];
		v += (unsigned int)buf[o+1];

		v += (unsigned int)buf[o+width-1];
		v += (unsigned int)buf[o+width];
		v += (unsigned int)buf[o+width+1];

		new[o] = v / 9;
		//new[o] = buf[o];
		}
	   }

	return new;
	}
	   
/*
 * Subsample a region starting at *ptr of size X size pixels.
 * Always return a 16bpp result
 */
int
SubSample(unsigned char *ptr, int img_wid, int bpp, int x_size, int y_size)
	{
	int x,y,val=0;
	unsigned short v,*uptr;
	double d;
	
	switch(bpp) {
	   case 1:
		for(y=0; y<y_size; ++y) {
	   	   for(x=0; x<x_size; ++x) val += ((unsigned short)ptr[x]) << 8;
		   ptr += img_wid;
	   	   }
		break;

	   case 2:
		uptr = (unsigned short *)ptr;
		for(y=0; y<y_size; ++y) {
	   	   for(x=0; x<x_size; ++x) val += uptr[x];
		   uptr += img_wid;
	   	   }
		break;

	   case 3:
		// 24bpp data, 8 bit BGR, convert to monochrome
		for(y=0; y<y_size; ++y) {
	   	   for(x=0; x<x_size*3; x+=3) {
			d = ((double)ptr[x])*0.114;
			d += ((double)ptr[x+1])*0.587;
			d += ((double)ptr[x+2])*0.299;
			val += ((unsigned short)(d))<<8;
			}
		   ptr += img_wid*3;
	   	   }
		break;

	   default:
		Print("SubSample: Unsupported depth %d\n",bpp);
		exit(1);
	   }

	return val / (x_size * y_size);
	}

#define MINIMUM(a,b) ((a)<(b)?(a):(b))
#define MAXIMUM(a,b) ((a)>(b)?(a):(b))

static int
MinimumFilter(unsigned short *in, unsigned short *out, int width, int height)
	{
	int x,y,d;
	unsigned short *ptr,*nptr;

	for(y=1; y<height-1; ++y) {
	   ptr = in + y * width;  	// start of row y
	   nptr = out + y * width;  	// start of row y

	   for(x=1; x<width-1; ++x) {
		d = ptr[x]; 

		d = MINIMUM(d,ptr[x-width-1]);
		d = MINIMUM(d,ptr[x-width]);
		d = MINIMUM(d,ptr[x-width+1]);

		d = MINIMUM(d,ptr[x-1]);
		d = MINIMUM(d,ptr[x+1]);

		d = MINIMUM(d,ptr[x+width-1]);
		d = MINIMUM(d,ptr[x+width]);
		d = MINIMUM(d,ptr[x+width+1]);

		nptr[x] = d;
		}
	   }
	return(1);
	}

static int
MaximumFilter(unsigned short *in, unsigned short *out, int width, int height)
	{
	int x,y,d;
	unsigned short *ptr,*nptr;

	for(y=1; y<height-1; ++y) {
	   ptr = in + y * width;  	// start of row y
	   nptr = out + y * width;  	// start of row y

	   for(x=1; x<width-1; ++x) {
		d = ptr[x]; 

		d = MAXIMUM(d,ptr[x-width-1]);
		d = MAXIMUM(d,ptr[x-width]);
		d = MAXIMUM(d,ptr[x-width+1]);

		d = MAXIMUM(d,ptr[x-1]);
		d = MAXIMUM(d,ptr[x+1]);

		d = MAXIMUM(d,ptr[x+width-1]);
		d = MAXIMUM(d,ptr[x+width]);
		d = MAXIMUM(d,ptr[x+width+1]);

		nptr[x] = d;
		}
	   }
	return(1);
	}

static double
Gradient(unsigned short *buf, int width, int height)
	{
	int pixels;
	int x,y;
	int yborder=height*QMargin + 1, xborder=width*QMargin + 1;
	double d1,d2;
	double val,avg=0;
	int threshhold = ThreshHold << 8;
	unsigned char *map = ZeroMalloc(width * height);

	// pass 1 locate all pixels > threshhold and flag the 3x3 region
	// around them for inclusion in the algorithm
	pixels=avg=0;
	for(y=yborder; y<height-yborder; ++y) {
	   int o = y * width+xborder;
	   for(x=xborder; x<width-xborder; ++x,++o) if (buf[o] >= threshhold) {
		map[o-width-1] = map[o-width] = map[o-width+1] = 1;
		map[o-1] = map[o] = map[o+1] = 1;
		map[o+width-1] = map[o+width] = map[o+width+1] = 1;
		++pixels;
		avg += buf[o];
		}
	   }

	// Average of the significant pixels
	if (!pixels) { val = -1.0; goto end;}
	avg /= (double)pixels;

	val=pixels=0;
	for(y=yborder; y<height-yborder; ++y) {
	   int o = y * width + xborder;  	// start of row y
	   for(x=xborder; x<width-xborder; ++x,++o) if (map[o]) {
		// Pixel differences
		d1 = d2 = buf[o];
		d1 -= (int)(buf[o+1]); 		if (d1<0) d1=-d1;
		d2 -= (int)(buf[o+width]);	if (d2<0) d2=-d2;

		val += d1 + d2;

		pixels++;
		}
	   }

	val /= (double) pixels; // normalise value to per-pixel

	val *= 50/avg;	// arbitrary increase + 
			// normalise differences to compensate for dim/bright images

	end:
	free(map);
	return val;
	}

static double
AnalyseImage16(unsigned short *buf, int width, int height)
	{
	unsigned short *newbuf;
	double total;

	//newbuf = (unsigned short *)ZeroMalloc(width * height * sizeof(unsigned short));
	//MinimumFilter(buf,newbuf,width,height);
	//total = Gradient(newbuf,width,height);
	//MaximumFilter(buf,newbuf,width,height);
	//total += Gradient(newbuf,width,height);

	//free(newbuf);

	total = Gradient(buf,width,height);

	return total;
	}

/*
 * Distance in colourspace between the pixels, normalised
 */
static double
PixDiff(unsigned char *p1, unsigned char *p2)
	{
	int r = *(p1++) - *(p2++);
	int g = *(p1++) - *(p2++);
	int b = *(p1++) - *(p2++);

	return (r*r) + (g*g) + (b*b) / 1966080.0;
	}

/*
 * Distance in colourspace between the pixels
 */
inline double
PixDiffLC1(unsigned char *p1, unsigned char *p2, int depth, double l, double c)
	{
	int b,g,r;

	if (depth==24) {
	   b = *(p1++) - *(p2++);	// 0.0 - 255.0
	   g = *(p1++) - *(p2++);	// 0.0 - 255.0
	   r = *(p1++) - *(p2++);	// 0.0 - 255.0
	   double L,C;

	   // square to exaggerate differences
	   r *= r; g *= g; b *= b;

	   // Luminance difference (square again)
	   L = r * 0.299 + g * 0.587 + b * 0.114;
	   L *= L;

	   // Colourspace difference (square again)
	   C = r + g + b;
	   C *= C;

	   // Weighted result, reduced
	   return (L * l + C * c) / (195075.0 * (l+c));
	   }

	if (depth==8) {
	   b = *(p1++) - *(p2++);
	   return ((double)b * (double)b) / 65025.0;
	   }
	}