proj/mag/harp/libs/matlab/mfile-mex/fits/savefitsa.c

#include "mex.h"  /* must appear first */
#include <math.h>
#include <string.h>
#include "fitsio.h"
#include "mexhead.h"
#include "Doc/savefitsa_docstring.h"  /* autogenerated from this file */

/****************************************************************

%savefitsa: save array as file in FITS format
% 
% error = savefitsa(array, file, bitpix, bscale, bzero, blank, mode)
% * array is a d-dimensional real-valued array; it is written in
% natural Matlab/Fortran order to filename.  (The transposed
% version is just "savefits".)
% * filename is silently overwritten if it exists.
% * bitpix must belong to {8,16,32,64,-32,-64}; negative numbers imply
% floating-point FITS output, 8 is *unsigned* character, and 16, 32,
% 64 are signed integers.
% * bscale and bzero are most useful when integers are written.
% Whenever they are supplied, the value actually written to the file, 
% which must belong to the corresponding range, is:
%    (array - bzero) / bscale
% Therefore, unsigned integers may be written using bscale = 1
% and bzero = 2^15, 2^31, etc.
% * blank is the value of the blank flag to be used in the FITS file. 
% A blank will be inserted when an IEEE NaN is encountered in
% the image input.  This will be turned back into a NaN by a FITS 
% reader. 
% * If the blank is supplied as NaN, and bitpix > 0, a standard 
% default is used:  bitpix=8 --> blank = 255; 16 --> -32768; 
% 32 --> -2^31; 64 -> -2^52.  
% (The last is the largest double-precision integer.)  
% * If the matrix contains NaNs, and bitpix belongs to {8,16,32,64},
% the blank argument must be given.  In this case, there
% could be a conflict between of the legitimate integer codings 
% and its use as a blank flag.  Of course, whatever value of blank
% is used, the image must not contain any value such that 
% (array - bzero)/bscale == blank, or it too will be interpreted 
% as a NaN.
% * Note that the blank argument is the value put into the FITS
% file, not in the original image.  It is not translated using
% bscale and bzero. 
% * The mode can contain +f or -f to turn on or off a comma-separated
% list of features. Currently we support date (date keyword), 
% history (history saying the save was from this program), and check 
% (FITS checksums).  So, '-date,-history' turns off date and history, 
% in case needed for compatibility.
% * Error is 1 if there was an error, 0 otherwise.  An error
% synopsis is printed if error = 1.
%
% Inputs:
%   real array image;
%   string filename;
%   opt int bitpix = -64;
%   opt real bscale = 1;
%   opt real bzero = 0;
%   opt real blank = {255,-32768,-2^31,-2^52};  -- unused if bitpix<0
%   opt string mode = '+date,+history,+check';
%
% Outputs:
%   int error;
%
% See Also: loadfitsa, savefits

% implemented as a MEX file.  
% turmon, rewritten, aug 2002
% turmon, updated for 64-bit matlab, 4/2010
% if the transposition routine is in an available library, 
% you can get the transposed behavior by defining FITS_TRANSPOSE_MODE
% at compile time.

****************************************************************/
/* Updated by -mex2pymex.py-ver1- on Wed Sep 23 16:55:50 2009 */

#define NARGIN_MIN  2
#define NARGIN_MAX  7
#define NARGOUT_MIN 0
#define NARGOUT_MAX 1

#define ARG_ARRAY    0
#define ARG_FILENAME 1
#define ARG_BITPIX   2
#define ARG_BSCALE   3
#define ARG_BZERO    4
#define ARG_BLANK    5
#define ARG_MODE     6

#define ARG_ERROR    0

#define PROGNAME savefitsa
static const char *progname = "savefitsa";
static const char *in_specs[NARGIN_MAX] = { "RA", "SV", "IS", "RS", "RS", "IS", "SV" };
static const char *in_names[NARGIN_MAX] = { "array", "filename", 
                        "bitpix", "bscale", "bzero", "blank", "mode" };
static const char *out_names[NARGOUT_MAX] = { "error" };

/* 
 * Print out cfitsio error messages, if any, and exit program 
 * Note that stdout and stderr don't exist for matlab.
 */
static
void
fits_printerror(int status)
{
  char errmsg[81]; // cfitsio caps messages at 80 chars

  mexPrintf("%s: Runtime error (%s).  Exiting.\n",
        progname, 
        (status == 0) ? "not a FITSIO error" : "FITSIO trace follows");
  if (status != 0) {
    fits_get_errstatus(status, errmsg);   /* get the error description */
    mexPrintf("status = %d: %s\n", status, errmsg);
    /* get first message; null if stack is empty */
    mexPrintf("Error message stack:\n");
    while (fits_read_errmsg(errmsg))
      mexPrintf("\t%s\n", errmsg);
  }
}

