#include "ninox.h"

static struct Image *Ref = NULL;

static void Deprotect8(struct Image *Img);
static int RoughAlign(struct Image *RealImage);
static void get_min_max(struct Image *img, int *MIN, int *MAX);
static double imagediff_8(struct Image *ref, struct Image *img, int x1, int y1, int x2, int y2, int x_off, int y_off, double err, int debug);
static inline int Pixel(unsigned char *data, int width, int depth, int x, int y);
static int LoadReference(char *fname);
static int morph_points(struct Image *ref, struct Image *img, struct morph_point *m, int n);
static int regenerate_image(struct Image *img, struct Image *RealImage, struct morph_point *m, int npoints);
static int smooth_image3x3(struct Image *img);
static int smooth_image5x5(struct Image *img);

static void Deprotect8(struct Image *Img);

#ifdef MSWIN32
unsigned random(void) { return 0; }
unsigned srandom(int x) { return 0; }
#endif

#define NBITS 0xf8

static int
LoadReference(char *fname)
	{
	int i,w,h,bpp;
	struct Image *r;

	// Unload any previous reference image
	DestroyImage(Ref);
	Ref=NULL;
	   
	if (access(fname,0)) {
	   Print("LoadReference: Image file '%s' does not exist\n",fname);
	   return(0);
	   }

	r = LoadImage(fname,"morph_ref.fit");
	if (! r) return 0;

	// convert to FITS 8bpp monochrome
	Ref = ConvertImage(r,IMG_FIT,8);
	if (!Ref) {
	   Print("LoadReference: Conversion to 8bpp failed\n");
	   exit(1);
	   }

	Deprotect8(Ref);
	DestroyImage(r);

	//Print("Loaded reference image %x (%d,%d,%d)\n",Ref,Ref->width,Ref->height,Ref->depth);

        if (UpScale > 1)
           upscale_image(Ref,UpScale);

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

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

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

	if (DisplayFrames) {
   	   ShowImage(Ref);
   	   Usleep(3000000);
   	   }

	return 1;
	}
	
#define MAXP 6

// image data guaranteed to be 16bpp monochrome
static void
get_min_max(struct Image *img, int *MIN, int *MAX)
	{
	unsigned char *buf = img->data;
	unsigned short *ubuf = (unsigned short *)img->data;
	int i,j,stop,min;
	int npixels = img->width * img->height;
	int max,maxp[MAXP];  // The brightest non-white pixels

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

	switch(img->depth) {
		case 8:
			max=0; min=255;
			for(i=0; i<npixels; ++i) {

           		   if (buf[i] > maxp[2] && buf[i] < 252) {
              			if (buf[i]>maxp[0]) stop=0;
              			else if (buf[i]>maxp[1]) stop=1;
              			else stop=2;

				for(j=MAXP-1; j>stop; --j) maxp[j] = maxp[j-1]; maxp[j]=buf[i];
              			}

	   		   if (buf[i]<min) min=buf[i];
	   		   }
			break;
		case 16:
			max=0; min=65535;
			for(i=0; i<npixels; ++i) {

           		   if (ubuf[i] > maxp[2] && ubuf[i] < 65530) {
              			if (ubuf[i]>maxp[0]) stop=0;
              			else if (ubuf[i]>maxp[1]) stop=1;
              			else stop=2;

				for(j=MAXP-1; j>stop; --j) maxp[j] = maxp[j-1]; maxp[j]=ubuf[i];
              			}

	   		   if (ubuf[i]<min) min=ubuf[i];
	   		   }
			break;
		default: Print("get_min_max: depth %d not supported\n",img->depth);
			exit(1);
			break;
		}

	j = MAXP/2;
        for(i=j,max=0; i<MAXP; ++i) max += maxp[i]; max /= (MAXP-j);

	*MIN = min;
	*MAX = max;
	}
	
