#include "ninox.h"

static struct Image *Ref = NULL;

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);
static int translate_image(struct Image *img, int x_off, int y_off);
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 void Deprotect8(struct Image *Img);

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

static int
LoadReference(char *fname)
	{
	int 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);
	   }

	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;

	smooth_image3x3(Ref);

	if (DisplayFrames) {
   	   struct Image *R = ConvertImage(Ref,IMG_FIT,16);
   	   ShowImage(R);
	   DestroyImage(R);
   	   Usleep(5000000);
   	   }

	return 1;
	}
	
// 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,min,max;
	int npixels = img->width * img->height;

	switch(img->depth) {
		case 8:
			max=0; min=255;
			for(i=0; i<npixels; ++i) {
	   		   if (buf[i]<252 && buf[i]>max) max=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]<65530 && ubuf[i]>max) max=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;
		}

	*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("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] = 0;
	   Img->data[1] = 0;
	   }
	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.
//
// We do this by downsampling our images and then working from smallest -> original size
// and aligning as best we can at each level.

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];
	static int last_x_off = 999;
	static int last_y_off = 999;
	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(Ref8);
	Deprotect8(Img);

	smooth_image3x3(Img);

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

	// Smallest image we bother with is 16 pixels across
	for(i=1; (1<<i) < Img->width && (1<<i) < Img->height; ++i); --i;
	i -= Morph_Min_Divisor;
	if (i<1) {
	   Print("RoughAlign: Can't find a good divisor!\n");
	   DestroyImage(Img);
	   return 0;
	   }

	div = 1<<i;

	// Generate the list of required image sizes
	n=0; while(div>1) {
	   x_sizes[n] = Img->width/div;
	   y_sizes[n] = Img->height/div;
//printf("divisor %d, w=%d, h=%d\n",div,x_sizes[n],y_sizes[n]);
	   ++n; div>>=1;
	   }

	reg.x1 = 0; reg.x2 = Img->width-1;
	reg.y1 = 0; reg.y2 = Img->height-1;

	// Get the downsampled set of reference images
	rlist = GenerateSubSampledImages(Ref8,&reg,x_sizes,y_sizes,8,n);
	if (! rlist) {
	   Print("Error generating subsampled images\n");
	   exit(1);
	   }

	// Get the downsampled set of src images
	list = GenerateSubSampledImages(Img,&reg,x_sizes,y_sizes,8,n);
	if (! list) {
	      Print("Error generating subsampled images\n");
	      exit(1);
	      }

	x_off = y_off = 0;
	for(l=0; l<=n; ++l) {
	   if (l<n) {
	      	ref = rlist[l];
	      	img = list[l];
		}
	   else { ref = Ref8; img = Img; }

	   stretch_to_ref(ref,img);

	   // Width and Height of subsampled images
	   w = ref->width;
	   h = ref->height;
	   npixels = w*h;
	   x_off *= 2; y_off *= 2;

	   // find the best alignment
	   if (first) d=4; else d=2;

	   // scan centre location first
	   X = x_off; Y=y_off;
	   min= imagediff_8(ref,img,0,0,w-1,h-1,X,Y,9999.0);

	   // generate the list of offsets and randomise them
	   for(y=y_off-d,N=0; y<=y_off+d; ++y) y_offsets[N++]=y;
	   for(x=x_off-d,N=0; x<=x_off+d; ++x) x_offsets[N++]=x;

	   for(i=0; i<N; ++i) {
		int j=Random(N);
		int k=Random(N);
		int tmp = y_offsets[j]; y_offsets[j]=y_offsets[k]; y_offsets[k]=tmp;

		j=Random(N); k=Random(N);
		tmp = x_offsets[j]; x_offsets[j]=x_offsets[k]; x_offsets[k]=tmp;
		}

	   for(i=0; i<N; ++i) {
		y=y_offsets[i];
		for(j=0; j<N; ++j) {
		   x = x_offsets[j];
		   if (y==y_off && x==x_off) continue;
		   v = imagediff_8(ref,img,10,10,w-10,h-10,x,y,min+1.0);
		   if (v<min) { min = v; X = x; Y = y; }
		   }
		}

	   // extend on edges
	   if (0)
	   if (X==x_off-d || Y==y_off-d || X==x_off+d || Y==y_off+d) {
	      x_off = X; y_off = Y;
	      for(y=y_off-1; y<=y_off+1; ++y) for(x=x_off-1; x<=x_off+1; ++x) {
		if (y==y_off && x==x_off) continue;
		v = imagediff_8(ref,img,10,10,w-10,h-10,x,y,min+1.0);
		if (v<min) { min = v; X = x; Y = y; }
		}
	      }

	   first=0;
	   x_off = X;
	   y_off = Y;
	   DestroyImage(ref);
	   DestroyImage(img);

//printf("%d: x_off=%d y_off=%d\n",l,x_off,y_off);
	   }

	if (last_x_off != 999) {
	   // Clamp to make sure we don't go more than 10 pixels from last frame
	   int d = last_x_off - x_off;
	   if (d>10) d=10; if (d<-10) d=-10; 
	   x_off = last_x_off - d;

	   d = last_y_off - y_off;
	   if (d>10) d=10; if (d<-10) d=-10; 
	   y_off = last_y_off - d;
	   }

	// 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;

	translate_image(RealImage,-x_off,-y_off);
	last_x_off = x_off;
	last_y_off = y_off;
	return 1;
	}