/*
 * extracts various mode parameters from a free-form string
 *
 * returns 0 if the string was OK, an explanatory tag otherwise
 */
static
char *
extract_mode(char *mode, int *Mdate, int *Mhist, int *Mcsum) 
{
  char *word;
  char *sep = ", \t";

  if (!mode) return "null"; /* should not happen */
  for (word = strtok(mode, sep); word; word = strtok(NULL, sep)) {
    /* unscaled */
    if      (strcasecmp(word, "+date"   ) == 0) { *Mdate = 1; }
    else if (strcasecmp(word, "-date"   ) == 0) { *Mdate = 0; }
    else if (strcasecmp(word, "+history") == 0) { *Mhist = 1; }
    else if (strcasecmp(word, "-history") == 0) { *Mhist = 0; }
    else if (strcasecmp(word, "+check"  ) == 0) { *Mcsum = 1; }
    else if (strcasecmp(word, "-check"  ) == 0) { *Mcsum = 0; }
    /* unrecognized word: not OK */
    else return word; 
  }
  return 0; /* OK */
}



static
int
fio_PutFITS(const mxArray *m, char *filename, int bitpix, 
        double bscale, double bzero,
        double blank, int use_blank,
        int put_date, int put_hist, int put_csum)
{
  fitsfile *fptr;  // pointer to the FITS file, defined in fitsio.h
  int status = 0;  // cfitsio status, starts cleared 
  long naxis;      // number of dimensions 
  const mwSize *naxesM;  // length of each dim (matlab mwSize)
  long *naxesF;    // length of each dim (fitsio)
  long nelem;      // total number of elements 
  long i;          // loops over array elements, so is long
  double *data;

  /* set up some constants */
  naxis  = mxGetNumberOfDimensions(m);
  nelem  = mxGetNumberOfElements(m);
  naxesM = mxGetDimensions(m);
  naxesF = mxCalloc(naxis, sizeof(long));
  for (i = 0; i < naxis; i++) 
    naxesF[i] = (long) naxesM[i];
  /* create new FITS file */
  unlink(filename); /* Delete old file if it already exists */
  fits_create_file(&fptr, filename, &status);
  /* initialize the image HDU */
  fits_create_img(fptr, bitpix, naxis, naxesF, &status);
  mxFree(naxesF);
  /* insert bscale and bzero into header */
  if (bscale != 1 || bzero != 0) {
    fits_write_key(fptr, TDOUBLE, "BSCALE", &bscale, "scale", &status);
    fits_write_key(fptr, TDOUBLE, "BZERO",  &bzero,  "offset",  &status);
  }
  /* insert BLANK keyword */
  /* NB: must be inserted as an int */
  if (use_blank) {
    int int_blank = blank;
    fits_write_key(fptr, TINT, "BLANK", &int_blank, "blank", &status);
  }
  /* check for any errors above */
  if (status) {
    fits_printerror(status);
    unlink(filename);
    return 1;
  }
  /* write the array of doubles to the FITS file */
  data = mxGetPr(m);
  if ((bitpix < 0) || (use_blank == 0)) {
    /* do not allow fitsio to muck with special values */
    fits_write_img(fptr, TDOUBLE, 1L, nelem, data, &status);
  } else /* use_blank == 1 */ {
    /* We need to translate all NaNs in data to the fits BLANK number.
     * Although cfitsio provides blank-aware writing routines, they
     * rely on passing in the exact value of the NaN number.  Because
     * IEEE nans can be several values, we must test for them outside
     * cfitsio, using isnan() */
    double *data2; /* data, but with translated NaNs */
    const double blank_value = bscale*blank+bzero; /* = BLANK when scaled */
    if ((data2 = mxCalloc(nelem, sizeof(double))) == NULL) {
      fits_printerror(status);
      unlink(filename); /* delete the partial filename */
      return 1;
    }
    for (i = 0; i < nelem; i++)
      data2[i] = isnan(data[i]) ? blank_value : data[i];
    fits_write_img(fptr, TDOUBLE, 1L, nelem, data2, &status);
    mxFree(data2);
  }
  // metadata
  if (put_date)
    fits_write_date(fptr, &status);
  if (put_hist)
    fits_write_history(fptr, "Saved from Matlab, savefitsa ver. 1.2", &status);
  if (put_csum)
    fits_write_chksum(fptr, &status);
  // check for any errors above
  if (status) {
    fits_printerror(status);
    unlink(filename);
    return 1;
  }
  /* close the file */
  if (fits_close_file(fptr, &status)) {
    /* close failed: bail */
    fits_printerror(status);
    unlink(filename); /* delete the partial filename */
    return 1;
  }
  /* success! */
  return 0;
}

