/*
 * xvfits.c - load/save routines for 'fits' format pictures
 *
 * LoadFITS(fname, pinfo, quick)  -  loads a FITS file
 * WriteFITS(fp,pic,ptype,w,h,r,g,b,numcols,style,raw,cmt,comment)
 */

/*
 * Copyright 1992, 1993 by David Robinson.
 *
 * The software may be modified for your own purposes.
 */

#include "xv.h"

#define NCARDS (36)
#define BLOCKSIZE (2880)

/* data types */
enum datatype { T_INT, T_LOG, T_STR, T_NOVAL };

typedef struct {
   FILE *fp;             /* file pointer */
   int bitpix;           /* number of bits per pixel */
   int size;             /* Size of each pixel, in bytes */
   int naxis;            /* number of axes */
   long int axes[3];     /* size of each axis */
   long int ndata;       /* number of elements in data */
   long int cpos;        /* current position in data file */
   char *comment;        /* malloc'ed comment string, or NULL if none */
} FITS;

/* static data */

/* block workspace of BLOCKSIZE bytes */
static char *fits_block=NULL;

/* end static data */


#ifdef __STDC__
#define __(x) x
#else
#define __(x) ()
#endif

static int splitfits __((byte *image, char *comment, int nx, int ny, int nz,
    char *basename));
static char *ftopen3d __((FITS *fs, char *file, int *nx, int *ny, int *nz,
    int *bitpix));
static void ftclose __((FITS *fs));
static int ftgbyte __((FITS *fs, unsigned char *buffer, int nelem));
static char *rdheader __((FITS *fs));
static char *wrheader __((FILE *fp, int nx, int ny, char *comment));
static char *rdcard __((char *card, char *name, enum datatype dtype,
    long int *kvalue));
static void wrcard __((char *card, char *name, enum datatype dtype,
    int kvalue, char *svalue));
static int ftgdata __((FITS *fs, void *buffer, int nelem));
static void ftfixdata __((FITS *fs, void *buffer, int nelem));
static void flip __((byte *buffer, int nx, int ny));

#undef __

/*******************************************/
int LoadFITS(fname, pinfo, quick)
     char    *fname;
     PICINFO *pinfo;
     int      quick;
/*******************************************/
{
  /* returns '1' on success */

  FITS fs;
  int i, nx, ny, nz, bitpix, np, nrd, ioerror;
  byte *image;
  char *error;
  char basename[64];

  if (fits_block == NULL)
  {
      fits_block = malloc(BLOCKSIZE);
      if (fits_block == NULL)
	  FatalError("Insufficient memory for FITS block buffer");
  }

  /* open the file */
  error = ftopen3d(&fs, fname, &nx, &ny, &nz, &bitpix);
  if (error != NULL)
  {
      SetISTR(ISTR_WARNING, "%s", error);
      return 0;
  }
  if (quick) nz = 1;  /* only load first plane */
  np = nx * ny * nz;
/* allocate memory for image and read it in */
  image = (byte *)malloc(np);
  if (image == NULL) FatalError("Insufficient memory for image");

/*
 * Read in image. For three dimensional images, read it in in one go
 * to ensure that we get that same scaling for all planes.
 */
  nrd = ftgbyte(&fs, image, np);
  ioerror = ferror(fs.fp);
  ftclose(&fs);

  if (nrd == 0)  /* didn't read any data at all */
  {
      if (ioerror)
	  SetISTR(ISTR_WARNING, "%s", "I/O error reading FITS file");
      else
	  SetISTR(ISTR_WARNING, "%s", "Unexpected eof reading FITS file");
      free(image);
      return 0;
  } else if (nrd < np)
  {       /* read partial image */
      if (ioerror)
	  SetISTR(ISTR_WARNING, "%s", "Truncated FITS file due to I/O error");
      else
	  SetISTR(ISTR_WARNING, "%s", "Truncated FITS file");
      memset(image+nrd, 0x80, np-nrd);  /* pad with grey */
  }

  if (nz > 1)
  {
/* how many planes do we actually have? */
      nz = (nrd-1)/(nx*ny) + 1;
/* returns how many sub-files created */
      nz = splitfits(image, fs.comment, nx, ny, nz, basename);
      np = nx * ny;
      image = (byte *)realloc(image, np);  /* throw away all but first */
  }

/* There seems to be a convention that fits files be displayed using
 * a cartesian coordinate system. Thus the first pixel is in the lower left
 * corner. Fix this by reflecting in the line y=ny/2.
 */
  flip(image, nx, ny);

/* Success! */
  pinfo->pic = image;
  pinfo->w = nx;
  pinfo->h = ny;
  for (i=0; i < 256; i++) pinfo->r[i] = pinfo->g[i] = pinfo->b[i] = i;
  pinfo->frmType = F_FITS;
  pinfo->colType = F_GREYSCALE;
  sprintf(pinfo->fullInfo, "FITS, bitpix: %d", bitpix);
  sprintf(pinfo->shrtInfo, "%dx%d FITS.", nx, ny);
  pinfo->comment = fs.comment;
  if (nz > 1)
  {
      pinfo->numpages = nz;
      strcpy(pinfo->pagebname, basename);
  }

  return 1;
}  
	  