static 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;
	}

// 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 width = ref->width;
	int height = ref->height;
	int x,y,count=0,diff;
	int o,off,p1,p2;
	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);
	   }

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

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

		   d += diff * diff;
		   ++count;
		}
	   }

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

	return d/(double)count;
	}

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

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

	for(j=0,y=y1; y<=y2; ++y) {
	   o = y*w + x1;
	   for(x=x1; x<=x2; ++x,++o,++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 = 1.0;
		wlist[i].ptr = pt;
		if (Morph_Debug&1) draw_point(img,x,y,(int)(wlist[i].wgt * 255));
		}
	   }
	return 1;
	}

// Add the influence of 2 control points. p1 decreases left->right, p2 increases left -> right
static int
add_dual_point_h(struct Image *img, struct weight *mesh,int w, struct morph_point *p1, struct morph_point *p2)
	{
	int o,x,y,i,j;
	struct weight *wlist;
	int x1 = p1->x - p1->w/2;
	int x2 = x1 + p1->w - 1;
	int y1 = p1->y - p1->h/2;
	int y2 = y1 + p1->h - 1;
	static double *wgt_buffer1 = NULL;
	static double *wgt_buffer2 = NULL;
	static int wgt_w=-1,wgt_h=-1;

	// adjust (x1,y1) and (x2,y2) to cover tile between p1 and p2
	x1 += p1->w;
	x2 += p1->w;

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

	if (wgt_w != p1->w || wgt_h != p1->h) {
	   wgt_w = p1->w;
	   wgt_h = p1->h;
	   if (wgt_buffer1) {free(wgt_buffer1); free(wgt_buffer2);}
	   wgt_buffer1 = create_weight_buffer_decrease_lr(wgt_w,wgt_h);
	   wgt_buffer2 = create_weight_buffer_increase_lr(wgt_w,wgt_h);
	   }

	for(j=0,y=y1; y<=y2; ++y) {
	   o = y*w + x1;
	   for(x=x1; x<=x2; ++x,++o,++j) {
		wlist = mesh + o*4;

		for(i=0; i<4 && wlist[i].active; ++i); if (i>2) {
		   Print("Mesh overfull at (%d,%d)\n",x,y);
		   exit(1);
		   }

		wlist[i].active=1;
		wlist[i].wgt = wgt_buffer1[j]; // wgt_buffer[j];
		wlist[i].ptr = p1;
		if (Morph_Debug&2) draw_point(img,x,y,(int)(wlist[i].wgt * 255));

		wlist[++i].active=1;
		wlist[i].wgt = wgt_buffer2[j]; // wgt_buffer[j];
		wlist[i].ptr = p2;
		if (Morph_Debug&4) draw_point(img,x,y,(int)(wlist[i].wgt * 255));
		}
	   }
	return 1;
	}

