#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 struct Image * ConvertToPPM(struct Image *img, int newdepth);

static int load_ppm_data(FILE *in, struct Image *img);
static int load_ppm_data_8(FILE *in, struct Image *img);
static int load_ppm_data_16(FILE *in, struct Image *img);

static int ArchiveNextFilename_fta(char *this, char *next);
static int ArchiveNextFilename_ser(char *this, char *next);

static int WriteSerHeader(struct ser_info *s, FILE *out);
static int ReadSerHeader(struct ser_info *s, FILE *in);

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

// When reading from video streams, this int tells us the
// channel number for the current video frame
static int ch_num = 0;

int ShowHist_PrintToStdout = 1;

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

//******************************************************************************
//
// Do our own buffering  for file I/O
//
//******************************************************************************

#define CACHE_DISABLED 1

#define CACHE_BUFFER_SIZE (128 * 1024 * 1024)

off64_t CurrentFilePos[128]; 		// Current file position for each file descriptior
					// taking into consideration the read caching
static int cache_fd = -1;
static off64_t cache_offset = -1;
static unsigned char *cache_buffer = NULL;
static off64_t cache_buffer_start = -1;
static off64_t cache_buffer_end = -1;
static int cache_valid = 0;

static int cache_debug=0;

static int
_lowlevel_read(int fd, char *buf, int nbytes)
	{
	int retry_count=0;
	off64_t fd_offset = Tell(fd);
	
	if (cache_debug) {
	   Print("__lowlevel_read fd=%d ",fd);
	   Print("offset="); PrintUInt64(fd_offset);
	   Print(" nbytes=%d\n",nbytes);
	   }

	// If we try to read something too big for the cache then bypass it
	if (CACHE_DISABLED || nbytes > CACHE_BUFFER_SIZE) {
		int n;
		cache_valid=0;
		n = read(fd,buf,nbytes);
		CurrentFilePos[fd] = Tell(fd);
		if (cache_debug) Print("_lowlevel_read: read %d bytes\n",n);
		return n;
		}

	if (fd != cache_fd)
		cache_valid=0;

again:
	if (retry_count>2) {
		Print("_lowlevel_read: error reading offset=%d nbytes=%d, cache_start=%d, cache_end=%d\n",
				fd_offset, nbytes, cache_buffer_start, cache_buffer_end);
		exit(1);
		}

	// Has the cache been allocated?
	if (! cache_valid) {
		int n;
		off64_t tmp;

		if (cache_buffer == NULL)
			cache_buffer = Malloc(CACHE_BUFFER_SIZE);

		// Read some data without moving the file offset
		tmp = Tell(fd);
		n = read(fd,cache_buffer,CACHE_BUFFER_SIZE);
		lseek64(fd,tmp,SEEK_SET);

		if (n <= 0) return n;

		if (cache_debug) {
			Print("Read %d bytes from fd=%d into cache, offset=",n,fd);
			PrintUInt64(fd_offset);
			Print("\n");
			}

		cache_buffer_start = fd_offset;
		cache_buffer_end = fd_offset+n-1;
		cache_fd = fd;
		cache_valid=1;
		if (cache_debug) Print("cache contains %d - %d, valid=%d\n",cache_buffer_start, cache_buffer_end, cache_valid);
		}

	// Is this in the buffer?
	if (cache_buffer && fd_offset >= cache_buffer_start && (fd_offset+nbytes-1) <= cache_buffer_end) {
	   memcpy(buf,cache_buffer+(fd_offset-cache_buffer_start), nbytes);

	   if (cache_debug) Print("returning data from buffer\n");
	   // seek forward to pretend we read some data from the file
	   lseek64(fd,nbytes,SEEK_CUR);
	   CurrentFilePos[fd] = Tell(fd);
	   return nbytes;
	   }
	
	cache_valid=0;
	++retry_count;
	goto again;

	return 0; // notreached
	}

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

	if (fd==0) 
	     err = read(fd,&ch,1);
	else err = _lowlevel_read(fd,&ch,1);

	if (err==1) {
	   ++TotalFreadBytes;
	   CurrentFilePos[fd] = Tell(fd);
	   return ch;
	   }

	// some other error
	return EOF;
	}

static off64_t
Fseek(FILE *in, off64_t nbytes, int whence)
	{
	int fd = fileno(in);

	return lseek64(fd,nbytes,whence);
	}

off64_t Tell(int fd) { return lseek64(fd,0L,SEEK_CUR); }
off64_t FTell(FILE *in) { return Tell(fileno(in)); }


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

	while(bytes_remaining && failures < 5) {
	   int n;

	   if (fd==0) {
		n = read(fd,buffer,bytes_remaining);
		CurrentFilePos[fd] = Tell(fd);
		}
	   else n = _lowlevel_read(fd,buffer,bytes_remaining);

	   if (cache_debug) printf("Fread: read %d bytes\n",n);

	   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 %ul\n",
		bytes_remaining,size*count,off);

	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;
	img->histo_stretch = 1.0;
	img->auto_black = 0;

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

	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_ppm(struct Image *src, int depth)
	{
	Print("convert_to_ppm: 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,0,-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;
	   case IMG_PPM: return ConvertToPPM(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);
		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, char *ppmtype, 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;
	img->depth = 0;
	img->histo_stretch = 1;
	incomment=0; i=0;

	// If reading from stdin then we have already got the first line of the PPM/PGM
	if (ppmtype && ppmtype[0]=='P') {
	   strncpy(pi->ppm_type,ppmtype,2); pi->ppm_type[2]=0;
	   switch(ppmtype[1]) {
	   	case '3': img->depth = 24; break;
	   	case '6': img->depth = 24; break;
	   	case '5': img->depth = 8; break;
	   	}
	   }

	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') switch (pi->ppm_type[1]) {
				   case '3': img->depth = 24; break;
				   case '6': img->depth = 24; break;
				   case '5': img->depth = 8; break;  // default, may be overridden by maxval below
				   }

			   if (! img->depth) {
				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;
				if (pi->maxval>255 && pi->maxval<65536) img->depth=16;
				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, char *tmpbuf, 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;
	img->histo_stretch = 1;

        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);
           hline[80]=0;

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

           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 && d != -32) {
		   Print("FITS error: depth %d not understood\n",d);
		   exit(1);
		   }
		img->depth = d;
		}

           if (!strcmp(arg,"NAXIS")) {
		int n = atoi(val);

		 // If we've specified -grey on the commandline and we are given a
		 // colour FITS image then let it through, it will be converted into
		 // greyscale when we load the data
		if (n == 3 && Grey) {
			Print("Loading 3-colour FITS image as greyscale\n");
			}
		else if (n != 2) {
                   Print("AlignFIT: Unsupported header '%s = %s'\n",arg,val);
                   exit(1);
		   }

		fi->naxis = n;
                }

           if (!strcmp(arg,"NAXIS1")) { int n=atoi(val); fi->naxis1=n; img->width = n;}
           if (!strcmp(arg,"NAXIS2")) { int n=atoi(val); fi->naxis2=n; img->height = n;}
           if (!strcmp(arg,"NAXIS3")) { int n=atoi(val); fi->naxis3=n; }

           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 year=0,month=0,day=0,hr=0,min=0,n;
		double secs=0;
		unsigned char *ptr = val;

		while(*ptr && ! isdigit(*ptr)) ++ptr;

		// Look for optional YYYY-MM-DD
		if (sscanf(ptr,"%u-%u-%u %n",&year,&month,&day,&n) == 3) {
		   img->tm_sec = (double)DateTime_to_Unix_ms(year,month,day,0,0,0.0) / 1000.0;
		   //Print("Consumed %d chars\n",n);
		   ptr += n;
		   }

		// Look for HH:MM:SS.SSS
		if (sscanf(ptr,"%u:%u:%lf%n",&hr,&min,&secs,&n) == 3) {
		   img->tm_sec += (hr * 3600 + min * 60 + secs);
		   //Print("Consumed %d chars\n",n);
		   ptr += n;
		   }

		while(*ptr && *ptr==' ') ++ptr;

#ifdef MSWIN32
		if (sscanf(ptr,"/%I64u",&tsc) == 1)
#else
		if (sscanf(ptr,"/%llu",&tsc) == 1)
#endif
		   img->tm_tsc = tsc;
		}

	      // Load the (x,y) position on the ccd if given for the top left corner of the
	      // image
	      if (!strcasecmp(arg,"COMMENT") && !strncasecmp(val,"'__FORMAT7ROI'",14)) {
		unsigned char *ptr = val + 13;
		int xpos,ypos,width,height;

		// seek past the keyword to the first value
		while(*ptr && ! isdigit(*ptr)) ++ptr;

		if (sscanf(ptr,"%d %d %d %d",&xpos,&ypos,&width,&height) == 4) {
		   img->xpos = xpos;
		   img->ypos = ypos;
		   }
		else
		   Print("FORMAT7ROI: string not parsed\n");
	      	}
	   }

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

	// read all remaining header lines until we fine the END
	while (! found_end) {
	   //Print("Warning: END not present in first 36 records, continuing...\n");
	   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);
	   }

	// make sure we've read a multiple of 36 lines of 80 chars (2880 bytes)
	while(i%36) {
	   if (Fread(hline,80,1,in) != 1) {
		Print("Error reading header lines past END\n");
		exit(1);
		}
	   ++i;
	   }

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

        return 1;
        }

