home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / netpbma.zip / pnm / fitstopnm.c < prev    next >
C/C++ Source or Header  |  1993-12-21  |  12KB  |  426 lines
  1.  /* fitstopnm.c - read a FITS file and produce a portable anymap
  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.  ** Hacked up version by Daniel Briggs  (dbriggs@nrao.edu)  20-Oct-92
  13.  **
  14.  ** Include floating point formats, more or less.  Will only work on
  15.  ** machines that understand IEEE-754.  Added -scanmax -printmax
  16.  ** -min -max and -noraw.  Ignore axes past 3, instead of error (many packages
  17.  ** use pseudo axes).  Use a finite scale when max=min.  NB: Min and max
  18.  ** are the real world FITS values (scaled), so watch out when bzer & bscale
  19.  ** are not 0 & 1.  Datamin & datamax interpreted correctly in scaled case,
  20.  ** and initialization changed to less likely values.  If datamin & max are
  21.  ** not present in the header, the a first pass is made to determine them
  22.  ** from the array values.
  23.  **
  24.  ** Modified by Alberto Accomazzi (alberto@cfa.harvard.edu), Dec 1, 1992.
  25.  **
  26.  ** Added support for 3-planes FITS files, the program is renamed fitstopnm.
  27.  ** Fixed some non-ansi declarations (DBL_MAX and FLT_MAX replace MAXDOUBLE
  28.  ** and MAXFLOAT), fixed some scaling parameters to map the full FITS data
  29.  ** resolution to the maximum PNM resolution, disabled min max scanning when
  30.  ** reading from stdin.
  31.  */
  32.  
  33. #include "pnm.h"
  34.  
  35. /* Do what you have to, to get a sensible value here.  This may not be so */
  36. /* portable as the rest of the program.  We want to use MAXFLOAT and */
  37. /* MAXDOUBLE, so you could define them manually if you have to. */
  38. /* #include <values.h> */
  39. #include <float.h>
  40.  
  41. static struct FITS_Header {
  42.   int simple;        /* basic format or not */
  43.   int bitpix;        /* number of bits per pixel */
  44.   int naxis;        /* number of axes */
  45.   int naxis1;        /* number of points on axis 1 */
  46.   int naxis2;        /* number of points on axis 2 */
  47.   int naxis3;        /* number of points on axis 3 */
  48.   double datamin;    /* min # (Physical value!) */
  49.   double datamax;    /* max #     "       "     */
  50.   double bzer;        /* Physical value = Array value*bscale + bzero */
  51.   double bscale;
  52. };
  53.  
  54. static void read_fits_header ARGS(( FILE* fp, struct FITS_Header* hP ));
  55. static void read_card ARGS(( FILE* fp, char* buf ));
  56. static void read_val ARGS(( FILE* fp, int bitpix, double* vp ));
  57.      
  58. int
  59. main( argc, argv )
  60. int argc;
  61. char* argv[];
  62. {
  63.   FILE* ifp;
  64.   register gray* gP;
  65.   int argn, imagenum, image, row;
  66.   register int col;
  67.   xelval maxval;
  68.   double val, fmaxval, frmin, frmax, dmax, dmin, scale, t;
  69.   int rows, cols, images, format;
  70.   struct FITS_Header h;
  71.   xel** pnmarray;
  72.   xelval tx, txv[4];
  73.   char* fits_name;
  74.   char* usage = "[-image N] [-noraw] [-scanmax] [-printmax] [-min f] [-max f] [FITSfile]";
  75.  
  76.   int doscan = 0;
  77.   int forceplain = 0;
  78.   int forcemin = 0;
  79.   int forcemax = 0;
  80.   int printmax = 0;
  81.   
  82.   pnm_init( &argc, argv );
  83.   
  84.   argn = 1;
  85.   imagenum = 0;
  86.   
  87.   while ( argn < argc && argv[argn][0] == '-' && argv[argn][1] != '\0' )
  88.     {
  89.       if ( pm_keymatch( argv[argn], "-image", 2 ) )
  90.     {
  91.       ++argn;
  92.       if ( argn == argc || sscanf( argv[argn], "%d", &imagenum ) != 1 )
  93.         pm_usage( usage );
  94.     }
  95.       else if ( pm_keymatch( argv[argn], "-max", 3 ) )
  96.     {
  97.       ++argn;
  98.       forcemax = 1;
  99.       if ( argn == argc || sscanf( argv[argn], "%lf", &frmax ) != 1 )
  100.         pm_usage( usage );
  101.     }
  102.       else if ( pm_keymatch( argv[argn], "-min", 3 ) )
  103.     {
  104.       ++argn;
  105.       forcemin = 1;
  106.       if ( argn == argc || sscanf( argv[argn], "%lf", &frmin ) != 1 )
  107.         pm_usage( usage );
  108.     }
  109.       else if ( pm_keymatch( argv[argn], "-scanmax", 2 ) )
  110.     doscan = 1;
  111.       else if ( pm_keymatch( argv[argn], "-noraw", 2 ) )
  112.     forceplain = 1;
  113.       else if ( pm_keymatch( argv[argn], "-printmax", 2 ) )
  114.     printmax = 1;
  115.       else
  116.     pm_usage( usage );
  117.       ++argn;
  118.     }
  119.   
  120.   if ( argn < argc )
  121.     {
  122.       fits_name = argv[argn];
  123.       ifp = pm_openr( fits_name );
  124.       ++argn;
  125.     }
  126.   else
  127.     ifp = stdin;
  128.  
  129.   if ( argn != argc )
  130.     pm_usage( usage );
  131.   
  132.   read_fits_header( ifp, &h );
  133.   
  134.   if ( ! h.simple )
  135.     pm_error( "FITS file is not in simple format, can't read" );
  136.   switch ( h.bitpix )
  137.     {
  138.     case 8:
  139.       fmaxval = 255.0;
  140.       break;
  141.       
  142.     case 16:
  143.       fmaxval = 65535.0;
  144.       break;
  145.       
  146.     case 32:
  147.       fmaxval = 4294967295.0;
  148.       break;
  149.       
  150.     case -32:
  151.       fmaxval = FLT_MAX;
  152.       break;
  153.       
  154.     case -64:
  155.       fmaxval = DBL_MAX;
  156.       break;
  157.       
  158.     default:
  159.       pm_error( "unusual bits per pixel (%d), can't read", h.bitpix );
  160.     }
  161.  
  162.   if ( h.naxis != 2 && h.naxis != 3 )
  163.     pm_message( "Warning: FITS file has %d axes", h.naxis );
  164.   cols = h.naxis1;
  165.   rows = h.naxis2;
  166.   if ( h.naxis == 2 )
  167.     images = imagenum = 1;
  168.   else
  169.     images = h.naxis3;
  170.   if ( imagenum > images )
  171.     pm_error( "only %d plane%s in this file", images, images > 1 ? "s" : "" );
  172.   if ( images != 3 && images > 1 && imagenum == 0 )
  173.     pm_error( "need to specify a plane using the -imagenum flag" );
  174.  
  175.   format = PGM_FORMAT;
  176.   if ( images == 3 && imagenum == 0 ) 
  177.     format = PPM_FORMAT;
  178.  
  179.   if (forcemin) h.datamin = frmin;
  180.   if (forcemax) h.datamax = frmax;
  181.  
  182.   if ( doscan || h.datamin == - DBL_MAX || h.datamax == DBL_MAX )
  183.     {
  184.     if ( ifp == stdin )
  185.       pm_message( "warning: cannot scan max and read from stdin" );
  186.     else 
  187.       {
  188.       /* OK, We have to scale this the hard way. */
  189.       pm_message( "Scanning file for scaling parameters" );
  190.       dmax = -fmaxval;
  191.       dmin = fmaxval;
  192.       for ( image = 1; image <= images; ++image )
  193.     {
  194.       for ( row = 0; row < rows; ++row)
  195.         {
  196.           for ( col = 0; col < cols; ++col)
  197.         {
  198.           read_val (ifp, h.bitpix, &val);
  199.           if ( image == imagenum || format == PPM_FORMAT ) {
  200.             if ( val > dmax ) dmax = val;
  201.             if ( val < dmin ) dmin = val;
  202.           }
  203.             }
  204.         }
  205.       pm_close( ifp );
  206.       if (h.bscale < 0.0) {
  207.         t = dmax;
  208.         dmax = dmin;
  209.         dmin = t;
  210.       }
  211.       ifp = pm_openr( fits_name );
  212.       read_fits_header( ifp, &h );
  213.       h.datamin = forcemin ? frmin : dmin * h.bscale + h.bzer;
  214.       h.datamax = forcemax ? frmax : dmax * h.bscale + h.bzer;
  215.     }
  216.       }
  217.     }
  218.  
  219.  
  220.   maxval = min( h.datamax - h.datamin, PNM_MAXMAXVAL );
  221.   scale = 1.0;
  222.   /* Just in case we have a file with constant value. */
  223.   if ( h.datamin != h.datamax )
  224.     scale = maxval / ( h.datamax - h.datamin );
  225.  
  226.   /* If printmax option is true, just print and exit. */
  227.   /* For use in shellscripts.  Ex:                    */
  228.   /* eval `fitstopnm -printmax $filename | \          */
  229.   /*   awk '{min = $1; max = $2}\                     */
  230.   /*         END {print "min=" min; " max=" max}'`    */
  231.   if (printmax) {
  232.     printf( "%lf %lf\n", h.datamin, h.datamax);
  233.     pm_close( ifp );
  234.     pm_close( stdout );
  235.     exit( 0 );
  236.     }
  237.  
  238.   pnmarray = pnm_allocarray( cols, rows );
  239.  
  240.   switch( PNM_FORMAT_TYPE( format ))
  241.     {
  242.     case PGM_TYPE:
  243.       pm_message( "writing PGM file" );
  244.       for ( image = 1; image <= imagenum; ++image )
  245.         {
  246.         if ( image != imagenum )
  247.           pm_message( "skipping image plane %d of %d", image, images );
  248.         else if ( images > 1 )
  249.           pm_message( "reading image plane %d of %d", image, images );
  250.         for ( row = 0; row < rows; ++row)
  251.           for ( col = 0; col < cols; ++col )
  252.             {
  253.             read_val (ifp, h.bitpix, &val);
  254.             t = scale * ( val * h.bscale + h.bzer - h.datamin );
  255.         tx = max( 0, min( t, maxval ) );
  256.         if ( image == imagenum )
  257.               PNM_ASSIGN1( pnmarray[row][col], tx );
  258.             }
  259.     }
  260.       break;
  261.     case PPM_TYPE:
  262.       pm_message( "writing PPM file" );
  263.       for ( image = 1; image <= images; image++ )
  264.     {
  265.     pm_message( "reading image plane %d of %d", image, images );
  266.     for ( row = 0; row < rows; row++ )
  267.       for ( col = 0; col < cols; col++ )
  268.             {
  269.             read_val (ifp, h.bitpix, &val);
  270.         txv[1] = PPM_GETR( pnmarray[row][col] );
  271.         txv[2] = PPM_GETG( pnmarray[row][col] );
  272.         txv[3] = PPM_GETB( pnmarray[row][col] );
  273.             t = scale * ( val * h.bscale + h.bzer - h.datamin );
  274.         txv[image] = max( 0, min( t, maxval ));
  275.             PPM_ASSIGN( pnmarray[row][col], txv[1], txv[2], txv[3] );
  276.             }
  277.     }
  278.       break;
  279.     default:
  280.       pm_error( "can't happen" );
  281.       break;
  282.     }
  283.  
  284.   pnm_writepnm( stdout, pnmarray, cols, rows, maxval, format, forceplain );
  285.   pnm_freearray( pnmarray, rows );
  286.  
  287.   pm_close( ifp );
  288.   pm_close( stdout );
  289.   
  290.   exit( 0 );
  291. }
  292.  
  293. /*
  294.  ** This code will deal properly with integers, no matter what the byte order
  295.  ** or integer size of the host machine.  Sign extension is handled manually
  296.  ** to prevent problems with signed/unsigned characters.  Floating point
  297.  ** values will only be read properly when the host architecture is IEEE-754
  298.  ** conformant.  If you need to tweak this code for other machines, you might
  299.  ** want to snag a copy of the FITS documentation from nssdca.gsfc.nasa.gov
  300.  */
  301.  
  302. static void
  303.   read_val (fp, bitpix, vp)
  304. FILE *fp;
  305. int bitpix;
  306. double *vp;
  307.  
  308. {
  309.   int i, ich, ival;
  310.   long int lval;
  311.   unsigned char c[8];
  312.   
  313.   switch ( bitpix )
  314.     {
  315.       /* 8 bit FITS integers are unsigned */
  316.     case 8:
  317.       ich = getc( fp );
  318.       if ( ich == EOF )
  319.     pm_error( "EOF / read error" );
  320.       *vp = ich;
  321.       break;
  322.       
  323.     case 16:
  324.       ich = getc( fp );
  325.       if ( ich == EOF )
  326.     pm_error( "EOF / read error" );
  327.       c[0] = ich;
  328.       ich = getc( fp );
  329.       if ( ich == EOF )
  330.     pm_error( "EOF / read error" );
  331.       c[1] = ich;
  332.       if (c[0] & 0x80)
  333.     ival = ~0xFFFF | c[0]<<8 | c[1];
  334.       else
  335.     ival = c[0]<<8 | c[1];
  336.       *vp = ival;
  337.       break;
  338.       
  339.     case 32:
  340.       for (i=0; i<4; i++) {
  341.     ich = getc( fp );
  342.     if ( ich == EOF )
  343.       pm_error( "EOF / read error" );
  344.     c[i] = ich;
  345.       }
  346.       if (c[0] & 0x80)
  347.     lval = ~0xFFFFFFFF |
  348.       c[0]<<24 | c[1]<<16 | c[2]<<8 | c[3];
  349.       else
  350.     lval = c[0]<<24 | c[1]<<16 | c[2]<<8 | c[3];
  351.       *vp = lval;
  352.       
  353.     case -32:
  354.       for (i=0; i<4; i++) {
  355.     ich = getc( fp );
  356.     if ( ich == EOF )
  357.       pm_error( "EOF / read error" );
  358.     c[i] = ich;
  359.       }
  360.       *vp = *( (float *) c);
  361.       break;
  362.       
  363.     case -64:
  364.       for (i=0; i<8; i++) {
  365.     ich = getc( fp );
  366.     if ( ich == EOF )
  367.       pm_error( "EOF / read error" );
  368.     c[i] = ich;
  369.       }
  370.       *vp = *( (double *) c);
  371.       break;
  372.       
  373.     default:
  374.       pm_error( "Strange bitpix in read_value" );
  375.     }
  376. }
  377.  
  378. static void
  379.   read_fits_header( fp, hP )
  380. FILE* fp;
  381. struct FITS_Header* hP;
  382. {
  383.   char buf[80];
  384.   int seen_end;
  385.   int i;
  386.   char c;
  387.   
  388.   seen_end = 0;
  389.   hP->simple = 0;
  390.   hP->bzer = 0.0;
  391.   hP->bscale = 1.0;
  392.   hP->datamin = - DBL_MAX;
  393.   hP->datamax = DBL_MAX;
  394.   
  395.   while ( ! seen_end )
  396.     for ( i = 0; i < 36; ++i )
  397.       {
  398.     read_card( fp, buf );
  399.     
  400.     if ( sscanf( buf, "SIMPLE = %c", &c ) == 1 )
  401.       {
  402.         if ( c == 'T' || c == 't' )
  403.           hP->simple = 1;
  404.       }
  405.     else if ( sscanf( buf, "BITPIX = %d", &(hP->bitpix) ) == 1 );
  406.     else if ( sscanf( buf, "NAXIS = %d", &(hP->naxis) ) == 1 );
  407.     else if ( sscanf( buf, "NAXIS1 = %d", &(hP->naxis1) ) == 1 );
  408.     else if ( sscanf( buf, "NAXIS2 = %d", &(hP->naxis2) ) == 1 );
  409.     else if ( sscanf( buf, "NAXIS3 = %d", &(hP->naxis3) ) == 1 );
  410.     else if ( sscanf( buf, "DATAMIN = %lf", &(hP->datamin) ) == 1 );
  411.     else if ( sscanf( buf, "DATAMAX = %lf", &(hP->datamax) ) == 1 );
  412.     else if ( sscanf( buf, "BZERO = %lf", &(hP->bzer) ) == 1 );
  413.     else if ( sscanf( buf, "BSCALE = %lf", &(hP->bscale) ) == 1 );
  414.     else if ( strncmp( buf, "END ", 4 ) == 0 ) seen_end = 1;
  415.       }
  416. }
  417.  
  418. static void
  419.   read_card( fp, buf )
  420. FILE* fp;
  421. char* buf;
  422. {
  423.   if ( fread( buf, 1, 80, fp ) == 0 )
  424.     pm_error( "error reading header" );
  425. }
  426.