// Add the influence of 2 control points. p1 decreases top->bottom, p2 decreases bottom -> top
static int
add_dual_point_v(struct Image *img, struct weight *mesh,int w, struct morph_point *p1, struct morph_point *p2)
	{
	int o,x,y,i,j;
	struct weight *wlist;
	int x1 = p1->x - p1->w/2;
	int x2 = x1 + p1->w - 1;
	int y1 = p1->y - p1->h/2;
	int y2 = y1 + p1->h - 1;
	static double *wgt_buffer1 = NULL;
	static double *wgt_buffer2 = NULL;
	static int wgt_w=-1,wgt_h=-1;

	// adjust (x1,y1) and (x2,y2) to cover tile between p1 and p2
	y1 += p1->h;
	y2 += p1->h;

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

	if (wgt_w != p1->w || wgt_h != p1->h) {
	   wgt_w = p1->w;
	   wgt_h = p1->h;
	   if (wgt_buffer1) {free(wgt_buffer1); free(wgt_buffer2);}
	   wgt_buffer1 = create_weight_buffer_decrease_ud(wgt_w,wgt_h);
	   wgt_buffer2 = create_weight_buffer_increase_ud(wgt_w,wgt_h);
	   }

	for(j=0,y=y1; y<=y2; ++y) {
	   o = y*w + x1;
	   for(x=x1; x<=x2; ++x,++o,++j) {
		wlist = mesh + o*4;

		for(i=0; i<4 && wlist[i].active; ++i); if (i>2) {
		   Print("Mesh overfull at (%d,%d)\n",x,y);
		   exit(1);
		   }

		wlist[i].active=1;
		wlist[i].wgt = wgt_buffer1[j]; // wgt_buffer[j];
		wlist[i].ptr = p1;
		if (Morph_Debug&8) draw_point(img,x,y,wlist[i].wgt * 255);

		wlist[++i].active=1;
		wlist[i].wgt = wgt_buffer2[j]; // wgt_buffer[j];
		wlist[i].ptr = p2;
		if (Morph_Debug&16) draw_point(img,x,y,wlist[i].wgt * 255);
		}
	   }
	return 1;
	}

// Add the influence of 2 control points. p1 decreases top->bottom, p2 decreases bottom -> top
static int
add_above(struct Image *img, struct weight *mesh,int w, struct morph_point *p1)
	{
	int o,x,y,i,j;
	struct weight *wlist;
	int x1 = p1->x - p1->w/2;
	int x2 = x1 + p1->w - 1;
	int y1 = p1->y - p1->h/2;
	int y2 = y1 + p1->h - 1;
	static double *wgt_buffer1 = NULL;
	static int wgt_w=-1,wgt_h=-1;

	// adjust (x1,y1) and (x2,y2) to cover tile between p1 and p2
	y1 -= p1->h;
	y2 -= p1->h;

	if (wgt_w != p1->w || wgt_h != p1->h) {
	   wgt_w = p1->w;
	   wgt_h = p1->h;
	   if (wgt_buffer1) free(wgt_buffer1);
	   wgt_buffer1 = create_weight_buffer_increase_ud(wgt_w,wgt_h);
	   }

	for(y=y1; y<=y2; ++y) {
	   if (y<0 || y>=img->height) continue;
	   o = y*w + x1; j=(y-y1)*wgt_w;
	   for(x=x1; x<=x2; ++x,++o,++j)
		if (x>=0 && x<img->width) {
		   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 = wgt_buffer1[j]; // wgt_buffer[j];
		   wlist[i].ptr = p1;
		   if (Morph_Debug&16) draw_point(img,x,y,wlist[i].wgt * 255);
		   }
	   }
	return 1;
	}

// Add the influence of 2 control points. p1 decreases top->bottom, p2 decreases bottom -> top
static int
add_below(struct Image *img, struct weight *mesh,int w, struct morph_point *p1)
	{
	int o,x,y,i,j;
	struct weight *wlist;
	int x1 = p1->x - p1->w/2;
	int x2 = x1 + p1->w - 1;
	int y1 = p1->y - p1->h/2;
	int y2 = y1 + p1->h - 1;
	static double *wgt_buffer1 = NULL;
	static int wgt_w=-1,wgt_h=-1;

	// adjust (x1,y1) and (x2,y2) to cover tile between p1 and p2
	y1 += p1->h;
	y2 += p1->h;

	if (wgt_w != p1->w || wgt_h != p1->h) {
	   wgt_w = p1->w;
	   wgt_h = p1->h;
	   if (wgt_buffer1) free(wgt_buffer1);
	   wgt_buffer1 = create_weight_buffer_decrease_ud(wgt_w,wgt_h);
	   }

	for(y=y1; y<=y2; ++y) {
	   if (y<0 || y>=img->height) continue;
	   o = y*w + x1; j=(y-y1)*wgt_w;
	   for(x=x1; x<=x2; ++x,++o,++j)
		if (x>=0 && x<img->width) {
		   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 = wgt_buffer1[j]; // wgt_buffer[j];
		   wlist[i].ptr = p1;
		   if (Morph_Debug&8) draw_point(img,x,y,wlist[i].wgt * 255);
		   }
	   }
	return 1;
	}