/*******************************************/
int WriteFITS(fp,pic,ptype,w,h,rmap,gmap,bmap,numcols,colorstyle,comment)
     FILE *fp;
     byte *pic;
     int   ptype, w,h;
     byte *rmap, *gmap, *bmap;
     int   numcols, colorstyle;
     char *comment;
{
  register int i, j, np, nend;
  register byte *ptr;
  char *error;
  byte rgb[256];

  if (fits_block == NULL)
  {
      fits_block = malloc(BLOCKSIZE);
      if (fits_block == NULL)
	  FatalError("Insufficient memory for FITS block buffer");
  }

  error = wrheader(fp, w, h, comment);
  if (error != NULL)
  {
    SetISTR(ISTR_WARNING, "%s", error);
    return -1;
  }

  if (ptype == PIC8)
  {
/* If ptype is PIC8, we need to calculate the greyscale colourmap */
    for (i=0; i < numcols; i++) rgb[i] = MONO(rmap[i], gmap[i], bmap[i]);
/* flip line ordering when writing out */
    for (i=h-1; i >= 0; i--)
    {
	ptr = &pic[i*w];
	for (j=0; j < w; j++, ptr++) putc(rgb[*ptr], fp);
    }
  } else
  {
/* flip line ordering when writing out */
    for (i=h-1; i >= 0; i--)
    {
	ptr = &pic[i*w*3];
	for (j=0; j < w; j++, ptr += 3) putc(MONO(ptr[0], ptr[1], ptr[2]), fp);
    }
     
  }
  np = w*h;
/* nend is the number of padding characters at the end of the last block */
  nend = ((np+BLOCKSIZE-1)/BLOCKSIZE)*BLOCKSIZE - np;
  if (nend) for (i=0; i < nend; i++) putc('\0', fp);
 
  return 0;
}


/*
 * Given a 3 dimensional FITS image, this splits it up into nz 2-d files.
 * It returns the number of files actually stored.
 * If only one file could be written, then no split files are created.
 * It returns the basename of the split files in bname.
 * If there was a problem writing files, then a error message will be set.
 */
static int splitfits(image, comment, nx, ny, nz, basename)
byte *image;
char *comment;
int nx, ny, nz;
char *basename;
{
    int i, np=nx * ny, ioerror, nwrt;
    FILE *fp;
    char *error;
    byte *work;
    char filename[70];

#ifndef VMS
    sprintf(basename, "%s/xvpgXXXXXX", tmpdir);
#else
    sprintf(basename, "Sys$Disk:[]xvpgXXXXXX");
#endif

    mktemp(basename);
    if (basename[0] == '\0')
    {   /* mktemp() blew up */
	SetISTR(ISTR_WARNING, "%s", "Unable to build temporary filename");
	return 1;
    }
    strcat(basename, ".");

    error = NULL; /* no error */
    for (i=0; i < nz && error == NULL; i++)
    {
	sprintf(filename, "%s%d", basename, i+1);
	fp = fopen(filename, "wb");
	if (fp == NULL)
	{
	    error = "Unable to open temporary file";
	    break;
	}
        if (wrheader(fp, nx, ny, comment) != NULL)
	{
	    error = "I/O error writing temporary file";
	    fclose(fp);
	    unlink(fp);
	    break;
	}
	nwrt = fwrite(image+i*np, sizeof(unsigned char), np, fp);
	fclose(fp);
	if (nwrt == 0)  /* failed to write any data */
	{
	    error = "I/O error writing temporary file";
	    unlink(fp);
	    break;
	} else if (nwrt < np)
	    error = "I/O error writing temporary file";
    }

/* at this point, i is the number of files created */
    if (i == 0) return 1;
    if (i == 1)
    {
	sprintf(filename, "%s1", basename);
	unlink(filename);
    }
    if (error != NULL) SetISTR(ISTR_WARNING, "%s", error);
    return i;
}
	
