/******************************************************************************
*    fitstohdf
*
*
* Purpose:  to convert image data stored in FITS files to hdf scientific 
*           data set (SDS) format, storing the results in an hdf file.  
*           Currently only works with simplei FITS files, ignoring all
*           keywords except SIMPLE, BITPIX, NAXIS, NAXISn, BZERO, BSCALE,
*           and END. 
*
*                                   -----------
*                                  |           |
*       FITS image data   -------> | fitstohdf | ----------> SDS
*       (simple FITS file)         |           | 
*                                   -----------
*
*
*  Synopsis:
*       fitstohdf infile outfile 
*
*       infile: input file, containing at least the following keywords
*               in the header:
*                   SIMPLE, BITPIX, NAXIS, NAXISn, BZERO, BSCALE, END. 
*
*       outfile: output file, in hdf format, containing a single SDS. 
*
*
*  Notes:
*      1. This is a very preliminary version, designed primarily to elicit
*  comments.  Please let the author know what works and what doesn't
*  work, and what features you would like to see added or changed.
*  In particular, what features would you like to see made optional.
*
*      2. The treatment of BLANK, for instance, almost certainly needs work.
*
*      3. This version does not handle the case BITPIX = -32.  That option
*  will be added.
*
*      4. This version writes both an SDS and a corresponding raster image.
*  The presence of the raster image makes debugging easier, because it
*  lets you look at the corresponding image immediately.  
*
*
*  National Center for Supercomputing Applications
*  University of Illinois, Urbana-Champaign
*
*  by Mike Folk (mfolk@ncsa.uiuc.edu; 217-244-0647)
*  Beta version: 02/21/90
*
*  Modified by Peter Webb, 10/15/90
*    - New error handling
*    - New interface
*
*  This program is in the public domain
*
*/

#include <ctype.h>
#include <stdio.h>

/* Package include files */

#include "fits2hdf.h"
#include "types.h"      /* Package-wide type definitions */
#include "error.h"      /* Error codes */
#include "extern.h"     /* External routines and global variables */
#include "reformat.h"   /* Package-wide contants */



/**********************************************************************
** get_required_keywords: read in required keyword records: SIMPLE, 
**                        BITPIX, NAXIS, NAXIS1,... 
**********************************************************************/
static ErrorCode get_required_keywords(ff)
FITSFILE *ff;
{
    TOKEN token;
    char c;
    int dim, i, oldrank;

    if ( 0 > check_next_card(ff, &token, SIMPLE))
       return(err_msg(Fits2HDF, CorruptedInputFile));

    ff->val.simple = c = token.str_val[0];          /* SIMPLE */
    if ( (c!='T') && (c!='t') ) {  
      return(err_msg(Fits2HDF, BadInputVersion));
    }

    if ( 0 > check_next_card(ff, &token, BITPIX))
       return(err_msg(Fits2HDF, CorruptedInputFile));

    ff->val.bitpix = atoi(token.str_val);           /* BITPIX */

    if ( 0 > check_next_card(ff, &token, NAXIS))
       return(err_msg(Fits2HDF, CorruptedInputFile));

    ff->val.rank = token.int_val;                   /* NAXIS */
    ff->val.dimsizes = (int *) malloc(ff->val.rank*sizeof(int32));

/* get dimensions */

    oldrank = ff->val.rank;  /* use oldrank in case ff->val.rank gets 
                               changed in loop. (ff->val.rank gets
                               decreased by 1 for each dim value of 1) */
    for (dim=0, i=0; i < oldrank; i++) {
        if ( 0 > check_next_card(ff, &token, DIMSIZE))
	  return(err_msg(Fits2HDF, CorruptedInputFile));
        if ( (token.subclass == (i+1)) ) {  /* make sure NAXIS card correct */
            if (token.int_val == 1)
                ff->val.rank--;              /* if dimension = 1, ignore it */
            else 
                ff->val.dimsizes[dim++] = token.int_val;     /* else add it */
        } else 
	  return(err_msg(Fits2HDF, CorruptedInputFile));
    }
    return(AllOk);
}