static int
add_above_corners(struct Image *img, struct weight *mesh,int w, struct morph_point *p1, struct morph_point *p2)
	{
	int o,x,y,i,j;
	struct weight *wlist;
	int x1 = p1->x - p1->w/2;
	int x2 = x1 + p1->w - 1;
	int y1 = p1->y - p1->h/2;
	int y2 = y1 + p1->h - 1;
	static double *wgt_buffer1 = NULL;
	static double *wgt_buffer2 = NULL;
	static int wgt_w=-1,wgt_h=-1;

	// adjust (x1,y1) and (x2,y2) to cover tile above centre of p1 and p2
	x1 += p1->w; x2 += p1->w;
	y1 -= p1->h; y2 -= p1->h;

	if (wgt_w != p1->w || wgt_h != p1->h) {
	   wgt_w = p1->w;
	   wgt_h = p1->h;
	   if (wgt_buffer1) { free(wgt_buffer1); free(wgt_buffer2);}
	   wgt_buffer1 = create_weight_buffer_decrease_bl(wgt_w,wgt_h);
	   wgt_buffer2 = create_weight_buffer_decrease_br(wgt_w,wgt_h);
	   }

	for(y=y1; y<=y2; ++y) {
	   if (y<0 || y>=img->height) continue;
	   o = y*w + x1; j=(y-y1)*wgt_w;
	   for(x=x1; x<=x2; ++x,++o,++j)
		if (x>=0 && x<img->width) {
		   wlist = mesh + o*4;

		   for(i=0; i<4 && wlist[i].active; ++i); if (i>2) {
		      Print("Mesh overfull at (%d,%d)\n",x,y);
		      exit(1);
		      }

		   wlist[i].active=1;
		   wlist[i].wgt = wgt_buffer1[j]; // wgt_buffer[j];
		   wlist[i].ptr = p1;
		   if (Morph_Debug&128) draw_point(img,x,y,wlist[i].wgt * 255);

		   wlist[++i].active=1;
		   wlist[i].wgt = wgt_buffer2[j]; // wgt_buffer[j];
		   wlist[i].ptr = p2;
		   if (Morph_Debug&256) draw_point(img,x,y,wlist[i].wgt * 255);
		   }
	   }
	return 1;
	}

static int
add_below_corners(struct Image *img, struct weight *mesh,int w, struct morph_point *p1, struct morph_point *p2)
	{
	int o,x,y,i,j;
	struct weight *wlist;
	int x1 = p1->x - p1->w/2;
	int x2 = x1 + p1->w - 1;
	int y1 = p1->y - p1->h/2;
	int y2 = y1 + p1->h - 1;
	static double *wgt_buffer1 = NULL;
	static double *wgt_buffer2 = NULL;
	static int wgt_w=-1,wgt_h=-1;

	// adjust (x1,y1) and (x2,y2) to cover tile above centre of p1 and p2
	x1 += p1->w; x2 += p1->w;
	y1 += p1->h; y2 += p1->h;

	if (wgt_w != p1->w || wgt_h != p1->h) {
	   wgt_w = p1->w;
	   wgt_h = p1->h;
	   if (wgt_buffer1) { free(wgt_buffer1); free(wgt_buffer2);}
	   wgt_buffer1 = create_weight_buffer_decrease_tl(wgt_w,wgt_h);
	   wgt_buffer2 = create_weight_buffer_decrease_tr(wgt_w,wgt_h);
	   }

	for(y=y1; y<=y2; ++y) {
	   if (y<0 || y>=img->height) continue;
	   o = y*w + x1; j=(y-y1)*wgt_w;
	   for(x=x1; x<=x2; ++x,++o,++j)
		if (x>=0 && x<img->width) {
		   wlist = mesh + o*4;

		   for(i=0; i<4 && wlist[i].active; ++i); if (i>2) {
		      Print("Mesh overfull at (%d,%d)\n",x,y);
		      exit(1);
		      }

		   wlist[i].active=1;
		   wlist[i].wgt = wgt_buffer1[j]; // wgt_buffer[j];
		   wlist[i].ptr = p1;
		   if (Morph_Debug&32) draw_point(img,x,y,wlist[i].wgt * 255);

		   wlist[++i].active=1;
		   wlist[i].wgt = wgt_buffer2[j]; // wgt_buffer[j];
		   wlist[i].ptr = p2;
		   if (Morph_Debug&64) draw_point(img,x,y,wlist[i].wgt * 255);
		   }
	   }
	return 1;
	}