// Adjust the intensity values in our image to match those in the reference
static int
stretch_to_ref(struct Image *ref, struct Image *img)
	{
	int rmin,rmax,bmin,bmax;
	int i,j,r_range,b_range,o;
	double scale;
	unsigned char *buf = img->data;
	unsigned short *ubuf = (unsigned short *)img->data;

	if (ref->depth != img->depth) {
	   Print("stretch_to_ref: ref and img must both be same depth, not %d/%d\n",ref->depth,img->depth);
	   exit(1);
	   }

	// Intensity balance
	get_min_max(ref,&rmin,&rmax);
	get_min_max(img,&bmin,&bmax);

	r_range = rmax - rmin;
	b_range = bmax - bmin;
	// Print("stretch_to_ref: ref range: %d-%d(%d), img range %d-%d(%d)\n",rmin,rmax,r_range,bmin,bmax,b_range);

	o = rmin - bmin;
	scale = (double)r_range / (double)b_range;
	switch(img->depth) {
	   case 8:
		i = img->width * img->height;
		while(i--) {
	   	   j = *buf; j+=o; j*= scale;
	   	   if (j<0) j=0; if (j>255) j=255;
	   	   *(buf++)=j;
	   	   }
		break;
	   case 16:
		i = img->width * img->height;
		while(i--) {
	   	   j = *ubuf; j+=o; j*= scale;
	   	   if (j<0) j=0; if (j>65535) j=65535;
	   	   *(ubuf++)=j;
	   	   }
		break;
	   default:
		Print("stretch_to_ref: depth %d not supported\n",img->depth);
		exit(1);
		break;
	   }

	//Print("Adjust image: o=%d, scale=%lf\n",o,scale);
	return 1;
	}

// remove any possible histogram protection on 8 bit images
static void
Deprotect8(struct Image *Img)
	{
	if (Img->depth == 8 && Img->data[0]==0 && Img->data[1]>0) {
	   Img->data[0] = Img->data[Img->width];
	   Img->data[1] = Img->data[Img->width+1];
	   }
	else {
	   //printf("not deprotecting: depth=%d, 0/1 = %d/%d\n",Img->depth,Img->data[0],Img->data[1]);
	   }
	}

// Roughly align the given image with our reference image.
//

static int
RoughAlign(struct Image *RealImage)
	{
	struct Image *Img;
	struct Region reg;
	unsigned short *rbuf;
	unsigned short *buf;
	unsigned char *ptr;
	int n,i,j,l,npixels,X,Y,x,y,div,x_off,y_off;
	int N,o,d;
	int first=1,b_range,r_range,w,h;
	double v,min,scale;
	struct Image **rlist,*Ref8,*ref,*img;
	struct Image **list;
	int x_sizes[32],y_sizes[32];
	int y_offsets[10],x_offsets[10];

	Ref8 = ConvertImage(Ref,IMG_FIT,8);
	if (!Ref8) {
	   Print("morph_image: cannot convert Ref image to 8bpp\n");
	   exit(1);
	   }
	Ref8->src_fname=strdup("ref");

	// Create an 8bpp copy
	Img = ConvertImage(RealImage,IMG_FIT,8);
	if (!Img) {
	   Print("morph_image: cannot convert to 8bpp\n");
	   exit(1);
	   }
	Img->src_fname=strdup("img");

	// Make sure these images are not protected - ie that pixels 0 and 1 are
	// both black.
	Deprotect8(Img);

	if (Img->depth != Ref->depth) {
	   Print("RoughAlign: Mismatched depth between ref and %s\n",Img->src_fname);
	   DestroyImage(Img);
	   DestroyImage(Ref8);
	   return 0;
	   }

	x_off=y_off=0;
	fft_register_pad(Ref,Img,FFT_LPF_Radius1,FFT_LPF_Radius2,FFT_LPF_Pow,Morph_MaxTranslate,&x_off,&y_off);

	// Clamp to make sure we don't go too far
	if (x_off >  Morph_MaxTranslate) x_off =  Morph_MaxTranslate;
	if (x_off < -Morph_MaxTranslate) x_off = -Morph_MaxTranslate;
	if (y_off >  Morph_MaxTranslate) y_off =  Morph_MaxTranslate;
	if (y_off < -Morph_MaxTranslate) y_off = -Morph_MaxTranslate;

Print("translate [%d,%d]\n",x_off,y_off);

	translate_image(RealImage,x_off,y_off);

	DestroyImage(Img);
	DestroyImage(Ref8);

	return 1;
	}