/*****************************************************************
** get_data_and_image: get the data; convert to float format AND image 
*****************************************************************/

static ErrorCode get_data_and_image(file, elementsize, tot_elements, outdata,
				    image)
FITSFILE *file;
int elementsize, tot_elements;
float *outdata;
unsigned char *image;
{
    int    count, i, j, x; 
    char  *data2, *data2ptr;          /* for input if bitpix = 16 */
    char  *data4, *data4ptr;          /* for input if bitpix = 32 */
    char  c, neg;
    float *outdataptr;                /* for float output */
    float  bzero, bscale, temp, spread;
    unsigned char *imptr;              /* for image output */
 
    outdataptr = outdata;
    bzero = file->val.bzero;
    bscale = file->val.bscale;

/* Compute value to store as "blank" in output.  */
/* NOTE: This procedure is especially designed to produce nice    */
/*       images using NCSA tools, but is not necessarily suited   */
/*       to other applications.  Other options should be provided */
/*       in future versions of this program.                      */
/* The procedure: if max and min are provided, make "blank a bit below */
/*       min, so it will take on the background color when displayed.  */
/*       Otherwise assign it a "NOTaNUMBER" value.                */

    if (file->opt.ismax && file->opt.ismin) {
        temp = (file->val.max - file->val.min)/IMAGE_SPREAD; 
        file->val.outblank =  file->val.min - temp;
    } else 
        file->val.outblank = NOTaNUMBER;


/* convert elements to floating point */
/* Two cases: elementsize==2 => bitpix==16          */
/*            elementsize==4 => bitpix==32          */
/* Not yet covered: bitpix = -32                    */    

    if (elementsize == 2) {
        data2 = data2ptr= (char *) malloc( tot_elements*elementsize);
        count = fread(data2, elementsize, tot_elements, file->fp); 
        c = *data2ptr++;          /* get high byte  */
        neg = 0x80 & c;           /* mask off sign bit */
        c = neg ? ~c : c;         /* if negative, take complement */
        x = c*256;                /* shift right 8 bits */
        c = *data2ptr++;          /* get low byte   */         
        c = neg ? ~c : c;         /* if neg, take complement */
        x += c;                   /* add to first half */
        x = neg ? -(x+1) : x;     /* if neg, taking2's complement, add 1 */

        if (x == file->val.blank && file->opt.isblank)  /* If BLANK, assign */
             *outdataptr++ = file->val.outblank;         /*    output blank  */
        else
             *outdataptr++ = bzero + bscale * (float) x;  /* else scale it */

/* same procedure for the rest of the numbers */

        for (i=1; i<tot_elements; i++) {
            c = *data2ptr++; 
            neg = 0x80 & c; 
            c = neg ? ~c : c;
            x = c*256;      
            c = *data2ptr++;
            c = neg ? ~c : c;
            x += c;         
            x = neg ? -(x+1) : x; 

            if (x == file->val.blank && file->opt.isblank) 
                 *outdataptr++ = file->val.outblank;    
            else
                 *outdataptr++ = bzero + bscale * (float) x;
        }

    } 
    else {   /* elementsize == 4 */
        data4 = data4ptr= (char *) malloc( tot_elements*elementsize);
        count = fread(data4, elementsize, tot_elements, file->fp); 
        c = *data4ptr++;          /* get high byte  */
        neg = 0x80 & c;           /* mask off sign bit */
        c = neg ? ~c : c;         /* if negative, take complement */
        x = c;                    /* store in x */
        for (j=0; j<3; j++) {
            c = *data4ptr++;      /* get next byte   */         
            c = neg ? ~c : c;     /* if neg, take complement */
            x = 256*x + c;        /* shift and add next byte */
        }
        x = neg ? -(x+1) : x;     /* if neg, taking 2's complement, add 1 */

        if (x == file->val.blank && file->opt.isblank)   /* if BLANK, assign */
             *outdataptr++ = file->val.outblank;         /*    output blank */
        else
             *outdataptr++ = bzero + bscale * (float) x;  /* else scale it */

/* same procedure for the rest of the numbers */

        for (i=1; i<tot_elements; i++) {  
            c = *data4ptr++;    
            neg = 0x80 & c;    
            c = neg ? ~c : c; 
            x = c;           
            for (j=0; j<3; j++) {
                c = *data4ptr++;    
                c = neg ? ~c : c;    
                x = 256*x + c;      
            }
            x = neg ? -(x+1) : x;  

            if (x == file->val.blank && file->opt.isblank) 
                 *outdataptr++ = file->val.outblank;    
            else
                 *outdataptr++ = bzero + bscale * (float) x;
        }
    }

    if ( count != tot_elements)
      return(err_msg(Fits2HDF, InternalError));

/* convert to image -- hard-wired for HDF Images */

    outdataptr = outdata;
    imptr = image;

/* If there is a max and a min, generate an image, else skip it */

    if (file->opt.ismax && file->opt.ismin) {
        spread = file->val.max - file->val.min;
        for (i=0; i<tot_elements; i++, outdataptr++, imptr++) {
            if (*outdataptr == file->val.outblank)
                *imptr = 0;
            else
                *imptr = ((*outdataptr - file->val.min) / spread) *
		  IMAGE_SPREAD;
        }
    } else
      return(err_msg(Fits2HDF, ConvertFailed));
    return(AllOk);
}


