home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 4 / DATAFILE_PDCD4.iso / utilities / utilsm / netpbmsca / pgm / c / fitstopgm < prev    next >
Encoding:
Text File  |  1993-10-04  |  5.9 KB  |  243 lines
  1. /* fitstopgm.c - read a FITS file and produce a portable graymap
  2. **
  3. ** Copyright (C) 1989 by Jef Poskanzer.
  4. **
  5. ** Permission to use, copy, modify, and distribute this software and its
  6. ** documentation for any purpose and without fee is hereby granted, provided
  7. ** that the above copyright notice appear in all copies and that both that
  8. ** copyright notice and this permission notice appear in supporting
  9. ** documentation.  This software is provided "as is" without express or
  10. ** implied warranty.
  11. */
  12.  
  13. #include "pgm.h"
  14.  
  15. struct FITS_Header {
  16.     int simple;        /* basic format or not */
  17.     int bitpix;        /* number of bits per pixel */
  18.     int naxis;        /* number of axes */
  19.     int naxis1;        /* number of points on axis 1 */
  20.     int naxis2;        /* number of points on axis 2 */
  21.     int naxis3;        /* number of points on axis 3 */
  22.     double datamin;    /* min # */
  23.     double datamax;    /* max # */
  24.     double bzer;    /* ??? */
  25.     double bscale;    /* ??? */
  26.     };
  27.  
  28. static void read_fits_header ARGS(( FILE* fp, struct FITS_Header* hP ));
  29. static void read_card ARGS(( FILE* fp, char* buf ));
  30.  
  31. int
  32. main( argc, argv )
  33. int argc;
  34. char* argv[];
  35.     {
  36.     FILE* ifp;
  37.     gray* grayrow;
  38.     register gray* gP;
  39.     int argn, imagenum, image, row;
  40.     register int col;
  41.     gray maxval;
  42.     double fmaxval, scale;
  43.     int rows, cols, images;
  44.     struct FITS_Header h;
  45.     char* usage = "[-image N] [FITSfile]";
  46.  
  47.  
  48.     pgm_init( &argc, argv );
  49.  
  50.     argn = 1;
  51.     imagenum = 1;
  52.  
  53.     if ( argn < argc && argv[argn][0] == '-' && argv[argn][1] != '\0' )
  54.     {
  55.     if ( pm_keymatch( argv[argn], "-image", 2 ) )
  56.         {
  57.         ++argn;
  58.         if ( argn == argc || sscanf( argv[argn], "%d", &imagenum ) != 1 )
  59.         pm_usage( usage );
  60.         }
  61.     else
  62.         pm_usage( usage );
  63.     ++argn;
  64.     }
  65.  
  66.     if ( argn < argc )
  67.     {
  68.     ifp = pm_openr( argv[argn] );
  69.     ++argn;
  70.     }
  71.     else
  72.     ifp = stdin;
  73.  
  74.     if ( argn != argc )
  75.     pm_usage( usage );
  76.  
  77.     read_fits_header( ifp, &h );
  78.  
  79.     if ( ! h.simple )
  80.     pm_error( "FITS file is not in simple format, can't read" );
  81.     switch ( h.bitpix )
  82.     {
  83.     case 8:
  84.     fmaxval = 255.0;
  85.     break;
  86.  
  87.     case 16:
  88.     fmaxval = 65535.0;
  89.     break;
  90.  
  91.     case 32:
  92.     fmaxval = 4294967295.0;
  93.     break;
  94.  
  95.     default:
  96.     pm_error( "unusual bits per pixel (%d), can't read", h.bitpix );
  97.     }
  98.     if ( h.naxis != 2 && h.naxis != 3 )
  99.     pm_error( "FITS file has %d axes, can't read", h.naxis );
  100.     cols = h.naxis1;
  101.     rows = h.naxis2;
  102.     if ( h.naxis == 2 )
  103.     images = 1;
  104.     else
  105.     images = h.naxis3;
  106.     if ( imagenum > images )
  107.     pm_error( "only %d image%s in this file",
  108.           images, images > 1 ? "s" : "" );
  109.     maxval = min( fmaxval, PGM_MAXMAXVAL );
  110.     if ( h.datamin != 0.0 || h.datamax != 1.0 )
  111.     scale = maxval / ( h.datamax - h.datamin );
  112.     else
  113.     {
  114.     /* FITS images are not required to contain DATAMIN and DATAMAX cards
  115.     ** In this case, it would be necessary to read through the image
  116.     ** twice to properly scale.  Take this shortcut instead, remembering
  117.     ** that all FITS integers are signed values.  This scheme works
  118.     ** on most images because most astronomical data reduction packages
  119.     ** scale the images when writing so as to make maximal use of the
  120.     ** dynamic range of the output format. */
  121.     h.bscale = 1.;
  122.     h.bzer = ( fmaxval + 1.0 ) * 0.5;
  123.     scale = maxval / fmaxval;
  124.     }
  125.  
  126.     pgm_writepgminit( stdout, cols, rows, maxval, 0 );
  127.     grayrow = pgm_allocrow( cols );
  128.     for ( image = 1; image <= imagenum; ++image )
  129.     {
  130.     if ( image != imagenum )
  131.         pm_message( "skipping image %d of %d", image, images );
  132.     else if ( images > 1 )
  133.         pm_message( "reading image %d of %d", image, images );
  134.     for ( row = 0; row < rows; ++row)
  135.         {
  136.         for ( col = 0, gP = grayrow; col < cols; ++col, ++gP )
  137.         {
  138.         int ich;
  139.         double val;
  140.  
  141.         switch ( h.bitpix )
  142.             {
  143.             case 8:
  144.             ich = getc( ifp );
  145.             if ( ich == EOF )
  146.             pm_error( "EOF / read error" );
  147.             val = ich;
  148.             break;
  149.  
  150.             case 16:
  151.             ich = getc( ifp );
  152.             if ( ich == EOF )
  153.             pm_error( "EOF / read error" );
  154.             val = ich << 8;
  155.             ich = getc( ifp );
  156.             if ( ich == EOF )
  157.             pm_error( "EOF / read error" );
  158.             val += ich;
  159.             break;
  160.  
  161.             case 32:
  162.             ich = getc( ifp );
  163.             if ( ich == EOF )
  164.             pm_error( "EOF / read error" );
  165.             val = ich << 24;
  166.             ich = getc( ifp );
  167.             if ( ich == EOF )
  168.             pm_error( "EOF / read error" );
  169.             val += ich << 16;
  170.             ich = getc( ifp );
  171.             if ( ich == EOF )
  172.             pm_error( "EOF / read error" );
  173.             val += ich << 8;
  174.             ich = getc( ifp );
  175.             if ( ich == EOF )
  176.             pm_error( "EOF / read error" );
  177.             val += ich;
  178.             break;
  179.  
  180.             default:
  181.             pm_error( "can't happen" );
  182.             }
  183.         *gP = (gray) ( scale * ( val * h.bscale + h.bzer - h.datamin) );
  184.         }
  185.         if ( image == imagenum )
  186.         pgm_writepgmrow( stdout, grayrow, cols, maxval, 0 );
  187.         }
  188.     }
  189.     pm_close( ifp );
  190.     pm_close( stdout );
  191.  
  192.     exit( 0 );
  193.     }
  194.  
  195. static void
  196. read_fits_header( fp, hP )
  197. FILE* fp;
  198. struct FITS_Header* hP;
  199.     {
  200.     char buf[80];
  201.     int seen_end;
  202.     int i;
  203.     char c;
  204.  
  205.     seen_end = 0;
  206.     hP->simple = 0;
  207.     hP->bzer = 0.0;
  208.     hP->bscale = 1.0;
  209.     hP->datamin = 0.0;
  210.     hP->datamax = 1.0;
  211.  
  212.     while ( ! seen_end )
  213.     for ( i = 0; i < 36; ++i )
  214.         {
  215.         read_card( fp, buf );
  216.  
  217.         if ( sscanf( buf, "SIMPLE = %c", &c ) == 1 )
  218.         {
  219.         if ( c == 'T' || c == 't' )
  220.             hP->simple = 1;
  221.         }
  222.         else if ( sscanf( buf, "BITPIX = %d", &(hP->bitpix) ) == 1 );
  223.         else if ( sscanf( buf, "NAXIS = %d", &(hP->naxis) ) == 1 );
  224.         else if ( sscanf( buf, "NAXIS1 = %d", &(hP->naxis1) ) == 1 );
  225.         else if ( sscanf( buf, "NAXIS2 = %d", &(hP->naxis2) ) == 1 );
  226.         else if ( sscanf( buf, "NAXIS3 = %d", &(hP->naxis3) ) == 1 );
  227.         else if ( sscanf( buf, "DATAMIN = %lf", &(hP->datamin) ) == 1 );
  228.         else if ( sscanf( buf, "DATAMAX = %lf", &(hP->datamax) ) == 1 );
  229.         else if ( sscanf( buf, "BZERO = %lf", &(hP->bzer) ) == 1 );
  230.         else if ( sscanf( buf, "BSCALE = %lf", &(hP->bscale) ) == 1 );
  231.         else if ( strncmp( buf, "END ", 4 ) == 0 ) seen_end = 1;
  232.         }
  233.     }
  234.  
  235. static void
  236. read_card( fp, buf )
  237. FILE* fp;
  238. char* buf;
  239.     {
  240.     if ( fread( buf, 1, 80, fp ) == 0 )
  241.     pm_error( "error reading header" );
  242.     }
  243.