#include "ninox.h"

static unsigned int max_histo_8(struct Image *img);
static unsigned int max_histo_16(struct Image *img);

static struct Image * ConvertToBMP(struct Image *img, int newdepth);
static struct Image * ConvertToFITS(struct Image *img, int newdepth);

static inline void
ByteSwapC(unsigned char *ptr)
	{
	unsigned char x = *ptr;
	*ptr = *(ptr+1);
	*(ptr+1)=x;
	}

static inline void
ByteSwapS(unsigned short *sptr)
	{
	unsigned char *ptr = (unsigned char *)sptr;
	unsigned char x = *ptr;
	*ptr = *(ptr+1);
	*(ptr+1)=x;
	}

static inline u32
ByteSwapLV(u32 in)
	{
	register u32 out= (in>>24) & 0x000000ff;

	out |= (in>>8)  & 0x0000ff00;
	out |= (in<<8)  & 0x00ff0000;
	out |= (in<<24) & 0xff000000;

	return out;
	}
	
static inline void
ByteSwapL(u32 *sptr)
	{
	unsigned char *ptr = (unsigned char *)sptr;
	int x;

	x = *ptr; *ptr = *(ptr+3); *(ptr+3)=x;
	x = *(ptr+1); *(ptr+1) = *(ptr+2); *(ptr+2) = x;
	}

//******************************************************************************
//
// wrap up file i/o where we try a few times in case the data is coming from an 
// unusual source that occasionally fails

static int TotalFreadBytes=0;

int ReadByte(FILE *in) 
	{
	int fd = fileno(in);
	int err;
	char ch;

	err = read(fd,&ch,1);

	if (err==1) {
	   ++TotalFreadBytes;
	   return ch;
	   }

	// some other error
	return EOF;
	}

int
Fread(unsigned char *buffer, int size, int count, FILE *in)
	{
        int failures=0;
	int fd = fileno(in);
	int bytes_remaining = size * count;

	//printf("reading count=%d size=%d from %0x offset %d total=%d\n",count,size,(unsigned int)in,ftell(in),TotalFreadBytes);
	//fflush(stdout);

	while(bytes_remaining && failures < 5) {
	   int n = read(fd,buffer,bytes_remaining);
//printf("Read %d of %d bytes\n",n,bytes_remaining); fflush(stdout);
	   if (n>0) {
	      buffer += n;
	      bytes_remaining -= n;
	      TotalFreadBytes += n;
	      }
	   else { failures++; }
	   }

	if (! bytes_remaining) return count;

	Print("read: error, %d bytes remaining unread from %d, could not read from input at offset %d\n",
		bytes_remaining,size*count,TotalFreadBytes);

	return 0;
	}

//******************************************************************************

#define NBUCKETS_32 	4096
#define NBUCKETS_16 	656
#define NBUCKETS_8	256

struct Image *
CreateImage(void)
	{
	struct Image *img = (struct Image *)ZeroMalloc(sizeof(struct Image));
	return img;
	}

// Create a FITS image from raw data
struct Image *
CreateFitsImage(char *dst_fname, int width, int height, int depth,
		unsigned char *data, u32 bzero, double bscale)
	{
	struct Image *img = CreateImage();
	struct fits_info *fi = ZeroMalloc(sizeof(struct fits_info));

	img->type=IMG_FIT;
	img->width  = img->dst_width  = width;
	img->height = img->dst_height = height;
	img->depth = depth;

	if (bzero>=0) fi->bzero = bzero;
	else {
	   if (depth==8) fi->bzero=0;
	   if (depth==16)fi->bzero=32768;
	   if (depth==32)fi->bzero=0;
	   }

	if (bscale>0) fi->bscale = bscale;
	else fi->bscale=1;

	fi->header = NULL;

	img->info = fi;
	img->data = data;

	img->src_fname = strdup("<internal>");
	img->dst_fname = strdup(dst_fname);

	img->cutout.x1 = img->cutout.y1 = 0;
	img->cutout.x2 = width-1;
	img->cutout.y2 = height-1;

	return img;
	}

int
SetCutoutRegion(struct Image *img, int x1, int y1, int x2, int y2)
	{
	img->cutout.x1 = x1;
	img->cutout.y1 = y1;
	img->cutout.x2 = x2;
	img->cutout.y2 = y2;
	return 1;
	}

static struct Image *
convert_to_bmp(struct Image *src, int depth)
	{
	Print("convert_to_bmp: unsupported\n");
	exit(1);
	}

static struct Image *
convert_to_fits(struct Image *img, int depth)
	{
	struct Image *dst;
	struct fits_info *fi = img->info;
	unsigned char *src,*data;
	unsigned short *udata,*usrc;
	unsigned int v,*uisrc;
	int i,o,npixels = img->width * img->height;
	int bpp = depth/8;
	char fname[128];
	
	src = img->data;
	usrc = (unsigned short *)src;
	uisrc = (unsigned int *)src;

	data = Malloc(img->width * img->height * bpp);
	udata = (unsigned short *)data;

	if (depth==img->depth) memcpy(data,src,npixels*bpp);
	else switch(depth) {
	   case 8:
		if (img->depth==16) for(i=0; i<npixels; ++i) data[i] = usrc[i]>>8;
		else if (img->depth==24)
		   for(i=o=0; i<npixels; ++i) 
			data[i] = (unsigned)(src[o++]*0.114 + src[o++]*0.587 + src[o++]*0.299);
		else {
		   Print("convert_to_fits: source image depth %d not supported\n",img->depth);
		   exit(1);
		   }
		break;
	   case 16:
		if (img->depth==8) for(i=0; i<npixels; ++i) udata[i] = (unsigned short)(src[i])<<8;
		else if (img->depth==24) // BGR from BMP format
		   for(i=o=0; i<npixels; ++i) {
			unsigned short B = src[o++]; B<<=8; B *= 0.114;
			unsigned short G = src[o++]; G<<=8; G *= 0.587;
			unsigned short R = src[o++]; R<<=8; R *= 0.299;
			udata[i] = R + G + B;
			}
		else if (img->depth == 32) // 32 bit FITS
		   for(i=o=0; i<npixels; ++i) udata[i] = uisrc[i] >> 16;
		else {
		   Print("convert_to_fits: source image depth %d not supported\n",img->depth);
		   exit(1);
		   }
		break;
	   case 24:
		if (img->depth==8)       for(i=o=0; i<npixels; ++i) {v=src[i];     data[o++]=v; data[o++]=v; data[o++]=v;}
		else if (img->depth==16) for(i=o=0; i<npixels; ++i) {v=usrc[i]>>8; data[o++]=v; data[o++]=v; data[o++]=v;}
		else {
		   Print("convert_to_fits: source image depth %d not supported\n",img->depth);
		   exit(1);
		   }
		break;
	   }

	sprintf(fname,"convert_%d_%d.fit",img->depth,depth);
	dst = CreateFitsImage(fname,img->width,img->height,depth,data,-1,-1.0);
	if (! dst) {
	   Print("convert_to_fits: Error\n");
	   exit(1);
	   }

	dst->cutout = img->cutout;
	return dst;
	}

struct Image *
ConvertImage(struct Image *src, int type, int depth)
	{
	switch(type) {
	   case IMG_FIT: return convert_to_fits(src,depth); break;
	   case IMG_BMP: return convert_to_bmp(src,depth); break;
	   default:
		Print("ConvertImage: Unsupported image type %d\n",type);
		exit(1);
	   }

	// notreached
	return 0;
	}

int
DestroyImage(struct Image *img)
	{
	struct bmp_info *bi;
	struct fits_info *fi;

	if (img==NULL) return(1);

	if (img->src_fname) free(img->src_fname);
	if (img->dst_fname) free(img->dst_fname);

        switch(img->type) {
	   case IMG_PPM:
		// Nothing to see here
		break;
	   case IMG_BMP:
		bi = (struct bmp_info *)img->info;
		if (bi && bi->cmap) free(bi->cmap);
		break;
	   case IMG_FIT:
		fi = (struct fits_info *)img->info;
		if (fi && fi->header)
		   free(fi->header);
		break;

	   default:
		Print("DestroyImage: Unknown Image type %d\n",img->type);
		exit(1);
		break;
	   }

        if (img->data) free(img->data);
        if (img->info) free(img->info);

        if (img->q) {
	   if (img->q->newname) free(img->q->newname);
	   free(img->q);
	   }

	free(img);
	return(1);
	}
		
//###################################################################################
//
// Image Reading Support Routines
//
//###################################################################################

static int
load_ppm_header(FILE *in, struct Image *img)
	{
	struct ppm_info *pi = (struct ppm_info *)ZeroMalloc(sizeof(struct ppm_info));
	int ch,i,incomment;
	char buf[32];

	pi->ppm_type[0]=0; img->width=-1; img->height=-1; pi->maxval=-1;
	incomment=0; i=0;
	while(1) {
	   switch(ch=ReadByte(in))
		{
		case EOF: return 1;
		case '#': incomment=1; break;
		case ' ': case '\t':
			if (incomment) break;
		case '\n': case '\r':
			if (incomment) {incomment=0; break;}
			if (pi->ppm_type[0]==0) {
			   strncpy(pi->ppm_type,buf,2); pi->ppm_type[2]=0;
			   if (pi->ppm_type[0]!='P' || (pi->ppm_type[1]!='3' && pi->ppm_type[1]!='6')) {
				Print("load_ppm_header: Unsupported format '%s'\n",buf);
				return 0;
			   	}
			   }
			else if (img->width<0) img->width=atoi(buf);
			else if (img->height<0) img->height=atoi(buf);
			else if (pi->maxval<0) {
				// Successfuly read the header
				pi->maxval=atoi(buf);
				img->info = pi;
				return(1);
				}
			i=0; buf[0]=0; break;
			break;
		default:
			if (i<31) buf[i++]=ch;
			buf[i]=0;
			break;
	   	}
	   }
	return 1;
	}