/******************************************************************
** check_next_card: get next card and check to see that it is in the
**                specified class; if not, report an error
*******************************************************************/
static int check_next_card(file, token, class)
FITSFILE *file;
TOKEN *token;
int class;
{
    next_card(file);
    if (0 > parse_keyword_card(file->card,token))
        return(-1);
    if (token->class != class) {          /* is class correct? */
      err_msg(Fits2HDF, BadFileData);
      return(-1);
    }
    return(1);
}


/**********************************************************************
** next_card: read next card image; if at end of block, get a new block
**********************************************************************/
static int next_card(file)
FITSFILE *file;
{
    int i;

    if (file->blockptr >= BLOCKLENGTH ) {
        file->blockptr = 0;
        if (0 == fread(file->block, BLOCKLENGTH, 1, file->fp) ) {
	  err_msg(Fits2HDF, FileSystemError);
	  return(-1);
        }
    }
    for (i=0; i<CARDLENGTH; )
        file->card[i++] = file->block[file->blockptr++];
    file->card[i] = '\0';

/*printf("%s\n",file->card);*/
    return(1);
}
    

/**********************************************************************
** clear_token: clear all token elements
**********************************************************************/
static void clear_token(token)
TOKEN *token;
{
    token->class = ILLEGAL;
    token->subclass = ILLEGAL;
    token->int_val = 0;
    token->float_val = 0.0;
    strcpy(token->str_val,"");
    strcpy(token->comment,"");
}

/**********************************************************************
** parse_keyword_card: interpret keyword card and store info in struct 'token'
**********************************************************************/

