#include "ninox.h"

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

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, int bzero, double bscale)
	{
	struct Image *img = CreateImage();
	struct fits_info *fi = ZeroMalloc(sizeof(struct fits_info));

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

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

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

	fi->header = NULL;

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

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

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

	return img;
	}

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

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

static struct Image *
convert_to_fits(struct Image *img, int depth)
	{
	struct Image *dst;
	struct fits_info *fi = img->info;
	unsigned char *src,*data;
	unsigned short *udata,*usrc;
	unsigned int v;
	int i,o,npixels = img->width * img->height;
	int bpp = depth/8;
	char fname[128];
	
	src = img->data;
	usrc = (unsigned short *)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 {
		   Print("convert_to_fits: source image depth %d not supported\n",img->depth);
		   exit(1);
		   }
		break;
	   case 24:
		if (img->depth==8)       for(i=o=0; i<npixels; ++i) {v=src[i];     data[o++]=v; data[o++]=v; data[o++]=v;}
		else if (img->depth==16) for(i=o=0; i<npixels; ++i) {v=usrc[i]>>8; data[o++]=v; data[o++]=v; data[o++]=v;}
		else {
		   Print("convert_to_fits: source image depth %d not supported\n",img->depth);
		   exit(1);
		   }
		break;
	   }

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

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

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

	// notreached
	return 0;
	}

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

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

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

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

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

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

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

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

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

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

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

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

	fi->bzero = 0;
	fi->bscale = 1.0;

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

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

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

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

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

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

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

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

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

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

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

		while(*ptr && ! isdigit(*ptr)) ++ptr;
#ifdef MSWIN32
		if (sscanf(ptr,"%u/%I64u",&secs,&tsc) == 2) {
#else
		if (sscanf(ptr,"%u/%llu",&secs,&tsc) == 2) {
#endif
		   img->tm_sec = secs;
		   img->tm_tsc = tsc;
		   }
		}
	   }

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

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

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

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

        return 1;
        }

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

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

	img->info = bi;

	/* Look for magic */
	if (! fread(type,2,1,in)) {
	   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)) {
	   Print("Short read\n");
	   exit(1);
	   }

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

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

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

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

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

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

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

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

	      /* Read colourmap */
	      if (H.cmapentries && fread(bi->cmap, H.cmapentries * 4,1, in)==0 ) {
		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);
	}

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

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

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

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

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

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

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

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

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

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

	      }
           free(rowbuf);

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

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

	img->data = buffer;

	return 1;
	}

// For 8 or 16 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 *uptr = (unsigned short *) 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 = uptr[o1];
			uptr[o1]=uptr[o2];
			uptr[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;
	   default:
		break; // nothing to do
	   }
	return 1;
	}


static int
load_fits_data(FILE *in, struct Image *img, int flags)
	{
	int bufsize = img->width * img->height * img->depth/8;
	struct fits_info *fi = img->info;
	int max=0,i,err;
	int nrecords, total, padding;
	int do_stretch = InputHistoStretch;

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

	img->data = NULL;

	if (flags == FITS_LOAD) {
           img->data = (unsigned char *)Malloc(bufsize);
           err = fread(img->data,1,bufsize,in);
           if (err != bufsize) {
              Print("Error in reading data for %s (%d x %d): only read %d bytes but expected %d bytes\n",
			img->src_fname,img->width,img->height,err,bufsize);
              exit(1);
              }
	   }
	else if (flags == FITS_SKIP){
	   // we are not going to use the data so just seek past it
	   err = fseek(in,bufsize,SEEK_CUR);
	   if (err) {
              Print("Error seeking on %s\n",img->src_fname);
              exit(1);
              }
	   }
	else {
	   Print("load_fits_data: Unknown flags value %d\n",flags);
	   exit(1);
	   }

       // Undo the endian-ness of the data for 16bpp FITS data
        if (img->depth == 16 && flags==FITS_LOAD) {
           unsigned char *tmp,*ptr;
           unsigned short *iptr;
	   int last_print=-1;

           ptr = img->data; iptr = (unsigned short *)ptr;
           i=img->width * img->height;
           while(i--) {
                int x = *ptr; *ptr = *(ptr+1); *(ptr+1) = x;

                if (fi->bzero) *iptr += fi->bzero;
                //if (fi->bscale != 1.0) *iptr *= fi->bscale;
		if (*iptr>max && *iptr < 65530) max=*iptr;
                ptr += 2; iptr++;
                }

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

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

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

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

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

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

	   // set the BZERO header so that images we write
	   // based on this data will be correct
	   fi->bzero = 32768;
           }
        else if (img->depth == 8 && flags==FITS_LOAD) {
           // Unsigned 8 bit data
           unsigned char *ptr = img->data;
           i = img->width * img->height;
           while(i--) {
                if (fi->bzero) *ptr += fi->bzero;
                //if (fi->bscale != 1.0) *ptr *= fi->bscale;
		if (*ptr>max && *ptr < 250) max=*ptr;
                ptr++;
                }

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

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

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

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

	   // flip to to bottom
	   InvertImage(img);

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

	// Read the padding data
	if ((bufsize % 2880) != 0) {
	   nrecords = bufsize / 2880 + 1;
	   total = nrecords * 2880;
	   padding = total - bufsize;
	   if (padding > 0) {
		err = fseek(in,padding,SEEK_CUR);
	         if (err) {
		   Print("Error seeking for %s\n",img->src_fname);
		   exit(1);
		   }
		}
	   }
	   
	return 1;
	}

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

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

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

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

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

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

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

struct Image *
LoadImage(char *src_fname,char *dst_fname)
	{
	FILE *in;
	struct Image *img = CreateImage();

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

	img->type = parse_file_type(src_fname);

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

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

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

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

        //if (!Quiet) Print("(%dx%dx%d) ",img->width,img->height,img->depth);

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

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

	fclose(in);
	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;
	struct Image *img;
	int ch;

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

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

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

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

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

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

	if (flags == ARCHIVE_LOAD) flags = FITS_LOAD;
	else flags = FITS_SKIP;

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

	*err = 1;
	return img;
	}

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

	int err;
	struct Image *img;

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

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

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

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

	   DestroyImage(img);
	   }

	fclose(in);
	}

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

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

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

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

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

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

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

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

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