static int
load_bmp_header(FILE *in, char *tmpbuf, 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;
	img->histo_stretch = 1;
	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);
	         if (! BmpIgnoreOffset) 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, char *buf, 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,buf,img);
		break;
	   case IMG_BMP:
		return load_bmp_header(in,buf,img);
		break;
	   case IMG_FIT:
		return load_fits_header(in,buf,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)
	{
	if (img->depth == 8)  return load_ppm_data_8(in,img);
	if (img->depth == 16) return load_ppm_data_16(in,img);

	Print("load_ppm_data: depth %d not supported\n",img->depth);
	return 0;
	}

static int
load_ppm_data_8(FILE *in, struct Image *img)
	{
	struct ppm_info *pi = img->info;
        unsigned char *ptr = NULL;
	int i,j,do_stretch = InputHistoStretch;
	int max,maxp[MAXP];  // The brightest non-white pixels
	int npix = img->width * img->height;
	int overflow=0;
	int underflow=0;

	switch(img->depth) {
	   case 8:  // P5 PGM

		// Allocate space for 8 bit data
		img->data = (unsigned char *)Malloc(npix);
		if (Fread(img->data,npix,1,in) != 1) {
		   Print("load_ppm_data_8: failed to read %d pixels\n",npix);
		   return 0;
		   }
		break;
	   default:
		Print("load_ppm_data_8: Format not supported\n");
		return 0;
	   }

	if (ShowHistogram) ShowHistogram_8(img->data, npix, ShowHist_PrintToStdout);

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

        for(i=0; i<npix; ++i,++ptr) {
	   // Apply bias value if given
	   if (Bias) {
	 	j=*ptr; j+=Bias;
		// clamp the result
		if (j>255) { ++overflow; j=255; }
		if (j<0) { ++underflow; j=0; }
		*ptr = j;
		}

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

           }

        if (overflow) { Print("load_ppm_data_8: Warning: bias %d causes 8 bit overflow in %d pixels\n",Bias,overflow); }
        if (underflow) { Print("load_ppm_data_8: Warning: bias %d causes 8 bit underflow in %d pixels\n",Bias,underflow); }

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

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

	// DarkFrame subtract
	if (DoDarkFrameSubtract) SubtractDarkFrame(img);

	// Flat Frame (gaincomp) adjust
	if (ApplyGainComp) ApplyGainCompensation(img);

	// Automatic black level adjustment, average top 3 rows, and bottom 3 rows, exclude left and right pixels
	// and choose whichever is brighter
	if (AutoBlack) {
	   double total_t=0,total_b=0;
	   int x,y,v,count_t=0,count_b=0;
	   int abl, abl_t, abl_b;
	   int width = img->width;
	   int height = img->height;

	   for(y=0; y<3; ++y) {
		x=2;
	   	ptr=img->data + y*width + x;
		while (x<width-2) {
			total_t += (double) *(ptr++);
			++x; ++count_t;
			}
		}

	   for(y=height-4; y<height-1; ++y) {
		x=2;
	   	ptr=img->data + y*width + x;
		while (x<width-2) {
			total_b += (double) *(ptr++);
			++x; ++count_b;
			}
		}

	   abl_t = total_t / count_t;
	   abl_b = total_b / count_b;
	   if (abl_t > abl_b) { abl = abl_t; }
	   else abl = abl_b;

	   //Print("Automatic Black Level adjustment: Black = %d\n",abl);
	   ptr = (unsigned char *)img->data;
	   x=npix; while(x--) {
		// Subtract with clamping to zero
		v = ptr[x]; v -= abl; if (v<0) v=0;
		ptr[x] = v;
	   	}

	   img->auto_black = abl;
	   if (max >= abl) max -= abl;
	   else max=0;
	   }

	if (do_stretch == 1) {
	   double scale = (255.0 * MAX_STRETCH) / max;

	   if (scale > 0.25) {
		int o = img->width * img->height - 1;
           	ptr = (unsigned char *)img->data;

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

	   if (!Quiet) Print("InputHistoStretch: ppm_8 max_stretch=%1.2f, max=%d, scale = %2.2f\n",MAX_STRETCH,max,img->histo_stretch);
	   }

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

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

	return 1;
	}

static int
load_ppm_data_16(FILE *in, struct Image *img)
	{
	struct ppm_info *pi = img->info;
        unsigned short *uptr = NULL;
	int i,j,do_stretch = InputHistoStretch;
	int max,maxp[MAXP];  // The brightest non-white pixels
	int npix = img->width * img->height;
	int overflow=0;
	int underflow=0;

	switch(img->depth) {
	   case 16:  // P5 16 bit PGM

		// Allocate space for 16 bit data
		img->data = (unsigned char *)Malloc(npix * 2);
		if (Fread(img->data,npix*2,1,in) != 1) {
		   Print("load_ppm_data_16: failed to read %d pixels\n",npix);
		   return 0;
		   }
		break;

	   default:
		Print("load_ppm_data_16: Format not supported\n");
		return 0;
	   }

        // byteswap
        i=npix;
        uptr = (unsigned short *) img->data;
        while(i--) ByteSwapS(uptr++);

	if (ShowHistogram) ShowHistogram_16((unsigned short *)img->data, npix, ShowHist_PrintToStdout);

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

        for(i=0; i<npix; ++i,++uptr) {
	   // Apply bias value if given
	   if (Bias) {
	 	j=*uptr; j+=Bias;
		// clamp the result
		if (j>65535) { ++overflow; j=65535; }
		if (j<0) { ++underflow; j=0; }
		*uptr = j;
		}

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

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

           }

        if (overflow) { Print("load_ppm_data_16: Warning: bias %d causes 8 bit overflow in %d pixels\n",Bias,overflow); }
        if (underflow) { Print("load_ppm_data_16: Warning: bias %d causes 8 bit underflow in %d pixels\n",Bias,underflow); }

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

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

	// DarkFrame subtract
	if (DoDarkFrameSubtract) SubtractDarkFrame(img);

	// Flat Frame (gaincomp) adjust
	if (ApplyGainComp) ApplyGainCompensation(img);

	// Automatic black level adjustment, average top 3 rows, and bottom 3 rows, exclude left and right pixels
	// and choose whichever is brighter
	if (AutoBlack) {
	   double total_t=0,total_b=0;
	   int x,y,v,count_t=0,count_b=0;
	   int abl, abl_t, abl_b;
	   int width = img->width;
	   int height = img->height;

	   for(y=0; y<3; ++y) {
		x=2;
	   	uptr=(unsigned short *)img->data + y*width + x;
		while (x<width-2) {
			total_t += (double) *(uptr++);
			++x; ++count_t;
			}
		}

	   for(y=height-4; y<height-1; ++y) {
		x=2;
	   	uptr=(unsigned short *)img->data + y*width + x;
		while (x<width-2) {
			total_b += (double) *(uptr++);
			++x; ++count_b;
			}
		}

	   abl_t = total_t / count_t;
	   abl_b = total_b / count_b;
	   if (abl_t > abl_b) { abl = abl_t; }
	   else abl = abl_b;

	   //Print("Automatic Black Level adjustment: Black = %d\n",abl);
	   uptr = (unsigned short *)img->data;
	   x=npix; while(x--) {
		// Subtract with clamping to zero
		v = uptr[x]; v -= abl; if (v<0) v=0;
		uptr[x] = v;
	   	}

	   img->auto_black = abl;
	   if (max >= abl) max -= abl;
	   else max=0;
	   }

	if (do_stretch == 1) {
	   double scale = (65635.0 * MAX_STRETCH) / max;

	   if (scale > 0.25) {
		int o = img->width * img->height - 1;
           	uptr = (unsigned short *)img->data;

		img->histo_stretch = scale;
		while(o>=0) {
		   int x = uptr[o]; x *= scale;
		   if (x<0) x=0; if (x>65535) x=65535;
		   uptr[o--]=x;
		   }
		}
	   else img->histo_stretch = 1;

	   if (!Quiet) Print("InputHistoStretch: ppm16 max_stretch=%1.2f, max=%d, scale = %2.2f\n",MAX_STRETCH,max,img->histo_stretch);
	   }

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

	return 1;
	}

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);
	int npix = img->width * img->height;

        // 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;
		      }
		   }
		if (ShowHistogram) ShowHistogram_8(buffer, npix, ShowHist_PrintToStdout);
		}
	   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;
	   }
	
	}