static int parse_keyword_card(card, token)
char  *card;
TOKEN *token;
{
    char string[CARDLENGTH+1];
    int  i, j;
    double atof();

    clear_token(token);

    /* get keyword */
    for ( i=0; i<MAXKEYWORDCOL  &&  card[i]==' '; i++)
        ;
    if (i == MAXKEYWORDCOL) {   /* if cols 1-8 blank, assume a comment */
        token->class = COMMENT;
    } else {                    /* otherwise get the keyword */
        for ( j=0; i<MAXDATACOL  &&  card[i]!=' ' && card[i]!='=' ; i++, j++)
            string[j] = card[i];
        string[j] = '\0';
    }   
    /* look up class of keyword */
    if (0 > look_up_class(string, token)) 
        return(-1);

    /* process rest of card */ 
    if (token->class == COMMENT) {   /* if a comment; take the whole card */   
        for ( j=0; i<CARDLENGTH ; i++, j++)
            token->str_val[j] = card[i];
        token->str_val[j] = '\0';
    }
    else if (token->class != END) {  /* if not a comment or END card... */
        /* get past '=' */
        for ( ; i<MAXDATACOL  &&  card[i]!='=' ; i++)
            ;
        if (card[i]!='=') {
	    err_msg(Fits2HDF, CorruptedInputFile);
            printf("Error: Missing \'=\' on keyword card:\n%s\n",card);
            return(-1);
        }
    
        /* get value part */
        for (i++ ; i<MAXDATACOL  &&  card[i]==' ' ; i++)
            ;
        for ( j=0; i<MAXDATACOL && card[i]!=' ' && card[i]!='/'; i++, j++)
            token->str_val[j] = card[i];
        token->str_val[j] = '\0';
    
        /* get comment */
        for ( ; i<MAXKEYWORDCOL  &&  card[i]!='/' ; i++)
            ;
        strcpy(token->comment, &card[i]);
    }

    /* assign token value */
    switch (token->class) { 
        case SIMPLE:  break;
        case BITPIX:  token->int_val = atoi(token->str_val); break;
        case NAXIS:   token->int_val = atoi(token->str_val); break;
        case DIMSIZE: token->int_val = atoi(token->str_val); break;
        case OTHER:   break;
        case BSCALE:  atof(token->str_val); 
        case BZERO:   token->float_val = atof(token->str_val); break;
        case DATAMAX: token->float_val = atof(token->str_val); break;
        case DATAMIN: token->float_val = atof(token->str_val); break;
        case BLANK:   token->int_val   = atoi(token->str_val); break;
        case COMMENT: break;
        case END:     break;
        default:
	  err_msg(Fits2HDF, CorruptedInputFile);
	  printf("Error: Invalid keyword card:\n%s\n",card);
	  return(-1);
    }
    return(1);
}


/**********************************************************************
** look_up_class: find class of keyword, storing it in token->class
**********************************************************************/

static int look_up_class(word,token)
char *word;
TOKEN *token;
{
    int i = 0;

    while ( i<NCLASSES && 0!=strcmp(word,classtab[i].word) )
        i++;

    if ( i < NCLASSES)                     /* keyword found in table */
        token->class = classtab[i].value;
     
    else if (0==strncmp(word, "NAXIS", 5)) {  /* keyword is NAXISn      */
        token->class = DIMSIZE;
        token->subclass = atoi(&word[5]);
    } 
    else
        token->class = OTHER;              /* unrecognized token */

    return(token->class);
}


/**********************************************************************
** clear_file_struct: clear struct with fits file info 
**********************************************************************/
static void clear_file_struct(file)
FITSFILE *file;
{
    file->val.simple='\0';
    file->val.bitpix = 0;
    file->val.rank   = 0;
    file->val.dimsizes = (int *) NULL;
    file->opt.issimple   =  FALSE;
    file->opt.isbzero    = 	FALSE;
    file->opt.isbscale   = 	FALSE;
    file->opt.ismax      = 	FALSE;
    file->opt.ismin      = 	FALSE;
    file->opt.isblank    = 	FALSE;
}
       



/**********************************************************************
** print_fits_info: for debugging, print info extracted from file
**********************************************************************/
static void print_fits_info(file)
FITSFILE *file;
{
    int i;

    printf("simple  = %c\n",file->val.simple);
    printf("bitpix  = %d\n",file->val.bitpix);
    printf("rank    = %d\n",file->val.rank);
    printf("dimsizes:\n");

    for (i=0; i<file->val.rank; i++)
        printf("         %d\n",file->val.dimsizes[i]);

    if (file->opt.isblank)
        printf("blank   = %d\n",file->val.blank);
    if (file->opt.isbscale)
        printf("bscale  = %f\n",file->val.bscale);
    if (file->opt.isbzero)
        printf("bzero   = %f\n",file->val.bzero);
    if (file->opt.ismax)
        printf("max     = %f\n",file->val.max);
    if (file->opt.ismin)
        printf("min     = %f\n",file->val.min);
}