#ifdef StaticP  /* undefined under mex */
StaticP
#endif
void 
mexFunction(int nlhs, 
        mxArray *plhs[], 
        int nrhs, 
        const mxArray *prhs[])
{
   char *filename;
   char *mode, *word;   // label for converted mode string and substring
   int error;           // error status
   char errstr[256];
   // default fits parameters
   int    bitpix = -64;
   double bscale = 1.0;
   double bzero  = 0.0;
   int use_blank = 0;
   double blank;
   // defaults for modes
   int mode_date = 1;
   int mode_hist = 1;
   int mode_csum = 1;

  /* Hook for introspection (function signature, docstring) */
  if (nrhs < 0) { 
    plhs[0] = mxt_PackSignature((mxt_Signature) (-nrhs), 
                NARGIN_MIN, NARGIN_MAX, 
                NARGOUT_MIN, NARGOUT_MAX, 
                in_names, in_specs, out_names, docstring);
    return;
  }
  /* check number of arguments */
  if ((nrhs < NARGIN_MIN) || (nrhs > NARGIN_MAX))
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: Expect %d <= input args <= %d",
               progname, NARGIN_MIN, NARGIN_MAX), errstr));
  if ((nlhs < NARGOUT_MIN) || (nlhs > NARGOUT_MAX))
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: Expect %d <= output args <= %d",
               progname, NARGOUT_MIN, NARGOUT_MAX), errstr));
  mexargparse(nrhs, prhs, in_names, in_specs, NULL, progname);

  /* arg checks */
  start_sizechecking();
  sizeinit(prhs[ARG_FILENAME]);
  sizeisM(1);
  sizecheck_msg(progname, in_names, ARG_FILENAME);

  /* read in filename */
  filename = mxArrayToString(prhs[ARG_FILENAME]);
  if (!filename)
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: Trouble converting filename string",
               progname), errstr));
  
  /* read in bitpix, if appropriate */
  if (nrhs > ARG_BITPIX)
    bitpix = (int) mxGetScalar(prhs[ARG_BITPIX]);
  if ((bitpix !=   8) && (bitpix !=  16) && (bitpix != 32) && (bitpix != 64) &&
      (bitpix != -32) && (bitpix != -64))
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
             "%s: can't do bitpix = %d.", 
              progname, bitpix), errstr));

  /* read in bscale and bzero if appropriate */
  if (nrhs > ARG_BSCALE)
    bscale = mxGetScalar(prhs[ARG_BSCALE]);
  if (nrhs > ARG_BZERO)
    bzero = mxGetScalar(prhs[ARG_BZERO]);

  /* read in blank; only allow its use for non-floatingpoint output */
  if (nrhs > ARG_BLANK && bitpix > 0) {
    use_blank = 1;
    blank = mxGetScalar(prhs[ARG_BLANK]);
    // turmon: fill in if NaN, but allow 0 to be a blank
    if (isnan(blank)) {
      /* substitute the "standard value" */
      if (bitpix == 8)       blank =  255.0; /* 2^8-1 */
      else if (bitpix == 16) blank = -32768.0; /* -2^15 */
      else if (bitpix == 32) blank = -2147483648.0; /* -2^31 */
      else if (bitpix == 64) blank = -4503599627370496.0; /* -2^52 */
    }
  }
  // default keyword modes
  if (nrhs > ARG_MODE) {
    // OK not to specify it
    if ((mode = mxArrayToString(prhs[ARG_MODE])) == NULL)
      mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
                             "%s: bad mode (non-string?). "
                             "Could not convert mode arg to string.", 
                             progname), errstr));
    if ((word = extract_mode(mode, &mode_date, &mode_hist, &mode_csum)) != NULL)
      mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
                             "%s: bad word <%s> in mode <%.80s>", 
                             progname, word, /* need a fresh copy of the mode */
                             mxArrayToString(prhs[ARG_MODE])), errstr));
    mxFree(mode); // done with the mode string
  }


  /* write image */
  error = fio_PutFITS(prhs[ARG_ARRAY], 
              filename, 
              bitpix, bscale, bzero, blank, use_blank,
              mode_date, mode_hist, mode_csum);
  mxFree(filename);
  /* plhs[0] always exists */
  plhs[ARG_ERROR] = mxCreateDoubleScalar((double) error);
}


/* Hook for generic tail matter */
#ifdef MEX2C_TAIL_HOOK
#include "mex2c_tail.h"
#endif

---