/* Writes a minimalist FITS file header */
static char *wrheader(fp, nx, ny, comment)
FILE *fp;
int nx, ny;
char *comment;
{
    char *block=fits_block;
    int i, j, wrotehist, lenhist;
    char *history="Written by xv 3.00";

    memset(block, ' ', BLOCKSIZE);

    lenhist = strlen(history);

    i = 0;
    wrcard(&block[80*i++], "SIMPLE", T_LOG, 1, NULL);   /* write SIMPLE card */
    wrcard(&block[80*i++], "BITPIX", T_INT, 8, NULL);   /* write BITPIX card */
    wrcard(&block[80*i++], "NAXIS", T_INT, 2, NULL);    /* write NAXIS card */
    wrcard(&block[80*i++], "NAXIS1", T_INT, nx, NULL);  /* write NAXIS1 card */
    wrcard(&block[80*i++], "NAXIS2", T_INT, ny, NULL);  /* write NAXIS2 card */
/* Write HISTORY keyword */
    wrcard(&block[80*i++], "HISTORY", T_STR, lenhist, history); 
    if (comment != NULL && *comment != '\0')
    {
	while (*comment == '\n') comment++;  /* Skip any blank lines */
	while (*comment != '\0')
	{
	    for (j=0; j < 72; j++)
		if (comment[j] == '\0' || comment[j] == '\n') break;
/* 
 * Check to see if it is an xv history record; if so, then avoid
 * duplicating it.
 */
	    if (j != lenhist || memcmp(comment, history, j) != 0)
		wrcard(&block[80*i++], "COMMENT", T_STR, j, comment);
	    if (i == NCARDS)  /* Filled up a block */
	    {
		i = fwrite(block, sizeof(char), BLOCKSIZE, fp);
		if (i != BLOCKSIZE) return "Error writing FITS file";
		i = 0;
		memset(block, ' ', BLOCKSIZE);
	    }
		
	    comment += j;
	    while (*comment == '\n') comment++;  /* Skip any blank lines */
	}
    }
    wrcard(&block[80*i++], "END", T_NOVAL, 0, NULL);    /* write END keyword */
    i = fwrite(block, sizeof(char), BLOCKSIZE, fp);
    if (i != BLOCKSIZE) return "Error writing FITS file";
    return NULL;
}

/* open a 2 or 3-dimensional fits file.
 * Stores the dimensions of the file in nx, ny and nz, and updates the FITS
 * structure passed in fs.
 * If successful, returns NULL otherwise returns an error message.
 * Will return an error message if the primary data unit is not a
 * 2 or 3-dimensional array.
 */
static char *ftopen3d(fs, file, nx, ny, nz, bitpix)
FITS *fs;
char *file;
int *nx, *ny, *nz, *bitpix;
{
    FILE *fp;
    int naxis, i;
    char *error;
    
    fp = fopen(file, "rb");
    if (fp == NULL) return "Unable to open FITS file";

    fs->fp = fp;
    fs->bitpix = 0;
    fs->naxis = 0;
    fs->cpos = 0;

/* read header */
    error = rdheader(fs);

    if (error != NULL)
    {
	ftclose(fs);
	return error;
    }

    naxis = fs->naxis;

/* get number of data */
    fs->ndata = 1;
    for (i=0; i < naxis; i++) fs->ndata = fs->ndata*fs->axes[i];
    
    *nx = fs->axes[0];
    *ny = fs->axes[1];
    if (naxis == 2) *nz = 1;
    else *nz = fs->axes[2];

    *bitpix = fs->bitpix;

    return NULL;
}

/* closes a fits file */
static void ftclose(fs)
FITS *fs;
{
    if (fs == NULL) return;
    if (fs->fp != NULL) fclose(fs->fp);
}

/* reads the fits header, and updates the FITS structure fs.
 * Returns NULL on success, or an error message otherwise.
 */