// 32bpp double precision floating point.
// load and convert to U32
// NOT IMPLEMENTED
static void
_load_fixup_fits_F32(struct Image *img)
	{
	Print("Load FITS -32: Not implemented\n");
	exit(1);
        }

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;
	s32 bias = Bias * 256*256*256;  // promote 8bit value to 32bit, assume bias is small enough to fit
	double bf;
	u32 umin,umax,bzero = fi->bzero;
	s32 smin,smax;
	int overflow=0;
	int underflow=0;

        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;

           // Apply bias value if given
           if (bias) {
                bf=v; bf+=bias;
                // clamp the result
                if (bf>UMAX32) { ++overflow; bf=UMAX32;}
                if (bf<0) { ++underflow; bf=0; }
                v=(u32)bf;
                }

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

        if (overflow) { Print("Warning: bias %d causes 32 bit overflow in %d pixels\n",Bias,overflow); }
        if (underflow) { Print("Warning: bias %d causes 32 bit underflow in %d pixels\n",Bias,underflow); }

	// 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);
	   img->histo_stretch = 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;
        }

static void
_load_fixup_fits_M32(struct Image *img)
	{
	struct fits_info *fi = img->info;
        u32 *uptr = (u32 *)img->data;
	int i,npix,do_stretch = InputHistoStretch;
	u32 umin,umax,bzero = fi->bzero;

        npix=img->width * img->height;
	umax=0; umin=UMAX32; 

	img->depth = 32;

       	// Undo the endian-ness of the data for 32bpp FITS data
	// and convert to unsigned 32 bit
        for(i=0; i<npix; ++i) {
	   u32 v;
	   union {
		u32 u;
		float f;
	   	} uf;

           uf.u = ByteSwapLV(uptr[i]);
           v = (u32) uf.f;
           uptr[i]= v;
	   if (i>2) {
	      if (v>umax) umax=v;
	      if (v<umin) umin=v;
	      }
           }

	Print("max=%u (%2.2f %%) min=%u\n",umax,(double)umax/(double)UMAX32 * 100.0,umin);

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

static int Hist[65536];

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,*uptr;
	int npix,i,j,do_stretch = InputHistoStretch;
	int max,maxp[MAXP];  // The brightest non-white pixels
	int overbright = 0;  // count of overbright pixels
	int bias = Bias * 256;
	int overflow=0;
	int underflow=0;

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


           // Apply bias value if given
           if (bias) {
                j=*iptr; j+=bias;
                // clamp the result
                if (j>65535) { ++overflow;  j=65535; }
                if (j<0) { ++underflow; j=0; }
                *iptr=j;
                }

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

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

	   if (i>20 && *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 (ShowHistogram) ShowHistogram_16((unsigned short *)img->data, npix, ShowHist_PrintToStdout);

        if (overflow) { Print("Warning: bias %d causes 16 bit overflow in %d pixels\n",Bias,overflow); }
        if (underflow) { Print("Warning: bias %d causes 16 bit underflow in %d pixels\n",Bias,underflow); }

	// 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 *)img->data;
	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;

	if (EnableHistoWarning && overbright > 100 && ! Quiet) {
	   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<1) max=1;

	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 (doRowBoost) ApplyRowBoost(img);

	// Apply antighost
	if (AntiGhostEnable)
		ApplyAntiGhost(img);

	// DarkFrame subtract
	if (DoDarkFrameSubtract) SubtractDarkFrame(img);

	// Flat Frame (gaincomp) adjust
	if (ApplyGainComp) ApplyGainCompensation(img);

	// Automatic black level adjustment, average top 3 rows, and bottom 3 rows, exclude left and right 20 pixels
	// and choose whichever is brighter
	if (AutoBlack) {
	   double total_t=0,total_b=0;
	   int x,y,v,count_t=0,count_b=0;
	   int abl, abl_t, abl_b;
	   int width = img->width;
	   int height = img->height;

	   for(y=1; y<4; ++y) {
		x=2;
	   	uptr=(unsigned short *)img->data + y*width + x;
		while (x<width-2) {
			if (*uptr < 65000) { total_t += (double) *uptr; ++count_t; }
			++x; ++uptr;
			}
		}

	   for(y=height-5; y<height-2; ++y) {
		x=2;
	   	uptr=(unsigned short *)img->data + y*width + x;
		while (x<width-2) {
			if (*uptr < 65000) { total_b += (double) *uptr; ++count_b; }
			++x; ++uptr;
			}
		}

	   abl_t = total_t / count_t;
	   abl_b = total_b / count_b;
	   if (abl_t > abl_b) { abl = abl_t; }
	   else abl = abl_b;

	   if (! Quiet) Print("Automatic Black Level adjustment: Black = %d, max=%d\n",abl,max);
	   uptr = (unsigned short *)img->data;
	   x=npix; while(x--) {
		// Subtract with clamping to zero
		v = uptr[x]; v -= abl; if (v<0) v=0;
		uptr[x] = v;
	   	}

	   img->auto_black = abl;
	   if (max >= abl) max -= abl;
	   else max=0;

	   if (! Quiet) Print("Autoblack: fits16 max=%d\n",max);
	   }

	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 = (65535.0 * MAX_STRETCH) / max;

	   if (scale > 0.25) {
	      int o = img->width * img->height - 1;
              unsigned short *iptr = (unsigned short *)img->data;

	      img->histo_stretch = scale;
	      while(o>=0) {
		int x = iptr[o]; x *= scale;
		if (x<0) x=0; if (x>65535) x=65535;
		iptr[o--]=x;
		}
	      }
	   else img->histo_stretch = 1;

	   if (!Quiet) Print("InputHistoStretch: fits16 max_stretch=%1.2f, max=%d, scale = %2.2f\n",MAX_STRETCH,max,img->histo_stretch);
	   }
		
	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;
		}
	   }
	   
	// Finished
        }

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
	int npix = img->width * img->height;
	int overflow=0;
	int underflow=0;

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

	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;

	   // Apply bias value if given
	   if (Bias) {
	 	j=*ptr; j+=Bias;
		// clamp the result
		if (j>255) { ++overflow; j=255; }
		if (j<0) { ++underflow; j=0; }
		*ptr=j;
		}

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

	   if (i>20 && *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++;
           }

	if (ShowHistogram) ShowHistogram_8(img->data, npix, ShowHist_PrintToStdout);

        if (overflow) { Print("Warning: bias %d causes 8 bit overflow in %d pixels\n",Bias,overflow); }
        if (underflow) { Print("Warning: bias %d causes 8 bit underflow in %d pixels\n",Bias,underflow); }

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

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

	//printf("do_stretch=%d max=%d\n",do_stretch,max);
	//fflush(stdout);

	if (doRowBoost) ApplyRowBoost(img);

	// DarkFrame subtract
	if (DoDarkFrameSubtract) SubtractDarkFrame(img);

	// Flat Frame (gaincomp) adjust
	if (ApplyGainComp) ApplyGainCompensation(img);

	// Automatic black level adjustment, average top 3 rows, and bottom 3 rows, exclude left and right pixels
	// and choose whichever is brighter
	if (AutoBlack) {
	   double total_t=0,total_b=0;
	   int x,y,v,count_t=0,count_b=0;
	   int abl, abl_t, abl_b;
	   int width = img->width;
	   int height = img->height;

	   for(y=0; y<3; ++y) {
		x=20;
	   	ptr=img->data + y*width + x;
		while (x<width-20) {
			total_t += (double) *(ptr++);
			++x; ++count_t;
			}
		}

	   for(y=height-4; y<height-1; ++y) {
		x=20;
	   	ptr=img->data + y*width + x;
		while (x<width-20) {
			total_b += (double) *(ptr++);
			++x; ++count_b;
			}
		}

	   abl_t = total_t / count_t;
	   abl_b = total_b / count_b;
	   if (abl_t > abl_b) { abl = abl_t; }
	   else abl = abl_b;

	   //Print("Automatic Black Level adjustment: Black = %d\n",abl);
	   ptr = (unsigned char *)img->data;
	   x=npix; while(x--) {
		// Subtract with clamping to zero
		v = ptr[x]; v -= abl; if (v<0) v=0;
		ptr[x] = v;
	   	}

	   img->auto_black = abl;
	   if (max >= abl) max -= abl;
	   else max=0;
	   }

	if (do_stretch == 1) {
	   double scale = (255.0 * MAX_STRETCH) / max;

	   if (scale > 0.25) {
		int o = img->width * img->height - 1;
           	ptr = img->data;

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

	    if (!Quiet) Print("InputHistoStretch: fits8 max_stretch=%1.2f, scale = %2.2f\n",MAX_STRETCH,img->histo_stretch);
	    }

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


	}