static int
add_left_above(struct Image *img, struct weight *mesh,int w, struct morph_point *p1)
	{
	int o,x,y,i,j;
	struct weight *wlist;
	int x1 = p1->x - p1->w/2;
	int x2 = x1 + p1->w - 1;
	int y1 = p1->y - p1->h/2;
	int y2 = y1 + p1->h - 1;
	static double *wgt_buffer1 = NULL;
	static int wgt_w=-1,wgt_h=-1;

	// adjust (x1,y1) and (x2,y2) to cover tile between p1 and p2
	x1 -= p1->w; x2 -= p1->w;
	y1 -= p1->h; y2 -= p1->h;

	if (wgt_w != p1->w || wgt_h != p1->h) {
	   wgt_w = p1->w;
	   wgt_h = p1->h;
	   if (wgt_buffer1) free(wgt_buffer1);
	   wgt_buffer1 = create_weight_buffer_decrease_br(wgt_w,wgt_h);
	   }

	for(y=y1; y<=y2; ++y) {
	   if (y<0 || y>=img->height) continue;
	   o = y*w + x1; j=(y-y1)*wgt_w;
	   for(x=x1; x<=x2; ++x,++o,++j)
		if (x>=0 && x<img->width) {
		   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 = wgt_buffer1[j]; // wgt_buffer[j];
		   wlist[i].ptr = p1;
		   if (Morph_Debug&4) draw_point(img,x,y,wlist[i].wgt * 255);
		   }
	   }

	return 1;
	}

static int
add_left_below(struct Image *img, struct weight *mesh,int w, struct morph_point *p1)
	{
	int o,x,y,i,j;
	struct weight *wlist;
	int x1 = p1->x - p1->w/2;
	int x2 = x1 + p1->w - 1;
	int y1 = p1->y - p1->h/2;
	int y2 = y1 + p1->h - 1;
	static double *wgt_buffer1 = NULL;
	static int wgt_w=-1,wgt_h=-1;

	// adjust (x1,y1) and (x2,y2) to cover tile between p1 and p2
	x1 -= p1->w; x2 -= p1->w;
	y1 += p1->h; y2 += p1->h;

	if (wgt_w != p1->w || wgt_h != p1->h) {
	   wgt_w = p1->w;
	   wgt_h = p1->h;
	   if (wgt_buffer1) free(wgt_buffer1);
	   wgt_buffer1 = create_weight_buffer_decrease_tr(wgt_w,wgt_h);
	   }

	for(y=y1; y<=y2; ++y) {
	   if (y<0 || y>=img->height) continue;
	   o = y*w + x1; j=(y-y1)*wgt_w;
	   for(x=x1; x<=x2; ++x,++o,++j)
		if (x>=0 && x<img->width) {
		   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 = wgt_buffer1[j]; // wgt_buffer[j];
		   wlist[i].ptr = p1;
		   if (Morph_Debug&4) draw_point(img,x,y,wlist[i].wgt * 255);
		   }
	   }

	return 1;
	}