static char *rdheader(fs)
FITS *fs;
{
    int i, j, res, commlen, commsize;
    char name[9];
    char *block=fits_block, *p;
    char *error;
    long int val;         /* the value */

    fs->comment = NULL;
    commlen = 0;
    commsize = 256;

    res = fread(block, sizeof(char), BLOCKSIZE, fs->fp);
    if (res != BLOCKSIZE) return "Error reading FITS file";
    i = 0;

/* read SIMPLE key */
    error = rdcard(block, "SIMPLE", T_LOG, &val);
    if (error != NULL) return error;
    if (val == 0) return "Not a SIMPLE FITS file";
    i++;

/* read BITPIX key */
    error = rdcard(&block[80], "BITPIX", T_INT, &val);
    if (error != NULL) return error;
    if (val != 8 && val != 16 && val != 32 && val != 64 && val != -32 &&
	val != -64)
	   return "Bad BITPIX value in FITS file";
    fs->bitpix = val;
    j = fs->bitpix;
    if (j < 0) j = -j;
    fs->size = j/8;
    i++;

/* read NAXIS key */
    error = rdcard(&block[2*80], "NAXIS", T_INT, &val);
    if (error != NULL) return error;
    if (val < 0 || val > 999)
	return "Bad NAXIS value in FITS file";
    if (val == 0) return "FITS file does not contain an image";
    if (val < 2) return "FITS file has fewer than two dimensions";
    fs->naxis = val;
    i++;

/* read NAXISnnn keys.
 * We allow NAXIS to be > 3 iff the dimensions of the extra axes are 1
 */
    for (j=0; j < fs->naxis; j++)
    {
	if (i == NCARDS)
	{
	    res = fread(block, sizeof(char), BLOCKSIZE, fs->fp);
	    if (res != BLOCKSIZE) return "Error reading FITS file";
	    i = 0;
	}

	sprintf(name, "NAXIS%d", j+1);
	error = rdcard(&block[i*80], name, T_INT, &val);
	if (error != NULL) return error;
	if (val < 0) return "Bad NAXISn value in FITS file";
	if (val == 0) return "FITS file does not contain an image";
	if (j < 3) fs->axes[j] = val;
	else if (val != 1) return "FITS file has more than three dimensions";
	i++;
    }
    if (fs->naxis > 3) fs->naxis = 3;

/* do remainder, looking for comment cards */
/* Section 5.2.2.4 of the NOST standard groups COMMENT and HISTORY keywords
 * under 'Commentary Keywords'. Many applications write HISTORY records in
 * preference to COMMENT records, so we recognise both types here.
 */
    for (;;)
    {
	if (i == NCARDS)
	{
	    res = fread(block, sizeof(char), BLOCKSIZE, fs->fp);
	    if (res != BLOCKSIZE) return "Unexpected eof in FITS file";
	    i = 0;
	}
	p = &block[i*80];
	if (strncmp(p, "END     ", 8) == 0) break;
	if (strncmp(p, "HISTORY ", 8) == 0 || strncmp(p, "COMMENT ", 8) == 0)
	{
	    p += 8;  /* skip keyword */
	    for (j=71; j >= 0; j--) if (p[j] != ' ') break;
	    j++;     /* make j length of comment */
	    if (j > 0) /* skip blank comment cards */
	    {
		if (fs->comment == NULL)
		{
		    fs->comment = malloc(commsize);
		    if (fs->comment == NULL)
			FatalError("Insufficient memory for comment buffer");
		}
		if (commlen + j + 2 > commsize) /* if too small */
		{
		    char *new;
		    commsize += commsize;      /* double size of array */
		    new = malloc(commsize);
		    if (new == NULL)
			FatalError("Insufficient memory for comment buffer");
		    if (commlen) memcpy(new, fs->comment, commlen);
		    free(fs->comment);
		    fs->comment = new;
		}
		memcpy(&fs->comment[commlen], p, j);  /* add string */
		commlen += j;
		fs->comment[commlen++] = '\n';       /* with trailing cr */
		fs->comment[commlen] = '\0';
	    }
	}
	i++;
    }
    return NULL;
}

