diff options
Diffstat (limited to 'converter/other/fitstopnm.c')
| -rw-r--r-- | converter/other/fitstopnm.c | 203 |
1 files changed, 128 insertions, 75 deletions
diff --git a/converter/other/fitstopnm.c b/converter/other/fitstopnm.c index 73564c4b..bdf5c78a 100644 --- a/converter/other/fitstopnm.c +++ b/converter/other/fitstopnm.c @@ -34,10 +34,16 @@ The official specification of FITS format (which is for more than just visual images) is at ftp://legacy.gsfc.nasa.gov/fits_info/fits_office/fits_standard.pdf + + An example FITS file is at + + http://fits.gsfc.nasa.gov/nrao_data/tests/incunabula/mndrll-8.fits + */ #include <string.h> #include <float.h> +#include <assert.h> #include "pm_config.h" #include "pm_c_util.h" @@ -48,7 +54,7 @@ -struct cmdlineInfo { +struct CmdlineInfo { const char * inputFileName; unsigned int image; /* zero if unspecified */ float max; @@ -69,8 +75,8 @@ struct cmdlineInfo { static void -parseCommandLine(int argc, char ** argv, - struct cmdlineInfo * const cmdlineP) { +parseCommandLine(int argc, const 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. @@ -82,7 +88,7 @@ parseCommandLine(int argc, char ** argv, was passed to us as the argv array. We also trash *argv. --------------------------------------------------------------------------*/ optEntry * option_def; - /* Instructions to optParseOptions3 on how to parse our options. */ + /* Instructions to pm_optParseOptions3 on how to parse our options. */ optStruct3 opt; unsigned int imageSpec; @@ -114,7 +120,7 @@ parseCommandLine(int argc, char ** argv, /* Set some defaults the lazy way (using multiple setting of variables) */ - optParseOptions3(&argc, argv, opt, sizeof(opt), 0); + pm_optParseOptions3(&argc, (char**)argv, opt, sizeof(opt), 0); /* Uses and sets argc, argv, and some of *cmdlineP and others. */ if (imageSpec) { @@ -139,13 +145,17 @@ parseCommandLine(int argc, char ** argv, pm_error("Too many arguments (%u). The only non-option argument " "is the input file name.", argc-1); } + free(option_def); } struct FITS_Header { int simple; /* basic format or not */ - int bitpix; /* number of bits per pixel */ + int bitpix; + /* number of bits per pixel, positive for integer, negative + for floating point + */ int naxis; /* number of axes */ int naxis1; /* number of points on axis 1 */ int naxis2; /* number of points on axis 2 */ @@ -157,6 +167,16 @@ struct FITS_Header { }; +typedef enum { + VF_CHAR, VF_SHORT, VF_LONG, VF_FLOAT, VF_DOUBLE +} valFmt; + +struct fitsRasterInfo { + valFmt valFmt; + double bzer; + double bscale; +}; + /* This code deals properly with integers, no matter what the byte order or integer size of the host machine. We handle sign extension manually to prevent problems with signed/unsigned characters. We read floating point @@ -276,34 +296,56 @@ readFitsDouble(FILE * const ifP, +static valFmt +valFmtFromBitpix(int const bitpix) { +/*---------------------------------------------------------------------------- + Return the format of a "value" in the FITS file, given the value + of the BITPIX header in the FITS file. + + BITPIX has a stupid format wherein it is fundamentally the number + of bits per value, but its sign indicates whether it is integer + or floating point. +-----------------------------------------------------------------------------*/ + switch (bitpix) { + case +8: return VF_CHAR; + case +16: return VF_SHORT; + case +32: return VF_LONG; + case -32: return VF_FLOAT; + case -64: return VF_DOUBLE; + default: + /* Every possibility is covered above. */ + assert(false); + return 0; /* quiet compiler warning */ + } +} + + + static void readVal(FILE * const ifP, - int const bitpix, + valFmt const fmt, double * const vP) { - switch (bitpix) { - case 8: + switch (fmt) { + case VF_CHAR: readFitsChar(ifP, vP); break; - case 16: + case VF_SHORT: readFitsShort(ifP, vP); break; - case 32: + case VF_LONG: readFitsLong(ifP, vP); break; - case -32: + case VF_FLOAT: readFitsFloat(ifP, vP); break; - case -64: + case VF_DOUBLE: readFitsDouble(ifP, vP); break; - - default: - pm_error("Strange bitpix value %d in readVal()", bitpix); } } @@ -419,7 +461,7 @@ scanImageForMinMax(FILE * const ifP, unsigned int const images, int const cols, int const rows, - unsigned int const bitpix, + valFmt const valFmt, double const bscale, double const bzer, unsigned int const imagenum, @@ -427,6 +469,9 @@ scanImageForMinMax(FILE * const ifP, double * const dataminP, double * const datamaxP) { + /* Note that a value in the file might be Not-A-Number. We ignore + such entries in computing the minimum and maximum for the image. + */ double dmax, dmin; unsigned int image; pm_filepos rasterPos; @@ -436,14 +481,12 @@ scanImageForMinMax(FILE * const ifP, pm_message("Scanning file for scaling parameters"); - switch (bitpix) { - case 8: fmaxval = 255.0; break; - case 16: fmaxval = 65535.0; break; - case 32: fmaxval = 4294967295.0; break; - case -32: fmaxval = FLT_MAX; break; - case -64: fmaxval = DBL_MAX; break; - default: - pm_error("unusual bits per pixel (%u), can't read", bitpix); + switch (valFmt) { + case VF_CHAR: fmaxval = 255.0; break; + case VF_SHORT: fmaxval = 65535.0; break; + case VF_LONG: fmaxval = 4294967295.0; break; + case VF_FLOAT: fmaxval = FLT_MAX; break; + case VF_DOUBLE: fmaxval = DBL_MAX; break; } dmax = -fmaxval; @@ -454,10 +497,11 @@ scanImageForMinMax(FILE * const ifP, unsigned int col; for (col = 0; col < cols; ++col) { double val; - readVal(ifP, bitpix, &val); + readVal(ifP, valFmt, &val); if (image == imagenum || multiplane ) { - dmax = MAX(dmax, val); - dmin = MIN(dmin, val); + /* Note: if 'val' is NaN, result is 2nd operand */ + dmax = MAX(val, dmax); + dmin = MIN(val, dmin); } } } @@ -509,7 +553,8 @@ computeMinMax(FILE * const ifP, if (datamin == -DBL_MAX || datamax == DBL_MAX) { double scannedDatamin, scannedDatamax; - scanImageForMinMax(ifP, images, cols, rows, h.bitpix, h.bscale, h.bzer, + scanImageForMinMax(ifP, images, cols, rows, + valFmtFromBitpix(h.bitpix), h.bscale, h.bzer, imagenum, multiplane, &scannedDatamin, &scannedDatamax); @@ -525,8 +570,8 @@ computeMinMax(FILE * const ifP, static xelval -determineMaxval(struct cmdlineInfo const cmdline, - struct FITS_Header const fitsHeader, +determineMaxval(struct CmdlineInfo const cmdline, + valFmt const valFmt, double const datamax, double const datamin) { @@ -535,7 +580,7 @@ determineMaxval(struct cmdlineInfo const cmdline, if (cmdline.omaxvalSpec) retval = cmdline.omaxval; else { - if (fitsHeader.bitpix < 0) { + if (valFmt == VF_FLOAT || valFmt == VF_DOUBLE) { /* 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 @@ -561,16 +606,16 @@ determineMaxval(struct cmdlineInfo const cmdline, 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) { +convertPgmRaster(FILE * const ifP, + unsigned int const cols, + unsigned int const rows, + xelval const maxval, + unsigned int const desiredImage, + unsigned int const imageCount, + struct fitsRasterInfo const rasterInfo, + double const scale, + double const datamin, + xel ** const xels) { /* Note: the FITS specification does not give the association between file position and image position (i.e. is the first pixel in the @@ -581,8 +626,7 @@ convertPgmRaster(FILE * const ifP, */ unsigned int image; - pm_message("Writing PPM file " - "(Probably not what you want - consider an -image option)"); + pm_message("writing PGM file"); for (image = 1; image <= desiredImage; ++image) { unsigned int row; @@ -594,10 +638,10 @@ convertPgmRaster(FILE * const ifP, unsigned int col; for (col = 0; col < cols; ++col) { double val; - readVal(ifP, fitsHdr.bitpix, &val); + readVal(ifP, rasterInfo.valFmt, &val); { double const t = scale * - (val * fitsHdr.bscale + fitsHdr.bzer - datamin); + (val * rasterInfo.bscale + rasterInfo.bzer - datamin); xelval const tx = MAX(0, MIN(t, maxval)); if (image == desiredImage) PNM_ASSIGN1(xels[row][col], tx); @@ -610,14 +654,14 @@ convertPgmRaster(FILE * const ifP, 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) { +convertPpmRaster(FILE * const ifP, + unsigned int const cols, + unsigned int const rows, + xelval const maxval, + struct fitsRasterInfo const rasterInfo, + 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 @@ -625,7 +669,8 @@ convertPpmRaster(FILE * const ifP, -----------------------------------------------------------------------------*/ unsigned int plane; - pm_message("writing PPM file"); + pm_message("Writing PPM file " + "(Probably not what you want - consider an -image option)"); for (plane = 0; plane < 3; ++plane) { unsigned int row; @@ -635,10 +680,10 @@ convertPpmRaster(FILE * const ifP, unsigned int col; for (col = 0; col < cols; ++col) { double val; - readVal(ifP, fitsHdr.bitpix, &val); + readVal(ifP, rasterInfo.valFmt, &val); { double const t = scale * - (val * fitsHdr.bscale + fitsHdr.bzer - datamin); + (val * rasterInfo.bscale + rasterInfo.bzer - datamin); xelval const sample = MAX(0, MIN(t, maxval)); switch (plane) { @@ -655,17 +700,17 @@ convertPpmRaster(FILE * const ifP, static void -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) { +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 fitsRasterInfo const rasterInfo, + double const scale, + double const datamin) { xel ** xels; int format; @@ -674,12 +719,12 @@ convertRaster(FILE * const ifP, if (multiplane) { format = PPM_FORMAT; - convertPpmRaster(ifP, cols, rows, maxval, fitsHdr, scale, datamin, + convertPpmRaster(ifP, cols, rows, maxval, rasterInfo, scale, datamin, xels); } else { format = PGM_FORMAT; convertPgmRaster(ifP, cols, rows, maxval, - desiredImage, imageCount, fitsHdr, scale, datamin, + desiredImage, imageCount, rasterInfo, scale, datamin, xels); } pnm_writepnm(stdout, xels, cols, rows, maxval, format, forceplain); @@ -689,15 +734,16 @@ convertRaster(FILE * const ifP, int -main(int argc, char * argv[]) { +main(int argc, const char * argv[]) { - struct cmdlineInfo cmdline; + struct CmdlineInfo cmdline; FILE * ifP; unsigned int cols, rows; xelval maxval; double scale; double datamin, datamax; struct FITS_Header fitsHeader; + struct fitsRasterInfo rasterInfo; unsigned int imageCount; unsigned int desiredImage; @@ -709,7 +755,7 @@ main(int argc, char * argv[]) { is undefined */ - pnm_init( &argc, argv ); + pm_proginit(&argc, argv); parseCommandLine(argc, argv, &cmdline); @@ -726,6 +772,10 @@ main(int argc, char * argv[]) { cols = fitsHeader.naxis1; rows = fitsHeader.naxis2; + rasterInfo.bscale = fitsHeader.bscale; + rasterInfo.bzer = fitsHeader.bzer; + rasterInfo.valFmt = valFmtFromBitpix(fitsHeader.bitpix); + interpretPlanes(fitsHeader, cmdline.image, cmdline.verbose, &imageCount, &multiplane, &desiredImage); @@ -735,7 +785,7 @@ main(int argc, char * argv[]) { cmdline.min, cmdline.max, &datamin, &datamax); - maxval = determineMaxval(cmdline, fitsHeader, datamax, datamin); + maxval = determineMaxval(cmdline, rasterInfo.valFmt, datamax, datamin); if (datamax - datamin == 0) scale = 1.0; @@ -747,10 +797,13 @@ main(int argc, char * argv[]) { else convertRaster(ifP, cols, rows, maxval, cmdline.noraw, multiplane, desiredImage, imageCount, - fitsHeader, scale, datamin); + rasterInfo, scale, datamin); pm_close(ifP); pm_close(stdout); return 0; } + + + |