static int
add_right_above(struct Image *img, struct weight *mesh,int w, struct morph_point *p1)
	{
	int o,x,y,i,j;
	struct weight *wlist;
	int x1 = p1->x - p1->w/2;
	int x2 = x1 + p1->w - 1;
	int y1 = p1->y - p1->h/2;
	int y2 = y1 + p1->h - 1;
	static double *wgt_buffer1 = NULL;
	static int wgt_w=-1,wgt_h=-1;

	// adjust (x1,y1) and (x2,y2) to cover tile between p1 and p2
	x1 += p1->w; x2 += p1->w;
	y1 -= p1->h; y2 -= p1->h;

	if (wgt_w != p1->w || wgt_h != p1->h) {
	   wgt_w = p1->w;
	   wgt_h = p1->h;
	   if (wgt_buffer1) free(wgt_buffer1);
	   wgt_buffer1 = create_weight_buffer_decrease_bl(wgt_w,wgt_h);
	   }

	for(y=y1; y<=y2; ++y) {
	   if (y<0 || y>=img->height) continue;
	   o = y*w + x1; j=(y-y1)*wgt_w;
	   for(x=x1; x<=x2; ++x,++o,++j)
		if (x>=0 && x<img->width) {
		   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 = wgt_buffer1[j]; // wgt_buffer[j];
		   wlist[i].ptr = p1;
		   if (Morph_Debug&4) draw_point(img,x,y,wlist[i].wgt * 255);
		   }
	   }

	return 1;
	}

static int
add_right_below(struct Image *img, struct weight *mesh,int w, struct morph_point *p1)
	{
	int o,x,y,i,j;
	struct weight *wlist;
	int x1 = p1->x - p1->w/2;
	int x2 = x1 + p1->w - 1;
	int y1 = p1->y - p1->h/2;
	int y2 = y1 + p1->h - 1;
	static double *wgt_buffer1 = NULL;
	static int wgt_w=-1,wgt_h=-1;

	// adjust (x1,y1) and (x2,y2) to cover tile between p1 and p2
	x1 += p1->w; x2 += p1->w;
	y1 += p1->h; y2 += p1->h;

	if (wgt_w != p1->w || wgt_h != p1->h) {
	   wgt_w = p1->w;
	   wgt_h = p1->h;
	   if (wgt_buffer1) free(wgt_buffer1);
	   wgt_buffer1 = create_weight_buffer_decrease_tl(wgt_w,wgt_h);
	   }

	for(y=y1; y<=y2; ++y) {
	   if (y<0 || y>=img->height) continue;
	   o = y*w + x1; j=(y-y1)*wgt_w;
	   for(x=x1; x<=x2; ++x,++o,++j)
		if (x>=0 && x<img->width) {
		   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 = wgt_buffer1[j]; // wgt_buffer[j];
		   wlist[i].ptr = p1;
		   if (Morph_Debug&4) draw_point(img,x,y,wlist[i].wgt * 255);
		   }
	   }

	return 1;
	}

static int
add_left(struct Image *img, struct weight *mesh,int w, struct morph_point *p1)
	{
	int o,x,y,i,j;
	struct weight *wlist;
	int x1 = p1->x - p1->w/2;
	int x2 = x1 + p1->w - 1;
	int y1 = p1->y - p1->h/2;
	int y2 = y1 + p1->h - 1;
	static double *wgt_buffer1 = NULL;
	static int wgt_w=-1,wgt_h=-1;

	// adjust (x1,y1) and (x2,y2) to cover tile between p1 and p2
	x1 -= p1->w;
	x2 -= p1->w;

	if (wgt_w != p1->w || wgt_h != p1->h) {
	   wgt_w = p1->w;
	   wgt_h = p1->h;
	   if (wgt_buffer1) free(wgt_buffer1);
	   wgt_buffer1 = create_weight_buffer_increase_lr(wgt_w,wgt_h);
	   }

	for(y=y1; y<=y2; ++y) {
	   if (y<0 || y>=img->height) continue;
	   o = y*w + x1; j=(y-y1)*wgt_w;
	   for(x=x1; x<=x2; ++x,++o,++j)
		if (x>=0 && x<img->width) {
		   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 = wgt_buffer1[j]; // wgt_buffer[j];
		   wlist[i].ptr = p1;
		   if (Morph_Debug&4) draw_point(img,x,y,wlist[i].wgt * 255);
		   }
	   }

	add_left_above(img,mesh,w,p1);
	add_left_below(img,mesh,w,p1);
	return 1;
	}