/* write a header record into the 80 byte buffer card.
 * The keyword name is passed in name. The value type is in dtype; this
 * can have the following values:
 *    dtype = T_NOVAL
 *         no keyword value is written
 *    dtype = T_LOG
 *         a logical value, either 'T' or 'F' in column 30 is written
 *    dtype = T_INT
 *         an integer is written, right justified in columns 11-30
 *    dtype = T_STR
 *         a string 'svalue' of length kvalue is written, in columns 9-80.
 */
static void wrcard(card, name, dtype, kvalue, svalue)
char *card, *name;
enum datatype dtype;   /* type of value */
int kvalue;
char *svalue;
{
    int l;
    memset(card, ' ', 80);
    l = strlen(name);
    if (l) memcpy(card, name, l);   /* copy name */

    if (dtype == T_NOVAL) return;

    if (dtype == T_STR)
    {
	l = kvalue;
	if (l <= 0) return;
	if (l > 72) l = 72;
	memcpy(&card[8], svalue, l);
	return;
    }

    card[8] = '=';

    if (dtype == T_LOG)
	card[29] = kvalue ? 'T' : 'F';
    else /* T_INT */
    {
	sprintf(&card[10], "%20d", kvalue);
	card[30] = ' ';
    }
}

/* Read a header record, from the 80 byte buffer card.
 * the keyword name must match 'name'; and parse its value according to
 * dtype. This can have the following values:
 *    dtype = T_LOG
 *        value is logical, either 'T' or 'F' in column 30.
 *    dtype = T_INT
 *        value is an integer, right justified in columns 11-30.
 *
 * The value is stored in kvalue.
 * It returns NULL on success, or an error message otherwise.
 */
static char *rdcard(card, name, dtype, kvalue)
char *card, *name;
enum datatype dtype;   /* type of value */
long int *kvalue;
{
    int i, ptr;
    char namestr[9];
    static char error[45];

    memcpy(namestr, card, 8);
    i = 8;
    do i--; while (i >= 0 && namestr[i] == ' ');
    namestr[i+1] = '\0';

    if (strcmp(namestr, name) != 0)
    {
	sprintf(error, "Keyword %s not found in FITS file", name);
	return error;
    }

/* get start of value */
    ptr = 10;
    while (ptr < 80 && card[ptr] == ' ') ptr++;
    if (ptr == 80) return "FITS file has missing keyword value"; /* no value */

    if (dtype == T_LOG)
    {
	if (ptr != 29 || (card[29] != 'T' && card[29] != 'F'))
	    return "Keyword has bad logical value in FITS file";
	*kvalue = (card[29] == 'T');
    } else /* an integer */
    {
	int j;
	long int ival;
	char num[21];

	if (ptr > 29) return "Keyword has bad integer value in FITS file";
	memcpy(num, &card[ptr], 30-ptr);
	num[30-ptr] = '\0';
	j = sscanf(num, "%ld", &ival);
	if (j != 1) return "Keyword has bad integer value in FITS file";
	*kvalue = ival;
    }

    return NULL;
}


/* reads nelem values into the buffer.
 * returns NULL for success or an error message.
 * Copes with the fact that the last 2880 byte record of the FITS file
 * may be truncated, and should be padded out with zeros.
 *  bitpix   type of data
 *    8        unsigned char
 *   16        short int       (NOT 2-byte integer)
 *   32        int             (NOT 4-byte integer)
 *  -32        float
 *  -64        double
 *
 * Returns the number of elements actually read.
 */
static int ftgdata(fs, buffer, nelem)
FITS *fs;
void *buffer;
int nelem;
{
    int res;

    if (nelem == 0) return 0;

    res = fread(buffer, fs->size, nelem, fs->fp);
/* if failed to read all the data because at end of file */
    if (res != nelem && feof(fs->fp))
    {
/* nblock is the number of elements in a record. size is always a factor of
   BLOCKSIZE */
	int loffs, nblock=BLOCKSIZE/fs->size;
/*
 * the last record might be short; check this.
 * loffs is the offset of the start of the last record from the current
 * position.
 */
	loffs = ((fs->ndata+nblock-1)/nblock - 1)*nblock - fs->cpos;
/* if we read to the end of the penultimate record */
	if (res >= loffs)
	{
/* pad with zeros */
	    memset((char *)buffer+res*fs->size, '\0', (nelem-res)*fs->size);
	    res = nelem;
	}
    }
    fs->cpos += res;
    if (res) ftfixdata(fs, buffer, res);
    return res;
}