int
translate_image(struct Image *img, int x1, int y1)
	{
	unsigned char *buf,*sptr,*dptr;
	int width=img->width;
	int height=img->height;
	int bpp = img->depth/8;
	int rowbytes = width * bpp;
	int copybytes;
	int offset,o;
	int i,x2,y2,pfx,sfx;

	// offset between input & output
	offset = (y1*width + x1) * bpp;

	// Determine the output rectangle
	x2 = x1 + width - 1;
	y2 = y1 + height - 1;
	Clip(width,height,&x1,&y1,&x2,&y2);

	buf = Malloc(width * height * bpp);

	// src and dst pointers for first row to copy
	o = (y1*width+x1)*bpp;
	sptr = img->data + o - offset;

	// zero all the rows before our output
	dptr=buf;
	i = y1*rowbytes; while(i--) *(dptr++)=0;

	copybytes = (x2-x1+1) * bpp;
	pfx=x1*bpp; sfx=rowbytes-copybytes-pfx;

	// active area in output image for next stage of morphing
	//img->morph->reg.x1 = x1;
	//img->morph->reg.y1 = y1;
	//img->morph->reg.x2 = x2;
	//img->morph->reg.y2 = y2;

	while(y1<=y2) {
	   i=pfx; while(i--) *(dptr++)=0;
	   memcpy(dptr,sptr,copybytes); sptr+=rowbytes;
	   dptr +=copybytes;
	   i=sfx; while(i--) *(dptr++)=0;
	   ++y1;
	   }

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

// Sample a pixel form the image, return as 8bpp monochrome
static inline int
Pixel(unsigned char *data, int width, int depth, int x, int y)
	{
	unsigned int v;
	int o = y*width + x;

	switch(depth) {
	   case 8:
		return data[o];
		break;
	   case 16:
	 	v = *((unsigned short *)data + o);
	 	return v>>8;
		break;
	   case 24:
		// BGR from BMP image format
		o*=3;
		return (int)((data[o]*0.114) + (data[o+1]*0.587) + (data[o+2]*0.299));
		break;
	   }

	Print("Pixel: depth %d not understood\n",depth);
	exit(1);
	return 0;
	}

/*
 * Create a map of double precision floats of dimension (w,h)
 * The central third of this map has weight 1.0 and weights decrease to
 * zero at the edges
 */

static double *
DiffWeightMap1(int w, int h)
	{
	double *map = ZeroMalloc(w*h*sizeof(double));
	int x1,y1,x2,y2;
	int o,x,y;
	double d,y_inf,x_inf;

	// Define the central region where influence is 1.0
	x1 = w/4; x2 = x1 + w/2;
	y1 = h/4; y2 = y1 + h/2;

	for(y=o=0; y<h; ++y) {
	   y_inf = 1.0;
	   if (y<y1) y_inf = (double)y/(double)y1;
	   else if (y>y2) y_inf = 1.0 - (double)(y-y2)/(double)(h-y2);

	   for(x=0; x<w; ++x,++o) {
		x_inf = 1.0;
	      	if (x<x1) x_inf = (double)x/(double)x1;
	      	else if (x>x2) x_inf = 1.0 - (double)(x-x2)/(double)(w-x2);

		d = x_inf * y_inf;
		if (d>1.0) {
		   printf("oops, d=%lf, x_inf=%lf y_inf=%lf\n",d,x_inf,y_inf);
		   exit(1);
		   }
		map[o] = d;
		}
	   }

	return map;
	}

/*
 * Create a map of double precision floats of dimension (w,h)
 * The central quarter of this map has weight 1.0 and weights decrease to
 * zero at the edges
 */

static double *
DiffWeightMap2(int w, int h)
	{
	double *map = ZeroMalloc(w*h*sizeof(double));
	int x1,y1,x2,y2;
	int o,x,y;
	double d,y_inf,x_inf;

	// Define the central region where influence is 1.0
	x1 = w/3; x2 = x1 + w/3;
	y1 = h/3; y2 = y1 + h/3;

	for(y=o=0; y<h; ++y) {
	   y_inf = 1.0;
	   if (y<y1) y_inf = (double)y/(double)y1;
	   else if (y>y2) y_inf = 1.0 - (double)(y-y2)/(double)(h-y2);

	   for(x=0; x<w; ++x,++o) {
		x_inf = 1.0;
	      	if (x<x1) x_inf = (double)x/(double)x1;
	      	else if (x>x2) x_inf = 1.0 - (double)(x-x2)/(double)(w-x2);

		d = x_inf * y_inf;
		if (d>1.0) {
		   printf("oops, d=%lf, x_inf=%lf y_inf=%lf\n",d,x_inf,y_inf);
		   exit(1);
		   }
		map[o] = d * d;
		}
	   }

	return map;
	}

// both images are guaranteed 8bpp monochrme
static double
imagediff_8(struct Image *ref, struct Image *img, int x1, int y1, int x2, int y2, int x_off, int y_off, double err, int debug)
	{
	int width = ref->width;
	int height = ref->height;
	int x,y,count=0,diff;
	int o,oo,off,p1,p2,count2;
	static double *wmap = NULL;
	static int wmap_w = -1,wmap_h = -1;
	double d=0;

	if (ref->depth != 8 || img->depth != 8) {
	   Print("imagediff_8: images must be 8bpp and not %d/%d\n",ref->depth,img->depth);
	   exit(1);
	   }

	if (wmap_w != (x2-x1+1) || wmap_h != (y2-y1+1) || wmap == NULL) {
	   wmap_w = x2-x1+1;
	   wmap_h = y2-y1+1;
	   if (wmap) free(wmap);
	   wmap = DiffWeightMap1(wmap_w,wmap_h);
	   }

	off = y_off * width + x_off;
	count2=0;
	for(y=y1; y<=y2; ++y) {
	   if (y+y_off < 0 || y+y_off >=height) continue;
	   o = y * width + x1;

	   oo = (y-y1) * wmap_w;

	   for(x=x1; x<=x2; ++x,++o,++oo)
		//if (x+x_off>=0 && x+x_off<width && (img->data[o]>=ThreshHold || ref->data[o+off]>=ThreshHold)) {

		// if (x+x_off>=0 && x+x_off<width && img->data[o]>=ThreshHold && ref->data[o+off]>=ThreshHold) {

		if (x+x_off>=0 && x+x_off<width) {

//if (debug)
//printf("o=%d off=%d (%d,%d) = %d / (%d,%d) = %d\n",o,off,x,y,img->data[o],x+x_off,y+y_off,ref->data[o+off]);

		   diff =  img->data[o];
		   diff -= ref->data[o+off];

		   d += diff * diff * wmap[oo];
//if (debug) img->data[o] = 255 * wmap[oo];
		   ++count;
		}
	   }

	if (count==0) {
	   return err;
	   }

	return d/(double)count;
	}

	
// NEW Add the influence of a control point to a region
static int
add_tile_influence(struct Image *img, struct weight *mesh,int w, struct morph_point *pt)
	{
	int o,oo,x,y,i,j;
	struct weight *wlist;
	int x1 = pt->x - pt->w;
	int x2 = x1 + pt->w*2 - 1;
	int y1 = pt->y - pt->h;
	int y2 = y1 + pt->h*2 - 1;
	static double *wmap = NULL;
	static int wmap_w = -1,wmap_h = -1;

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

	if (wmap_w != (x2-x1+1) || wmap_h != (y2-y1+1) || wmap == NULL) {
	   wmap_w = x2-x1+1;
	   wmap_h = y2-y1+1;
	   if (wmap) free(wmap);
	   wmap = DiffWeightMap2(wmap_w,wmap_h);
	   }

	for(j=0,y=y1; y<=y2; ++y) {
	   o = y*w + x1;
	   oo = (y-y1) * wmap_w;
	   for(x=x1; x<=x2; ++x,++o,++oo,++j) {
		wlist = mesh + o*4;
		for(i=0; i<4 && wlist[i].active; ++i); if (i>3) {
		   Print("Mesh overfull at (%d,%d)\n",x,y);
		   exit(1);
		   }
		wlist[i].active=1;
		wlist[i].wgt = wmap[oo];
		wlist[i].ptr = pt;
		if (Morph_Debug) draw_point(img,x,y,(int)(wlist[i].wgt * 255));
		}
	   }
	if (Morph_Debug) { printf("morph_debug: add_tile_influence\n"); fflush(stdout); ShowImage(img); Usleep(200*1000); }
	return 1;
	}

// ref and img and guaranteed to be 8bpp monochrome
static int
morph_points(struct Image *ref, struct Image *img, struct morph_point *m, int n)
	{
	struct morph_point *pt;
	int i,j,x1,y1,x2,y2,x,y,d,X,Y,changed=0;
	int acount,p,nx,ny;
	double min,v;

	struct Point {
	   int x,y;
	   } *plist;

	int npoints;

	d = Morph_Drift;
	if (UpScale>1) d *= UpScale;
	if (DownScale>1) d /= DownScale;
	if (d<1) d = 1;

	// set up a mapping to randomly scan the region (-d,-d) .. (+d,+d)
	// make sure the first point we scan is the centre to set a good default value
	// (ie no movement).

	npoints = d*2+1; npoints *= npoints;

	plist = (struct Point *)ZeroMalloc(sizeof(struct Point) * npoints);
	for(i=0,y=-d; y<=+d; ++y)
	   for(x=-d; x<=d; ++x) {
		plist[i].x = x;
		plist[i].y = y;
		i++;
		}

	acount=0;
	for(i=0; i<n; ++i) {
	   pt = m+i;

	   if (pt->w==0 && pt->h == 0) continue;
	   if (pt->active == 0) continue;

	   x1 = pt->x - pt->w/2; x2 = x1 + pt->w - 1;
	   y1 = pt->y - pt->h/2; y2 = y1 + pt->h - 1;

	   // randomise
	   for(j=0; j<npoints; ++j) {
	      int r1 = Random(npoints);
	      int r2 = Random(npoints);
	      if (r1 != r2) {
		int x = plist[r1].x;
		int y = plist[r1].y;

		plist[r1].y = plist[r2].y; plist[r2].y = y;
		plist[r1].x = plist[r2].x; plist[r2].x = x;
		}
	      }

	   // Move the (0,0) point to the front
	   for(j=1; j<npoints; ++j) if (plist[j].x==0 && plist[j].y==0) {
	      int x = plist[0].x;
	      int y = plist[0].y;

	      plist[0].y = plist[j].y; plist[j].y = y;
	      plist[0].x = plist[j].x; plist[j].x = x;
	      break;
	      }
	
	   min = 999999999999.0;
	   X=pt->x_offset;
	   Y=pt->y_offset;
	   nx=X; ny=Y;

	   for(p=0; p<npoints; ++p) {
		struct Image *I1,*R1;
		unsigned short *id;
		unsigned short *rd;
		int xx,yy;
		x = plist[p].x + X;
		y = plist[p].y + Y;

//I1 = ConvertImage(img,IMG_FIT,16);
//R1 = ConvertImage(ref,IMG_FIT,16);

//printf("img:\n"); fflush(stdout);
//ShowImage(I1); Usleep(5 * 1000 * 1000);

//printf("img:\n"); fflush(stdout);
//ShowImage(R1); Usleep(5 * 1000 * 1000);

//printf("scanning (%d,%d) - (%d,%d) + [%d,%d] diff=%lf\n",x1,y1,x2,y2,x,y);
//fflush(stdout);
		v = imagediff_8(ref,img,x1,y1,x2,y2,x,y,min,1);
		if (v< (min-min/20)) { min = v; nx = x; ny = y;}

//printf("img: (%d,%d) - (%d,%d) + [%d,%d] diff=%lf\n",x1,y1,x2,y2,x,y,v);
#if 0
if (p==0) {
   struct Image *I = ConvertImage(img,IMG_FIT,16);
   ShowImage(I);
   Usleep(5*1000*1000);
   }
#endif
		}

//printf("point %d, offsets %d,%d\n",i,nx,ny);

	   if (nx != X || ny != Y) {
		changed ++;
	   	pt->x_offset = nx;
	   	pt->y_offset = ny;
		}
	   }

	free(plist);
	return changed;
	}

int MAX_OFFSET_X = 0;
int MAX_OFFSET_Y = 0;

static int
read_pixel(struct Image *img, int x, int y, struct weight *wlist, unsigned char *pixel)
	{
	int i,o;
	double x_offset=0;
	double y_offset=0;
	unsigned char *data = img->data;
	unsigned short *udata = (unsigned short *)img->data;
	unsigned short *up = (unsigned short *)pixel;

	// up to 4 control points contribute to the offset
	for(i=0; i<4 && wlist[i].active; ++i) {
	   x_offset += wlist[i].wgt * wlist[i].ptr->x_offset;
	   y_offset += wlist[i].wgt * wlist[i].ptr->y_offset;
	   }

	x -= x_offset;
	y -= y_offset;

	if (x<0 || x >= img->width || y<0 || y >= img->height) {
	   //Print("merge_pixel: (%d,%d) out of bounds: (%d,%d) + (%d,%d)\n",x,y,X,Y,pt->x_offset,pt->y_offset);
	   return 0;
	   }

if (abs(x_offset)>MAX_OFFSET_X) MAX_OFFSET_X=x_offset;
if (abs(y_offset)>MAX_OFFSET_Y) MAX_OFFSET_Y=y_offset;

	o = y*img->width + x;

	switch(img->depth) {
	   case 8:
		*pixel = data[o];
		break;
	   case 16:
		*up= udata[o];
		break;
	   case 24:
		o*=3;
		*(pixel++)= data[o++];
		*(pixel++)= data[o++];
		*(pixel)= data[o];
		break;
	   }

	return 1;
	}

// Apply the morphing to generate a new image
static int
regenerate_image(struct Image *img, struct Image *RealImage, struct morph_point *m, int npoints)
	{
	unsigned char *buf;
	struct weight *mesh,*wlist;
	int width = img->width;
	int height = img->height;
	int bpp = img->depth/8;
	struct morph_point *pt;
	int x,y,i,j,n,x1,y1,x2,y2;
	int o,mo,v,missed;
	unsigned char pixel[3];

	n = width * height;
	mesh = (struct weight *) ZeroMalloc(sizeof(struct weight) * n * 4);

	// Add the influence of each tile to its surroundings
	for(i=0; i<npoints; ++i) {
	   pt = m + i;
	   if (pt->w || pt->h)
	      add_tile_influence(img,mesh,img->width,pt);
	   }

#if 0
if (Morph_Debug) {
   printf("morph_debug: showing mesh\n");
   ShowImage(img);
   while(1) Usleep(1000*1000);
   return (1);
   }
#endif

	// pass 2: generate output image
	buf = ZeroMalloc(RealImage->width * RealImage->height * RealImage->depth/8);

MAX_OFFSET_X=0;
MAX_OFFSET_Y=0;

	wlist = mesh; o=missed=0;
	for(y=0; y<height; ++y) {
	   switch(RealImage->depth) {
		case 8:
	   		for(x=0; x<width; ++x,wlist+=4)
			   if (read_pixel(RealImage,x,y,wlist,pixel))
				buf[o++]=*pixel;
			   else { buf[o] = RealImage->data[o++]; missed++;}
			break;
		case 16:
	   		for(x=0; x<width; ++x,wlist+=4)
			   if (read_pixel(RealImage,x,y,wlist,pixel)) {
				buf[o++]=*pixel;
				buf[o++]=*(pixel+1);
				}
			   else { buf[o]=RealImage->data[o++]; buf[o]=RealImage->data[o++]; missed++;}
			break;
		case 24:
	   		for(x=0; x<width; ++x,wlist+=4)
			   if (read_pixel(RealImage,x,y,wlist,pixel)){
				buf[o++]=*pixel;
				buf[o++]=*(pixel+1);
				buf[o++]=*(pixel+2);
				}
			   else { buf[o]=RealImage->data[o++]; buf[o]=RealImage->data[o++]; buf[o]=RealImage->data[o++]; missed++;}
			break;
		default:
			Print("regenerate: invalid depth %d\n",RealImage->depth);
			exit(1);
			break;
		}
	   }

	free(RealImage->data);
	RealImage->data = buf;

	// free the mesh
	free(mesh);

printf("Max Offsets: %d,%d\n",MAX_OFFSET_X,MAX_OFFSET_Y);

	return 1;
	}

// Given an image of type 8-bit greyscale, and a region in that image,
// return the number of pixels above the lower threshhold brightness
static int
u8_pixels_above_threshhold(struct Image *img, int x1, int y1, int w, int h)
	{
	int x2 = x1 + w -1;
	int y2 = y1 + h -1;
	int x,y,count;
	unsigned char *ptr;

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

	count=0; 
	for(y=y1; y<=y2; ++y) {
	   ptr = img->data + y * img->width + x1;
	   for(x=x1; x<=x2; ++x,++ptr)
		if (*ptr > ThreshHold) ++count;
	   }
		
	return count;
	}

static void
allocate_points(struct Image *img,struct morph_point *p, int npoints, int x1, int y1, int x2, int y2, int w, int h)
	{
	int x,y,n,W,H,X1,Y1,X2,Y2;
	int x_offset=0,y_offset=0;
	int cutoff;

	// If the tiles are an odd size then make an adjustment otherwise they will not
	// pack properly and we'll get gaps...
	if (w&1) w++; if (h&1) h++;

	W = x2-x1+1; H=y2-y1+1;

	x_offset = Random(w);
	y_offset = Random(h);

	X1 = x1 + x_offset; X2 = X1 + W - 1;
	Y1 = y1 + y_offset; Y2 = Y1 + H - 1;

	if (X1<0) X1=0; if (Y1<0) Y1=0;
	if (X2>x2) X2=x2; if (Y2>y2) Y2=y2;

	// set a cutoff for active pixel count in each tile. Tiles that do
	// not meet this count are not scanned.
	cutoff = (w * h) * 0.05 ;

	for(n=0,y=Y1; y<=Y2-h; y+=h)
	   for(x=X1; x<=X2-w; x+=w,++n) {
	   	p[n].x = x + w/2;
		p[n].y = y + h/2;
		p[n].w = w;
		p[n].h = h;
		p[n].x_offset=0;
		p[n].y_offset=0;
		p[n].scale=1;
		if (u8_pixels_above_threshhold(img,x,y,w,h) > cutoff) p[n].active = 1;
		else p[n].active = 0;
	   	}
		
	while(n<npoints) {
	   p[n].x=0; p[n].y=0;
	   p[n].w=0; p[n].h=0;
	   p[n].active = 0;
	   ++n;
	   }

	if (0 && DisplayFrames) {
   	   struct Image *I = ConvertImage(img,IMG_FIT,16);
	   int i;

	   for(i=0; i<npoints; ++i) if (p[i].active) {
		int x = p[i].x;
		int y = p[i].y;
		int x1 = x - p[i].w/2;
		int y1 = y - p[i].h/2;
	 	int x2 = x1 + p[i].w-1;
	 	int y2 = y1 + p[i].h-1;
		int X,Y;
		int T = ThreshHold << 8; T *= 2;
		unsigned short *ptr;

		for(Y=y1; Y<=y2; ++Y) {
		   ptr = (unsigned short *)I->data + Y*I->width + x1;
		   for(X=x1; X<=x2; ++X) {
			int v = *ptr; if (v<T) {v += 16*255; if (v>65535) v=65535; *ptr = v;}
			++ptr;
			}
		   }
		}

   	   ShowImage(I);
   	   //Usleep(2000000);
	   DestroyImage(I);
   	   }
	if (n>npoints) 
	   Print("Allocate_1: n=%d\n",n);

	return;
	}

// Morph the image *img to approximately the same shape as the
// reference image *Ref
static int
MorphAlign(struct Image *ref, struct Image *RealImage)
	{
	struct morph_point *points;
	struct Image *img;
	int n,x,y,x1,y1,x2,y2,w,h,i,j;
	char tmpstr[32];

	if (Morph_Across <= 0) {
	   Print("Missing morph_across parameter\n");
	   return 0;
	   }

	x1=0; x2 = RealImage->width-1;
	y1=0; y2 = RealImage->height-1;

	// Distance between control points
	w = (x2-x1+1)/Morph_Across;

	// make tile width into the nearest power of 2
	i=1; while((i<<1) < w) i<<=1;
	w = (w + (i>>1)) & (i<<1) ? i<<1 : i;

	h = w;  // square tiles

	// Given only one param, calculate the other
	Morph_Down = RealImage->height / h;

	if (Morph_Across <= 0 || Morph_Down <= 0) {
	   Print("Invalid: morph_across=%d and morph_down=%d\n",Morph_Across,Morph_Down);
	   return 0;
	   }

	n = Morph_Across * Morph_Down;

	Print("calculated %d morph control points (%d,%d)\n",n,Morph_Across,Morph_Down);

	points = (struct morph_point *)ZeroMalloc(n * sizeof(struct morph_point));
	
	Print("tiles are %d x %d\n",w,h);

	img = ConvertImage(RealImage,IMG_FIT,8);
	if (!img) {
	      Print("morph_image: cannot convert to 8bpp\n");
	      exit(1);
	      }

	Deprotect8(img);
	stretch_to_ref(ref,img);
	smooth_image5x5(img);

	if (img->depth != ref->depth) {
	   Print("MorphAlign: Image %s has depth %d, doesn't match reference image depth of %d\n",img->src_fname,img->depth,ref->depth);
	   DestroyImage(img); free(points);
	   return 0;
	   }

	if (0 && DisplayFrames) {
   	   ShowImage(img);
   	   //Usleep(1000000);
   	   }

	// Create our grid of control points, move them around to reduce
	// artifacts
	allocate_points(img,points,n,x1,y1,x2,y2,w,h);
	morph_points(ref,img,points,n);
	regenerate_image(img,RealImage,points,n);

	DestroyImage(img);
	free(points);
	return 1;
	}

static int
smooth_image3x3(struct Image *img)
	{
	int width = img->width;
	int height = img->height;
	int nbytes = width * height * img->depth/8;
	unsigned char *buf = Malloc(nbytes);
	unsigned char *src = img->data;
	int x,y,o;

	if (img->depth != 8) {
	   Print("smooth_image3x3: unsupported depth %d\n",img->depth);
	   exit(1);
	   }

	memcpy(buf,img->data,nbytes);

	for(y=1; y<height-1; ++y) {
	   o = y*width + 1;
	   for(x=1; x<width-1; ++x,++o)
	   	buf[o] = (src[o] + src[o-1] + src[o+1] + src[o-width] + src[o+width])/5;
	    }

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

// 5x5 smoothing, discard brightest and dimmest pixel in each block, average remaining 23
static int
smooth_image5x5(struct Image *img)
	{
	int width = img->width;
	int height = img->height;
	int nbytes = width * height * img->depth/8;
	unsigned char *buf = Malloc(nbytes);
	unsigned char *src = img->data;
	int x,y,o;

	if (img->depth != 8) {
	   Print("smooth_image3x3: unsupported depth %d\n",img->depth);
	   exit(1);
	   }

	memcpy(buf,img->data,nbytes);

	for(y=2; y<height-2; ++y) {
	   o = y*width + 2;
	   for(x=2; x<width-2; ++x,++o) {
	   	unsigned int v;
		int oo = o-width*2;

		v = src[oo-2] + src[oo-1] + src[oo] + src[oo+1] + src[oo+2];
		oo += width; v += src[oo-2] + src[oo-1] + src[oo] + src[oo+1] + src[oo+2];
		oo += width; v += src[oo-2] + src[oo-1] + src[oo] + src[oo+1] + src[oo+2];
		oo += width; v += src[oo-2] + src[oo-1] + src[oo] + src[oo+1] + src[oo+2];
		oo += width; v += src[oo-2] + src[oo-1] + src[oo] + src[oo+1] + src[oo+2];
	
		buf[o] = v / 25;
		}
	    }

	free(img->data);
	img->data = buf;
	return 1;
	}
int
morph_image(struct Image *I, int level)
	{
	int i;

	if (!Ref) LoadReference(Morph_Ref);
	if (! Ref) {
	   Print("Error: Could not load Morph reference image '%s'\n",Morph_Ref);
	   return 0;
	   }

	if (0 && DisplayFrames) ShowImage(I);

	// Check that we're not aligning the reference image with itself!
	if (!strcmp(I->src_fname,Ref->src_fname)) {
	   Print("Aligning reference frame with itself, ignored\n");
	   return 0;
	   }

	if (I->width != Ref->width || I->height != Ref->height) {
	   Print("reference frame '%s' is %dx%d, image '%s' is %dx%d, mismatched!\n",
		Ref->src_fname,Ref->width,Ref->height,I->src_fname,I->width,I->height);
	   return 0;
	   }

	for(i=0; i<Morph_Iterations; ++i) {
	   // Step 1 is an overall alignment
	   if (! RoughAlign(I)) {
	      Print("Rough Alignment failed\n");
	      return 0;
	      }

	   if (level>1 && ! MorphAlign(Ref,I)) {
	      Print("Morph Alignment failed\n");
	      return 0;
	      }
	   if (0 && DisplayFrames) ShowImage(I);
	   }

	if (level==3 && ! RoughAlign(I)) {
	   Print("Final Alignment failed\n");
	   return 0;
	   }

	return 1;
	}