static int
add_right(struct Image *img, struct weight *mesh,int w, struct morph_point *p1)
	{
	int o,x,y,i,j;
	struct weight *wlist;
	int x1 = p1->x - p1->w/2;
	int x2 = x1 + p1->w - 1;
	int y1 = p1->y - p1->h/2;
	int y2 = y1 + p1->h - 1;
	static double *wgt_buffer1 = NULL;
	static int wgt_w=-1,wgt_h=-1;

	// adjust (x1,y1) and (x2,y2) to cover tile to the right of p1
	x1 += p1->w; x2 += p1->w;

	if (wgt_w != p1->w || wgt_h != p1->h) {
	   wgt_w = p1->w; wgt_h = p1->h;
	   if (wgt_buffer1) free(wgt_buffer1);
	   wgt_buffer1 = create_weight_buffer_decrease_lr(wgt_w,wgt_h);
	   }

	for(y=y1; y<=y2; ++y) {
	   if (y<0 || y>=img->height) continue;
	   o = y*w + x1; j=(y-y1)*wgt_w;
	   for(x=x1; x<=x2; ++x,++o,++j)
		if (x>=0 && x<img->width) {
		   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 = wgt_buffer1[j]; // wgt_buffer[j];
		   wlist[i].ptr = p1;
		   if (Morph_Debug&8) draw_point(img,x,y,wlist[i].wgt * 255);
		   }
	   }

	add_right_above(img,mesh,w,p1);
	add_right_below(img,mesh,w,p1);
	return 1;
	}