static int
write_fits_header(FILE *out, struct Image *img)
	{
	static double last_bscale = -1;
	static int last_depth = -1;
	static int last_bzero = -1;
	static int last_w=-1, last_h=-1;
	static unsigned int 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;

	// 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=%21d%50s","BZERO",(int)fi->bzero,""); ++i;
	sprintf(buffer+i*80,"%-8.8s=%21d%50s","DATAMIN",0,""); ++i;
	sprintf(buffer+i*80,"%-8.8s=%21d%50s","DATAMAX",((img->depth==8)?255:65535),""); ++i;
	if (img->tm_sec || img->tm_tsc) {
#ifdef MSWIN32
	   sprintf(buffer+i*80,"%-8.8s= '__TIMESTAMP' %10u/%-20I64u%26s","COMMENT",img->tm_sec,img->tm_tsc,""); ++i;
#else
	   sprintf(buffer+i*80,"%-8.8s= '__TIMESTAMP' %10u/%-20llu%26s","COMMENT",img->tm_sec,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;
	}

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

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

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

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

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

        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;

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

#define NBUCKETS_16 	656
#define NBUCKETS_8	256

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 int
write_fits_data(FILE *out, struct Image *img)
	{
	int x1,y1,x2,y2,i,y;
	unsigned char *outbuffer;
	unsigned char *sptr,x;
	struct fits_info *fi = img->info;
	int max,bzero = fi->bzero;
	int ndatabytes;

	if (img->depth != 8 && img->depth != 16) {
	   Print("write_fits_data: Invalid depth %d\n",img->depth);
	   return 0;
	   }

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

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

        // Convert the buffer back to original format
        if (img->depth == 8) {
           unsigned char *ptr = img->data;
           unsigned char *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) {
		   unsigned int v = ptr[o];
		   v *= scale; if (v<0) v=0; if (v>255) v=255;
		   ptr[o] = v;
		   }
	 	}

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

           if (HistoProtect) {
              double scale = (double)White/255.0;
              unsigned int val = 255.0 * scale;

	      if (val<0) val=0; if (val>255) val=255;
              *(dptr++) = 0; *(dptr++) = val;
	      ptr += 2; i-=2;
              }

           while(i--) {
                //if (fi->bscale != 1.0) *(dptr++) /= fi->bscale;
                *(dptr++) = *(ptr++) - bzero;
                }
           }
        else if (img->depth == 16) {
           unsigned short *iptr =  (unsigned short *) img->data;
           unsigned short *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) {
		   unsigned int v = iptr[o];
		   v *= scale; if (v<0) v=0; if (v>65535) v=65535;
		   iptr[o] = v;
		   }
	 	}

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

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

           }

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

	free(outbuffer);

	// pad to multiple of 2880
	if (ndatabytes % 2880) {
	   char ch=0;
	   int padbytes = 2880 - (ndatabytes%2880);
	   while(padbytes--) fwrite(&ch,1,1,out);
	   }

        return(1);
	}

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

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

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

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

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

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

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

	fclose(out);
	return(1);
	}

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

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

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

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

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

	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("ConvertToBMP: Buffer conversion failed\n");
	   return NULL;
	   }

	fi = ZeroMalloc(sizeof(struct fits_info));

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

	fi->bscale = 1;

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

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

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

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

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

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

	if (newtype == img->type && newdepth == img->depth) return img;

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