// Column correction, average the top and bottom 20 rows, keep averages for each column and
// determine correction value per column
void
ApplyColumnCorrection(struct Image *img)
	{
	static double *col_avg = NULL;
	double avg=0;
	int nrows = ColumnCorrection;
	int x,y,cols_count=0;
	int avg_count=0;
	int width = img->width;
	int height = img->height;
	unsigned short *uptr;
	unsigned char *ptr;

	if (! ColumnCorrection) return;

	if (col_avg == NULL) col_avg = ZeroMalloc(width * sizeof(double));
	else for(x=0; x<width; ++x) col_avg[x]=0;

	avg_count=0;
	for(y=1; y<=nrows; ++y,++cols_count) {
	   x=0;
	   uptr=(unsigned short *)img->data + y*width + x;
	   ptr=(unsigned char *)img->data + y*width + x;
	   if (img->depth==8) while (x<width) {
		avg += *ptr; avg_count++;
		col_avg[x++] += (double) *(ptr++);
		}
	   else if (img->depth==16) while (x<width) {
		avg += *uptr; avg_count++;
		col_avg[x++] += (double) *(uptr++);
		}
	   }

	for(y=height-nrows-1; y<height-1; ++y,++cols_count) {
	   x=0;
	   uptr=(unsigned short *)img->data + y*width + x;
	   ptr=(unsigned char *)img->data + y*width + x;
	   if (img->depth==8) while (x<width) {
		avg += *ptr; avg_count++;
		col_avg[x++] += (double) *(ptr++);
		}
	   else if (img->depth==16) while (x<width) {
		avg += *uptr; avg_count++;
		col_avg[x++] += (double) *(uptr++);
		}
	   }

	// overall  average of all pixels counted
	avg /= avg_count;
	//Print("Column Correction: average = %lf, %d pixels in %d cols\n",avg,avg_count,cols_count);

	// Adjust all the column multipliers to compensate
	for(x=0; x<width; ++x) {
	   col_avg[x] /= cols_count;	// get column average
	   col_avg[x] = avg - col_avg[x];	// generate diff
	   }

	// Apply the compensation
	for(y=0; y<height; ++y) {
	   uptr = (unsigned short *)img->data + y*width;
	   ptr = (unsigned char *)img->data + y*width;

	   if (img->depth==8) for(x=0; x<width; ++x, ++ptr) {
		int v = *ptr + col_avg[x];
		if (v>255) v=255; if (v<0) v=0;
		*ptr = v;
		}
	   else if (img->depth==16) for(x=0; x<width; ++x, ++uptr) {
		int v = *uptr + col_avg[x];
		if (v>255*255) v=255*255; if (v<0) v=0;
		*uptr = v;
		}
	   }
	}

// Row correction, average the left and right ncols columns, keep averages for each row and
// determine correction value per row
void
ApplyRowCorrection(struct Image *img)
	{
	static double *row_avg = NULL;
	double avg=0;
	int ncols = RowCorrection;
	int x,y,rows_count=0;
	int avg_count=0;
	int width = img->width;
	int height = img->height;
	unsigned short *uptr;
	unsigned char *ptr;

	if (! RowCorrection) return;

	if (row_avg == NULL) row_avg = ZeroMalloc(height * sizeof(double));
	else for(y=0; y<height; ++y) row_avg[y]=0;

	avg_count=0;
	for(x=1; x<=ncols; ++x,++rows_count) {
	   y=0;
	   uptr=(unsigned short *)img->data + y*width + x;
	   ptr=(unsigned char *)img->data + y*width + x;
	   if (img->depth==8) while (y<height) {
		avg += *ptr; avg_count++;
		row_avg[y++] += (double) *ptr; ptr += width;
		}
	   else if (img->depth==16) while (y<height) {
		avg += *uptr; avg_count++;
		row_avg[y++] += (double) *uptr; uptr += width;
		}
	   }

	for(x=width-ncols-1; x<width-1; ++x,++rows_count) {
	   y=0;
	   uptr=(unsigned short *)img->data + y*width + x;
	   ptr=(unsigned char *)img->data + y*width + x;
	   if (img->depth==8) while (y<height) {
		avg += *ptr; avg_count++;
		row_avg[y++] += (double) *ptr; ptr += width;
		}
	   else if (img->depth==16) while (y<height) {
		avg += *uptr; avg_count++;
		row_avg[y++] += (double) *uptr; uptr += width;
		}
	   }

	// overall  average of all pixels counted
	avg /= avg_count;
	//Print("Row Correction: average = %lf, %d pixels in %d rows\n",avg,avg_count,rows_count);

	// Adjust all the row averages to be differences from the global average
	for(y=0; y<height; ++y) {
	   row_avg[y] /= rows_count;
	   row_avg[y] = avg - row_avg[y];	// generate diff
	   //Print("row_avg %d = %lf\n",y,row_avg[y]);
	   }

	// Apply the compensation
	for(y=0; y<height; ++y) {
	   uptr = (unsigned short *)img->data + y*width;
	   ptr = (unsigned char *)img->data + y*width;

	   if (img->depth==8) for(x=0; x<width; ++x, ++ptr) {
		int v = *ptr; v += row_avg[y];
		if (v>255) v=255; if (v<0) v=0;
		*ptr = v;
		}
	   else if (img->depth==16) for(x=0; x<width; ++x, ++uptr) {
		int v = *uptr; v += row_avg[y];
		if (v>65535) v=65535; if (v<0) v=0;
		*uptr = v;
		}
	   }
	}

static int
load_fits_data(FILE *in, struct Image *img, int flags)
	{
	struct fits_info *fi = img->info;
	int bufsize = img->width * img->height * abs(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;

	   // If we have 3 RGB channels and "Grey" commandline option then read
	   // the data and convert to greyscale
	   if (fi->naxis == 3 && Grey) bufsize *= 3;

	   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:  
			if (fi->naxis == 3 && Grey) {
			   unsigned char *dst = img->data;
			   unsigned char *src = img->data;
			   int npix = img->width * img->height;
			   int n;

			   n=npix;
			   while(n-- > 0) {
				   int R = *src;
				   int G = *(src+npix);
				   int B = *(src+npix*2);

				   *(dst++) = (double)R * 0.299 + (double)G * 0.587 + (double)B * 0.114;
				   src++;
			   	   }
			   }
			   

			_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;
		case -32: _load_fixup_fits_M32(img); break;
		default: Print("load_fits_data: Unknown depth %d\n",img->depth);
			 exit(1);
			 break;
		}
	   }
	else if (flags & FITS_SKIP){
	   off64_t err;

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

	   if (err <= 0) {
              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;

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

// When reading from stdin, ch_num must be set manually to the current 
// filter channel.
void SetChannel(int v) { ch_num = v; }

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[256],tmpbuf[256];
	static char stream_name[256] = {0};
	static int frame_number = 0;	// For AVI streams we keep a count of frames

	tmpbuf[0]=0;
	if (dst_fname == NULL) dst_fname = "";

	img->src_fname = strdup(src_fname);
	img->dst_fname = strdup(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,"-")) {
	   // FTA format: read null-terminated filename from standard input to replace 
	   // the parameter "-"
	   // PGM format: read "P5" on the first line, use our internal frame counter to generate a pseudo-filename

	   // Optionally recognise 'STREAM-NAME:xxx' to be the full path to the file providing the data stream
	   // Also recognise 'STREAM-INFO:name=value' as extra info about the stream

	   // remove the default names
	   if (img->src_fname) free(img->src_fname); img->src_fname = NULL;
	   if (img->dst_fname) free(img->dst_fname); img->dst_fname = NULL;

	   // If we've seen a sequence of 5 blank images in a row then assume this is a
	   // filter change
	   if (BlankImageCount == 5 && ch_num>0) {
		   ++ch_num;
		   Printf("Implied channel change => %d\n",ch_num);
		   BlankImageCount = -1;
		   }

	   int i=0;
	   while(i<255) {
		int ch = ReadByte(in);

		if (feof(in)) {
		   Printf("End of input stream\n");
		   return NULL;
		   }

		if (ch == '\r') continue;

		if (ch=='\n' && !strncmp(fname,"STREAM-NAME:",12)) {
			char *sfx = NULL;

			// Get the path component, use this as a filename prefix
			// for the generated names below

			// *Note* Only remember the first stream_name if there is a 
			// series, so the channels stay together in the output
			//
			if (! stream_name[0])
			   strcpy(stream_name,fname+12);

			// Try to parse the channel number
			if (strstr(fname+12,"_R")) 	sfx = "R";
			if (strstr(fname+12,"_G")) 	sfx = "G";
			if (strstr(fname+12,"_B")) 	sfx = "B";
			if (strstr(fname+12,"_IR")) 	sfx = "IR";
			if (strstr(fname+12,"_CH4")) 	sfx = "CH4";
			if (strstr(fname+12,"_UV")) 	sfx = "UV";
	   		if (strstr(fname+12,"IR610")) 	sfx = "IR610";
	   		if (strstr(fname+12,"IR750")) 	sfx = "IR750";
	   		if (strstr(fname+12,"IR800")) 	sfx = "IR800";
	   		if (strstr(fname+12,"IR807")) 	sfx = "IR807";
	   		if (strstr(fname+12,"IR850")) 	sfx = "IR850";
	   		if (strstr(fname+12,"IR900")) 	sfx = "IR900";
	   		if (strstr(fname+12,"IR950")) 	sfx = "IR950";
	   		if (strstr(fname+12,"IR1000")) 	sfx = "IR1000";

			if (sfx) {
			   ch = SfxToChannel(sfx);
			   if (ch>0) SetChannel(ch);
			   }

			fname[0]=0; i=0; continue;
			}
		if (ch=='\n' && !strncmp(fname,"STREAM-INFO:",12)) {
			char *name = fname+12; while(*name && *name==' ') ++name;
			char *sep = strstr(name,"=");
			char *value = NULL;

			if (sep) { *sep=0; value=sep+1;}

			if (name && sep) {
			   Printf("Meta: [%s] = [%s]\n",name,value);
			   }

			fname[0]=0; i=0; continue;
			}

		if (ch==0 || ch==EOF) {
		   if (fname[0]) {
		   	char buf[256];
		   	if (stream_name) 	sprintf(buf,"%s/%s",stream_name,fname);
		   	else			strcpy(buf,fname);

			img->src_fname = strdup(buf);
			img->dst_fname = strdup(fname);
		   	}
		   break;
		   }

		fname[i++] = ch; fname[i]=0;

		if (i==3 && fname[0]=='P' && fname[1]=='5' && fname[2]=='\n') {
		   char buf[256];
		   // Detected PGM frames on standard input

		   strcpy(tmpbuf,fname);  // preserve to pass to load_header below
		   sprintf(buf,"%05d-%d.pgm",++frame_number,ch_num);
	   	   img->dst_fname = strdup(buf);

		   if (stream_name) 	sprintf(buf,"%s/%05d-%d.pgm",stream_name,frame_number,ch_num);
		   else			sprintf(buf,"%05d-%d.pgm",frame_number,ch_num);
		   img->src_fname = strdup(buf);

		   i=strlen(fname)+1;

		   if (strcasestr(stream_name,".avi")) doInvertImage=0;

		   break;
		   }
		}
		
	   if (! img->src_fname || ! img->dst_fname) {
		printf("LoadImage: Could not set src and dst filename from data stream, abort.\n");
		return NULL;
		}
	   }

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

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

	if (!load_image_header(in,tmpbuf,img)) {
	   Print("cannot read fits file '%s'\n",img->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);

		img->histo_stretch = 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");
	   }

	// Store the centre of the image from last FindCentre() call
	img->sub_x = Sub_x;
	img->sub_y = Sub_y;

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

	if (doMirrorImage) MirrorImage(img);

	return img;
	}