/* convert the raw data, as stored in the FITS file, to the format
 * appropiate for the data representation of the host computer.
 * Assumes that
 *  short int = 2 or more byte integer
 *  int       = 4 or more byte integer
 *  float     = 4 byte floating point, not necessarily IEEE.
 *  double    = 8 byte floating point.
 * This can exit for lack of memory on the first call with float or
 * double data.
 */
static void ftfixdata(fs, buffer, nelem)
FITS *fs;
void *buffer;
int nelem;
{
    register int i, n=nelem;
    register unsigned char *ptr=buffer;
/*
 * conversions. Although the data may be signed, reverse using unsigned 
 * variables.
 * Because the native int types may be larger than the types in the file,
 * we start from the end and work backwards to avoid overwriting data
 * prematurely.
 */
/* convert from big-endian two-byte signed integer to native form */
    if (fs->bitpix == 16)
    {
	unsigned short int *iptr=(unsigned short int *)ptr;
	iptr += n-1;    /* last short int */
	ptr += (n-1)*2; /* last pair of bytes */
	for (i=0; i < n; i++, ptr-=2, iptr--)
	    *iptr = (((int)*ptr) << 8) | (int)(ptr[1]);
    }
/* convert from big-endian four-byte signed integer to native form */
    else if (fs->bitpix == 32)
    {
	unsigned int *iptr = (unsigned int *)ptr;
	iptr += n-1;    /* last integer */
	ptr += (n-1)*4; /* last 4 bytes */
	for (i=0; i < n; i++, ptr-=4, iptr--)
	    *iptr = (((unsigned int)*ptr) << 24) |
		    ((unsigned int)ptr[1] << 16) |
		    ((unsigned int)ptr[2] << 8)  |
		    (unsigned int)ptr[3];
    }
/* convert from IEE 754 single precision to native form */
    else if (fs->bitpix == -32)
    {
	register int j, k, expo;
	static float *exps=NULL;

	if (exps == NULL)
	{
	    exps = (float *)malloc(256 * sizeof(float));
	    if (exps == NULL) FatalError("Insufficient memory for exps store");
	    exps[150] = 1.;
	    for (i=151; i < 256; i++) exps[i] = 2.*exps[i-1];
	    for (i=149; i >= 0; i--) exps[i] = 0.5*exps[i+1];
	}
	      

	for (i=0; i < n; i++, ptr+=4)
	{
	    k = (int)*ptr;
	    j = ((int)ptr[1] << 16) | ((int)ptr[2] << 8) | (int)ptr[3];
	    expo = ((k & 127) << 1) | (j >> 23);
	    if ((expo | j) == 0) *(float *)ptr = 0.;
	    else 
		*(float *)ptr = exps[expo]*(float)(j | 0x800000);
	    if (k & 128) *(float *)ptr = - *(float *)ptr;
	}
/* convert from IEE 754 double precision to native form */
    } else if (fs->bitpix == -64)
    {
	register int expo, k, l;
	register unsigned int j;
	static double *exps=NULL;

	if (exps == NULL)
	{
	    exps = (double *)malloc(2048 * sizeof(double));
	    if (exps == NULL) FatalError("Insufficient memory for exps store");
	    exps[1075] = 1.;
	    for (i=1076; i < 2048; i++) exps[i] = 2.*exps[i-1];
	    for (i=1074; i >= 0; i--) exps[i] = 0.5*exps[i+1];
	}
	      
	for (i=0; i < n; i++, ptr+=8)
	{
	    k = (int)*ptr;
	    j = ((unsigned int)ptr[1] << 24) | ((unsigned int)ptr[2] << 16) |
		((unsigned int)ptr[3] << 8) | (unsigned int)ptr[4];
	    l = ((int)ptr[5] << 16) | ((int)ptr[6] << 8) | (int)ptr[7];
	    expo = ((k & 127) << 4) | (j >> 28);
	    if ((expo | j | l) == 0) *(double *)ptr = 0.;
	    else 
		*(double *)ptr = exps[expo] * (16777216. *
			    (double)((j&0x0FFFFFFF)|0x10000000) + (double)l);
	    if (k & 128) *(double *)ptr = - *(double *)ptr;
	}
    }
}

#define maxmin(x, max, min) {\
 maxmin_t=(x); if(maxmin_t > max) max=maxmin_t; if (maxmin_t<min) min=maxmin_t;}

