/* fitstopnm.c - read a FITS file and produce a PNM.
**
** Copyright (C) 1989 by Jef Poskanzer.
**
** Permission to use, copy, modify, and distribute this software and its
** documentation for any purpose and without fee is hereby granted, provided
** that the above copyright notice appear in all copies and that both that
** copyright notice and this permission notice appear in supporting
** documentation. This software is provided "as is" without express or
** implied warranty.
**
** Hacked up version by Daniel Briggs (dbriggs@nrao.edu) 20-Oct-92
**
** Include floating point formats, more or less. Will only work on
** machines that understand IEEE-754. Added -scanmax -printmax
** -min -max and -noraw. Ignore axes past 3, instead of error (many packages
** use pseudo axes). Use a finite scale when max=min. NB: Min and max
** are the real world FITS values (scaled), so watch out when bzer & bscale
** are not 0 & 1. Datamin & datamax interpreted correctly in scaled case,
** and initialization changed to less likely values. If datamin & max are
** not present in the header, the a first pass is made to determine them
** from the array values.
**
** Modified by Alberto Accomazzi (alberto@cfa.harvard.edu), Dec 1, 1992.
**
** Added understanding of 3-plane FITS files, the program is renamed
** fitstopnm. Fixed some non-ansi declarations (DBL_MAX and FLT_MAX
** replace MAXDOUBLE and MAXFLOAT), fixed some scaling parameters to
** map the full FITS data resolution to the maximum PNM resolution,
** disabled min max scanning when reading from stdin.
*/
/*
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"
#include "mallocvar.h"
#include "floatcode.h"
#include "shhopt.h"
#include "pnm.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 pm_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) */
pm_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 */
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 */
int naxis3; /* number of points on axis 3 */
double datamin; /* min # (Physical value!) */
double datamax; /* max # " " */
double bzer; /* Physical value = Array value*bscale + bzero */
double bscale;
};
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
values properly only when the host architecture conforms to IEEE-754. If
you need to tweak this code for other machines, you might want to get a
copy of the FITS documentation from nssdca.gsfc.nasa.gov
*/
static void
readFitsChar(FILE * const ifP,
double * const vP) {
/* 8 bit FITS integers are unsigned */
int const ich = getc(ifP);
if (ich == EOF)
pm_error("EOF / read error");
else
*vP = ich;
}
static void
readFitsShort(FILE * const ifP,
double * const vP) {
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;
}
static void
readFitsLong(FILE * const ifP,
double * const vP) {
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;
}
static void
readFitsFloat(FILE * const ifP,
double * const vP) {
unsigned int i;
pm_bigendFloat bigend;
for (i = 0; i < 4; ++i) {
int const ich = getc(ifP);
if (ich == EOF)
pm_error("EOF / read error");
bigend.bytes[i] = ich;
}
*vP = pm_floatFromBigendFloat(bigend);
}
static void
readFitsDouble(FILE * const ifP,
double * const vP) {
unsigned int i;
pm_bigendDouble bigend;
for (i = 0; i < 8; ++i) {
int const ich = getc(ifP);
if (ich == EOF)
pm_error("EOF / read error");
bigend.bytes[i] = ich;
}
*vP = pm_doubleFromBigendDouble(bigend);
}
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,
valFmt const fmt,
double * const vP) {
switch (fmt) {
case VF_CHAR:
readFitsChar(ifP, vP);
break;
case VF_SHORT:
readFitsShort(ifP, vP);
break;
case VF_LONG:
readFitsLong(ifP, vP);
break;
case VF_FLOAT:
readFitsFloat(ifP, vP);
break;
case VF_DOUBLE:
readFitsDouble(ifP, vP);
break;
}
}
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
scanImageForMinMax(FILE * const ifP,
unsigned int const images,
int const cols,
int const rows,
valFmt const valFmt,
double const bscale,
double const bzer,
unsigned int const imagenum,
bool const multiplane,
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;
double fmaxval;
pm_tell2(ifP, &rasterPos, sizeof(rasterPos));
pm_message("Scanning file for scaling parameters");
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;
dmin = fmaxval;
for (image = 1; image <= images; ++image) {
unsigned int row;
for (row = 0; row < rows; ++row) {
unsigned int col;
for (col = 0; col < cols; ++col) {
double val;
readVal(ifP, valFmt, &val);
if (image == imagenum || multiplane ) {
/* Note: if 'val' is NaN, result is 2nd operand */
dmax = MAX(val, dmax);
dmin = MIN(val, dmin);
}
}
}
if (bscale < 0.0) {
double const origDmax = dmax;
dmax = dmin;
dmin = origDmax;
}
}
*dataminP = dmin * bscale + bzer;
*datamaxP = dmax * bscale + bzer;
pm_message("Scan results: min=%f max=%f", *dataminP, *datamaxP);
pm_seek2(ifP, &rasterPos, sizeof(rasterPos));
}
static void
computeMinMax(FILE * const ifP,
unsigned int const images,
int const cols,
int const rows,
struct FITS_Header const h,
unsigned int const imagenum,
bool const multiplane,
bool const forcemin,
bool const forcemax,
double const frmin,
double const frmax,
double * const dataminP,
double * const datamaxP) {
double datamin, datamax;
datamin = -DBL_MAX; /* initial assumption */
datamax = DBL_MAX; /* initial assumption */
if (forcemin)
datamin = frmin;
if (forcemax)
datamax = frmax;
if (datamin == -DBL_MAX)
datamin = h.datamin;
if (datamax == DBL_MAX)
datamax = h.datamax;
if (datamin == -DBL_MAX || datamax == DBL_MAX) {
double scannedDatamin, scannedDatamax;
scanImageForMinMax(ifP, images, cols, rows,
valFmtFromBitpix(h.bitpix), h.bscale, h.bzer,
imagenum, multiplane,
&scannedDatamin, &scannedDatamax);
if (datamin == -DBL_MAX)
datamin = scannedDatamin;
if (datamax == DBL_MAX)
datamax = scannedDatamax;
}
*dataminP = datamin;
*datamaxP = datamax;
}
static xelval
determineMaxval(struct cmdlineInfo const cmdline,
valFmt const valFmt,
double const datamax,
double const datamin) {
xelval retval;
if (cmdline.omaxvalSpec)
retval = cmdline.omaxval;
else {
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
maxval = max - min for floating point as well as
integer samples.
*/
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);
}
}
return retval;
}
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 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
file the top left, bottom left, etc.). We use the common sense,
popular order of row major, top to bottom, left to right. This
has been the case and accepted since 1989, but in 2008, we discovered
that Gimp and ImageMagick do bottom to top.
*/
unsigned int image;
pm_message("writing PGM file");
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, rasterInfo.valFmt, &val);
{
double const t = scale *
(val * rasterInfo.bscale + rasterInfo.bzer - datamin);
xelval const tx = MAX(0, MIN(t, maxval));
if (image == desiredImage)
PNM_ASSIGN1(xels[row][col], tx);
}
}
}
}
}
static void
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
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, rasterInfo.valFmt, &val);
{
double const t = scale *
(val * rasterInfo.bscale + rasterInfo.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;
}
}
}
}
}
}
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 fitsRasterInfo const rasterInfo,
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, rasterInfo, scale, datamin,
xels);
} else {
format = PGM_FORMAT;
convertPgmRaster(ifP, cols, rows, maxval,
desiredImage, imageCount, rasterInfo, scale, datamin,
xels);
}
pnm_writepnm(stdout, xels, cols, rows, maxval, format, forceplain);
pnm_freearray(xels, rows);
}
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;
struct fitsRasterInfo rasterInfo;
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
*/
pnm_init( &argc, argv );
parseCommandLine(argc, argv, &cmdline);
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;
rasterInfo.bscale = fitsHeader.bscale;
rasterInfo.bzer = fitsHeader.bzer;
rasterInfo.valFmt = valFmtFromBitpix(fitsHeader.bitpix);
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, rasterInfo.valFmt, 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,
rasterInfo, scale, datamin);
pm_close(ifP);
pm_close(stdout);
return 0;
}