static void
ByteSwap(unsigned char *data, int nbytes)
        {
        int i;

        while(nbytes) {
           unsigned char tmp = *data;
           *data = *(data+1);
           *(data+1) = tmp;
           nbytes -= 2;
           }
        }

int
SfxToChannel(char *sfx)
	{
	if (!strcasecmp(sfx,"R")) 	return CHANNEL_R;
	if (!strcasecmp(sfx,"G")) 	return CHANNEL_G;
	if (!strcasecmp(sfx,"B")) 	return CHANNEL_B;
	if (!strcasecmp(sfx,"IR")) 	return CHANNEL_IR;
	if (!strcasecmp(sfx,"CH4")) 	return CHANNEL_CH4;
	if (!strcasecmp(sfx,"UV")) 	return CHANNEL_UV;
	if (!strcasecmp(sfx,"IR610")) 	return CHANNEL_IR610;
	if (!strcasecmp(sfx,"IR750")) 	return CHANNEL_IR750;
	if (!strcasecmp(sfx,"IR800")) 	return CHANNEL_IR800;
	if (!strcasecmp(sfx,"IR807")) 	return CHANNEL_IR807;
	if (!strcasecmp(sfx,"IR850")) 	return CHANNEL_IR850;
	if (!strcasecmp(sfx,"IR900")) 	return CHANNEL_IR900;
	if (!strcasecmp(sfx,"IR950")) 	return CHANNEL_IR950;
	if (!strcasecmp(sfx,"IR1000")) 	return CHANNEL_IR1000;

	return 0;
	}

struct ser_info si_current;
unsigned si_framecount = 0;	// Frame number in current SER archive
unsigned si_frametotal = 0;	// Total frame count across all SER archives
int si_channel;
char si_fname[512];

FILE *
OpenArchive(char *archive_name)
	{
	char *ptr;

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

	if (! strcmp(archive_name+strlen(archive_name)-4,".ser")) {
	   // SER file, read the header
	   off64_t fsize = FileSize(archive_name);;
	   off64_t ts_offset;
	   unsigned long long st_tm;
	   char str[15],*ptr;
	   int c,ts_size,supported=0;

	   if (! ReadSerHeader(&si_current,in)) {
		Print("OpenArchive: Failed to read SER header from '%s'\n",archive_name);
		return NULL;
		}

	   Print("Loaded SER header: LuID=%d, ColorID=%d, %s, %d X %d X %d, %d frames\n",
		si_current.LuID,
		si_current.ColorID,
		si_current.LittleEndian==1 ? "Little Endian" : "Big Endian",
		si_current.ImageWidth,
		si_current.ImageHeight,
		si_current.PixelDepth,
		si_current.FrameCount
		);

	   Print("Observer   = [%-40.40s]\n",si_current.Observer);
	   Print("Instrument = [%-40.40s]\n",si_current.Instrument);
	   Print("Telescope  = [%-40.40s]\n",si_current.Telescope);

	   // Check the start time for the data
	   //Print("Start Time      : %s\n", SerTime_utc(si_current.dt));
	   //Print("Start Time (utc): %s\n", SerTime_utc(si_current.dt_utc));

	   // Automatically set the output file depth to 16 bit
	   // since Registax6 can't handle 8 bit fits.
	   // Don't override commandline choice
	   if (OutputFileDepth < 0) OutputFileDepth = 16;

	   if (si_current.PixelDepth != 8 && si_current.PixelDepth != 16) {
	 	Print("Pixel depth not supported\n");
		exit(1);
		}

	   if (si_current.ColorID != SER_COLORID_MONO) {
		Print("Only SER_COLORID_MONO data type supported\n");
		exit(1);
	  	}

	   // Are there timestamps provided?
	   ts_offset = (unsigned long long)SER_HEADER_SIZE + (unsigned long long)si_current.FrameCount 
							* (unsigned long long)si_current.ImageWidth 
							* (unsigned long long)si_current.ImageHeight 
							* (unsigned long long)si_current.PixelDepth/8;

	   ts_size = fsize - ts_offset;

	   if (ts_size/8 >= si_current.FrameCount) {
		int i,lower_n,upper_n;
		unsigned long long lower, upper;

		Print("Found timestamp array for %d frames (%d bytes)\n",si_current.FrameCount,ts_size);

		si_current.ts = (unsigned long long *)ZeroMalloc(ts_size);
		// Position and read the array of timestamps
		Fseek(in,(off64_t)ts_offset,SEEK_SET);
		if (Fread((unsigned char *)si_current.ts, ts_size, 1, in) != 1) {
		   Print("Error reading timestamp trailer, size should be %d bytes\n",ts_size);
		   exit(1);
		   }

		// Find the lowest / highest timestamps
		lower = upper = si_current.ts[0];
		lower_n = upper_n = 0;
		for(i=0; i< si_current.FrameCount; ++i) {
		   if (si_current.ts[i] < lower) { lower = si_current.ts[i]; lower_n = i; }
		   if (si_current.ts[i] > upper) { upper = si_current.ts[i];  upper_n = i; }
		   }
		Print("Lowest ts = %s (frame %d), Highest ts = %s (frame %d)\n",
			SerTime_utc(lower),lower_n,SerTime_utc(upper),upper_n);
		}
	   else {
		si_current.ts = NULL;
		Print("No timestamp array provided\n");
	 	}

	   // position the stream ready to read data
	   Fseek(in,(off64_t)SER_HEADER_SIZE,SEEK_SET);
	   si_framecount = 0;

	   // look for known channel identifier in filename
	   ptr = archive_name + strlen(archive_name)-4;
	   si_channel=1;
	   c=0;
	   while((ptr != archive_name) && *(ptr-1) != '_' && *(ptr-1) != '/' ) { --ptr; ++c;}
	   if (c>0) {
		char buf[32];
		int ch;

		strncpy(buf,ptr,c);
		buf[c]=0;
		ch = SfxToChannel(buf);
		if (ch>0) si_channel = ch;
		}

	   // Find the last path component and the filename
	   ptr = archive_name + strlen(archive_name);
	   while(ptr != archive_name && *(ptr-1) != '/' && *(ptr-1) != '\\') --ptr;
	   //if (ptr != archive_name) --ptr;
	   //while(ptr != archive_name && *(ptr-1) != '/' && *(ptr-1) != '\\') --ptr;
	   strcpy(si_fname,ptr);
	   }

	return in;
	}