/* Reads a byte image from the FITS file fs. The image contains nelem pixels.
 * If bitpix = 8, then the image is loaded as stored in the file.
 * Otherwise, it is rescaled so that the minimum value is stored as 0, and
 * the maximum is stored as 255.
 * Returns the number of pixels read.
 */
static int ftgbyte(fs, cbuff, nelem)
FITS *fs;
unsigned char *cbuff;
int nelem;
{
    void *voidbuff;
    register int i, n, nrd;

/* if the data is unsigned char, then read it directly */
    if (fs->bitpix == 8)
	return ftgdata(fs, cbuff, nelem);

/* allocate a buffer to store the image */
    if (fs->bitpix == 16)
	voidbuff = (void *)malloc(nelem * sizeof(short int));
    else if (fs->bitpix == 32)
	voidbuff = (void *)malloc(nelem * sizeof(int));
    else
	voidbuff = (void *)malloc(nelem*fs->size);  /* float, double */
    if (voidbuff == NULL)
    {
	char emess[60];
	sprintf(emess, "Insufficient memory for raw image of %d bytes", 
		nelem*fs->size);
	FatalError(emess);
    }
    nrd = ftgdata(fs, voidbuff, nelem);
    if (nrd == 0) return 0;
    n = nrd;
/* convert short int to unsigned char */
    if (fs->bitpix == 16)
    {
	register short int *buffer=voidbuff;
	register int max, min, maxmin_t;
	register float scale;

	min = max = buffer[0];
	for (i=1; i < n; i++, buffer++) maxmin(*buffer, max, min);
	scale = (max == min) ? 0. : 255./(float)(max-min);

/* rescale and convert */
	for (i=0, buffer=voidbuff; i < n; i++)
	    cbuff[i] = (unsigned char)(scale*(float)((int)buffer[i]-min));

/* convert long int to unsigned char */
    } else if (fs->bitpix == 32)
    {
	register int *buffer=voidbuff;
	register int max, min, maxmin_t;
	register float scale, fmin;

	min = max = buffer[0];
	for (i=1; i < n; i++, buffer++) maxmin(*buffer, max, min);
	scale = (max == min) ? 1. : 255./((double)max-(double)min);
	fmin = (float)min;

/* rescale and convert */
	if (scale < 255./2.1e9) /* is max-min too big for an int ? */
	    for (i=0, buffer=voidbuff; i < n; i++)
		cbuff[i] = (unsigned char)(scale*((float)buffer[i]-fmin));
	else /* use integer subtraction */
	    for (i=0, buffer=voidbuff; i < n; i++)
		cbuff[i] = (unsigned char)(scale*(float)(buffer[i]-min));

/* convert float to unsigned char */
    } else if (fs->bitpix == -32)
    {
	register float *buffer=voidbuff;
	register float max, min, maxmin_t, scale;

	min = max = buffer[0];
	for (i=1; i < n; i++, buffer++) maxmin(*buffer, max, min);
	scale = (max == min) ? 0. : 255./(max-min);

/* rescale and convert */
	for (i=0, buffer=voidbuff; i < n; i++)
	    cbuff[i] = (unsigned char)(scale*(buffer[i]-min));

/* convert double to unsigned char */
    } else if (fs->bitpix == -64)
    {
	register double *buffer=voidbuff;
	register double max, min, maxmin_t, scale;

	min = max = buffer[0];
	for (i=1; i < n; i++, buffer++) maxmin(*buffer, max, min);
	scale = (max == min) ? 0. : 255./(max-min);

/* rescale and convert */
	for (i=0, buffer=voidbuff; i < n; i++)
	    cbuff[i] = (unsigned char)(scale*(buffer[i]-min));
    }

    free(voidbuff);
    return n;
}

#undef maxmin

/* reverse order of lines in image */
static void flip(buffer, nx, ny)
byte *buffer;
int nx;
int ny;
{
    int i;
    register int j, v;
    register byte *buff1, *buff2;

    for (i=0; i < ny/2; i++)
    {
	buff1 = &buffer[i*nx];
	buff2 = &buffer[(ny-1-i)*nx];
	for (j=0; j < nx; j++)
	{
	    v = *buff1;
	    *(buff1++) = *buff2;
	    *(buff2++) = v;
	}
    }
}