// Add the influence of 4 control points in a square, all decreasing to the centre
static int
add_quad_point(struct Image *img, struct weight *mesh,int w, struct morph_point *p1, struct morph_point *p2,
		struct morph_point *p3, struct morph_point *p4)
	{
	int o,x,y,i,j;
	struct weight *wlist;
	int x1 = p1->x - p1->w/2;
	int x2 = x1 + p1->w - 1;
	int y1 = p1->y - p1->h/2;
	int y2 = y1 + p1->h - 1;
	static double *wgt_buffer1 = NULL;
	static double *wgt_buffer2 = NULL;
	static double *wgt_buffer3 = NULL;
	static double *wgt_buffer4 = NULL;
	static int wgt_w=-1,wgt_h=-1;

	// adjust (x1,y1) and (x2,y2) to cover tile in the centre of (p1,p2,p3,p4)
	x1 += p1->w; x2 += p1->w;
	y1 += p1->h; y2 += p1->h;

	if (! Clip(img->width, img->height,&x1,&y1,&x2,&y2))
	   return(0);

	if (wgt_w != p1->w || wgt_h != p1->h) {
	   wgt_w = p1->w;
	   wgt_h = p1->h;
	   if (wgt_buffer1) {free(wgt_buffer1); free(wgt_buffer2); free(wgt_buffer3); free(wgt_buffer4);}
	   wgt_buffer1 = create_weight_buffer_decrease_tl(wgt_w,wgt_h);
	   wgt_buffer2 = create_weight_buffer_decrease_tr(wgt_w,wgt_h);
	   wgt_buffer3 = create_weight_buffer_decrease_bl(wgt_w,wgt_h);
	   wgt_buffer4 = create_weight_buffer_decrease_br(wgt_w,wgt_h);
	   }

	for(j=0,y=y1; y<=y2; ++y) {
	   o = y*w + x1;
	   for(x=x1; x<=x2; ++x,++o,++j) {
		wlist = mesh + o*4;

		for(i=0; i<4 && wlist[i].active; ++i); if (i) {
		   Print("Mesh overfull at (%d,%d)\n",x,y);
		   exit(1);
		   }

		wlist[i].active=1;
		wlist[i].wgt = wgt_buffer1[j]; // wgt_buffer[j];
		wlist[i].ptr = p1;
		if (Morph_Debug&32) draw_point(img,x,y,(int)(wlist[i].wgt * 255));

		wlist[++i].active=1;
		wlist[i].wgt = wgt_buffer2[j]; // wgt_buffer[j];
		wlist[i].ptr = p2;
		if (Morph_Debug&64) draw_point(img,x,y,(int)(wlist[i].wgt * 255));

		wlist[++i].active=1;
		wlist[i].wgt = wgt_buffer3[j]; // wgt_buffer[j];
		wlist[i].ptr = p3;
		if (Morph_Debug&128) draw_point(img,x,y,(int)(wlist[i].wgt * 255));

		wlist[++i].active=1;
		wlist[i].wgt = wgt_buffer4[j]; // wgt_buffer[j];
		wlist[i].ptr = p4;
		if (Morph_Debug&256) draw_point(img,x,y,(int)(wlist[i].wgt * 255));
		}
	   }
	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 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++;
		}

	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) {
		x = plist[p].x + X;
		y = plist[p].y + Y;

		v = imagediff_8(ref,img,x1,y1,x2,y2,x,y,min);
		if (v<min) { min = v; nx = x; ny = y;}
		}

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

	free(plist);
	return changed;
	}

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 += (int)(x_offset+0.5);
	y += (int)(y_offset+0.5);

	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;
	   }

	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;
	int minx,miny,maxx,maxy;
	unsigned char pixel[3];

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

	minx=m->x; maxx=minx;
	miny=m->y; maxy=miny;

	// Reduce the width and height of all tiles to half their original value
	for(i=0; i<npoints; ++i) {
	   pt = m + i;
	   pt->w /=2;
	   pt->h /=2;
	   if (pt->x > maxx) maxx=pt->x;
	   if (pt->y > maxy) maxy=pt->y;
	   }

	// each control point is authoritative for the central half of its original 
	// area.
	for(i=0; i<npoints; ++i) {
	   pt = m + i;
	   if (pt->w || pt->h)
	      add_single_point(img,mesh,img->width,pt);
	   }

	// determine the increment between vertically adjacent control points
	for(v=1; v<npoints; ++v) {
	   pt = m + v;
	   if (m->w && pt->w && m->x == pt->x) break;
	   }

	// Add the dual points (between control points)
	for(i=0; i<npoints-1; ++i) {
	   pt = m + i;
	   // adjacent horizontal points
	   if (pt->w && pt->y == (pt+1)->y) {
		if (pt->y == miny) 
		   // top row, so also add top corners
		   add_above_corners(img,mesh,img->width,pt,pt+1);
		if (pt->y == maxy) 
		   // top row, so also add top corners
		   add_below_corners(img,mesh,img->width,pt,pt+1);

		add_dual_point_h(img,mesh,img->width,pt,pt+1);
		}
	   // adjacent vertical points
	   if (i+v<npoints && pt->w && (pt+v)->w && pt->x == (pt+v)->x) {
		add_dual_point_v(img,mesh,img->width,pt,pt+v);
		}

	   // top row - complete above
	   if (pt->y == miny) add_above(img,mesh,img->width,pt);
	   if (pt->y == maxy) add_below(img,mesh,img->width,pt);

	   if (pt->x == minx) add_left(img,mesh,img->width,pt);
	   if (pt->x == maxx) add_right(img,mesh,img->width,pt);
	   }

	// Add the quad points the fall between 4 control points
	for(i=0;i<npoints-v-1; ++i) {
	   pt = m + i;
	   // points that form a square
	   int x1 = pt->x, y1=pt->y;
	   int x2 = (pt+1)->x, y2=(pt+1)->y;
	   int x3 = (pt+v)->x, y3=(pt+v)->y;
	   int x4 = (pt+v+1)->x, y4=(pt+v+1)->y;

	   //Print("(%d,%d) (%d,%d) (%d,%d) (%d,%d)\n",x1,y1,x2,y2,x3,y3,x4,y4);
	   if (x1 && x2 && x3 && x4 && (x1 == x3) && (y1 == y2))
		add_quad_point(img,mesh,img->width,pt,pt+1,pt+v,pt+v+1);
	   }

if (Morph_Debug) return (1);

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

	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);

	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. Set to 20%
	cutoff = (w * h) * 0.1 ;

	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;

//printf("(%d,%d) = %d\n",p[n].x,p[n].y,p[n].active);
	   	}
		
	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 (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;
	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);
	      }

	// convert input image to 6 bit for alignment purposes
	for(i=0; i<img->width * img->height; ++i) img->data[i] &= 0xf4;

	stretch_to_ref(ref,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 (DisplayFrames) {
   	   struct Image *I = ConvertImage(img,IMG_FIT,16);
   	   ShowImage(I);
   	   //Usleep(1000000);
	   DestroyImage(I);
   	   }

	// 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;
	}

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;
	   }

	// 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 (level==3 && ! RoughAlign(I)) {
	   Print("Final Alignment failed\n");
	   return 0;
	   }

	return 1;
	}