// Load next image from SER format file
// return image or NULL, set err to 1 on success, 0 on end of file, -1 on error
struct Image *
ArchiveLoadNextImage_SER(FILE *in, int *err, int flags)
	{
	struct Image *img=NULL;
	off64_t fd_offset;
	int depth = si_current.PixelDepth / 8;
	int nbytes = si_current.ImageWidth * si_current.ImageHeight * depth;
	unsigned char *data = Malloc(nbytes);
	int need_byteswap = 0;
	int need_invert = 0;

	*err=-1;

	// Maybe we're at the end of the file
	if (si_framecount>= si_current.FrameCount) {
	   Print("All done\n"); fflush(stdout);
	   *err=0;
	   return NULL;
	   }

	// Read the data
	if (! Fread(data,nbytes,1,in)) {
	   Print("Error reading SER data\n"); fflush(stdout);
	   *err = 0;
	   return NULL;
	   }

#if 0
	if (depth==2 && ! si_current.LittleEndian) {
	   Print("ByteSwap\n");
           int i=si_current.ImageWidth * si_current.ImageHeight;
           unsigned short *uptr = (unsigned short *) img->data;
           while(i--) ByteSwapS(uptr++);
	   need_byteswap = 0;
	   }
#endif

	img = CreateFitsImage(
		"unknown",
		si_current.ImageWidth,si_current.ImageHeight, si_current.PixelDepth,
		data, 0, -1
		);

	// work out which frame we have just read
	// based on the current offset
	fd_offset = FTell(in);
	si_framecount = (unsigned int) ((fd_offset - SER_HEADER_SIZE) / nbytes);

	img->src_fname = Malloc(strlen(si_fname) + 32);
	sprintf(img->src_fname,"%05d-%d.fit",si_frametotal++,si_channel);
	img->dst_fname = strdup(img->src_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;

	if (doInvertImage) need_invert = 1;
	else need_invert=0;

	if (depth==1) _load_fixup_fits_8(img);
	else _load_fixup_fits_16(img, need_byteswap, need_invert);

	if (doMirrorImage) MirrorImage(img);

	// Was there a timestamp?
	if (si_current.ts != NULL) {
	   unsigned long long t = si_current.ts[si_framecount - 1];
	   img->tm_sec = (double)SerTimeStamp_Unix_ms(t) / 1000.0;
	   }

	*err = 1;
	return img;
	}

// 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;
	int fd = fileno(in);

	// SER file format
	if (flags & ARCHIVE_SER)
	   return ArchiveLoadNextImage_SER(in,err,flags);

	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) {
		   Print("ArchiveLoadNextImage(compressed): ReadByte returned %d\n",ch);
		   *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;
		Print("ArchiveLoadNextImage: ReadByte returned %d at offset %ld\n",ch,ftell(in));
		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);
	   }

	// Check for channel overrides to change the filename
	if (OverrideChannel) {
	   int frnum,ch;
	   char sfx[256];

	   if (sscanf(fname,"%5d-%d.%s",&frnum,&ch,sfx)==3) {
		if (frnum >= OverrideChannel_Min && frnum <= OverrideChannel_Max) {
			sprintf(fname,"%05d-%d.%s",frnum,OverrideChannel_Value,sfx);
			Print("OverrideChannel changes filename to '%s'\n",fname);
			}
		}
	   }

	// load FITS header
	img = CreateImage();
	if (! load_fits_header(in,(char *)NULL,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) && ! NoSave) 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 1;
	}

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


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) {
	   int year=0,month=0,day=0;
	   unsigned int hr=0,min=0;
	   double secs=0;

	   if (img->tm_sec < 100000) { // seconds in the current day only
	      	hr = img->tm_sec / 3600;
	      	min = (img->tm_sec - hr*3600) / 60;
	      	secs = img->tm_sec - hr*3600 - min * 60;
		}
	   else if (img->tm_sec > 1000000) { // most likely a UNIX time 
		UnixTime_ms_utc(img->tm_sec * 1000, &year,&month,&day,&hr,&min,&secs);
		}

#ifdef MSWIN32
	   sprintf(buffer+i*80,"%-8.8s= '__TIMESTAMP' %04d-%02d-%02d %02u:%02u:%2.6lf/%-20I64u%9s","COMMENT",
		year,month,day,hr,min,secs,img->tm_tsc,""); ++i;
#else
	   sprintf(buffer+i*80,"%-8.8s= '__TIMESTAMP' %04d-%02d-%02d %02u:%02u:%2.6lf/%-20llu%9s","COMMENT",
		year,month,day,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:
		if (img->depth == 24) return write_ppm_header(out,img);
		else if (img->depth == 8) return write_pgm_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)
	{
	int npix = img->width * img->height;

	switch(img->depth) {
	   case 8:
		if (fwrite(img->data,npix,1,out) != 1) {
		   Print("write_ppm_data: short write on output\n");
		   return 0;
		   }
		return 1;
		break;

	   default:
		Print("write_ppm_data: format not supported\n");
		return 0;
	 	break;
	   }

	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;

	   if (rowbytes&3) {
	 	Print("_write_fits_8: rowbytes %d not multiple of 4\n",rowbytes);
		return 0;
		}

           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,rowpixels;

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

	rowpixels = x2-x1+1;
	rowbytes = rowpixels * 2;

	if (rowpixels&3) {
	   Print("_write_fits_16: width %d not multiple of 4\n",rowpixels);
	   return 0;
	   }

	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;
	int rowbytes,rowpixels;
	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) {
	   u32 val = UMAX32;

	   if (White != 255) {
              	double scale = (double)White/255.0;
	      	double dval = (double)UMAX32 * scale; if (dval > UMAX32) dval=UMAX32;
              	val = (u32)dval;
	  	}

	   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;

	rowpixels = x2-x1+1;
	rowbytes = rowpixels * 4;

	if (rowpixels&3) {
	   Print("_write_fits_32: width %d not multiple of 4\n",rowpixels);
	   return 0;
	   }

        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;
	int rowbytes,rowpixels;
	unsigned char *outbuffer;
        u32 *iuptr = (u32 *) img->data;
        float *fptr = (float *) img->data;
	u32 *ouptr;
	struct fits_info *fi = img->info;
	float fmin,fmax;
	u32 max;

        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

#if 0
	fmin = FLT_MAX;
	fmax = FLT_MIN;
        for(i=0;i<npixels; ++i) {
	   if (fptr[i] > fmax) fmax = fptr[i];
	   if (fptr[i] < fmin) fmin = fptr[i];
	   }
	printf("Min/Max val = %f/%f\n",fmin,fmax);
#endif

        // Swap high and low words (big-endian)
        ouptr = (u32 *) outbuffer;
        for(i=0;i<npixels; ++i) {
	   u32 v = iuptr[i];
	   v = ByteSwapLV(v);
	   ouptr[i] = v;
           }

        // Set pixel0 to black, pixel1 to white to prevent scaling
        if (HistoProtect) {
	   union {
              float val;
	      u32 uval;
	      } u;

	   // The data is really 32 bit unsigned.
	   u.val = (double)65535.0 * 65536.0;

	   if (White < 255) {
              	double scale = (double)White/255.0;
              	u.val *= scale;
		}

	   u32 *dp = ((u32 *)outbuffer) + y2*img->width + x1;
	   u32 v;

	   v = ByteSwapLV(0); *dp = v;
	   v = ByteSwapLV(u.uval); *(dp+1) = v;
           }

	ndatabytes=0;

	rowpixels = x2-x1+1;
	rowbytes = rowpixels * 4;

	if (rowpixels&3) {
	   Print("_write_fits_M32: width %d not multiple of 4\n",rowpixels);
	   return 0;
	   }

        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;

	// Enforce a width that is a multiple of 4 to save
	// later embarrasment when creating movies with some codecs
	if (img->dst_width&3) {
	   Print("Image width (%d) must be multiple of 4\n",img->dst_width);
	   return 0;
	   }

	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 n,npix = width * height;
	int absdepth = newdepth;

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

	if (newdepth == -32) absdepth = 32;
	else if (newdepth<=8) absdepth = 8;
	else absdepth = abs(newdepth);

	n = width * height * absdepth/8;
	//Print("read %d pixels, write %d bytes\n",npix,n);

	buf = Malloc(n);

	switch(newdepth) {
	   case 4: case 5: case 6: case 7: case 8:
		if (depth==16) {
		   unsigned short mask = 0xff00 << (8-newdepth);
		   ubuf = (unsigned short *)in;
		   i=width * height; while(i--) {
			   unsigned int u = (ubuf[i] & mask) >> 8;
			   buf[i] = u;
			   }
		   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--) {
			//Print("Convert %d\n",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 if (newdepth==32) fi->bzero=(u32)(1<<31);
	else if (newdepth==-32) fi->bzero=0;
	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;

	if (newdepth>=0 && newdepth<=8) img->depth = 8;
	else img->depth = newdepth;
	return img;
	}

static struct Image *
ConvertToPPM(struct Image *img, int newdepth)
	{
	int depth = img->depth;
	int width = img->width;
	int height = img->height;
	struct Image *new = NULL;
	unsigned char *buf;
	struct ppm_info *pi;
	int i,o;

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

	pi = ZeroMalloc(sizeof(struct ppm_info));

	if (newdepth==8) {
		sprintf(pi->ppm_type,"P5");
		pi->maxval = 255;
		}
	else if (newdepth==16) {
		sprintf(pi->ppm_type,"P5");
		pi->maxval = 65535;
		}
	else if (newdepth==24) {
		sprintf(pi->ppm_type,"P6");
		pi->maxval = 255;
		}

	// Create a copy of the original image with the new values
	// remember to clone strings
	new = ZeroMalloc(sizeof(struct Image));
	memcpy(new,img,sizeof(struct Image));

	if (img->src_fname) new->src_fname = strdup(img->src_fname);
	if (img->dst_fname) new->dst_fname = strdup(img->dst_fname);
	new->data = buf;
	new->info = pi;
	new->depth = newdepth;
	new->type = IMG_PPM;

	return new;
	}

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

static int
ArchiveNextFilename_fta(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<NUM_CHANNELS; ++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;
	}

/*
 * Given an SER filename like "yyyy_hhmmss_F.ser", try to determine the
 * next archive name in the sequence. 
 * - Increment filter (R -> G, G -> B , B = end)
 * - Look for yyyy_hhmmss_F1.ser as next file in directory with no intervening files
 */


/* NOT IMPLEMENTED */
static int
ArchiveNextFilename_ser(char *this, char *next)
	{
	int y,m,d,h,min,s;
	char *ptr,sfx[10];
	
	next[0]=0;

	// Locate the final component
	ptr = this + strlen(this) - 1;
	while(ptr!=this && *(ptr-1) != '/' && *(ptr-1) != '\\') --ptr;

	if (sscanf(ptr,"%04d%02d%02d_%02d%02d%02d_%s.ser",&y,&m,&d,&h,&min,&s,sfx) == 7) {
	   char *next_sfx;
	   if (!strcasecmp(sfx,"R")) next_sfx = "G";
	   else if (!strcasecmp(sfx,"G")) next_sfx = "B";
	   else if (!strcasecmp(sfx,"B")) return 0;
	   }

	return 0;
	}

int
ArchiveNextFilename(char *this, char *next)
	{
	int len = strlen(this);
	if (! strcasecmp(this+len-4,".fta")) return ArchiveNextFilename_fta(this,next);
	if (! strcasecmp(this+len-4,".ser")) return ArchiveNextFilename_ser(this,next);

	return 0;
	}

static struct ser_info s;
static unsigned long long *s_timestamp = NULL;
static int s_timestamp_size = 0;
static FILE *ser_last = NULL;

//*********************************************
// We need a specific routine for this thanks to stupid
// alignments on the in-memory struct
//*********************************************

static int
WriteSerHeader(struct ser_info *s, FILE *out)
	{
	int nbytes;

	fflush(out);
	Fseek(out,(off64_t)0,SEEK_SET);

	// Write the magic string
	if (! fprintf(out,"LUCAM-RECORDER")) {
	   Print("WriteSerHeader: Cannot write magic string 'LUCAM-RECORDER'\n");
	   return 0;
	   }
	   
	// Write all the int members of the struct and the 
	// 3 x 40 char strings. They will all be aligned ok
	nbytes = sizeof(int) * 7 + 120;
	if (fwrite((char *)s, 1, nbytes, out) != nbytes) {
	   Print("WriteSerHeader: Short write\n");
	   return 0;
	   }

	// Write the 2 x unsigned long long members,
	nbytes = sizeof(unsigned long long) * 2;
	if (fwrite((char *)&s->dt, 1, nbytes, out) != nbytes) {
	   Print("WriteSerHeader: Short write\n");
	   return 0;
	   }

	// Offset must be 178 bytes
	if (ftell(out) != SER_HEADER_SIZE) {
	   Print("Error: Offset after header is %ld, should be %ld\n",(long)ftell(out),(long)SER_HEADER_SIZE);
	   exit(1);
	   }

	return 1;
	}

static int
ReadSerHeader(struct ser_info *s, FILE *in)
	{
	int nbytes;
	int supported;
	char str[15];

	Fseek(in,(off64_t)0,SEEK_SET);

	// read the 14 byte identifier
	if (Fread(str,1,14,in) != 14) {
	   Print("ReadSerHeadr: Error reading SER string id\n");
	   return 0;
	   }

	supported=0;
	if (! strncasecmp(str,"LUCAM-RECORDER",14)) supported=1;
	if (! strncasecmp(str,"FirecaptureV24",14)) supported=1;

	if (! supported) {
	   Print("ReadSerHeader: Invalid SER header - unsupported signature %-14.14s\n",str);
	   return 0;
	   }

	// Read the 7 x int params and 3 x char strings
	nbytes = sizeof(int) * 7 + 120;
	
	if (Fread((unsigned char *)s, 1, nbytes, in) != nbytes) {
	   Print("ReadSerHeader: short read 1\n");
	   return 0;
	   }

	// Read the 2 x unsigned long long params
	nbytes = sizeof(unsigned long long) * 2;

	if (Fread((unsigned char *)&s->dt, 1, nbytes, in) != nbytes) {
	   Print("ReadSerHeader: short read 2\n");
	   return 0;
	   }

	return 1;
	}

// Open and truncate the given file, write a SER header
FILE *
OpenSerForWriting(char *filename, int width, int height, int depth, char *Obs, char *Inst, char *Tel)
	{
	FILE *out;
	int i;

	if (depth != 8 && depth != 16) {
	   Print("OpenSerForWriting: Unsupported depth %d\n",depth);
	   return NULL;
	   }

	s_timestamp_size = 0;
	s_timestamp = NULL;

	s.ImageWidth = width;
	s.ImageHeight = height;
	s.PixelDepth = depth;	// 8 or 16
	s.ColorID = SER_COLORID_MONO;
	s.LittleEndian = 0;
	s.FrameCount = 0;
	s.dt = 0;
	s.dt_utc = 0;

	strncpy(s.Observer,Obs,39);
	strncpy(s.Instrument,Inst,39);
	strncpy(s.Telescope,Tel,39);

	out = fopen(filename,"wb+");
	if (! out) {
	   Print("OpenSerForWriting: Cannot open file '%s'\n",filename);
	   return NULL;
	   }

	if (! WriteSerHeader(&s,out)) {
	   Print("OpenSerForWriting: Failed\n");
	   exit(1);
	   }

	Print("Written header to '%s'\n",filename);
	ser_last = out;
	return out;
	}

int
WriteSerFrame(struct Image *img, FILE *out)
	{
	struct ser_info ss;
	char *ptr;
	int nbytes = s.ImageWidth * s.ImageHeight * (s.PixelDepth/8);
	int i,hb;

	if (out == NULL) out=ser_last;

	// Check that the specs match
	if (s.ImageWidth != img->width) {
	   Print("WriteSerFrame: Given width %d differs from SER header value %d\n",img->width,s.ImageWidth);
	   return 0;
	   }

	if (s.ImageHeight != img->height) {
	   Print("WriteSerFrame: Given height %d differs from SER header value %d\n",img->height,s.ImageHeight);
	   return 0;
	   }

	if (s.PixelDepth != img->depth) {
	   Print("WriteSerFrame: Given depth %d differs from SER header value %d\n",img->depth,s.PixelDepth);
	   return 0;
	   }

	if (s_timestamp==NULL || (s.FrameCount >= s_timestamp_size)) {
	   s_timestamp_size += 1000;
	   s_timestamp = Realloc(s_timestamp, s_timestamp_size * sizeof (unsigned long long));
	   }

	s_timestamp[s.FrameCount] = SerTimeStamp_FromUnix_ms((unsigned long long)(img->tm_sec * 1000));
	if (s_timestamp[s.FrameCount] < s.dt_utc || s.FrameCount==0) {
	   // s.dt_utc is the lowest timestamp in the data stream, notionally the timestamp for the
	   // first frame (but not if they've been resorted on quality)
	   s.dt_utc = s_timestamp[s.FrameCount];
	   s.dt = s_timestamp[s.FrameCount];
	   }

	// Write the data
	s.FrameCount++;
	fflush(out);
	Fseek(out,(off64_t)0,SEEK_END);
	if (! fwrite((void *)img->data, nbytes, 1, out)) {
	   Print("WriteSerFrame: Error writing frame data (frame %d)\n",s.FrameCount);
	   return 0;
	   }

	if (! WriteSerHeader(&s,out)) {
	   Print("WriteSerFrame: Failed on frame %d\n",s.FrameCount);
	   fclose(out); return 0;
	   }

	return 1;
	}

int
CloseSerFile(FILE *out)
	{
	if (out == NULL) out=ser_last;
	if (out == NULL) {
	   Print("CloseSerFile: No file open\n");
	   return 0;
	   }

	if (! WriteSerHeader(&s,out)) {
	   Print("CloseSerFile: Error writing header\n");
	   fclose(out);
	   return 0;
	   }

	// Write the timestamp array
	fflush(out);
	Fseek(out,(off64_t)0,SEEK_END);
	if (fwrite((char *)s_timestamp,sizeof(unsigned long long),s.FrameCount,out) != s.FrameCount) {
		Print("Error writing timestamp trailer for SER file\n");
		exit(1);
		}

	fclose(out);
	ser_last = NULL;
	return 1;
	}

int ShowHistogram_img(struct Image *img, int print_to_stdout)
	{
	if (img->depth == 8)  return ShowHistogram_8(img->data, img->width * img->height, print_to_stdout);
	if (img->depth == 16) return ShowHistogram_16((unsigned short *)img->data, img->width * img->height, print_to_stdout);

	Print("ShowHistogram: depth %d unsupported\n",img->depth);
	return 0;
	}

int ShowHistogram_16(unsigned short *ptr, int npix, int print_to_stdout)
	{
	int i,j,u_count;
	unsigned short *Hist = (unsigned short *)ZeroMalloc(65536 * sizeof(short));
	int min=65535, max=0;

	if (print_to_stdout) { printf("histogram 16bpp\n"); fflush(stdout); }
	for(i=0; i<npix; ++i) { 
	   Hist[ptr[i]]++;
	   if (ptr[i]<min) min=ptr[i];
	   if (ptr[i]>max) max=ptr[i];
	   }

	for(i=u_count=0; i<65536; ++i) if (Hist[i]) u_count++;
	if (print_to_stdout) printf("%d unique values detected: min=%d max=%d\n",u_count,min,max);

	for(i=0; print_to_stdout && i<65536; i+=8) {
   	   int tot=0; for(j=i; j<i+8; ++j) tot += Hist[j];
   	   if (tot) {
		for(j=i; j<i+8; ++j) { printf("%04d:%04d ",j,Hist[j]); }
   	   	printf("\n");
		}
	   }
	   
	if (print_to_stdout) {
	   printf("%d unique values detected: min=%d max=%d\n",u_count,min,max);
	   printf("\n\n");
	   fflush(stdout);
	   }

	free(Hist);
	return u_count;
	}
	
int ShowHistogram_8(unsigned char *ptr, int npix, int print_to_stdout)
	{
	int i,j,u_count;
	unsigned short Hist[256];
	int min=255, max=0;

	if (print_to_stdout) { printf("histogram 8bpp\n"); fflush(stdout); }
	for(i=0; i<256; ++i) Hist[i]=0;
	for(i=0; i<npix; ++i) {
 	   Hist[ptr[i]]++;
	   if (ptr[i]<min) min=ptr[i];
	   if (ptr[i]>max) max=ptr[i];
	   }

	for(i=u_count=0; i<256; ++i) if (Hist[i]) u_count++;
	if (print_to_stdout) printf("%d unique values detected: min=%d, max=%d\n",u_count,min,max);

	for(i=0; print_to_stdout && i<256; i+=8) {
   	   for(j=i; j<i+8; ++j) { printf("%04d:%04d ",j,Hist[j]);}
   	   printf("\n");
	   }

	if (print_to_stdout) {
	   printf("%d unique values detected: min=%d max=%d\n",u_count,min,max);
	   printf("\n\n");
	   fflush(stdout);
	   }

	return u_count;
	}

int
ApplyRowBoost(struct Image *img)
	{
	int width = img->width;
	int height = img->height;
	int bpp = img->depth/8;
	unsigned short *udata = (unsigned short *)img->data;
	unsigned char *data = img->data;
	int x,y,o;
	register int v,ra8 = RowBoostAdd/8;

	if (! doRowBoost) return 1;

	if (bpp != 1 && bpp != 2) {
	   Print("RowBoost not supported for %d bpp\n",bpp);
	   exit(1);
	   }

	if (RowBoostType == RB_EVEN_MULT || RowBoostType == RB_EVEN_ADD) y=0;
	else if (RowBoostType == RB_ODD_MULT || RowBoostType == RB_ODD_ADD) y=1;

	while(y < height) {
	   o = y * width;
	   if (bpp==1) for(x=0; x<width; ++x) {
		v=data[o]; v = v * RowBoostMult + ra8;
		if (v>255) v=255; data[o++]=v;
		}
	   else if (bpp==2) for(x=0; x<width; ++x) {
		v=udata[o]; v = v * RowBoostMult + RowBoostAdd;
		if (v>65535) v=65535; udata[o++]=v;
		}

	   y += 2;
	   }

	return 1;
	}

// For the G3-U3 camera, automatically balance odd/even rows
int
AutoRowBalance(struct Image *img)
	{
	int width = img->width;
	int height = img->height;
	int bpp = img->depth/8;
	unsigned short *udata = (unsigned short *)img->data;
	unsigned char *data = img->data;
	int y;
	unsigned int thresh = ThreshHold/2;
	register double d,avg_e,avg_o;  // 0=even rows, 1=odd rows
	register int x,o,count_e,count_o;	// counts for even & odd

	if (bpp != 1 && bpp != 2) {
	   Print("AutoRowBalance not supported for %d bpp\n",bpp);
	   exit(1);
	   }

	// add up totals for odd & even rows for all pixels above threshhold

	avg_e=avg_o=0;
	count_e=count_o=0;

	if (bpp==2) {
	   thresh <<= 8;
	   for(y=0; y<height; y+=2) { o = y*width; for(x=0; x<width; ++x,++o) if (udata[o]>=thresh) { avg_e += udata[o]; count_e++; } }
	   for(y=1; y<height; y+=2) { o = y*width; for(x=0; x<width; ++x,++o) if (udata[o]>=thresh) { avg_o += udata[o]; count_o++; } }

	   if (count_e==0 || count_o==0) {
		Print("AutoRowBalance: Error: No data\n");
		return 0;
		}

	   // Print("AutoRowBalance: Counted %d Even, %d Odd, total = %d pixels\n",count_e,count_o,count_e+count_o);
	   avg_e /= count_e; avg_o /= count_o;

	   if (avg_e==0) { Print("AutoRowBalance: No data in even rows\n"); return 1; }
	   if (avg_o==0) { Print("AutoRowBalance: No data in odd rows\n"); return 1; }
	   d = avg_o/avg_e;
	   
	   Print("AutoRowBalance: Even avg=%lf, Odd avg=%lf, d = %lf\n",avg_e,avg_o,d);

	   // Alter the even rows
	   for(y=0; y<height; y+=2) {
		register int v;
		o = y*width;
		for(x=0; x<width; ++x,++o) { v = udata[o]; v *= d; if (v>65535)v=65535; udata[o]=v; }
		}
	   }
	else if (bpp==1) {
	   for(y=0; y<height; y+=2) { o = y*width; for(x=0; x<width; ++x,++o) if (data[o]>=thresh) { avg_e += data[o]; count_e++; } }
	   for(y=1; y<height; y+=2) { o = y*width; for(x=0; x<width; ++x,++o) if (data[o]>=thresh) { avg_o += data[o]; count_o++; } }

	   if (count_e==0 || count_o==0) {
		Print("AutoRowBalance: Error: No data\n");
		return 0;
		}

	   // Print("AutoRowBalance: Counted %d Even, %d Odd, total = %d pixels\n",count_e,count_o,count_e+count_o);
	   avg_e /= count_e; avg_o /= count_o;

	   if (avg_e==0) { Print("AutoRowBalance: No data in even rows\n"); return 1; }
	   if (avg_o==0) { Print("AutoRowBalance: No data in odd rows\n"); return 1; }
	   d = avg_o/avg_e;
	   
	   Print("AutoRowBalance: Even avg=%lf, Odd avg=%lf, d = %lf\n",avg_e,avg_o,d);

	   // Alter the even rows
	   if (d<0.999 || d > 1.001) for(y=0; y<height; y+=2) {
		register int v;
		o = y*width;
		for(x=0; x<width; ++x,++o) { v = data[o]; v *= d; if (v>255)v=255; data[o]=v; }
		}
	   }

	return 1;
	}