ErrorCode FitsToHDF(name, info)
     char *name;
     FileInfo *info;
{
    ErrorCode error;
    unsigned char *image;
    int   i, elementsize, tot_elements;
    float *fpdata;
    TOKEN token;
    FITSFILE ff;
    
    (void) clear_file_struct(&ff); /* init. struct that will hold file info */

    ff.name = name;
    if ( NULL == (ff.fp = fopen(ff.name, "rb")) )
      return(err_msg(Fits2HDF, InputOpenFailed));

    ff.block = (char *) malloc( (unsigned int) (BLOCKLENGTH+CARDLENGTH) );
    ff.blockptr = BLOCKLENGTH;  /* prime blockptr to cause first block read */
    
/* read in required keyword records: SIMPLE, BITPIX, NAXIS, NAXIS1,... */

    if ((error = get_required_keywords(&ff)) != AllOk) return(error);

/* process rest of cards */

    if (0 > next_card(&ff))                     /* get first card image */
      return(err_msg(Fits2HDF, CorruptedInputFile));

    if (0 > parse_keyword_card(ff.card,&token)) /* parse it             */
      return(ConvertFailed);

    while (token.class != END) {                /* loop for all header cards */
        switch (token.class) {
            case BSCALE:  ff.val.bscale = token.float_val; 
                          ff.opt.isbscale = TRUE; 
                          break;
            case BZERO:   ff.val.bzero = token.float_val; 
                          ff.opt.isbzero = TRUE; 
                          break;
            case DATAMAX: ff.val.max = token.float_val; 
                          ff.opt.ismax = TRUE; 
                          break;
            case DATAMIN: ff.val.min = token.float_val; 
                          ff.opt.ismin = TRUE; 
                          break;
            case BLANK:   ff.val.blank = token.int_val; 
                          ff.opt.isblank = TRUE; 
                          break;
            default: break;
        }
        if (0 > next_card(&ff))                /* get next card image */
	  return(err_msg(Fits2HDF, CorruptedInputFile));

        if (0 > parse_keyword_card(ff.card,&token)) /* parse it       */
	  return(ConvertFailed);
    }

    if (Debug) print_fits_info(&ff);       /* for debugging */

/* get element size and compute space needed for elements */

    elementsize = (ff.val.bitpix==16) ? 2 : 4;  
    for (i=0, tot_elements=1; i<ff.val.rank; i++)     /* # of elements */
        tot_elements *= ff.val.dimsizes[i];

/* allocate space for output data and image */

    fpdata = (float *) malloc(tot_elements*sizeof(float));
    image = (unsigned char *) malloc(tot_elements);      

/* read file, and generate the data and the image */

    if ((error = get_data_and_image(&ff, elementsize, tot_elements, fpdata,
				    image)) != AllOk)
      return(error);

    info->image = image;
    info->data = fpdata;
    info->width = ff.val.dimsizes[0];
    info->height = ff.val.dimsizes[1];
    info->dimension = ff.val.rank;
    info->ranks = (long *)malloc(ff.val.rank * sizeof(long));
    for (i=0; i<ff.val.rank; i++)
      *((info->ranks)+i) = *((ff.val.dimsizes)+i);

    fclose(ff.fp);
    free (ff.block);

/* Make a nice gray scale palette */

    for (i=0; i<(MAX_PAL/3); i++)
      info->palette[i*3] = info->palette[i*3+1] =
	info->palette[i*3+2] = (unsigned char) 255 - i;

    info->values |= (BitMask)CvtPalette;

    if (info->image && info->data)
      {
	info->values |= (BitMask)(CvtImage | CvtData);
	return AllOk;
      }
    else
      return (err_msg(Fits2HDF, ConvertFailed));
}