Rebase advanced branch to 10.39

git-svn-id: http://svn.code.sf.net/p/netpbm/code/advanced@330 9d0c8265-081b-0410-96cb-a4ca84ce46f8

1 files changed, 497 insertions, 309 deletions

diff --git a/converter/other/fitstopnm.c b/converter/other/fitstopnm.c
index b143882b..86bdfbb8 100644
--- a/converter/other/fitstopnm.c
+++ b/converter/other/fitstopnm.c
@@ -1,4 +1,4 @@
- /* fitstopnm.c - read a FITS file and produce a portable anymap
+ /* fitstopnm.c - read a FITS file and produce a PNM.
  **
  ** Copyright (C) 1989 by Jef Poskanzer.
  **
@@ -31,15 +31,109 @@
  */
 
 #include <string.h>
+#include <float.h>
 
 #include "pm_c_util.h"
+#include "mallocvar.h"
+#include "shhopt.h"
 #include "pnm.h"
 
-/* Do what you have to, to get a sensible value here.  This may not be so */
-/* portable as the rest of the program.  We want to use MAXFLOAT and */
-/* MAXDOUBLE, so you could define them manually if you have to. */
-/* #include <values.h> */
-#include <float.h>
+
+
+struct cmdlineInfo {
+    const char * inputFileName;
+    unsigned int image;  /* zero if unspecified */
+    double max;
+    unsigned int maxSpec;
+    double min;
+    unsigned int minSpec;
+    unsigned int scanmax;
+    unsigned int printmax;
+    unsigned int noraw;
+        /* This is for backward compatibility only.  Use the common option
+           -plain now.  (pnm_init() processes -plain).
+        */
+    unsigned int verbose;
+    unsigned int omaxval;
+    unsigned int omaxvalSpec;
+};
+
+
+
+static void 
+parseCommandLine(int argc, char ** argv, 
+                 struct cmdlineInfo * const cmdlineP) {
+/* --------------------------------------------------------------------------
+   Parse program command line described in Unix standard form by argc
+   and argv.  Return the information in the options as *cmdlineP.  
+
+   If command line is internally inconsistent (invalid options, etc.),
+   issue error message to stderr and abort program.
+
+   Note that the strings we return are stored in the storage that
+   was passed to us as the argv array.  We also trash *argv.
+--------------------------------------------------------------------------*/
+    optEntry * option_def;
+        /* Instructions to optParseOptions3 on how to parse our options. */
+    optStruct3 opt;
+
+    unsigned int imageSpec;
+    unsigned int option_def_index;
+
+    MALLOCARRAY_NOFAIL(option_def, 100);
+
+    option_def_index = 0;   /* incremented by OPTENT3 */
+    OPTENT3(0, "image",    OPT_UINT,
+            &cmdlineP->image,   &imageSpec,                            0);
+    OPTENT3(0, "min",      OPT_FLOAT,
+            &cmdlineP->min,     &cmdlineP->minSpec,                    0);
+    OPTENT3(0, "max",      OPT_FLOAT,
+            &cmdlineP->max,     &cmdlineP->maxSpec,                    0);
+    OPTENT3(0, "scanmax",  OPT_FLAG,
+            NULL,               &cmdlineP->scanmax,                    0);
+    OPTENT3(0, "printmax", OPT_FLAG,
+            NULL,               &cmdlineP->printmax,                   0);
+    OPTENT3(0, "noraw",    OPT_FLAG,
+            NULL,               &cmdlineP->noraw,                      0);
+    OPTENT3(0, "verbose",  OPT_FLAG,
+            NULL,               &cmdlineP->verbose,                    0);
+    OPTENT3(0, "omaxval",  OPT_UINT,
+            &cmdlineP->omaxval, &cmdlineP->omaxvalSpec,                0);
+
+    opt.opt_table = option_def;
+    opt.short_allowed = FALSE;  /* We have no short (old-fashioned) options */
+    opt.allowNegNum = FALSE;   /* We have no parms that are negative numbers */
+
+    /* Set some defaults the lazy way (using multiple setting of variables) */
+
+    optParseOptions3(&argc, argv, opt, sizeof(opt), 0);
+        /* Uses and sets argc, argv, and some of *cmdlineP and others. */
+
+    if (imageSpec) {
+        if (cmdlineP->image == 0)
+            pm_error("You may not specify zero for the image number.  "
+                     "Images are numbered starting at 1.");
+    } else
+        cmdlineP->image = 0;
+
+    if (cmdlineP->maxSpec && cmdlineP->minSpec) {
+        if (cmdlineP->max <= cmdlineP->min)
+            pm_error("-max must be greater than -min.  You specified "
+                     "-max=%f, -min=%f", cmdlineP->max, cmdlineP->min);
+    }
+
+    if (argc-1 < 1)
+        cmdlineP->inputFileName = "-";
+    else {
+        cmdlineP->inputFileName = argv[1];
+        
+        if (argc-1 > 1)
+            pm_error("Too many arguments (%u).  The only non-option argument "
+                     "is the input file name.", argc-1);
+    }
+}
+
+
 
 struct FITS_Header {
   int simple;       /* basic format or not */
@@ -54,10 +148,205 @@ struct FITS_Header {
   double bscale;
 };
 
-static void read_fits_header ARGS(( FILE* fp, struct FITS_Header* hP ));
-static void read_card ARGS(( FILE* fp, char* buf ));
-static void read_val ARGS(( FILE* fp, int bitpix, double* vp ));
-     
+
+
+/*
+ ** This code will deal properly with integers, no matter what the byte order
+ ** or integer size of the host machine.  Sign extension is handled manually
+ ** to prevent problems with signed/unsigned characters.  Floating point
+ ** values will only be read properly when the host architecture is IEEE-754
+ ** conformant.  If you need to tweak this code for other machines, you might
+ ** want to snag a copy of the FITS documentation from nssdca.gsfc.nasa.gov
+ */
+
+static void
+readVal(FILE *   const ifP,
+        int      const bitpix,
+        double * const vP) {
+
+    switch (bitpix) {
+        /* 8 bit FITS integers are unsigned */
+    case 8: {
+        int const ich = getc(ifP);
+        if (ich == EOF)
+            pm_error("EOF / read error");
+        *vP = ich;
+    } break;
+
+    case 16: {
+        int ich;
+        int ival;
+        unsigned char c[8];
+
+        ich = getc(ifP);
+        if (ich == EOF)
+            pm_error("EOF / read error");
+        c[0] = ich;
+        ich = getc(ifP);
+        if (ich == EOF)
+            pm_error("EOF / read error");
+        c[1] = ich;
+        if (c[0] & 0x80)
+            ival = ~0xFFFF | c[0] << 8 | c[1];
+        else
+            ival = c[0] << 8 | c[1];
+        *vP = ival;
+    } break;
+      
+    case 32: {
+        unsigned int i;
+        long int lval;
+        unsigned char c[4];
+
+        for (i = 0; i < 4; ++i) {
+            int const ich = getc(ifP);
+            if (ich == EOF)
+                pm_error("EOF / read error");
+            c[i] = ich;
+        }
+        if (c[0] & 0x80)
+            lval = ~0xFFFFFFFF | c[0] << 24 | c[1] << 16 | c[2] << 8 | c[3];
+        else
+            lval = c[0] << 24 | c[1] << 16 | c[2] << 8 | c[3] << 0;
+        *vP = lval;
+    } break;
+      
+    case -32: {
+        unsigned int i;
+        unsigned char c[4];
+
+        for (i = 0; i < 4; ++i) {
+            int const ich = getc(ifP);
+            if (ich == EOF)
+                pm_error("EOF / read error");
+            c[i] = ich;
+        }
+        *vP = *((float *)c);
+    } break;
+      
+    case -64: {
+        unsigned int i;
+        unsigned char c[8];
+
+        for (i = 0; i < 8; ++i) {
+            int const ich = getc(ifP);
+            if (ich == EOF)
+                pm_error("EOF / read error");
+            c[i] = ich;
+        }
+        *vP = *((double *)c);
+    } break;
+      
+    default:
+        pm_error("Strange bitpix value %d in readVal()", bitpix);
+    }
+}
+
+
+
+static void
+readCard(FILE * const ifP,
+         char * const buf) {
+
+    size_t bytesRead;
+
+    bytesRead = fread(buf, 1, 80, ifP);
+    if (bytesRead == 0)
+        pm_error("error reading header");
+}
+
+
+
+static void
+readFitsHeader(FILE *               const ifP,
+               struct FITS_Header * const hP) {
+
+    int seenEnd;
+  
+    seenEnd = 0;
+    /* Set defaults */
+    hP->simple  = 0;
+    hP->bzer    = 0.0;
+    hP->bscale  = 1.0;
+    hP->datamin = - DBL_MAX;
+    hP->datamax = DBL_MAX;
+  
+    while (!seenEnd) {
+        unsigned int i;
+        for (i = 0; i < 36; ++i) {
+            char buf[80];
+            char c;
+
+            readCard(ifP, buf);
+    
+            if (sscanf(buf, "SIMPLE = %c", &c) == 1) {
+                if (c == 'T' || c == 't')
+                    hP->simple = 1;
+            } else if (sscanf(buf, "BITPIX = %d", &(hP->bitpix)) == 1);
+            else if (sscanf(buf, "NAXIS = %d", &(hP->naxis)) == 1);
+            else if (sscanf(buf, "NAXIS1 = %d", &(hP->naxis1)) == 1);
+            else if (sscanf(buf, "NAXIS2 = %d", &(hP->naxis2)) == 1);
+            else if (sscanf(buf, "NAXIS3 = %d", &(hP->naxis3)) == 1);
+            else if (sscanf(buf, "DATAMIN = %lf", &(hP->datamin)) == 1);
+            else if (sscanf(buf, "DATAMAX = %lf", &(hP->datamax)) == 1);
+            else if (sscanf(buf, "BZERO = %lf", &(hP->bzer)) == 1);
+            else if (sscanf(buf, "BSCALE = %lf", &(hP->bscale)) == 1);
+            else if (strncmp(buf, "END ", 4 ) == 0) seenEnd = 1;
+        }
+    }
+}
+
+
+
+static void
+interpretPlanes(struct FITS_Header const fitsHeader,
+                unsigned int       const imageRequest,
+                bool               const verbose,
+                unsigned int *     const imageCountP,
+                bool *             const multiplaneP,
+                unsigned int *     const desiredImageP) {
+
+    if (fitsHeader.naxis == 2) {
+        *imageCountP   = 1;
+        *multiplaneP   = FALSE;
+        *desiredImageP = 1;
+    } else {
+        if (imageRequest) {
+            if (imageRequest > fitsHeader.naxis3)
+                pm_error("Only %u plane%s in this file.  "
+                         "You requested image %u", 
+                         fitsHeader.naxis3, fitsHeader.naxis3 > 1 ? "s" : "",
+                         imageRequest);
+            else {
+                *imageCountP   = fitsHeader.naxis3;
+                *multiplaneP   = FALSE;
+                *desiredImageP = imageRequest;
+            }
+        } else {
+            if (fitsHeader.naxis3 == 3) {
+                *imageCountP   = 1;
+                *multiplaneP   = TRUE;
+                *desiredImageP = 1;
+            } else if (fitsHeader.naxis3 > 1)
+                pm_error("This FITS file contains multiple (%u) images.  "
+                         "You must specify which one you want with a "
+                         "-image option.", fitsHeader.naxis3);
+            else {
+                *imageCountP   = fitsHeader.naxis3;
+                *multiplaneP   = FALSE;
+                *desiredImageP = 1;
+            }
+        }
+    }
+    if (verbose) {
+        
+        pm_message("FITS stream is %smultiplane", *multiplaneP ? "" : "not ");
+        pm_message("We will take image %u (1 is first) of %u "
+                   "in the FITS stream",
+                   *desiredImageP, *imageCountP);
+    }
+}
+
 
 
 static void
@@ -100,7 +389,7 @@ scanImageForMinMax(FILE *       const ifP,
             unsigned int col;
             for (col = 0; col < cols; ++col) {
                 double val;
-                read_val(ifP, bitpix, &val);
+                readVal(ifP, bitpix, &val);
                 if (image == imagenum || multiplane ) {
                     dmax = MAX(dmax, val);
                     dmin = MIN(dmin, val);
@@ -170,326 +459,225 @@ computeMinMax(FILE *             const ifP,
 
 
 
-int
-main(int argc, char * argv[]) {
-
-    FILE* ifp;
-    int argn, imagenum, image, row;
-    register int col;
-    xelval maxval;
-    double val, frmin, frmax, scale, t;
-    double datamin, datamax;
-    int rows, cols, images, format;
-    struct FITS_Header h;
-    xel** pnmarray;
-    xelval tx, txv[4];
-    const char* fits_name;
-    const char* const usage = "[-image N] [-scanmax] [-printmax] [-min f] [-max f] [FITSfile]";
-
-    int doscan = 0;
-    int forceplain = 0;
-    int forcemin = 0;
-    int forcemax = 0;
-    int printmax = 0;
-    bool multiplane;
-  
-    pnm_init( &argc, argv );
-  
-    argn = 1;
-    imagenum = 0;
-  
-    while ( argn < argc && argv[argn][0] == '-' && argv[argn][1] != '\0' )
-    {
-        if ( pm_keymatch( argv[argn], "-image", 2 ) )
-        {
-            ++argn;
-            if ( argn == argc || sscanf( argv[argn], "%d", &imagenum ) != 1 )
-                pm_usage( usage );
-        }
-        else if ( pm_keymatch( argv[argn], "-max", 3 ) )
-        {
-            ++argn;
-            forcemax = 1;
-            if ( argn == argc || sscanf( argv[argn], "%lf", &frmax ) != 1 )
-                pm_usage( usage );
-        }
-        else if ( pm_keymatch( argv[argn], "-min", 3 ) )
-        {
-            ++argn;
-            forcemin = 1;
-            if ( argn == argc || sscanf( argv[argn], "%lf", &frmin ) != 1 )
-                pm_usage( usage );
-        }
-        else if ( pm_keymatch( argv[argn], "-scanmax", 2 ) )
-            doscan = 1;
-        else if ( pm_keymatch( argv[argn], "-noraw", 2 ) )
-            /* This is for backward compatibility only.  Use the common option
-               -plain now.  (pnm_init() processes -plain).
+static xelval
+determineMaxval(struct cmdlineInfo const cmdline,
+                struct FITS_Header const fitsHeader,
+                double             const datamax,
+                double             const datamin) {
+
+    xelval retval;
+                
+    if (cmdline.omaxvalSpec)
+        retval = cmdline.omaxval;
+    else {
+        if (fitsHeader.bitpix < 0) {
+            /* samples are floating point, which means the resolution
+               could be anything.  So we just pick a convenient maxval
+               of 255.  Before Netpbm 10.20 (January 2004), we did
+               maxval = max - min for floating point as well as
+               integer samples.
             */
-            forceplain = 1;
-        else if ( pm_keymatch( argv[argn], "-printmax", 2 ) )
-            printmax = 1;
-        else
-            pm_usage( usage );
-        ++argn;
-    }
-  
-    if ( argn < argc )
-    {
-        fits_name = argv[argn];
-        ++argn;
+            retval = 255;
+            if (cmdline.verbose)
+                pm_message("FITS image has floating point samples.  "
+                           "Using maxval = %u.", (unsigned int)retval);
+        } else {
+            retval = MAX(1, MIN(PNM_OVERALLMAXVAL, datamax - datamin));
+            if (cmdline.verbose)
+                pm_message("FITS image has samples in the range %d-%d.  "
+                           "Using maxval %u.",
+                           (int)(datamin+0.5), (int)(datamax+0.5),
+                           (unsigned int)retval);
+        }
     }
-    else
-        fits_name = "-";
-
-    if ( argn != argc )
-        pm_usage( usage );
-
-    ifp = pm_openr_seekable(fits_name);
-  
-    read_fits_header( ifp, &h );
-  
-    if ( ! h.simple )
-        pm_error( "FITS file is not in simple format, can't read" );
-    if ( h.naxis != 2 && h.naxis != 3 )
-        pm_message( "Warning: FITS file has %d axes", h.naxis );
-    cols = h.naxis1;
-    rows = h.naxis2;
-    if ( h.naxis == 2 )
-        images = imagenum = 1;
-    else
-        images = h.naxis3;
-    if ( imagenum > images )
-        pm_error( "only %d plane%s in this file", 
-                  images, images > 1 ? "s" : "" );
-    if ( images != 3 && images > 1 && imagenum == 0 )
-        pm_error( "need to specify a plane using the -imagenum flag" );
+    return retval;
+}
 
-    multiplane = ( images == 3 && imagenum == 0 );
 
-    computeMinMax(ifp, images, cols, rows, h, imagenum, multiplane,
-                  forcemin, forcemax, frmin, frmax,
-                  &datamin, &datamax);
 
-    if (h.bitpix < 0) {
-        /* samples are floating point, which means the resolution could be
-           anything.  So we just pick a convenient maxval of 255.  We should
-           have a program option to choose the maxval.  Before Netpbm 10.20
-           (January 2004), we did maxval = max - min for floating point as
-           well as integer samples.
-        */
-        maxval = 255;
-    } else
-        maxval = MAX(1, MIN(PNM_OVERALLMAXVAL, datamax - datamin));
-
-    if ( datamax - datamin == 0 )
-        scale = 1.0;
-    else
-        scale = maxval / ( datamax - datamin );
-
-    /* If printmax option is true, just print and exit. */
-    /* For use in shellscripts.  Ex:                    */
-    /* eval `fitstopnm -printmax $filename | \          */
-    /*   awk '{min = $1; max = $2}\                     */
-    /*         END {print "min=" min; " max=" max}'`    */
-    if (printmax) {
-        printf( "%f %f\n", datamin, datamax);
-        pm_close( ifp );
-        pm_close( stdout );
-        exit( 0 );
-    }
+static void
+convertPgmRaster(FILE *             const ifP,
+                 unsigned int       const cols,
+                 unsigned int       const rows,
+                 xelval             const maxval,
+                 unsigned int       const desiredImage,
+                 unsigned int       const imageCount,
+                 struct FITS_Header const fitsHdr,
+                 double             const scale,
+                 double             const datamin,
+                 xel **             const xels) {
+        
+    unsigned int image;
 
-    if (multiplane)
-        format = PPM_FORMAT;
-    else
-        format = PGM_FORMAT;
+    pm_message("writing PGM file");
 
-    pnmarray = pnm_allocarray( cols, rows );
-
-    switch( PNM_FORMAT_TYPE( format ))
-    {
-    case PGM_TYPE:
-        pm_message( "writing PGM file" );
-        for ( image = 1; image <= imagenum; ++image )
-        {
-            if ( image != imagenum )
-                pm_message( "skipping image plane %d of %d", image, images );
-            else if ( images > 1 )
-                pm_message( "reading image plane %d of %d", image, images );
-            for ( row = 0; row < rows; ++row)
-                for ( col = 0; col < cols; ++col )
+    for (image = 1; image <= desiredImage; ++image) {
+        unsigned int row;
+        if (image != desiredImage)
+            pm_message("skipping image plane %u of %u", image, imageCount);
+        else if (imageCount > 1)
+            pm_message("reading image plane %u of %u", image, imageCount);
+        for (row = 0; row < rows; ++row) {
+            unsigned int col;
+            for (col = 0; col < cols; ++col) {
+                double val;
+                readVal(ifP, fitsHdr.bitpix, &val);
                 {
-                    read_val (ifp, h.bitpix, &val);
-                    t = scale * ( val * h.bscale + h.bzer - datamin );
-                    tx = MAX( 0, MIN( t, maxval ) );
-                    if ( image == imagenum )
-                        PNM_ASSIGN1( pnmarray[row][col], tx );
+                    double const t = scale *
+                        (val * fitsHdr.bscale + fitsHdr.bzer - datamin);
+                    xelval const tx = MAX(0, MIN(t, maxval));
+                    if (image == desiredImage)
+                        PNM_ASSIGN1(xels[row][col], tx);
                 }
+            }
         }
-        break;
-    case PPM_TYPE:
-        pm_message( "writing PPM file" );
-        for ( image = 1; image <= images; image++ )
-        {
-            pm_message( "reading image plane %d of %d", image, images );
-            for ( row = 0; row < rows; row++ )
-                for ( col = 0; col < cols; col++ )
+    } 
+}
+
+
+
+static void
+convertPpmRaster(FILE *             const ifP,
+                 unsigned int       const cols,
+                 unsigned int       const rows,
+                 xelval             const maxval,
+                 struct FITS_Header const fitsHdr,
+                 double             const scale,
+                 double             const datamin,
+                 xel **             const xels) {
+/*----------------------------------------------------------------------------
+   Read the FITS raster from file *ifP into xels[][].  Image dimensions
+   are 'cols' by 'rows'.  The FITS raster is 3 planes composing one
+   image: a red plane followed by a green plane followed by a blue plane.
+-----------------------------------------------------------------------------*/
+    unsigned int plane;
+
+    pm_message("writing PPM file");
+
+    for (plane = 0; plane < 3; ++plane) {
+        unsigned int row;
+        pm_message("reading image plane %u (%s)",
+                   plane, plane == 0 ? "red" : plane == 1 ? "green" : "blue");
+        for (row = 0; row < rows; ++row) {
+            unsigned int col;
+            for (col = 0; col < cols; ++col) {
+                double val;
+                readVal(ifP, fitsHdr.bitpix, &val);
                 {
-                    read_val (ifp, h.bitpix, &val);
-                    txv[1] = PPM_GETR( pnmarray[row][col] );
-                    txv[2] = PPM_GETG( pnmarray[row][col] );
-                    txv[3] = PPM_GETB( pnmarray[row][col] );
-                    t = scale * ( val * h.bscale + h.bzer - datamin );
-                    txv[image] = MAX( 0, MIN( t, maxval ));
-                    PPM_ASSIGN( pnmarray[row][col], txv[1], txv[2], txv[3] );
+                    double const t = scale *
+                        (val * fitsHdr.bscale + fitsHdr.bzer - datamin);
+                    xelval const sample = MAX(0, MIN(t, maxval));
+
+                    switch (plane) {
+                    case 0: PPM_PUTR(xels[row][col], sample); break;
+                    case 1: PPM_PUTG(xels[row][col], sample); break;
+                    case 2: PPM_PUTB(xels[row][col], sample); break;
+                    }
                 }
+            }
         }
-        break;
-    default:
-        pm_error( "can't happen" );
-        break;
     }
-
-    pnm_writepnm( stdout, pnmarray, cols, rows, maxval, format, forceplain );
-    pnm_freearray( pnmarray, rows );
-
-    pm_close( ifp );
-    pm_close( stdout );
-  
-    exit( 0 );
 }
 
-/*
- ** This code will deal properly with integers, no matter what the byte order
- ** or integer size of the host machine.  Sign extension is handled manually
- ** to prevent problems with signed/unsigned characters.  Floating point
- ** values will only be read properly when the host architecture is IEEE-754
- ** conformant.  If you need to tweak this code for other machines, you might
- ** want to snag a copy of the FITS documentation from nssdca.gsfc.nasa.gov
- */
+
 
 static void
-read_val (fp, bitpix, vp)
-    FILE *fp;
-    int bitpix;
-    double *vp;
-
-{
-    int i, ich, ival;
-    long int lval;
-    unsigned char c[8];
-  
-    switch ( bitpix )
-    {
-        /* 8 bit FITS integers are unsigned */
-    case 8:
-        ich = getc( fp );
-        if ( ich == EOF )
-            pm_error( "EOF / read error" );
-        *vp = ich;
-        break;
-      
-    case 16:
-        ich = getc( fp );
-        if ( ich == EOF )
-            pm_error( "EOF / read error" );
-        c[0] = ich;
-        ich = getc( fp );
-        if ( ich == EOF )
-            pm_error( "EOF / read error" );
-        c[1] = ich;
-        if (c[0] & 0x80)
-            ival = ~0xFFFF | c[0]<<8 | c[1];
-        else
-            ival = c[0]<<8 | c[1];
-        *vp = ival;
-        break;
-      
-    case 32:
-        for (i=0; i<4; i++) {
-            ich = getc( fp );
-            if ( ich == EOF )
-                pm_error( "EOF / read error" );
-            c[i] = ich;
-        }
-        if (c[0] & 0x80)
-            lval = ~0xFFFFFFFF |
-                c[0]<<24 | c[1]<<16 | c[2]<<8 | c[3];
-        else
-            lval = c[0]<<24 | c[1]<<16 | c[2]<<8 | c[3];
-        *vp = lval;
-        break;
-      
-    case -32:
-        for (i=0; i<4; i++) {
-            ich = getc( fp );
-            if ( ich == EOF )
-                pm_error( "EOF / read error" );
-            c[i] = ich;
-        }
-        *vp = *( (float *) c);
-        break;
-      
-    case -64:
-        for (i=0; i<8; i++) {
-            ich = getc( fp );
-            if ( ich == EOF )
-                pm_error( "EOF / read error" );
-            c[i] = ich;
-        }
-        *vp = *( (double *) c);
-        break;
-      
-    default:
-        pm_error( "Strange bitpix in read_value" );
+convertRaster(FILE *             const ifP,
+              unsigned int       const cols,
+              unsigned int       const rows,
+              xelval             const maxval,
+              bool               const forceplain,
+              bool               const multiplane,
+              unsigned int       const desiredImage,
+              unsigned int       const imageCount,
+              struct FITS_Header const fitsHdr,
+              double             const scale,
+              double             const datamin) {
+
+    xel ** xels;
+    int format;
+
+    xels = pnm_allocarray(cols, rows);
+
+    if (multiplane) {
+        format = PPM_FORMAT;
+        convertPpmRaster(ifP, cols, rows, maxval, fitsHdr, scale, datamin,
+                         xels);
+    } else {
+        format = PGM_FORMAT;
+        convertPgmRaster(ifP, cols, rows, maxval,
+                         desiredImage, imageCount, fitsHdr, scale, datamin,
+                         xels);
     }
+    pnm_writepnm(stdout, xels, cols, rows, maxval, format, forceplain);
+    pnm_freearray(xels, rows);
 }
 
-static void
-read_fits_header( fp, hP )
-    FILE* fp;
-    struct FITS_Header* hP;
-{
-    char buf[80];
-    int seen_end;
-    int i;
-    char c;
+
+
+int
+main(int argc, char * argv[]) {
+
+    struct cmdlineInfo cmdline;
+    FILE * ifP;
+    unsigned int cols, rows;
+    xelval maxval;
+    double scale;
+    double datamin, datamax;
+    struct FITS_Header fitsHeader;
+
+    unsigned int imageCount;
+    unsigned int desiredImage;
+        /* Plane number (starting at one) of plane that contains the image
+           we want.
+        */
+    bool multiplane;
+        /* This is a one-image multiplane stream; 'desiredImage'
+           is undefined
+        */
   
-    seen_end = 0;
-    hP->simple = 0;
-    hP->bzer = 0.0;
-    hP->bscale = 1.0;
-    hP->datamin = - DBL_MAX;
-    hP->datamax = DBL_MAX;
+    pnm_init( &argc, argv );
   
-    while ( ! seen_end )
-        for ( i = 0; i < 36; ++i )
-        {
-            read_card( fp, buf );
-    
-            if ( sscanf( buf, "SIMPLE = %c", &c ) == 1 )
-            {
-                if ( c == 'T' || c == 't' )
-                    hP->simple = 1;
-            }
-            else if ( sscanf( buf, "BITPIX = %d", &(hP->bitpix) ) == 1 );
-            else if ( sscanf( buf, "NAXIS = %d", &(hP->naxis) ) == 1 );
-            else if ( sscanf( buf, "NAXIS1 = %d", &(hP->naxis1) ) == 1 );
-            else if ( sscanf( buf, "NAXIS2 = %d", &(hP->naxis2) ) == 1 );
-            else if ( sscanf( buf, "NAXIS3 = %d", &(hP->naxis3) ) == 1 );
-            else if ( sscanf( buf, "DATAMIN = %lf", &(hP->datamin) ) == 1 );
-            else if ( sscanf( buf, "DATAMAX = %lf", &(hP->datamax) ) == 1 );
-            else if ( sscanf( buf, "BZERO = %lf", &(hP->bzer) ) == 1 );
-            else if ( sscanf( buf, "BSCALE = %lf", &(hP->bscale) ) == 1 );
-            else if ( strncmp( buf, "END ", 4 ) == 0 ) seen_end = 1;
-        }
-}
+    parseCommandLine(argc, argv, &cmdline);
 
-static void
-read_card( fp, buf )
-    FILE* fp;
-    char* buf;
-{
-    if ( fread( buf, 1, 80, fp ) == 0 )
-        pm_error( "error reading header" );
+    ifP = pm_openr(cmdline.inputFileName);
+
+    readFitsHeader(ifP, &fitsHeader);
+  
+    if (!fitsHeader.simple)
+        pm_error("FITS file is not in simple format, can't read");
+
+    if (fitsHeader.naxis != 2 && fitsHeader.naxis != 3)
+        pm_message("Warning: FITS file has %u axes", fitsHeader.naxis);
+
+    cols = fitsHeader.naxis1;
+    rows = fitsHeader.naxis2;
+
+    interpretPlanes(fitsHeader, cmdline.image, cmdline.verbose,
+                    &imageCount, &multiplane, &desiredImage);
+
+    computeMinMax(ifP, imageCount, cols, rows, fitsHeader,
+                  desiredImage, multiplane,
+                  cmdline.minSpec, cmdline.maxSpec,
+                  cmdline.min, cmdline.max,
+                  &datamin, &datamax);
+
+    maxval = determineMaxval(cmdline, fitsHeader, datamax, datamin);
+
+    if (datamax - datamin == 0)
+        scale = 1.0;
+    else
+        scale = maxval / (datamax - datamin);
+
+    if (cmdline.printmax)
+        printf("%f %f\n", datamin, datamax);
+    else
+        convertRaster(ifP, cols, rows, maxval, cmdline.noraw,
+                      multiplane, desiredImage, imageCount,
+                      fitsHeader, scale, datamin);
+
+    pm_close(ifP);
+    pm_close(stdout);
+
+    return 0;
 }