static int
load_fits_header(FILE *in, struct Image *img)
	{
	struct fits_info *fi;
        char hline[81],*arg,*val;
	int found_end=0;
        int i;

	fi = (struct fits_info *)ZeroMalloc(sizeof(struct fits_info));
	fi->header = (unsigned char *)ZeroMalloc(2880);

	fi->bzero = 0;
	fi->bscale = 1.0;
	img->tm_sec = 0;
	img->tm_tsc = 0;

	// FITS specific header and information
	img->info = fi;

        hline[80]=0;
        for(i=0; i<36; ++i) {
           if (Fread(hline,80,1,in) != 1) {
                Print("AlignFIT: Short read on header line %d\n",i+1);
                exit(1);
                }

	   // Save into our header
	   memcpy(fi->header+i*80, hline, 80);

           // Format is ARG = VAL
           arg=hline; while(*arg != ' ' && *arg != '=') ++arg; *arg=0;
           val = arg+1; arg=hline;

           if (!strcmp(arg,"END")) { ++i; found_end=1; break; }

           while(*val == ' ' || *val == '=') {
                ++val;
                }

           if (!strcmp(arg,"SIMPLE") && val[0] != 'T') {
                Print("AlignFIT: Only handle SIMPLE format\n");
                exit(1);
                }

           if (!strcmp(arg,"BITPIX")) {
		int d = atoi(val);
		if (d !=8 && d !=16 && d != 32) {
		   Print("FITS error: depth %d not understood\n",d);
		   exit(1);
		   }
		img->depth = d;
		}

           if (!strcmp(arg,"NAXIS") && atoi(val) != 2) {
                Print("AlignFIT: Unsupported header '%s = %s'\n",arg,val);
                exit(1);
                }
           if (!strcmp(arg,"NAXIS1")) img->width = atoi(val);
           if (!strcmp(arg,"NAXIS2")) img->height = atoi(val);
           if (!strcmp(arg,"BSCALE")) {
		double bscale = atof(val);
		if (bscale < 0.99 || bscale > 2.0) {
		   Print("FITS: BSCALE of %f is out of accepted range\n",bscale);
		   exit(1);
		   }
		fi->bscale = bscale;
		}

           if (!strcmp(arg,"BZERO")) {
		int v = (int)atof(val);
		if (v < 0 || v > 32768) {
		   Print("FITS: BZERO value of %d out of accepted range\n",v);
		   exit(1);
		   }
		fi->bzero = v;
		}

	   // Load the timestamp field
	   if (!strcasecmp(arg,"COMMENT") && !strncasecmp(val,"'__TIMESTAMP'",13)) {
		unsigned long long tsc;
		unsigned int hr,min;
		double secs;
		unsigned char *ptr = val;

		while(*ptr && ! isdigit(*ptr)) ++ptr;
#ifdef MSWIN32
		if (sscanf(ptr,"%u:%u:%lf/%I64u",&hr,&min,&secs,&tsc) == 4) {
#else
		if (sscanf(ptr,"%u:%u:%lf/%llu",&hr,&min,&secs,&tsc) == 4) {
#endif
		   img->tm_sec = hr * 3600 + min * 60 + secs;
		   img->tm_tsc = tsc;
		   // Print("[%u][%u][%lf] => %u (%lu)\n",hr,min,secs,img->tm_sec,tsc);
		   }
		}
	   }

        while(i < 36) {
	   if (Fread(hline,80,1,in) != 1) {
		Print("Error reading extra FITS header lines\n");
		exit(1);
		}

	   memcpy(fi->header+i*80, hline, 80);
           if (!strncmp(arg,"END ",4)) found_end=1;
	   ++i;
	   }

	if (! found_end) {
	   //Print("Warning: END not present in first 36 records, continuing...\n");
	   while(i < 72) {
	      if (Fread(hline,80,1,in) != 1) break;
              if (!strncmp(arg,"END ",4)) {
		found_end=1;
		//Print("END marker found on line %d of header\n",i+1);
		}
	      ++i;
	      }
	   }

	if (! found_end) {
	   Print("Error: END marker not found in first 72 lines, failed\n");
	   exit(1);
	   }

	// hack
	if (img->depth==8) fi->bzero=0;

        return 1;
        }

static int
load_bmp_header(FILE *in, struct Image *img)
	{
	struct bmp_info *bi = (struct bmp_info *)ZeroMalloc(sizeof(struct bmp_info));
	bmp_header H;
	char type[2];
	unsigned char *ptr;
	unsigned int R,G,B;
	int i,n,colour_data;

	img->width = 0;
	img->height = 0;
	bi->cmap = NULL;
	bi->cmap_entries = 0;

	img->info = bi;

	/* Look for magic */
	if (Fread(type,2,1,in) != 1) {
	   fprintf(stderr,"Short read\n");
	   exit(1);
	   }

	if (strncmp(type,"BM",2)) {
	   Print("Missing BM header - '%s' not a BMP file?\n",img->src_fname);
	   exit(1);
	   }

	/* Read the rest of the header */
	if (Fread((char *)&H,sizeof(H),1,in) != 1) {
	   Print("Short read\n");
	   exit(1);
	   }

	/* Header remaining size shoulf be 40 bytes */
	if (H.size != 40) {
	   Print("Header size %d not 40 bytes\n",H.size);
	   exit(1);
	   }

	img->depth = H.bitcount;
	bi->cmap_entries = H.cmapentries;

	if (img->depth == 24) {
	   /* Offset to image data should be 54 */
	   if (H.offset != 54) {
	      Print("24bpp: Image data offset %d not 54 bytes\n",H.offset);
	      exit(1);
	      }

	   /* sanity check */
	   if (H.planes == 1 && H.bitcount == 24) {
	      img->width = H.width;
	      img->height = H.height;
	      return 1;
	      }
	   else {
	      	Print("Error in header: planes=%d, bitcount=%d\n",H.planes,H.bitcount);
		exit(1);
	      	}
	   }

	if (img->depth == 8) {
	      // monochrome OR colour palette data

	      bi->cmap = (unsigned char *)ZeroMalloc(1024);  // 256 entries X 4 bytes per entry

	      // If there's no colourmap then the data is 8bpp monochrome
	      // provide a fake colourmap for that
	      if (H.cmapentries == 0) {
	         bi->cmap_entries = 256;
		 for(i=0; i<256; ++i) {
			bi->cmap[i*4] = i;
			bi->cmap[i*4+1] = i;
			bi->cmap[i*4+2] = i;
			}
		 }

	      /* Offset to image data should be 54  + colourmap(1024)*/
	      if (H.offset != H.cmapentries*4+54) {
	         Print("8bpp: Image data offset %d not %d bytes\n",H.offset,H.cmapentries*4+54);
	         exit(1);
	         }

	      /* Read colourmap */
	      if (H.cmapentries && Fread(bi->cmap, H.cmapentries * 4,1, in) != 1 ) {
		Print("Short read on colour map\n");
		exit(1);
	 	}

	      // to distinguish between RGB and mono indexed modes we have to
	      // look at the colourmap.
	      ptr = bi->cmap;
	      colour_data = 0;
	      for(i=0; i<256 && colour_data==0; ++i,ptr+=4) {
		   B = *(ptr);
		   G = *(ptr+1);
		   R = *(ptr+2);

		   if (R != G || R != B) colour_data=1;
		   }
	      if (colour_data) img->depth = 24;
		
	      img->info = bi;
	      img->width = H.width;
	      img->height = H.height;

	      return 1;
	      }

	Print("load_bmp_header: end of function reached\n");
	exit(1);
	}

int
load_image_header(FILE *in, struct Image *img)
	{
	if (img==NULL) {
	   Print("oops: load_image_header called with NULL ptr\n");
	   return 0;
	   }

	switch(img->type) {
	   case IMG_PPM:
		return load_ppm_header(in,img);
		break;
	   case IMG_BMP:
		return load_bmp_header(in,img);
		break;
	   case IMG_FIT:
		return load_fits_header(in,img);
		break;
	   default:
		Print("load_image_header: unknown image type %d\n",img->type);
		return 0;
		break;
	   }
	return 0;
	}

static int
load_ppm_data(FILE *in, struct Image *img)
	{
	Print("Not Implemented\n");
	return 0;
	}

static int
load_bmp_data(FILE *in, struct Image *img)
	{
	int i,n,err;
	int width = img->width;
	int height = img->height;
	int bpp = img->depth/8;
	struct bmp_info *bi = img->info;
	int rowbytes = width * bpp;		// number of active bytes in row
	int padding = ((rowbytes+3)/4)*4 - rowbytes;
        int bufsize = width * height * bpp;
        unsigned char *buffer = (unsigned char *)Malloc(bufsize);

        // If we have a colourmap, load a row at a time of 1 byte indexes and
	// fill from the colourmap. Might be 8 or 24 bpp data
        if (bi->cmap_entries > 0) {
           int x,y;
           unsigned char *rowbuf = (unsigned char *)Malloc(width);

	   for(y=height-1; y>=0; --y) {
              unsigned char *sptr=rowbuf;
	      unsigned char *dptr=buffer + y*width * bpp;

              err = Fread(rowbuf,1,width,in);
              if (err != width) {
                Print("Short read on image data\n");
                return 0;
                }

              for(x=0; x<width; ++x) {
                  n = *(sptr++);
                  if (n>=bi->cmap_entries) {
                      Print("Colourmap entry %d > colourmap size of %d\n",n,bi->cmap_entries);
                      return 0;
                      }

                  n *= 4;
                  *(dptr++) = bi->cmap[n];
		  if (bpp==3) {
                     *(dptr++) = bi->cmap[n+1];
                     *(dptr++) = bi->cmap[n+2];
		     }
                  }

	      if (padding) {
                 err = Fread(rowbuf,1,padding,in);
                 if (err != padding) {
                   Print("load_bmp_data: warning: short read on row padding\n");
                   return 0;
                   }
		}

	      }
           free(rowbuf);

	   // If the original data was 24bpp indexed then make sure
	   // we re-mark it as 24bpp RGB
	   if (bpp==3) bi->cmap_entries=0;
           }
        else {
	   unsigned char padbuffer[4];
	   int y;

	   if (bpp==1) {
		// 8bit monochrome, not indexed
	        for(y=height-1; y>=0; y--) {
		   int o = y*width;
           	   if (Fread(buffer+o,1,rowbytes,in) != rowbytes) {
              	      Print("Short read on image data\n");
              	      return 0;
              	      }
           	   if (Fread(padbuffer,1,padding,in) != padding) {
              	      Print("Short read on image data\n");
              	      return 0;
		      }
		   }
		}
	   else {
              // Data is BGR, rows padded to 4byte boundary
	      for(y=height-1; y>=0; y--) {
		int o = y*width*3;
           	if (Fread(buffer+o,1,rowbytes,in) != rowbytes) {
              	   Print("Short read on image data\n");
              	   return 0;
              	   }
           	if (Fread(padbuffer,1,padding,in) != padding) {
              	   Print("Short read on image data\n");
              	   return 0;
		   }
		}
              }
	   }

	img->data = buffer;

	return 1;
	}

// For 8/16/32 bit FITS images, invert the data after loading
// since that seems to be the way other programs do it.
int
InvertImage(struct Image *img)
	{
	int x,y;
	int w = img->width;
	int h = img->height;
	unsigned char *ptr = (unsigned char *)img->data;
	unsigned short *iptr = (unsigned short *) img->data;
	u32 *uptr = (u32 *) img->data;

	switch(img->depth) {
	   case 8:
		for(y=0; y<h/2; ++y) {
		   int o1 = y * w;
		   int o2 = (h-y-1) * w;
		   for(x=0; x<w; ++x,++o1,++o2) {
			int tmp = ptr[o1]; ptr[o1]=ptr[o2]; ptr[o2]=tmp; }
		   }
		break;
	   case 16:
		for(y=0; y<h/2; ++y) {
		   int o1 = y * w;
		   int o2 = (h-y-1) * w;
		   for(x=0; x<w; ++x,++o1,++o2) {
			int tmp = iptr[o1];
			iptr[o1]=iptr[o2];
			iptr[o2]=tmp;
			}
		   }
		break;
	   case 24:
		for(y=0; y<h/2; ++y) {
		   int o1 = y * w * 3;
		   int o2 = (h-y-1) * w * 3;
		   for(x=0; x<w; ++x,o1+=3,o2+=3) {
			int tmp;
			tmp = ptr[o1]; ptr[o1]=ptr[o2]; ptr[o2]=tmp; 
			tmp = ptr[o1+1]; ptr[o1+1]=ptr[o2+1]; ptr[o2+1]=tmp; 
			tmp = ptr[o1+2]; ptr[o1+2]=ptr[o2+2]; ptr[o2+2]=tmp; 
			}
		   }
		break;
	   case 32:
		for(y=0; y<h/2; ++y) {
		   int o1 = y * w;
		   int o2 = (h-y-1) * w;
		   for(x=0; x<w; ++x,++o1,++o2) {
			u32 tmp = uptr[o1];
			uptr[o1]=uptr[o2];
			uptr[o2]=tmp;
			}
		   }
		break;
	   default:
		break; // nothing to do
	   }
	return 1;
	}

// Flip the image left-right
int
MirrorImage(struct Image *img)
	{
	int x,y;
	int w = img->width;
	int h = img->height;
	unsigned char *ptr = (unsigned char *)img->data;
	unsigned short *iptr = (unsigned short *) img->data;
	u32 *uptr = (u32 *) img->data;

	switch(img->depth) {
	   case 8:
		for(y=0; y<h; ++y) {
		   int o1 = y * w;
		   int o2 = o1 + w-1;
		   for(x=0; x<=w/2; ++x,++o1,--o2) {
			int tmp = ptr[o1]; ptr[o1]=ptr[o2]; ptr[o2]=tmp; 
			}
		   }
		break;

	   case 16:
		for(y=0; y<h; ++y) {
		   int o1 = y * w;
		   int o2 = o1 + w-1;
		   for(x=0; x<=w/2; ++x,++o1,--o2) {
			int tmp = iptr[o1];
			iptr[o1]=iptr[o2];
			iptr[o2]=tmp;
			}
		   }
		break;

	   case 24:
		for(y=0; y<h; ++y) {
		   int o1 = y * w * 3;
		   int o2 = o1 + (w-1)*3;
		   for(x=0; x<=w/2; ++x,o1+=3,o2-=3) {
			int tmp;
			tmp = ptr[o1]; ptr[o1]=ptr[o2]; ptr[o2]=tmp; 
			tmp = ptr[o1+1]; ptr[o1+1]=ptr[o2+1]; ptr[o2+1]=tmp; 
			tmp = ptr[o1+2]; ptr[o1+2]=ptr[o2+2]; ptr[o2+2]=tmp; 
			}
		   }
		break;

	   case 32:
		for(y=0; y<h; ++y) {
		   int o1 = y * w;
		   int o2 = o1 + w - 1;
		   for(x=0; x<=w/2; ++x,++o1,--o2) {
			u32 tmp = uptr[o1];
			uptr[o1]=uptr[o2];
			uptr[o2]=tmp;
			}
		   }
		break;

	   default:
		Print("MirrorImage: Unsupported depth %d\n",img->depth);
	  	exit(0);
		break;
	   }

	return 1;
	}

// Convert a fits image of type "32" (unsigned int) to "-32" (float)
// take advantage of the fact that they are the same size elements to
// do an in-place conversion

static void
_convert_fits_32(struct Image *img)
	{
	int w = img->width;
	int h = img->height;
	int npix = w * h;
	float *fptr = (float *)img->data;
	u32 *uptr = (u32 *)img->data;
	int i;

	if (img->depth != 32) {
	   Print("convert_fits_32: Depth must be 32\n");
	   exit(1);
	   }

	img->depth = -32;

	for(i=0; i<npix; ++i) {
	   double d = uptr[i]; d /= (double)UMAX32;
	   fptr[i] = d;
	   }
	
	}

static void
_load_fixup_fits_U32(struct Image *img)
	{
	struct fits_info *fi = img->info;
        u32 *ptr,*uptr = (u32 *)img->data;
	int i,npix,do_stretch = InputHistoStretch;
	unsigned short *sbuffer;
	u32 umin,umax,bzero = fi->bzero;
	s32 smin,smax;

        npix=img->width * img->height;
	ptr=uptr; 
	umax=0; umin=UMAX32; 
	smax = -SMAX32; smin=SMAX32;

       	// Undo the endian-ness of the data for 32bpp FITS data
	// also check if the data is unsigned or signed as stored
        for(i=0; i<npix; ++i) {
           u32 v = ByteSwapLV(ptr[i]);

           v += bzero;
	   if ((s32)v>smax) smax=(s32)v;
	   if ((s32)v<smin) smin=(s32)v;
           ptr[i]=v;
           }

	// Adjust the bottom of our data to have 0 (unsigned) as the lowest value
	if (smin < 0) {
	   // Print("Data is signed 32 bit, subtracting offset %d to make unsigned\n",smin);
	   for(i=0; i<npix; ++i) ptr[i] -= smin;
	   smax -= smin;
	   smin=0;
	   }

	// Now switch to unsigned arithmetic
	umax=smax; umin=smin;

	if (do_stretch ==-1) {
	   // AUTO, so look at max to decide. Whatever we decide
	   // to do applies to all the images
	   if (umax < SMAX32) do_stretch = 1;
	   else do_stretch = 0;
	   }

	if (do_stretch == 1) {
	   double scale = (double)(UMAX32 * 0.9) / (double)umax;  // brightest pixel is 90% point

	   // if (! Quiet) Print("InputHistoStretch: max=%u, scale=%lf\n",umax,scale);

	   for(i=0; i<npix; ++i) {
		double x = uptr[i]; x *= scale;
		if (x<0) x=0; if (x>UMAX32) x=UMAX32;
		uptr[i]=(u32)x;
		}
	   }
		
	// If we want the image correct, then the FITS data must be
	// inverted.
	if (doInvertImage) InvertImage(img);

	// remove possible histprotection pixels
	uptr[0] = uptr[1] = 0;
        }

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

static void
_load_fixup_fits_16(struct Image *img, int need_byteswap, int need_invert)
	{
	struct fits_info *fi = img->info;
        unsigned char *tmp,*ptr;
        unsigned short *iptr;
	int npix,i,j,do_stretch = InputHistoStretch;
	int max,maxp[MAXP];  // The brightest non-white pixels
	int overbright = 0;  // count of overbright pixels

        ptr = img->data; iptr = (unsigned short *)ptr;
	for(i=0; i<MAXP; ++i) maxp[i] = 0;

        npix=img->width * img->height;
	i=0;

       	// Undo the endian-ness of the data for 16bpp FITS data
        while(i<npix) {

	   if (need_byteswap) {
           	int x = *ptr; *ptr = *(ptr+1); *(ptr+1) = x;
                if (fi->bzero) *iptr += fi->bzero;
		}

	   // Keep track of the brightest N pixels as part of a queue 
	   // MAXP long

	   if (*iptr >= 65472) ++overbright;

	   if (*iptr > maxp[2] && *iptr < 65530) {
	      if (*iptr>maxp[0]) {
		for(j=MAXP-1; j>0; --j) maxp[j] = maxp[j-1]; maxp[j]=*iptr;
		}
	      else if (*iptr>maxp[1]) {
		for(j=MAXP-1; j>1; --j) maxp[j] = maxp[j-1]; maxp[j]=*iptr;
		}
	      else {
		for(j=MAXP-1; j>2; --j) maxp[j] = maxp[j-1]; maxp[j]=*iptr;
		}
	      }

           ptr += 2; iptr++; i++;
           }

	if (EnableHistoWarning && overbright > 100) {
	   Print("Overbright warning - %d pixels \n",overbright);
	   img->histo_warning = 1;
	   }

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

	if (max<256) {
	   // Annoying bugs in other programs 
   	   Print("possible 8-bit unsigned data, scaling to 16bpp...\n");
   	   iptr=(unsigned short *)img->data;
   	   i=npix;
   	   while(i--) *(iptr++)<<=8;
	   max <<= 8;
   	   }

	if (do_stretch ==-1) {
	   // AUTO, so look at max to decide. Whatever we decide
	   // to do applies to all the images
	   if (max < 32000) do_stretch = 1;
	   else do_stretch = 0;
	    }

	if (do_stretch == 1) {
	   double scale = 59000.0 / max;  // brightest pixel is 90% point
	   int o = img->width * img->height - 1;
           unsigned short *iptr = (unsigned short *)img->data;

	   if (! Quiet) Print("InputHistoStretch: max=%d, scale=%lf\n",max,scale);

	   while(o>=0) {
		int x = iptr[o]; x *= scale;
		if (x<0) x=0; if (x>65535) x=65535;
		iptr[o--]=x;
		}
	   }
		
	if (IQuant>0 && IQuant<16) {
	   int o = img->width * img->height - 1;
           unsigned short *iptr = (unsigned short *)img->data;
	   unsigned short mask;

	   for(i=1,mask=0x8000; i<IQuant; ++i, mask>>=1);
	   mask = ~(mask-1);

	   while(o>=0) {
	 	int x = iptr[o];  x &= mask;
		if (x>65535) x=65535;
		iptr[o--]=x;
		}
	   }
	   
	// If we want the image correct, then the FITS data must be
	// inverted.
	if (doInvertImage && need_invert) InvertImage(img);

	// remove possible histprotection pixels
        iptr = (unsigned short *)ptr;
	iptr[0] = iptr[1] = 0;

	// set the BZERO header so that images we write
	// based on this data will be correct and unambiguous
	fi->bzero = 32768;
        }

static void
_load_fixup_fits_8(struct Image *img)
	{
	struct fits_info *fi = img->info;
        unsigned char *ptr = img->data;
	int i,j,do_stretch = InputHistoStretch;
	int max,maxp[MAXP];  // The brightest non-white pixels

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

        i = img->width * img->height;
        while(i--) {
           if (fi->bzero) *ptr += fi->bzero;
           //if (fi->bscale != 1.0) *ptr *= fi->bscale;

	   // Keep track of the brightest N pixels as part of a queue 
	   // MAXP long

	   if (*ptr > maxp[2] && *ptr < 250) {
	      if (*ptr>maxp[0]) {
		for(j=MAXP-1; j>0; --j) maxp[j] = maxp[j-1]; maxp[j]=*ptr;
		}
	      else if (*ptr>maxp[1]) {
		for(j=MAXP-1; j>1; --j) maxp[j] = maxp[j-1]; maxp[j]=*ptr;
		}
	      else {
		for(j=MAXP-1; j>2; --j) maxp[j] = maxp[j-1]; maxp[j]=*ptr;
		}
	      }

           ptr++;
           }

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

	if (do_stretch ==-1) {
		// AUTO, so look at max to decide. Whatever we decide
		// to do applies to all the images
		if (max < 128) do_stretch = 1;
		else do_stretch = 0;
	 	}

	if (do_stretch == 1) {
	  	double scale = 250.0 / max;  // brightest pixel is 90% point
		int o = img->width * img->height - 1;
           	ptr = img->data;

		if (! Quiet) Print("InputHistoStretch: max=%d, scale=%lf\n",max,scale);

		while(o>=0) {
		   int x = ptr[o]; x *= scale;
		   if (x<0) x=0; if (x>255) x=255;
		   ptr[o--]=x;
		   }
		}

	if (IQuant<8 && IQuant>0) {
	   int o = img->width * img->height - 1;
           ptr = img->data;
	   unsigned char mask;

	   for(i=1,mask=0x80; i<IQuant; ++i, mask>>=1);
	   mask = ~(mask-1);

	   while(o>=0) {
	 	int x = ptr[o];  x &= mask;
		if (x>255) x=255;
		ptr[o--]=x;
		}
	   }

	// flip to to bottom
	if (doInvertImage) InvertImage(img);

#if 0
	// remove possible hist protection pixels
	img->data[0] = img->data[1] = 0;
#endif
	}

static int
load_fits_data(FILE *in, struct Image *img, int flags)
	{
	int bufsize = img->width * img->height * img->depth/8;
	int nrecords, total, padding;
	int zbytes,nbytes;
	int need_byteswap = 1;
	int need_invert = 1;

        /*
         * Allocate the buffer and load the data
         */

	img->data = NULL;

	if (flags & ARCHIVE_COMPRESSED) {
		if (Fread((unsigned char *)&zbytes,sizeof(zbytes),1,in) != 1) {
		   Print("load_fits_data: Error reading zbytes\n");
		   exit(1);
		   }
		}

	if (flags & FITS_LOAD) {
	   int err;

           img->data = (unsigned char *)Malloc(bufsize);

	   if (flags & ARCHIVE_COMPRESSED) {
                unsigned char *buf = (unsigned char *)Malloc(zbytes);
           	if (Fread(buf,zbytes,1,in) != 1) {
              	   Print("Error in reading compressed data for %s (%d x %d)\n", img->src_fname,img->width,img->height);
              	   exit(1);
		   }

		nbytes = _ftz_decode_frame(img,buf,zbytes);
		free(buf);
		need_byteswap=0;
		need_invert=0;
		}
	   else {
		if (Fread(img->data,bufsize,1,in) != 1) {
              	   Print("Error in reading data for %s (%d x %d)\n",img->src_fname,img->width,img->height);
              	   exit(1);
		   }
		need_byteswap=1;
		need_invert = 1;
                }

	   switch(img->depth) {
		case 8:  _load_fixup_fits_8(img); break;
		case 16: _load_fixup_fits_16(img,need_byteswap,need_invert); break;
		case 32: _load_fixup_fits_U32(img); break;
		default: Print("load_fits_data: Unknown depth %d\n",img->depth);
			 exit(1);
			 break;
		}
	   }
	else if (flags & FITS_SKIP){
	   int err;

	   // we are not going to use the data so just seek past it
	   if (flags & ARCHIVE_COMPRESSED) err = fseek(in,zbytes,SEEK_CUR);
	   else err = fseek(in,bufsize,SEEK_CUR);

	   if (err) {
              Print("Error seeking on %s\n",img->src_fname);
              exit(1);
              }
	   }
	else {
	   Print("load_fits_data: Unknown flags value %d\n",flags);
	   exit(1);
	   }

	// Read the padding data
	if (! (flags & ARCHIVE_COMPRESSED) && (bufsize % 2880) != 0) {
	   nrecords = bufsize / 2880 + 1;
	   total = nrecords * 2880;
	   padding = total - bufsize;
	   if (padding > 0) {
		char buffer[3000];
		int err = Fread(buffer,1,padding,in);
	        if (err != padding) {
		   Print("load_fits_data: warning: error reading fits padding (%d bytes) for %s\n",padding,img->src_fname);
		   }
		}
	   }
	   
	return 1;
	}

static int
load_image_data(FILE *in,struct Image *img)
	{
	switch(img->type) {
	   case IMG_PPM:
		return load_ppm_data(in,img);
		break;
	   case IMG_BMP:
		return load_bmp_data(in,img);
		break;
	   case IMG_FIT:
		return load_fits_data(in,img,FITS_LOAD);
		break;
	   }

	Print("load_data: Unknown filetype %d for %s\n",img->type,img->src_fname);
	return 0;
	}

static int
parse_file_type(char *fname)
	{
	char *ptr = fname + strlen(fname)-1;

	while(ptr != fname && *ptr != '.') --ptr;
        if (! strcasecmp(ptr,".bmp")) return IMG_BMP;

        if (! strcasecmp(ptr,".fit")) return IMG_FIT;
        if (! strcasecmp(ptr,".fts")) return IMG_FIT;
        if (! strcasecmp(ptr,".fits")) return IMG_FIT;

        if (! strcasecmp(ptr,".ppm")) return IMG_PPM;
        if (! strcasecmp(ptr,".pgm")) return IMG_PPM; // Not yet supported

	Print("Unsupported filetype: '%s'\n",fname);
	return 0;
	}

struct Image *
LoadImage(char *src_fname,char *dst_fname)
	{
	static FILE *in = NULL;
	static using_stdin=0;
	struct Image *img = CreateImage();
	int need_close = 1;
	char fname[32];

	if (dst_fname == NULL) dst_fname = "";

	if (!strcmp(src_fname,"-")) {
	   need_close=0; 
	   if (! using_stdin) {
		using_stdin=1;
#ifdef MSWIN32
	 	// silly win32 has text mode stdio by default
		_setmode(fileno(stdin),O_BINARY);
#endif
	   	in = stdin;
		}
	   }
	else in = fopen(src_fname,"rb");

	if (in==NULL) {
	   Print("Cannot open file '%s' for reading\n",src_fname);
	   DestroyImage(img);
	   return NULL;
	   }

	if (!strcmp(src_fname,"-")) {
	   // read null-terminated filename from standard input to replace 
	   // the parameter "-"
	   int i=0;
	   while(i<31) {
		int ch = ReadByte(stdin);
		fname[i++] = ch;
		if (ch==0 || ch==EOF) break;
		}
		
	   if (i==32 || fname[i-1] != 0) {
		printf("LoadImage: Could not read filename from start of stream, abort\n");
		exit(1);
		}
	   src_fname = fname;
	   }

	if (InputFileType == -1)
	     img->type = parse_file_type(src_fname);
	else img->type = InputFileType;

	if (! img->type) {
	   Print("LoadImage: Image source '%s' format not recognised\n",src_fname);
	   DestroyImage(img);
	   return NULL;
	   }

	img->src_fname = strdup(src_fname);
	img->dst_fname = strdup(dst_fname);

	if (!load_image_header(in,img)) {
	   Print("cannot read fits file '%s'\n",src_fname);
	   if (need_close) {fclose(in); in=NULL;}
	   DestroyImage(img);
	   return NULL;
	   }

        //if (!Quiet) Print("(%dx%dx%d) ",img->width,img->height,img->depth);
	// Set defaults
	img->dst_width = img->width;
	img->dst_height = img->height;

	// Default cutout is whole image
	img->cutout.x1 = 0;
	img->cutout.y1 = 0;
	img->cutout.x2 = img->width-1;
	img->cutout.y2 = img->height-1;

	if (!load_image_data(in,img)) {
	   Print("read_image: failed\n");
	   if (need_close) { fclose(in); in=NULL;}
	   DestroyImage(img);
	   return NULL;
	   }

	// If we are asked for greyscale then convert to 16bpp monochrome
	if (img->depth == 24 && Grey) {
           int j,i = img->width * img->height;
	   int max = 0;
	   int do_stretch = InputHistoStretch;
	   unsigned short *udata;
	   char *ptr;

	   if (! ConvertToFITS(img, 16)) {
		Print("Convert BMP  24 bpp colour -> 16 bpp FITS failed\n");
		exit(1);
		}
	   Print("Convert 24bpp -> 16bpp FITS monochrome\n");

	   udata = (unsigned short *)img->data;

	   // Find the brightest pixel, ignore top and bottom rows
	   for(j=img->width+1; j<i-img->width-1; ++j) if (udata[j]>max) max=udata[j];

	   // Maybe histo-stretch the data
	   if (do_stretch ==-1) {
		// AUTO, so look at max to decide. Whatever we decide
		// to do applies to all the images
		if (max < 32767) do_stretch = 1;
		else do_stretch = 0;
	 	}

	   if (do_stretch == 1 && max < 60000) {
	  	double scale = 60000.0 / max;  // brightest pixel is 90% point
		int o = img->width * img->height - 1;

		if (! Quiet) Print("InputHistoStretch: max=%d, scale=%lf\n",max,scale);

		while(o>=0) {
		   unsigned int x = udata[o]; x *= scale;
		   if (x<0) x=0; if (x>65535) x=65535;
		   udata[o--]=x;
		   }
		}

	   // Change filename suffix to .FIT
	   ptr = img->dst_fname + strlen(img->dst_fname)-4;
	   strcpy(ptr,".fit");
	   }

	// If we are tracking then remember the last centre.
	// If this is the very first image then initialise the
	// subregion centre to the commandline values
	if (Sub_x && Sub_y) {
	   if (cur_Sub_x && cur_Sub_y) {
	      	img->sub_x = cur_Sub_x;
	      	img->sub_y = cur_Sub_y;
		}
	   else {
	      	img->sub_x = Sub_x;
	      	img->sub_y = Sub_y;
		}
	   }

	if (need_close) {fclose(in); in=NULL;}

	if (doMirrorImage) MirrorImage(img);

	return img;
	}

FILE *
OpenArchive(char *archive_name)
	{
#ifdef MSWIN32
        FILE *in = fopen(archive_name,"rb");
#else
        extern FILE * fopen64();
        FILE *in = fopen64(archive_name,"rb");
#endif

	return in;
	}


// return image or NULL, set err to 1 on success, 0 on end of file, -1 on error
struct Image *
ArchiveLoadNextImage(FILE *in, int *err, int flags)
	{
	char fname[64], *ptr;
	char archive_fname[1024];
	struct Image *img;
	int ch,nflags;

	if (flags & ARCHIVE_COMPRESSED) {
	   // read null-terminated filename for archive
	    ptr = archive_fname; *ptr=0;
	    while(1) {
	   	ch = ReadByte(in);
	   	if (ch==EOF) { *err=0; return NULL;}
	   	if (ch<0) { *err=-1; return(NULL);}

		if (strlen(archive_fname) >= 1023) {
		   Print("ArchiveLoadNextImage(compressed): filename too long\n");
		   exit(1);
		   }

 	   	*(ptr++) = ch; *ptr=0; 
	   	if (ch==0) break;
	   	}
	   }

	 // read null-terminated filename
	 ptr = fname; *ptr=0;
	 while(strlen(fname)<32) {
	   ch = ReadByte(in);
	   if (ch==EOF) { *err=0; return NULL;}
	   if (ch<0) { *err=-1; return(NULL);}
 	   *(ptr++) = ch; *ptr=0; 
	   if (ch==0) break;
	   }

	if (strlen(fname) >= 32) {
	   Print("Archive entry '%s' error: filename too long\n",fname);
	   *err=-1;
	   return(NULL);
	   }

	// load FITS header
	img = CreateImage();
	if (! load_fits_header(in,img)) {
	   Print("Failed to load FITS header for archive entry '%s'\n",fname);
	   *err = -1;
	   return NULL;
	   }

	img->type = parse_file_type(fname);
	if (! img->type) {
	   Print("Unknown image type for entry '%s'\n",fname);
	   *err = -1;
	   return NULL;
	   }

	img->src_fname = strdup(fname);
	img->dst_fname = strdup(fname);
	img->dst_width = img->width;
	img->dst_height = img->height;

	// Default cutout is whole image
	img->cutout.x1 = 0;
	img->cutout.y1 = 0;
	img->cutout.x2 = img->width-1;
	img->cutout.y2 = img->height-1;

	nflags=0;
	if (flags & ARCHIVE_LOAD) nflags |= FITS_LOAD;
	if (flags & ARCHIVE_SKIP) nflags |= FITS_SKIP;
	if (flags & ARCHIVE_COMPRESSED) nflags |= FITS_COMPRESSED;

	if (! load_fits_data(in,img,nflags)) {
	   Print("Failed to load FITS data for entry '%s'\n",fname);
	   *err = -1;
	   return NULL;
	   }

	if (doMirrorImage) MirrorImage(img);

	*err = 1;
	return img;
	}

int
ExtractArchive(char *archive_name, char *dstdir)
	{
#ifdef MSWIN32
	FILE *in = fopen(archive_name,"rb");
#else
	extern FILE * fopen64();
	FILE *in = fopen64(archive_name,"rb");
#endif

	int err;
	struct Image *img;

	if (in == NULL) {
	   Print("Cannot open archive '%s'\n",archive_name);
	   perror(NULL);
	   exit(1);
	   }

	if (! isDirectory(dstdir)) Mkdir(dstdir);

	while(1) {
	   img = ArchiveLoadNextImage(in,&err,ARCHIVE_LOAD);

	   if (err<0) {
		Print("Fatal error reading archive '%s'\n",archive_name);
		exit(1);
		}
	   else if (err==0) {
		Print("Finished reading archive\n");
		break;
		}
	   else if (err==1) {
		char *ptr = Malloc(strlen(dstdir) + strlen(img->dst_fname) + 2);
		sprintf(ptr,"%s/%s",dstdir,img->dst_fname);
		if (img->dst_fname) free(img->dst_fname);
		img->dst_fname = ptr;
	      	WriteImage(img);
	      	Print("%s:%s\n",archive_name,img->dst_fname);
		}
	   else {
		Print("Unknown return value %d\n",err);
		exit(1);
		}

	   DestroyImage(img);
	   }

	fclose(in);
	}

//###################################################################################
//
// Image Writing Support Routines
//
//###################################################################################
static int
write_bmp_header(FILE *out, struct Image *img)
	{
	bmp_header H;
	int i;
	int nw = img->dst_width;
	int nh = img->dst_height;
	struct bmp_info *bi = img->info;

        if (! fwrite("BM",2,1,out)) {
           Print("Short write on %s\n",img->dst_fname);
           return(0);
           }

        if (img->depth == 8) {
           H.filesz = 54 + (nw * nh) + 1024;
           H.offset = 54 + 1024;
           H.bitcount = 8;
           H.cmapentries = bi->cmap_entries;
           }
        else {
           H.filesz = 54 + (nw * nh * 3);
           H.offset = 54;
           H.bitcount = 24;
           H.cmapentries = 0;
           }

        H.r1 = H.r2 = 0;
        H.size = 40;
        H.width = nw;
        H.height = nh;
        H.planes = 1;
        H.compression = 0;
        H.sizeimage = 0;
        H.xppm = 0;
        H.yppm = 0;
        H.clrimportant = 0;

        if (! fwrite((char *)&H,sizeof(H),1,out)) {
           Print("cannot write to output file %s\n",img->dst_fname);
           return 0;
           }

        /* If we have a colourmap then write it out */
        if (bi->cmap_entries>0) {
	   if (img->depth==8) {
	      // Generate a new colourmap
	      for(i=0; i<bi->cmap_entries; ++i) {
		 bi->cmap[i*4]   = i;
		 bi->cmap[i*4+1] = i;
		 bi->cmap[i*4+2] = i;
		 bi->cmap[i*4+3] = 0;
	 	 }
	      }
	   else {
	      Print("write_bmp_header: a colourmap requires a depth of 8\n");
	      return 0;
	      }

           if (! fwrite(bi->cmap,bi->cmap_entries * 4,1,out)) {
              Print("Cannot write colourmap to output file %s\n",img->dst_fname);
              return 0;
              }
           }

	return 1;
	}
	
static int
write_ppm_header(FILE *out, struct Image *img)
	{
        fprintf(out,"P6\n%d %d\n255\n",img->width,img->height);
	return 0;
	}

static int
write_pgm_header(FILE *out, struct Image *img)
        {
	if (img->depth==8) fprintf(out,"P5\n%d %d\n255\n",img->width,img->height);
	else if (img->depth==16) fprintf(out,"P5\n%d %d\n65535\n",img->width,img->height);
        return(0);
        }


static int
write_fits_header(FILE *out, struct Image *img)
	{
	static double last_bscale = -1;
	static int last_depth = -1;
	static int last_bzero = -1;
	static int last_w=-1, last_h=-1;
	static double tm_sec = 0;
	static unsigned long long tm_tsc = 0;
	static char *buffer = NULL;
	struct fits_info *fi = (struct fits_info *)img->info;
	int i=0;
	u32 dmax;

	// If no change to params then there is no change to the header
	if (last_bscale == fi->bscale && last_bzero == fi->bzero && 
		last_depth == img->depth && last_w == img->dst_width && last_h == img->dst_height
		&& img->tm_sec == tm_sec && img->tm_tsc == tm_tsc && buffer)
	   goto writeout;

	if (buffer==NULL) buffer = (char *)Malloc(2881);
	if (! buffer) {
	   Print("make_fits_header: out of memory\n");
	   return 0;
	   }

	i=0;
	sprintf(buffer+i*80,"%-8.8s=%21.21s%50s","SIMPLE","T",""); ++i;
	sprintf(buffer+i*80,"%-8.8s=%21d%50s","BITPIX",img->depth,""); ++i;
	sprintf(buffer+i*80,"%-8.8s=%21s%50s","NAXIS","2",""); ++i;
	sprintf(buffer+i*80,"%-8.8s=%21d%50s","NAXIS1",img->dst_width,""); ++i;
	sprintf(buffer+i*80,"%-8.8s=%21d%50s","NAXIS2",img->dst_height,""); ++i;
	sprintf(buffer+i*80,"%-8.8s=%21d%50s","BSCALE",(int)fi->bscale,""); ++i;
	sprintf(buffer+i*80,"%-8.8s=%21u%50s","BZERO",(u32)fi->bzero,""); ++i;

#if 0
	sprintf(buffer+i*80,"%-8.8s=%21d%50s","DATAMIN",0,""); ++i;

	dmax=0;
	switch(img->depth) {
	   case 8: dmax = 255; break;
	   case 16: dmax = 65535; break;
	   case 32: dmax = UMAX32; break;
	   }

	sprintf(buffer+i*80,"%-8.8s=%21u%50s","DATAMAX",dmax,""); ++i;
#endif

	if (img->depth == 32 || img->depth == -32) {
	   sprintf(buffer+i*80,"%-8.8s=%21s%50s","BYTESWAP","T",""); ++i;
	   }

	if (img->tm_sec || img->tm_tsc) {
	   unsigned int hr = img->tm_sec / 3600;
	   unsigned int min = (img->tm_sec - hr*3600) / 60;
	   double secs = img->tm_sec - hr*3600 - min * 60;

#ifdef MSWIN32
	   sprintf(buffer+i*80,"%-8.8s= '__TIMESTAMP' %02u:%02u:%2.6lf/%-20I64u%21s","COMMENT",
		hr,min,secs,img->tm_tsc,""); ++i;
#else
	   sprintf(buffer+i*80,"%-8.8s= '__TIMESTAMP' %02u:%02u:%2.6lf/%-20llu%21s","COMMENT",
		hr,min,secs,img->tm_tsc,""); ++i;
#endif
	   }
	sprintf(buffer+i*80,"%-80.80s","END",""); ++i;
	while(i<36) {sprintf(buffer+i*80,"%80s",""); ++i;}

	writeout:
	if (fwrite(buffer,2880,1,out) != 1) {
	   Print("Short write on FITS header\n");
	   return 0;
	   }

	last_depth = img->depth;
	last_bzero = fi->bzero;
	last_bscale = fi->bscale;
	last_w = img->dst_width;
	last_h = img->dst_height;

	tm_sec = img->tm_sec;
	tm_tsc = img->tm_tsc;

	return 1;
	}

int
write_image_header(FILE *out, struct Image *img)
	{
	if (img==NULL) {
	   Print("oops: write_image_header called with NULL ptr\n");
	   return 0;
	   }

	switch(img->type) {
	   case IMG_PPM:
		return write_ppm_header(out,img);
		break;
	   case IMG_BMP:
		return write_bmp_header(out,img);
		break;
	   case IMG_FIT:
		return write_fits_header(out,img);
		break;
	   default:
		Print("write_image_header: unknown image type %d\n",img->type);
		return 0;
		break;
	   }
	return 0;
	}

static int
write_ppm_data(FILE *out, struct Image *img)
	{
	Print("Not Implemented\n");
	return 0;
	}

static int
write_bmp_data(FILE *out, struct Image *img)
	{
	int nw = img->dst_width;
	int nh = img->dst_height;
	int width = img->width;
	int height = img->height;
	int depth = img->depth;
	struct bmp_info *bi = img->info;
	int x,y,x1,y1,x2,y2;
	int bpp = img->depth/8;
	unsigned char *sptr,pad[4];
	int padding;

	/*
	 * If the WHITE point has changed then scale all the data
	 */
	if (White != 255 && White > 0) {
	   int o, npix = width * height;
	   double scale = (double)White / 255.0;
	   unsigned char *data = img->data;

	   if (depth==8) for(o=0; o<npix; ++o) {
		double v = (double)data[o]; 
		v /= scale;
		if (v<0) v=0; if (v>255) v=255;
		data[o] = (unsigned char)v;
		}
	   else if (depth == 24) for(o=0; o<npix*3; o+=3) {
		double R = (double)data[o];
		double G = (double)data[o+1];
		double B = (double)data[o+2];

		R /= scale; G /= scale; B /= scale;

		if (R<0) R=0; if (R>255) R=255;
		if (G<0) G=0; if (G>255) G=255;
		if (B<0) B=0; if (B>255) B=255;

		data[o]   = (unsigned char)R;
		data[o+1] = (unsigned char)G;
		data[o+2] = (unsigned char)B;
	      	}
	   }
	
        /* Now write the buffer according to our reverse colourmap */
        /* Write a window with dimensions (newWidth, newHeight) */
        y1 = (height - nh)/2; y2 = y1 + nh-1;
        x1 = (width - nw)/2; x2 = x1 + nw-1;

	padding = 4 - (nw*bpp)&3; if (padding==4) padding=0;
	pad[0]=pad[1]=pad[2]=pad[3]=0;

        x1 *= bpp;
        for(sptr = img->data + y2*width*bpp,y = y2; y>=y1; y--, sptr -= width*bpp) {
	   if (! fwrite(sptr+x1,nw*bpp,1,out)) return 0;
	   if (padding && ! fwrite(pad,padding,1,out)) return 0;
	   }

	return 1;
	}

static unsigned int
max_histo_8(struct Image *img)
	{
	int w = img->width;
	int h = img->height;
	int npix = w * h;
	int count,o,total;
	unsigned int pixval[NBUCKETS_8];
	unsigned char *data = (unsigned char *)img->data;
	int threshhold = ThreshHold;

	// init all buckets to empty
	for(o=0; o<NBUCKETS_8; ++o) pixval[o]=0;

	// generate histogram, counting only pixels > ThreshHold
	for(o=count=0; o<npix; ++o) {
	   unsigned int v = data[o];
	   if (v >= threshhold) {
	      	pixval[v]++;
		count++;
		}
	   }

	// Count back from brightest -> dimmest, find the bucket where we have 1% of the
	// total pixels on the right hand side. This rejects noise
	total=0; o=NBUCKETS_8-1;
	while(o) {
	   total += pixval[o];
	   if (total >= count/100)
		return o;
		
	   --o;
	   }

	return 0;
	}

static unsigned int
max_histo_16(struct Image *img)
	{
	int w = img->width;
	int h = img->height;
	int npix = w * h;
	int count,o,total;
	unsigned int pixval[NBUCKETS_16];   // 1% buckets
	unsigned short *udata = (unsigned short *)img->data;
	int threshhold = ThreshHold << 8;
	int div = (int)(65535 / (NBUCKETS_16-1));

	// init all buckets to empty
	for(o=0; o<NBUCKETS_16; ++o) pixval[o]=0;

	// generate histogram, counting only pixels > ThreshHold
	for(o=count=0; o<npix; ++o) {
	   unsigned short v = udata[o];
	   if (v >= threshhold) {
	      	pixval[v/div]++;
		count++;
		}
	   }

	// Count back from brightest -> dimmest, find the bucket where we have 1% of the
	// total pixels on the right hand side. This rejects noise
	total=0; o=NBUCKETS_16-1;
	while(o) {
	   total += pixval[o];
	   if (total >= count/100)
		return o * div;
		
	   --o;
	   }

	return 0;
	}
	
static u32
max_histo_32(struct Image *img)
	{
	int w = img->width;
	int h = img->height;
	int npix = w * h;
	int count,o,total;
	int pixval[NBUCKETS_32];   // 1% buckets
	u32 *udata = (u32 *)img->data;
	int threshhold = ThreshHold << 24;
	int div = (int)(UMAX32 / (NBUCKETS_32-1));

	// init all buckets to empty
	for(o=0; o<NBUCKETS_32; ++o) pixval[o]=0;

	// generate histogram, counting only pixels > ThreshHold
	for(o=count=0; o<npix; ++o) {
	   u32 v = udata[o];
	   if (v >= threshhold) {
	      	pixval[v/div]++;
		count++;
		}
	   }

	// Count back from brightest -> dimmest, find the bucket where we have 1% of the
	// total pixels on the right hand side. This rejects noise
	total=0; o=NBUCKETS_32-1;
	while(o) {
	   total += pixval[o];
	   if (total >= count/100)
		return o * div;
		
	   --o;
	   }

	return 0;
	}

static int
_write_fits_8(FILE *out,struct Image *img)
	{
	int x1,y1,x2,y2,i,y;
	int ndatabytes=0;
	unsigned char *outbuffer;
        unsigned char *ptr = img->data;
        unsigned char *dptr = outbuffer;
	struct fits_info *fi = img->info;
	int max,bzero = fi->bzero;
	int rowbytes;

        y1 = (img->height - img->dst_height)/2; y2 = y1 + img->dst_height-1;
        x1 = (img->width - img->dst_width)/2; x2 = x1 + img->dst_width-1;

	// make a copy of the data for mangling
	ndatabytes = img->width * img->height * img->depth/8;
	outbuffer = Malloc(ndatabytes);
        dptr = outbuffer;

	   if (HistoStretch) {
	   	unsigned int maxh = max_histo_8(img);
		unsigned int h = HistoStretch;
		double scale = (double)h / (double)maxh;
		int npix = img->width * img->height;
		int o;

		for(o=0; o<npix; ++o) {
		   int v = ptr[o];
		   v *= scale; if (v<0) v=0; if (v>255) v=255;
		   ptr[o] = v;
		   }
	 	}

           i=img->width * img->height;
           while(i--) {
                //if (fi->bscale != 1.0) *(dptr++) /= fi->bscale;
                *(dptr++) = *(ptr++) - bzero;
                }

           if (HistoProtect) {
              double scale = (double)White/255.0;
              unsigned int val = 255.0 * scale;
	      unsigned char *dp = outbuffer + y1*img->width + x1;

	      if (val<0) val=0; if (val>255) val=255;
              *dp = -bzero; *(dp+1) = val-bzero;
              }

	   ndatabytes=0;
	   rowbytes = x2-x1+1;
           dptr = (unsigned char *)outbuffer + y2 * img->width + x1;
           for(y=y2; y>=y1; dptr -= img->width, --y) {
                if (fwrite(dptr,rowbytes,1,out) != 1) {
                   Print("Short write on output to %s\n",img->dst_fname);
		   free(outbuffer);
                   return 0;
                   }
		ndatabytes += rowbytes;
		}

	free(outbuffer);
	return ndatabytes;
        }
	
static int
_write_fits_16(FILE *out, struct Image *img)
	{
	int x1,y1,x2,y2,i,y;
	int ndatabytes=0;
	unsigned char *outbuffer;
        unsigned short *iptr = (unsigned short *) img->data;
        unsigned short *diptr;
	struct fits_info *fi = img->info;
	int max,bzero = fi->bzero;
	int rowbytes;

        y1 = (img->height - img->dst_height)/2; y2 = y1 + img->dst_height-1;
        x1 = (img->width - img->dst_width)/2; x2 = x1 + img->dst_width-1;

	// make a copy of the data for mangling
	ndatabytes = img->width * img->height * abs(img->depth)/8;
	outbuffer = Malloc(ndatabytes);

        diptr = (unsigned short *) outbuffer;

	if (HistoStretch) {
	   unsigned int maxh = max_histo_16(img);
	   unsigned int h = HistoStretch << 8;
	   double scale = (double)h / (double)maxh;
	   int npix = img->width * img->height;
	   int o;

	   for(o=0; o<npix; ++o) {
		int v = iptr[o];
		v *= scale; if (v<0) v=0; if (v>65535) v=65535;
		iptr[o] = v;
		}
	   }

        // byteswap
        i=img->width * img->height;
        iptr = (unsigned short *) img->data;
        while(i--) {
	   unsigned short v = *(iptr++);
           //if (fi->bscale != 1.0) v /= fi->bscale;
	   *diptr = v - bzero; ByteSwapS(diptr++);
           }

        // Set pixel0 to black, pixel1 to white to prevent scaling
        if (HistoProtect) {
           double scale = (double)White/255.0;
           unsigned int val = 65535.0 * scale;
	   unsigned short v,*dp = (unsigned short *)outbuffer + y1*img->width + x1;

	   if (val<0) val=0; if (val>65535) val=65535;
           // byteswap
	   v = -bzero; *dp = v; ByteSwapS(dp);
	   v = val - bzero; *(dp+1) = v; ByteSwapS(dp+1);
           }

	rowbytes = (x2-x1+1) * 2;
	ndatabytes = 0;
        diptr = ((unsigned short *)outbuffer) + y2 * img->width + x1;
        for(y=y2; y>=y1; diptr -= img->width, --y) {
                if (fwrite((unsigned char *)diptr,rowbytes,1,out) != 1) {
                   Print("Short write on output to %s\n",img->dst_fname);
		   free(outbuffer);
                   return 0;
                   }
		ndatabytes += rowbytes;
		}

	free(outbuffer);
	return ndatabytes;
	}

static int
_write_fits_32(FILE *out, struct Image *img)
	{
	int x1,y1,x2,y2,i,y;
	int npixels,ndatabytes=0,rowbytes;
	unsigned char *outbuffer;
        u32 *uptr = (u32 *) img->data;
        u32 *duptr;
	struct fits_info *fi = img->info;
	u32 max,bzero = 1<<31;

        y1 = (img->height - img->dst_height)/2; y2 = y1 + img->dst_height-1;
        x1 = (img->width - img->dst_width)/2; x2 = x1 + img->dst_width-1;

	npixels = img->width * img->height;

	// make a copy of the data for mangling
	ndatabytes = npixels * abs(img->depth)/8;
	outbuffer = Malloc(ndatabytes);

        duptr = (u32 *) outbuffer;

	if (HistoStretch) {
	   u32 maxh = max_histo_32(img);
	   u32 h = HistoStretch << 24;
	   double scale = (double)h / (double)maxh;
	   int o;

	   for(o=0; o<npixels; ++o) {
		double v = uptr[o];
		v *= scale; if (v<0) v=0; if (v>UMAX32) v=UMAX32;
		uptr[o] = (u32)v;
		}
	   }

        // byteswap and convert to signed 32 bit
        uptr = (u32 *) img->data;
        i=npixels;
        while(i--) {
	   u32 v = *(uptr++);
           //if (fi->bscale != 1.0) v /= fi->bscale;
	   v -= bzero; v = ByteSwapLV(v); *duptr = v; duptr++;
           }

        // Set pixel0 to black, pixel1 to white to prevent scaling
        if (HistoProtect) {
           double scale = (double)White/255.0;
           u32 val = UMAX32 * scale;
	   u32 v,*dp = (u32 *)outbuffer + y1*img->width + x1;

           // byteswap
	   v = -bzero; *dp = v; ByteSwapL(dp);
	   v = val - bzero; *(dp+1) = v; ByteSwapL(dp+1);
           }

	ndatabytes=0;
	rowbytes = (x2-x1+1) * 4;
        duptr = ((u32 *)outbuffer) + y2 * img->width + x1;
        for(y=y2; y>=y1; duptr -= img->width, --y) {
                if (fwrite((unsigned char *)duptr,rowbytes,1,out) != 1) {
                   Print("Short write on output to %s\n",img->dst_fname);
		   free(outbuffer);
                   return 0;
                   }
		ndatabytes += rowbytes;
		}

	free(outbuffer);
	return rowbytes;
	}

static int
_write_fits_M32(FILE *out, struct Image *img)
	{
	int x1,y1,x2,y2,i,y;
	int npixels,ndatabytes=0,rowbytes;
	unsigned char *outbuffer;
        u32 *iuptr = (u32 *) img->data;
	u32 *ouptr;
	struct fits_info *fi = img->info;
	u32 max,bzero = 0;

        y1 = (img->height - img->dst_height)/2; y2 = y1 + img->dst_height-1;
        x1 = (img->width - img->dst_width)/2; x2 = x1 + img->dst_width-1;

	npixels = img->width * img->height;

	// make a copy of the data for mangling
	ndatabytes = npixels * abs(img->depth)/8;
	outbuffer = Malloc(ndatabytes);

        ouptr = (u32 *) outbuffer;

#if 0
	if (HistoStretch) {
	   u32 maxh = max_histo_32(img);
	   u32 h = HistoStretch << 24;
	   double scale = (double)h / (double)maxh;
	   int o;

	   for(o=0; o<npixels; ++o) {
		double v = uptr[o];
		v *= scale; if (v<0) v=0; if (v>UMAX32) v=UMAX32;
		uptr[o] = (u32)v;
		}
	   }
#endif

        // byteswap and convert to 32 bit floating point
        ouptr = (u32 *) outbuffer;
        for(i=0; i<npixels; ++i) {
	   u32 v = iuptr[i];
	   v = ByteSwapLV(v);
	   ouptr[i] = v;
           }

#if 0
        // Set pixel0 to black, pixel1 to white to prevent scaling
        if (HistoProtect) {
           double scale = (double)White/255.0;
           u32 val = UMAX32 * scale;
	   u32 v,*dp = (u32 *)outbuffer + y1*img->width + x1;

           // byteswap
	   v = -bzero; *dp = v; ByteSwapL(dp);
	   v = val - bzero; *(dp+1) = v; ByteSwapL(dp+1);
           }
#endif

	ndatabytes=0;
	rowbytes = (x2-x1+1) * 4;
        ouptr = ((u32 *)outbuffer) + y2 * img->width + x1;
        for(y=y2; y>=y1; ouptr -= img->width, --y) {
                if (fwrite((unsigned char *)ouptr,rowbytes,1,out) != 1) {
                   Print("Short write on output to %s\n",img->dst_fname);
		   free(outbuffer);
                   return 0;
                   }
		ndatabytes += rowbytes;
		}

	return ndatabytes;
	}

static int
write_fits_data(FILE *out, struct Image *img)
	{
	int ndatabytes = 0;

	switch(img->depth) {
	   case 8:  ndatabytes = _write_fits_8(out,img); break;
	   case 16: ndatabytes = _write_fits_16(out,img); break;
	   case 32: ndatabytes = _write_fits_32(out,img); break;
	   case -32: ndatabytes = _write_fits_M32(out,img); break;
	   default: 
		Print("Unknown depth %d in write_fits_data\n",img->depth);
		exit(1);
		break;
	   }

	// pad to multiple of 2880
	if (ndatabytes % 2880) {
	   char ch=0;
	   int padbytes = 2880 - (ndatabytes%2880);
	   ndatabytes += padbytes;
	   while(padbytes--) 
		if (fwrite(&ch,1,1,out) != 1)
		   Print("oops - write padding data failed\n");
	   }

        return(ndatabytes);
	}

int
write_image_data(FILE *out,struct Image *img)
	{
	switch(img->type) {
	   case IMG_PPM:
		return write_ppm_data(out,img);
		break;
	   case IMG_BMP:
		return write_bmp_data(out,img);
		break;
	   case IMG_FIT:
		return write_fits_data(out,img);
		break;
	   }

	Print("write_image_data: Unknown filetype %d for %s\n",img->type,img->src_fname);
	return 0;
	}

int
WriteImage(struct Image *img)
	{
	FILE *out;

	if (!img || ! img->dst_fname) {
	   Print("WriteImage: NULL parameter\n");
	   return 0;
	   }
	
	out = fopen(img->dst_fname,"wb");
	if (out==NULL) {
	   Print("Error writing file '%s'\n",img->dst_fname);
	   return 0;
	   }

	// Unless told otherwise we write the image out at the
	// same size as the source
	if (img->dst_width==0)  img->dst_width = img->width;
	if (img->dst_height==0) img->dst_height = img->height;

	if (! write_image_header(out,img)) {
	   Print("Error writing image header\n");
	   return 0;
	   }

	if (! write_image_data(out,img)) {
	   Print("Error writing image data\n");
	   return 0;
	   }

	fclose(out);
	return(1);
	}

//####################################################################################

// convert a buffer from one depth to another.
static unsigned char *
ConvertBufferDepth(unsigned char *in, int width, int height, int depth, int newdepth)
	{
	int i;
	unsigned char *buf;
	unsigned short *ubuf;
	int npix = width * height;

	// nothing to do if the depths are the same
	if (depth == newdepth) return in;

	buf = Malloc(width * height * abs(newdepth)/8);

	switch(newdepth) {
	   case 8:
		if (depth==16) {
		   ubuf = (unsigned short *)in;
		   i=width * height; while(i--) buf[i] = ubuf[i]>>8;
		   return buf;
		   }
		if (depth==24) {
		   i=width * height;
		   while(i--)
			buf[i] = in[i*3]   * 0.114 +   // Blue
				 in[i*3+1] * 0.587 +   // Green
				 in[i*3+2] * 0.299;    // red
		   return buf;
		   }
		break;
	   case 16:
		ubuf = (unsigned short *)buf;

		if (depth==8) {
		   i=width * height; while(i--) ubuf[i] = (unsigned)(in[i])<<8;
		   return buf;
		   }
		if (depth==24) {
		   for(i=0; i<npix; ++i)
			ubuf[i] = (unsigned)in[i*3]   * 29.184 +  // Blue
				  (unsigned)in[i*3+1] * 150.272 + // Green
				  (unsigned)in[i*3+2] * 76.544;   // Red
		   return buf;
		   }
		break;
	   case 24:
		if (depth==8) {
		   i=width * height; while(i--) {
			unsigned v = in[i];
			buf[i*3+0] = v;
			buf[i*3+1] = v;
			buf[i*3+2] = v;
			}
		   return buf;
		   }
		if (depth==16) {
		   ubuf = (unsigned short *)in;
		   i=width * height; while(i--) {
			unsigned v = ubuf[i]>>8;
			buf[i*3+0] = v;
			buf[i*3+1] = v;
			buf[i*3+2] = v;
			}
		   return buf;
		   }
		break;
	   case -32:
		// Convert between  32bpp unsigned int and 32bpp floating point
		// both monochrome
		if (depth == 32) {
		   u32 *sptr = (u32 *) in;
		   f32 *dptr = (f32 *) buf;
		   i = width * height;
		   while(i--) *(dptr++) = *(sptr++);
		   }
		else {
		   Print("ConvertBufferDepth: Cannot convert from %d to %d\n",depth,newdepth);
		   exit(0);
		   }

		return buf;
	   }

	Print("ConvertBufferDepth: Unknown depth converting %d -> %d\n",depth,newdepth);
	exit(1);
	}

static struct Image *
ConvertToBMP(struct Image *img, int newdepth)
	{
	int depth = img->depth;
	int width = img->width;
	int height = img->height;
	unsigned char *buf;
	struct bmp_info *bi;
	int i,o;


	buf = ConvertBufferDepth(img->data,width,height,depth,newdepth);

	if (buf == NULL) {
	   Print("ConvertToBMP: Buffer conversion failed\n");
	   return NULL;
	   }

	if (newdepth == 8) {
	   // generate a colourmap
	   bi = ZeroMalloc(sizeof(struct bmp_info));
	   bi->cmap = ZeroMalloc(1024);
	   for(i=o=0; i<256; ++i) {
	      bi->cmap[o++]=i;
	      bi->cmap[o++]=i;
	      bi->cmap[o++]=i;
	      bi->cmap[o++]=i;
	      }
	   bi->cmap_entries = 256;
	   if (img->info) free(img->info);
	   img->info = bi;
	   }

	if (buf != img->data) {
		free(img->data);
		img->data = buf;
		}

	img->type = IMG_BMP;
	img->depth = newdepth;
	return img;
	}

static struct Image *
ConvertToFITS(struct Image *img, int newdepth)
	{
	int depth = img->depth;
	int width = img->width;
	int height = img->height;
	unsigned char *buf;
	struct fits_info *fi;
	int i,o;

	buf = ConvertBufferDepth(img->data,width,height,depth,newdepth);
	if (buf == NULL) {
	   Print("ConvertToFITS: Buffer conversion failed\n");
	   return NULL;
	   }

	fi = ZeroMalloc(sizeof(struct fits_info));

	if (newdepth==16) fi->bzero=32768;
	else fi->bzero = 0;

	fi->bscale = 1;

	if (img->info) free(img->info);
	img->info = fi;

	if (buf != img->data) {
		free(img->data);
		img->data = buf;
		}

	img->type = IMG_FIT;
	img->depth = newdepth;
	return img;
	}

static struct Image *
ConvertToPPM(struct Image *img, int newdepth)
	{
	Print("ConvertToPPM: Not Implemented\n");
	exit(1);
	return img;
	}

struct Image *
ConvertToType(struct Image *img, int newtype, int newdepth)
	{
	char *ptr,fname[1024];
	int count=0;

	// maybe nothing to do
	if (newtype<0) {
	   Print("ConvertToType: invalid type %d\n",newtype);
	   return NULL;
	   }

	if (newtype != img->type || newdepth != img->depth) switch(newtype) {
	   case IMG_BMP: img = ConvertToBMP(img,newdepth); break;
	   case IMG_FIT: img = ConvertToFITS(img,newdepth); break;
	   case IMG_PPM: img = ConvertToPPM(img,newdepth); break;
	   default:
		Print("ConvertToType: new type %d not supported\n",newtype);
		return NULL;
		 break;
	   }

	// Make sure the file suffix is appropriate for the file type
	ptr = img->dst_fname + strlen(img->dst_fname)-1;
	count=0;
	while(ptr != img->dst_fname && *ptr != '.' && *ptr != '/' && *ptr != '\\') {--ptr; ++count;}
	if (*ptr != '.')
	   Print("Warning! Filename does not end with '.sfx'");
	else {
	   if (count<3) {
	      // bizzarre suffix like ".x" not enough space to expand
	      *ptr=0;
	      sprintf(fname,"%s.%s",img->dst_fname,".XXX");
	      free(img->dst_fname);
	      img->dst_fname = strdup(fname);
	      }

	    if (img->type == IMG_FIT) strcpy(ptr,".fit");
	    else if (img->type == IMG_BMP) strcpy(ptr,".bmp");
	    else if (img->type == IMG_PPM) strcpy(ptr,".ppm");
	    else {
		Print("Image type %d not understood\n",img->type);
		exit(1);
		}
	   }

	return img;
	}

/*
 * Given an archive name like "data0.fta" return "data1.fta" (old format)
 *
 * Given an archive name like "dataN-M-O.fta" return the next name using
 * this precedence:
 *	(N,M,O+1) 	Next part for this archive
 *	(N+1,M+1,0)	Changed wheel +1
 *	(N+1,M-1,0)	Changed wheel -1
 *	(N+1,*,0)	Changed wheel non-contiguous (maybe wrapped 6 -> 1)
 */

int
ArchiveNextFilename(char *this, char *next)
	{
	int count,m,n,o,i;
	char *ptr;
	
	// Locate the final component
	ptr = this + strlen(this) - 1;
	while(ptr!=this && *(ptr-1) != '/' && *(ptr-1) != '\\') --ptr;

	// count of how many prefix bytes to copy
	count = ptr-this;
	strncpy(next,this,count);

	if (sscanf(ptr,"data%d-%d-%d.fta",&n,&m,&o) == 3) {
	   // new format
	   ptr = next+count;
	   sprintf(ptr,"data%d-%d-%d.fta",n,m,o+1);
	   if (! access(next,0)) return m;

	   sprintf(ptr,"data%d-%d-%d.fta",n+1,m+1,0);
	   if (! access(next,0)) return m+1;

	   sprintf(ptr,"data%d-%d-%d.fta",n+1,m-1,0);
	   if (! access(next,0)) return m-1;

	   for(i=1; i<10; ++i) {
	      	sprintf(ptr,"data%d-%d-%d.fta",n+1,i,0);
	      	if (! access(next,0)) return i;
	  	}
	   }
	else if (sscanf(ptr,"data%d.fta",&n) == 1) {
	   // Old format, just increment
	   ptr = next+count;
	   sprintf(ptr,"data%d.fta",n+1);
	   if (! access(next,0)) return -1; // file present with unknown channel
	   return 0;
	   }

	next[0]=0;
	return 0;
	}