home *** CD-ROM | disk | FTP | other *** search
/ Encyclopedia of Graphics File Formats Companion / GFF_CD.ISO / software / unix / saoimage / sao1_07.tar / readfits.c < prev    next >
C/C++ Source or Header  |  1990-05-02  |  7KB  |  230 lines
  1. #ifndef lint
  2. static char SccsId[] = "%W%  %G%";
  3. #endif
  4.  
  5. /* Module:    readfits.c (Read FITS)
  6.  * Purpose:    Read a FITS file header and prepare file for basic array read
  7.  * Subroutine:    int init_fits()
  8.  * Note:    Only SIMPLE = T is permitted, but will handle BITPIX = 8, 16,
  9.  *        32, -16(u_short), -32(float), -64(double).  If img->nimage
  10.  *        is set to a value greater than 1, file will be lseeked to the
  11.  *        nimage'th image (assuming there are several images stored as
  12.  *        a stack (on 3rd dimension).
  13.  * Copyright:    1989 Smithsonian Astrophysical Observatory
  14.  *        You may do anything you like with this file except remove
  15.  *        this copyright.  The Smithsonian Astrophysical Observatory
  16.  *        makes no representations about the suitability of this
  17.  *        software for any purpose.  It is provided "as is" without
  18.  *        express or implied warranty.
  19.  * Modified:    {0} Michael VanHilst    initial version         23 January 1989
  20.  *        {n} <who> -- <does what> -- <when>
  21.  */
  22.  
  23. #include <stdio.h>        /* define stderr, NULL, etc. */
  24. #include <X11/Xlib.h>        /* needed for control.h */
  25. #include "hfiles/constant.h"
  26. #include "hfiles/control.h"    /* define IOP codes */
  27. #include "hfiles/image.h"
  28.  
  29. #define FITSBLOCK 2880
  30.  
  31. /*
  32.  * Subroutine:    init_fits
  33.  * Purpose:    Open fits file, get parameters and get ready to read data.
  34.  * Returns:    -1 if failure, else fd of file, positioned at start of data
  35.  * Note:    Data read will be handled as standard array read.
  36.  */
  37. int init_fits ( img )
  38.      struct imageRec *img;
  39. {
  40.   float scale, bias;
  41.   int fd;
  42.   int headlen;
  43.   int bitpix, naxis, naxes[3];
  44.   char *header;
  45.   static char *load_fitsheader();
  46.   int read_fitsheader(), open_disk(), lseek_disk();
  47.   void close_disk();
  48.  
  49.   /* open the image file */
  50.   if( (fd = open_disk(img->filename, IOP_Read, 0)) <= 0 )
  51.     return( -1 );
  52.   headlen = img->headersize;
  53.   /* read and interpret the header */
  54.   if( (header = load_fitsheader(fd, &headlen, img->filename)) == NULL ) {
  55.     close_disk(fd, img->filename);
  56.     return( -1 );
  57.   }
  58.   if( read_fitsheader(header, headlen, &bitpix, &naxis, naxes, &scale, &bias)
  59.       != 0 ) {
  60.     free(header);
  61.     switch( bitpix ) {
  62.     case 8:
  63.       img->storage_type = ARR_U1;
  64.       img->bytepix = 1;
  65.       break;
  66.     case 16:
  67.       img->storage_type = ARR_I2;
  68.       img->bytepix = 2;
  69.       break;
  70.     case 32:
  71.       img->storage_type = ARR_I4;
  72.       img->bytepix = 4;
  73.       break;
  74.     case -16:
  75.       img->storage_type = ARR_U2;
  76.       img->bytepix = 2;
  77.       break;
  78.     case -32:
  79.       img->storage_type = ARR_R4;
  80.       img->bytepix = 4;
  81.       break;
  82.     case -64:
  83.       img->storage_type = ARR_R8;
  84.       img->bytepix = 8;
  85.       break;
  86.     default:
  87.       (void)fprintf(stderr,"Illegal FITS BITPIX: %d\n",bitpix);
  88.       close_disk(fd, img->filename);
  89.       return( -1 );
  90.     }
  91.     if( img->nimage > 1 ) {
  92.       if( (naxis <= 2) || (naxes[2] < img->nimage) ) {
  93.     (void)fprintf(stderr,
  94.               "Only %d images in file %s\n", naxis, img->filename);
  95.     close_disk(fd, img->filename);
  96.     return( -1 );
  97.       }
  98.       img->headersize = headlen + (naxes[0] * naxes[1] * img->bytepix);
  99.       if( lseek_disk(fd, img->headersize, img->filename) < 0 ) {
  100.     img->headersize = 0;
  101.     close_disk(fd, img->filename);
  102.     return( -1 );
  103.       }
  104.     } else {
  105.       img->headersize = headlen;
  106.     }
  107.     img->filecols = naxes[0];
  108.     img->filerows = naxes[1];
  109.     if( (scale != 1.0) || (bias != 0.0) ) {
  110.       img->fscale = scale;
  111.       img->fbias = bias;
  112.       img->fscaled = 1;
  113.     } else
  114.       img->fscaled = 0;
  115.     return( fd );
  116.   } else {
  117.     free(header);
  118.     close_disk(fd, img->filename);
  119.     return( -1 );
  120.   }
  121. }
  122.  
  123.  
  124. /*
  125.  * Subroutine:    load_fitsheader
  126.  * Purpose:    Load the FITS header from the disk file into memory.
  127.  * Notes:    FITS headers are made in block units (2880 bytes each).
  128.  *        If length is not initially zero, it is taken as a directive
  129.  *        to override FITS standard (header of given size is loaded).
  130.  *        FITS header must have "SIMPLE =" in first 10 bytes (and 'T'
  131.  *        in next 20 for ximage use).
  132.  */
  133. static char *load_fitsheader ( fd, length, filename )
  134.      int fd;
  135.      int *length;
  136.      char *filename;
  137. {
  138.   int headlen, nbytes;
  139.   char *header;
  140.   char TorF[4];
  141.   static char *get_keyfield();
  142.   char *realloc(), *calloc_errchk();
  143.   int read_disk();
  144.   void no_fitscomment();
  145.  
  146.   if( *length > 0 )
  147.     headlen = *length;
  148.   else
  149.     headlen = FITSBLOCK;
  150.   /* allocate initial sized buffer */
  151.   header = calloc_errchk(headlen, 1, "FITS header");
  152.   /* read in first block */
  153.   if( (nbytes = read_disk(fd, header, headlen, 1, filename, "FITS header"))
  154.      != headlen ) {
  155.     free (header);
  156.     return( NULL );
  157.   }
  158.   /* consistency check for FITS header */
  159.   if( strncmp(header, "SIMPLE  =", 9) != 0) {
  160.     (void)fprintf(stderr,"No SIMPLE keyword in FITS header");
  161.     free(header);
  162.     return( NULL );
  163.   }
  164.   /* we only support SIMPLE = T */
  165.   no_fitscomment(header+10, 20);
  166.   if( (sscanf(header+10, "%1s", TorF) <= 0) || (TorF[0] != 'T') ) {
  167.     (void)fprintf(stderr,"SIMPLE = %c not supported\n", TorF[0]);
  168.     free(header);
  169.     return( NULL );
  170.   }
  171.   if( *length != 0 ) {
  172.     /* if header length was given, it must check out with an END */
  173.     if( get_keyfield(header, "END     ", headlen, 0) == NULL ) {
  174.       (void)fprintf(stderr,"%d byte FITS header has no END card!\n", headlen);
  175.       free(header);
  176.       return( NULL );
  177.     }
  178.   } else {
  179.     /* if END keyword not found, read in more header until END is found */
  180.     while( get_keyfield(header, "END     ", headlen, 0) == NULL ) {
  181.       if( (header = realloc(header, (unsigned)(headlen+FITSBLOCK))) == NULL ) {
  182.     fputs("Reallocation failure reading header\n", stderr);
  183.     return( NULL );
  184.       }
  185.       if( (nbytes = read_disk(fd, header+headlen, FITSBLOCK, 1,
  186.                   filename, "extended header")) != FITSBLOCK ) {
  187.     (void)fprintf(stderr, "Read %d bytes of header with no END card!\n",
  188.               headlen + nbytes);
  189.     (void)fflush(stderr);
  190.     free( header );
  191.     return( NULL );
  192.       }
  193.       headlen += FITSBLOCK;
  194.     }
  195.     *length = headlen;
  196.   }
  197.   return( header );
  198. }
  199.  
  200.  
  201. /*
  202.  * Subroutine:    get_keyfield
  203.  * Purpose:    Return the data field for a given FITS header keyword
  204.  * Returns:    If key not found, return NULL (0).
  205.  */
  206. static char *get_keyfield ( header, keyword, length, report_error )
  207.      char *header;    /* buffer start */
  208.      char *keyword;    /* keyword to match */
  209.      int length;    /* if zero, search up to "END" keyword */
  210.      int report_error;    /* if > 0, fatal if key not found */
  211. {
  212.   int key_not_end, i;
  213.   void no_fitscomment();
  214.  
  215.   key_not_end = (strncmp(keyword, "END     ", 8) != 0);
  216.   for( i=0; i<length; i+=80 ) {
  217.     /* check for END keyword marking end of header, unless END is the key */
  218.     if( key_not_end && strncmp(header+i,"END     ",8) == 0 )
  219.       break;
  220.     /* check for desired keyword */
  221.     if( strncmp(header+i,keyword,8) == 0 ) {
  222.       no_fitscomment (header+i+10, 20);
  223.       return( header+i+10 );
  224.     }
  225.   }
  226.   if( report_error )
  227.     (void)fprintf(stderr, "No `%s' keyword in FITS header\n", keyword);
  228.   return( (char *)NULL );
  229. }
  230.