/** @file psFits.c * * @brief Contains Fits I/O routines * * @ingroup FileIO * * @author Robert DeSonia, MHPCC * * @version $Revision: 1.25 $ $Name: $ * @date $Date: 2005/04/05 19:33:30 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii */ #include #include "psFits.h" #include "string.h" #include "psError.h" #include "psFileUtilsErrors.h" #include "psImageExtraction.h" #include "psMemory.h" #include "psString.h" #include "psLogMsg.h" #include "psTrace.h" #define MAX_STRING_LENGTH 256 // maximum length string for FITS routines static char* standardFitsKeys[] = { NULL }; static psElemType convertFitsToPsType(int datatype) { switch (datatype) { case TBYTE: return PS_TYPE_U8; case TSBYTE: return PS_TYPE_S8; case TSHORT: return PS_TYPE_S16; case TUSHORT: return PS_TYPE_U16; case TLONG: if (sizeof(long) == 8) { return PS_TYPE_S64; } // no break case TINT: return PS_TYPE_S32; case TULONG: if (sizeof(unsigned long) == 8) { return PS_TYPE_U64; } // no break case TUINT: return PS_TYPE_U32; case TLONGLONG: return PS_TYPE_S64; case TFLOAT: return PS_TYPE_F32; case TDOUBLE: return PS_TYPE_F64; case TCOMPLEX: return PS_TYPE_C32; case TDBLCOMPLEX: return PS_TYPE_C64; case TLOGICAL: return PS_TYPE_BOOL; default: psError(PS_ERR_IO, true, "Unknown FITS datatype, %d.", datatype); return PS_TYPE_PTR; } } static bool convertPsTypeToFits(psElemType type, int* bitPix, double* bZero, int* dataType) { int bitpix; int datatype; double bzero = 0.0; switch (type) { case PS_TYPE_U8: bitpix = BYTE_IMG; datatype = TBYTE; break; case PS_TYPE_BOOL: case PS_TYPE_S8: bitpix = BYTE_IMG; bzero = INT8_MIN; datatype = TSBYTE; break; case PS_TYPE_U16: bitpix = SHORT_IMG; bzero = -1.0 * INT16_MIN; datatype = TUSHORT; break; case PS_TYPE_S16: bitpix = SHORT_IMG; datatype = TSHORT; break; case PS_TYPE_U32: bitpix = LONG_IMG; bzero = -1.0 * INT32_MIN; datatype = TUINT; break; case PS_TYPE_S32: bitpix = LONG_IMG; datatype = TINT; break; case PS_TYPE_F32: bitpix = FLOAT_IMG; datatype = TFLOAT; break; case PS_TYPE_F64: bitpix = DOUBLE_IMG; datatype = TDOUBLE; break; default: { char* typeStr; PS_TYPE_NAME(typeStr,type); psError(PS_ERR_BAD_PARAMETER_TYPE, true, PS_ERRORTEXT_psFits_TYPE_UNSUPPORTED, typeStr); return false; } } // pass the requested parameters (NULL parameters are not set, of course). if (bitPix != NULL) { *bitPix = bitpix; } if (dataType != NULL) { *dataType = datatype; } if (bZero != NULL) { *bZero = bzero; } return true; } static bool convertMetadataTypeToBinaryTForm(psMetadataType type, char** fitsType) { switch (type) { case PS_META_BOOL: *fitsType = psStringCopy("1L"); break; case PS_META_S32: *fitsType = psStringCopy("1J"); break; case PS_META_F32: *fitsType = psStringCopy("1E"); break; case PS_META_F64: *fitsType = psStringCopy("1D"); break; // XXX: Handle other types, e.g., Vectors, etc. default: return false; } return true; } static bool isHDUEmpty(const psFits* fits) { /* check for keys - no keys means this is really an empty HDU */ int keysexist = -1; int morekeys; int status = 0; fits_get_hdrspace(fits->p_fd, &keysexist, &morekeys, &status); // if no keys exist and not primary HDU, this really is an empty HDU if (keysexist == 0) { return true; } return false; } static void fitsFree(psFits* fits) { int status = 0; if (fits != NULL) { (void)fits_close_file(fits->p_fd, &status); psFree((void*)fits->filename); } } psFits* psFitsAlloc(const char* name) { int status = 0; fitsfile *fptr = NULL; /* Pointer to the FITS file */ if (name == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_FILENAME_NULL); return NULL; } /* Open/Create the FITS file */ if (access(name, F_OK) == 0) { // file exists (void)fits_open_file(&fptr, name, READWRITE, &status); if (fptr == NULL) { // if failed, try openning as just read-only status = 0; (void)fits_open_file(&fptr, name, READONLY, &status); } if (fptr == NULL || status != 0) { char fitsErr[MAX_STRING_LENGTH]; fits_get_errstatus(status, fitsErr); psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psFits_FILENAME_INVALID, name, fitsErr); return NULL; } } else { // file does not exist, so create. (void)fits_create_file(&fptr, name, &status); if (fptr == NULL || status != 0) { char fitsErr[MAX_STRING_LENGTH]; fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_FILENAME_CREATE_FAILED, name, fitsErr); return NULL; } } psFits* fits = psAlloc(sizeof(psFits)); fits->filename = psAlloc(strlen(name)+1); fits->p_fd = fptr; strcpy((char*)fits->filename,name); p_psMemSetDeallocator(fits,(psFreeFcn)fitsFree); return fits; } bool psFitsMoveExtName(const psFits* fits, const char* extname) { int status = 0; if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); return false; } if (extname == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_EXTNAME_NULL); return false; } if (fits_movnam_hdu(fits->p_fd, ANY_HDU, (char*)extname, 0, &status) != 0) { char fitsErr[MAX_STRING_LENGTH]; fits_get_errstatus(status, fitsErr); psError(PS_ERR_LOCATION_INVALID, true, PS_ERRORTEXT_psFits_EXTNAME_INVALID, extname, fits->filename, fitsErr); return false; } return true; } bool psFitsMoveExtNum(const psFits* fits, int extnum, bool relative) { if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); return false; } int status = 0; int hdutype = 0; if (relative) { fits_movrel_hdu(fits->p_fd, extnum, &hdutype, &status); if (status != 0) { char fitsErr[MAX_STRING_LENGTH]; fits_get_errstatus(status, fitsErr); psError(PS_ERR_LOCATION_INVALID, true, PS_ERRORTEXT_psFits_EXTNUM_REL_MOVE_FAILED, extnum, fits->filename, fitsErr); return false; } } else { fits_movabs_hdu(fits->p_fd, extnum+1, &hdutype, &status); if (status != 0) { char fitsErr[MAX_STRING_LENGTH]; fits_get_errstatus(status, fitsErr); psError(PS_ERR_LOCATION_INVALID, true, PS_ERRORTEXT_psFits_EXTNUM_ABS_MOVE_FAILED, extnum, fits->filename, fitsErr); return false; } } return true; } int psFitsGetExtNum(const psFits* fits) { int hdunum; if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); return PS_FITS_TYPE_NONE; } return fits_get_hdu_num(fits->p_fd,&hdunum) - 1; } char* psFitsGetExtName(const psFits* fits) { if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); return NULL; } int status = 0; char name[MAX_STRING_LENGTH]; if (fits_read_key_str(fits->p_fd, "EXTNAME", name, NULL, &status) != 0) { status = 0; if (fits_read_key_str(fits->p_fd, "HDUNAME", name, NULL, &status) != 0) { int num = psFitsGetExtNum(fits); snprintf(name, MAX_STRING_LENGTH, "EXT-%3d",num); } } return psStringCopy(name); } bool psFitsSetExtName(const psFits* fits, const char* name) { if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); return false; } if (name == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_EXTNAME_NULL); return false; } int status = 0; if (fits_update_key_str(fits->p_fd, "EXTNAME", (char*)name, NULL, &status) != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_WRITE_FAILED, fits->filename, fitsErr); return false; } return true; } int psFitsGetSize(const psFits* fits) { if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); return 0; } int num = 0; int status = 0; if (fits_get_num_hdus(fits->p_fd, &num, &status) != 0) { char fitsErr[MAX_STRING_LENGTH]; fits_get_errstatus(status, fitsErr); psError(PS_ERR_LOCATION_INVALID, true, PS_ERRORTEXT_psFits_GETNUMHDUS_FAILED, fits->filename, fitsErr); return 0; } return num; } psFitsType psFitsGetExtType(const psFits* fits) { if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); return PS_FITS_TYPE_NONE; } int status = 0; int hdutype = PS_FITS_TYPE_NONE; if (fits_get_hdu_type(fits->p_fd, &hdutype, &status) != 0) { char fitsErr[MAX_STRING_LENGTH]; fits_get_errstatus(status, fitsErr); psError(PS_ERR_LOCATION_INVALID, true, PS_ERRORTEXT_psFits_GETHDUTYPE_FAILED, fits->filename, fitsErr); return PS_FITS_TYPE_NONE; } if (hdutype == PS_FITS_TYPE_IMAGE && psFitsGetExtNum(fits) > 0 && isHDUEmpty(fits)) { return PS_FITS_TYPE_ANY; } return hdutype; } psMetadata* psFitsReadHeader(psMetadata* out, const psFits* fits) { if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); return NULL; } if (out == NULL) { out = psMetadataAlloc(); if (out == NULL) { psError(PS_ERR_UNKNOWN, false, "Failed to allocate a new psMetadata container."); return NULL; } } // Get number of key names int numKeys = 0; int keyNum = 0; int status = 0; fits_get_hdrpos(fits->p_fd, &numKeys, &keyNum, &status); // Get each key name. Keywords start at one. char keyType; char keyName[MAX_STRING_LENGTH]; char keyValue[MAX_STRING_LENGTH]; char keyComment[MAX_STRING_LENGTH]; psBool tempBool; psBool success; psBool stdKey; for (int i = 1; i <= numKeys; i++) { fits_read_keyn(fits->p_fd, i, keyName, keyValue, keyComment, &status); stdKey = false; int stdKeyIdx = 0; while (standardFitsKeys[stdKeyIdx] != NULL && ! stdKey) { if (strcmp(keyName,standardFitsKeys[stdKeyIdx++]) == 0) { stdKey = true; } } if (keyValue[0] != 0) { // blank values are not handled by fits_get_keytype fits_get_keytype(keyValue, &keyType, &status); } else { keyType = 'C'; } if (status != 0) { break; } if (! stdKey) { switch (keyType) { case 'X': // bit case 'I': // short int. case 'J': // int. case 'B': // byte success = psMetadataAdd(out, PS_LIST_TAIL, keyName, PS_META_S32 | PS_META_DUPLICATE_OK, keyComment, atoi(keyValue)); break; case 'U': // unsigned int. may not fit in a psS32 case 'K': // long int. can't all fit in a psS32 case 'F': success = psMetadataAdd(out, PS_LIST_TAIL, keyName, PS_META_F64 | PS_META_DUPLICATE_OK, keyComment, atof(keyValue)); break; case 'C': // remove the single-quotes at front/end if (keyValue[0] == '\'' && keyValue[strlen(keyValue)-1] == '\'') { keyValue[strlen(keyValue)-1] = '\0'; success = psMetadataAdd(out, PS_LIST_TAIL, keyName, PS_META_STR | PS_META_DUPLICATE_OK, keyComment, keyValue+1); } else { success = psMetadataAdd(out, PS_LIST_TAIL, keyName, PS_META_STR | PS_META_DUPLICATE_OK, keyComment, keyValue); } break; case 'L': tempBool = (keyValue[0] == 'T') ? 1 : 0; success = psMetadataAdd(out, PS_LIST_TAIL, keyName, PS_META_BOOL | PS_META_DUPLICATE_OK, keyComment, tempBool); break; default: psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_METATYPE_INVALID, keyType); return out; } if (!success) { psError(PS_ERR_UNKNOWN, false, PS_ERRORTEXT_psFits_METADATA_ADD_FAILED, keyName); return out; } } } if ( status != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_METADATA_ADD_FAILED, fitsErr); return false; } return out; } psHash* psFitsReadHeaderSet(psHash* out, const psFits* fits) { if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); psFree(out); return NULL; } if (out == NULL) { out = psHashAlloc(10); // XXX: what is the appropriate number of buckets? Got 10 from psMetadataAlloc. if (out == NULL) { psError(PS_ERR_UNKNOWN, false, "Failed to allocate a new psHash container."); return NULL; } } int size = psFitsGetSize(fits); int origPosition = psFitsGetExtNum(fits); for (int lcv=0; lcv < size; lcv++) { psFitsMoveExtNum(fits, lcv, false); char* name = NULL; if (lcv == 0) { name = psStringCopy("PHU"); } else { name = psFitsGetExtName(fits); } psMetadata* header = psFitsReadHeader(NULL, fits); if (name != NULL && header != NULL) { psHashAdd(out, name, header); } else { // XXX: is this a warning or error? psLogMsg(__func__, PS_LOG_WARN, "Failed to read HDU#%d header data.", lcv); } psFree(name); psFree(header); } // reposition to the original position psFitsMoveExtNum(fits, origPosition, false); return out; } psImage* psFitsReadImage(psImage* output, // a psImage to recycle. const psFits* fits, // the psFits object psRegion region, // the region in the FITS image to read int z) // the z-plane in the FITS image cube to read { psS32 status = 0; /* CFITSIO file vars */ psS32 nAxis = 0; psS32 anynull = 0; psS32 bitPix = 0; /* Pixel type */ long nAxes[3]; long firstPixel[3]; /* lower-left corner of image subset */ long lastPixel[3]; /* upper-right corner of image subset */ long increment[3]; /* increment for image subset */ char fitsErr[80] = ""; /* CFITSIO error message string */ psS32 fitsDatatype = 0; psS32 datatype = 0; if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); psFree(output); return NULL; } // check to see if we even are positioned on an image HDU int hdutype; if ( fits_get_hdu_type(fits->p_fd,&hdutype, &status) != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED, fitsErr); return NULL; } if (hdutype != IMAGE_HDU) { psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_NOT_IMAGE_TYPE); return NULL; } /* Get the data type 'bitPix' from the FITS image */ if (fits_get_img_equivtype(fits->p_fd, &bitPix, &status) != 0) { fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_DATATYPE_UNKNOWN, fitsErr); psFree(output); return NULL; } /* Get the dimensions 'nAxis' from the FITS image */ if (fits_get_img_dim(fits->p_fd, &nAxis, &status) != 0) { (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_IMAGE_DIM_UNKNOWN, fitsErr); psFree(output); return NULL; } /* Validate the number of axis */ if ((nAxis < 2) || (nAxis > 3)) { psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_IMAGE_DIMENSION_UNSUPPORTED, nAxis); psFree(output); return NULL; } /* Get the Image size from the FITS file */ if (fits_get_img_size(fits->p_fd, nAxis, nAxes, &status) != 0) { (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_IMAGE_SIZE_UNKNOWN, fitsErr); psFree(output); return NULL; } firstPixel[0] = region.x0 + 1; firstPixel[1] = region.y0 + 1; firstPixel[2] = z + 1; if (region.x1 > 0) { lastPixel[0] = region.x1; } else { lastPixel[0] = nAxes[0] + region.x1; // n.b., region.x1 < 0 } if (region.y1 > 0) { lastPixel[1] = region.y1; } else { lastPixel[1] = nAxes[1] + region.y1; // n.b., region.y1 < 0 } lastPixel[2] = z + 1; increment[0] = 1; increment[1] = 1; increment[2] = 1; switch (bitPix) { case BYTE_IMG: datatype = PS_TYPE_U8; fitsDatatype = TBYTE; break; case SBYTE_IMG: datatype = PS_TYPE_S8; fitsDatatype = TSBYTE; break; case USHORT_IMG: datatype = PS_TYPE_U16; fitsDatatype = TUSHORT; break; case SHORT_IMG: datatype = PS_TYPE_S16; fitsDatatype = TSHORT; break; case ULONG_IMG: datatype = PS_TYPE_U32; fitsDatatype = TUINT; break; case LONG_IMG: datatype = PS_TYPE_S32; fitsDatatype = TINT; break; case LONGLONG_IMG: datatype = PS_TYPE_S64; fitsDatatype = TLONGLONG; break; case FLOAT_IMG: datatype = PS_TYPE_F32; fitsDatatype = TFLOAT; break; case DOUBLE_IMG: datatype = PS_TYPE_F64; fitsDatatype = TDOUBLE; break; default: psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_FITS_TYPE_UNSUPPORTED, bitPix); psFree(output); return NULL; } output = psImageRecycle(output, lastPixel[0]-firstPixel[0]+1, lastPixel[1]-firstPixel[1]+1, datatype); if (output == NULL) { psError(PS_ERR_UNKNOWN, false, "Failed to allocate a properly sized image."); return false; } // n.b., this assumes contiguous image buffer if (fits_read_subset(fits->p_fd, fitsDatatype, firstPixel, lastPixel, increment, NULL, output->data.V[0], &anynull, &status) != 0) { psFree(output); (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_READ_FAILED, fitsErr); return NULL; } return output; } bool psFitsWriteImage(const psFits* fits, const psMetadata* header, const psImage* input, int numZPlanes, char* extname) { if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); return false; } if (input == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_IMAGE_NULL); return false; } int numCols = input->numCols; int numRows = input->numRows; int status = 0; // determine the FITS-equivalent parameters int bitPix; double bZero; int dataType; if (! convertPsTypeToFits(input->type.type, &bitPix, &bZero, &dataType) ) { return false; } int naxis = 3; long naxes[3]; naxes[0] = numCols; naxes[1] = numRows; naxes[2] = numZPlanes; if (numZPlanes < 2) { naxis = 2; } fits_create_img(fits->p_fd, bitPix, naxis, naxes, &status); if (extname != NULL) { fits_update_key_str(fits->p_fd, "EXTNAME", (char*)extname, NULL, &status); } if (bZero != 0) { // set the bscale/bzero fits_write_key_dbl(fits->p_fd, "BZERO", bZero, 12, "Pixel Value Offset", &status); fits_write_key_dbl(fits->p_fd, "BSCALE", 1.0, 12, "Pixel Value Scale", &status); fits_set_bscale(fits->p_fd, 1.0, bZero, &status); } // write the header, if any. if (header != NULL) { psFitsWriteHeader(header, fits); } if (input->parent == NULL) { // if no parent, assume that the image data is contiguous fits_write_img(fits->p_fd, dataType, // datatype 1, // writing to the first z-plane numCols*numRows, // number of elements to write, i.e., the whole image input->data.V[0], // the data &status); } else { // image data may not be contiguous; write one row at a time int firstPixel = 1; for (int row = 0; row < numRows; row++) { fits_write_img(fits->p_fd, dataType, // datatype firstPixel, numCols, // number of elements to write, i.e., one row's worth input->data.V[row],// the raw row data &status); firstPixel += numCols; // move to next row } } if ( status != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_WRITE_FAILED, fits->filename, fitsErr); return false; } return true; } bool psFitsUpdateImage(const psFits* fits, const psImage* input, psRegion region, int z) { int status = 0; if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); return false; } if (input == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_IMAGE_NULL); return false; } // check to see if we are positioned on an image HDU int hdutype; if ( fits_get_hdu_type(fits->p_fd,&hdutype, &status) != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED, fitsErr); return NULL; } if (hdutype != IMAGE_HDU) { psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_NOT_IMAGE_TYPE); return NULL; } int numCols = input->numCols; int numRows = input->numRows; // determine the FITS-equivalent parameters int bitPix; double bZero; int dataType; if (! convertPsTypeToFits(input->type.type, &bitPix, &bZero, &dataType) ) { return false; } //check to see if the HDU has the same datatype int fileBitpix; int naxis; long nAxes[3]; nAxes[2] = 1; fits_get_img_param(fits->p_fd, 3, &fileBitpix, &naxis, nAxes, &status); //check to see if the HDU has the same datatype if (bitPix != fileBitpix) { char* fitsTypeStr; char* imageTypeStr; PS_TYPE_NAME(fitsTypeStr,fileBitpix); PS_TYPE_NAME(imageTypeStr,input->type.type); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_IMAGE_UPDATE_TYPE_MISMATCH, fitsTypeStr, imageTypeStr); return false; } //check if the HDU has the z-plane requested if (z >= nAxes[2]) { psError(PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_psFits_FITS_Z_SMALL, nAxes[2],z); return false; } // determine the region in the FITS file domain long firstPixel[3]; long lastPixel[3]; firstPixel[0] = region.x0 + 1; firstPixel[1] = region.y0 + 1; firstPixel[2] = z + 1; if (region.x1 > 0) { lastPixel[0] = region.x1; } else { lastPixel[0] = nAxes[0] + region.x1; // n.b., region.x1 < 0 } if (region.y1 > 0) { lastPixel[1] = region.y1; } else { lastPixel[1] = nAxes[1] + region.y1; // n.b., region.y1 < 0 } lastPixel[2] = z + 1; if (firstPixel[0] < 1 || firstPixel[0] > nAxes[0] || firstPixel[1] < 1 || firstPixel[1] > nAxes[1] || lastPixel[0] < 1 || lastPixel[0] > nAxes[0] || lastPixel[1] < 1 || lastPixel[1] > nAxes[1]) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Specified region [%d:%d,%d:%d], is not valid given the %dx%d FITS image.", region.y0,region.y1-1,region.x0,region.x1-1); return false; } int dx = lastPixel[0] - firstPixel[0]; int dy = lastPixel[1] - firstPixel[1]; if (dx > numCols || dy > numRows) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, "The region [%d:%d,%d:%d], is not valid given the input %dx%d image.", firstPixel[1]-1,lastPixel[1]-1, firstPixel[0]-1,lastPixel[0]-1, numCols, numRows); return false; } psImage* subset; if (dx != numCols || dy != numRows) { // the input image needs to be subsetted subset = psImageSubset((psImage*)input,0,0,dx+1,dy+1); } else { subset = psMemIncrRefCounter((psImage*)input); } fits_write_subset(fits->p_fd, dataType, firstPixel, lastPixel, subset->data.V[0], &status); if ( status != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_WRITE_FAILED, fits->filename, fitsErr); return false; } return true; } bool psFitsWriteHeader(const psMetadata* header, const psFits* fits) { if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); return false; } if (header == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_METADATA_NULL); return false; } int status = 0; //transverse the metadata list and add each key. psListIterator* iter = psListIteratorAlloc(header->list,PS_LIST_HEAD,true); psMetadataItem* item; while ( (item=psListGetAndIncrement(iter)) != NULL ) { switch (item->type) { case PS_META_BOOL: { int value = item->data.B; fits_update_key(fits->p_fd, TLOGICAL, item->name, &value, item->comment, &status); break; } case PS_META_S32: fits_update_key(fits->p_fd, TINT, item->name, &item->data.S32, item->comment, &status); break; case PS_META_F32: fits_update_key(fits->p_fd, TFLOAT, item->name, &item->data.F32, item->comment, &status); break; case PS_META_F64: fits_update_key(fits->p_fd, TDOUBLE, item->name, &item->data.F64, item->comment, &status); break; case PS_META_STR: fits_update_key(fits->p_fd, TSTRING, item->name, item->data.V, item->comment, &status); break; default: // all other META types are ignored break; } if ( status != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_WRITE_FAILED, fits->filename, fitsErr); return false; } } return true; } psMetadata* psFitsReadTableRow(const psFits* fits, int row) { long numRows; int numCols; int status = 0; if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); return NULL; } // check to see if we even are positioned on a table HDU int hdutype; fits_get_hdu_type(fits->p_fd,&hdutype, &status); if ( status != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED, fitsErr); return NULL; } if (hdutype != ASCII_TBL && hdutype != BINARY_TBL) { psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_NOT_TABLE_TYPE); return NULL; } // get the size of the FITS table fits_get_num_rows(fits->p_fd, &numRows, &status); fits_get_num_cols(fits->p_fd, &numCols, &status); if ( status != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_GET_TABLE_SIZE_FAILED, fitsErr); return NULL; } // the row parameter in the proper range? if (row < 0 || row >= numRows) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psFits_ROW_INVALID, row,numRows); return NULL; } psMetadata* data = psMetadataAlloc(); int typecode; long repeat; long width; char name[60]; for (int col = 1; col <= numCols; col++) { // get the column name if (hdutype == BINARY_TBL) { fits_get_bcolparms(fits->p_fd, col, name, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &status); } else { fits_get_acolparms(fits->p_fd, col, name, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &status); } // get the column type fits_get_coltype(fits->p_fd, col, &typecode, &repeat, &width, &status); if (status == 0) { #define READ_TABLE_ROW_CASE(FITSTYPE, NATIVETYPE, TYPE) \ case FITSTYPE: { \ NATIVETYPE value = 0; \ int anynul = 0; \ fits_read_col(fits->p_fd, FITSTYPE, col,row+1, \ 1, 1, NULL, &value, &anynul, &status); \ psMetadataAdd(data,PS_LIST_TAIL, name, \ PS_META_##TYPE, \ "", (ps##TYPE)value); \ break; \ } switch (typecode) { case TBYTE: case TSHORT: case TLONGLONG: READ_TABLE_ROW_CASE(TLONG, long, S32) READ_TABLE_ROW_CASE(TFLOAT, float, F32) READ_TABLE_ROW_CASE(TDOUBLE, double, F64) READ_TABLE_ROW_CASE(TLOGICAL, bool, BOOL); case TSTRING: { char* value; int anynul = 0; fits_read_col(fits->p_fd, TSTRING, col,row+1, 1, 1, NULL, &value, &anynul, &status); if (anynul == 0) { psMetadataAdd(data,PS_LIST_TAIL, name, PS_META_STR, "", value); } break; } default: psTrace("psFits.psFitsReadTableRow", 2, "Column %d or row %d was of a non primitive type, %d", col, row, typecode); } } if ( status != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_GET_TABLE_ELEMENT, col,row,fitsErr); psFree(data); return NULL; } } return data; } psArray* psFitsReadTableColumn(const psFits* fits, const char* colname) { int colnum = 0; int status = 0; if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); return NULL; } // check to see if we even are positioned on a table HDU int hdutype; if ( fits_get_hdu_type(fits->p_fd,&hdutype, &status) != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED, fitsErr); return NULL; } if (hdutype != ASCII_TBL && hdutype != BINARY_TBL) { psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_NOT_TABLE_TYPE); return NULL; } // find the column by name if ( fits_get_colnum(fits->p_fd, CASESEN, (char*)colname, &colnum, &status) != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_FIND_COLUMN, colname, fitsErr); return NULL; } // get the number of rows long numRows = 0; fits_get_num_rows(fits->p_fd, &numRows, &status); // get the column length. int width; if ( fits_get_col_display_width(fits->p_fd, colnum, &width, &status) != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_GET_COLTYPE, fitsErr); return NULL; } // allocate the buffers psArray* result = psArrayAlloc(numRows); for (int row = 0; row < numRows; row++) { result->data[row] = psAlloc((width+1)*sizeof(char)); } result->n = numRows; fits_read_col_str(fits->p_fd, colnum, 1, // firstrow 1, // firestelem numRows, "", // nulstr (char**)result->data, NULL, &status); if ( status != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_TABLE_READ_COL, fitsErr); return NULL; } return result; } psVector* psFitsReadTableColumnNum(const psFits* fits, const char* colname) { int status = 0; int colnum = 0; if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); return NULL; } // check to see if we even are positioned on a table HDU int hdutype; if ( fits_get_hdu_type(fits->p_fd,&hdutype, &status) != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED, fitsErr); return NULL; } if (hdutype != ASCII_TBL && hdutype != BINARY_TBL) { psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_NOT_TABLE_TYPE); return NULL; } // find the column by name if ( fits_get_colnum(fits->p_fd, CASESEN, (char*)colname, &colnum, &status) != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_FIND_COLUMN, colname, fitsErr); return NULL; } // get the number of rows long numRows = 0; fits_get_num_rows(fits->p_fd, &numRows, &status); // get the column datatype. int typecode; long repeat; long width; if ( fits_get_eqcoltype(fits->p_fd, colnum, &typecode, &repeat, &width, &status) != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_GET_COLTYPE, fitsErr); return NULL; } psVector* result = psVectorAlloc(numRows, convertFitsToPsType(typecode)); fits_read_col(fits->p_fd, typecode, colnum, 1 /* firstrow */, 1 /* firstelem */, numRows, NULL, (psPtr)(result->data.U8), NULL, &status); if ( status != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_TABLE_READ_COL, fitsErr); return NULL; } return result; } psArray* psFitsReadTable(const psFits* fits) { int status = 0; if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); return NULL; } // check to see if we even are positioned on a table HDU int hdutype; if ( fits_get_hdu_type(fits->p_fd,&hdutype, &status) != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED, fitsErr); return NULL; } if (hdutype != ASCII_TBL && hdutype != BINARY_TBL) { psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_NOT_TABLE_TYPE); return NULL; } // get the size of the FITS table long numRows = 0; fits_get_num_rows(fits->p_fd, &numRows, &status); if ( status != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_GET_TABLE_SIZE_FAILED, fitsErr); return NULL; } psArray* table = psArrayAlloc(numRows); for (int row = 0; row < numRows; row++) { table->data[row] = psFitsReadTableRow(fits,row); } return table; } bool psFitsWriteTable(const psFits* fits, psMetadata* header, psArray* table, char* extname) { int status = 0; psMetadataItem* item; if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); return false; } if (table == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_IMAGE_NULL); return false; } int rows = table->n; if (rows < 1) { // no table data, what can I do? psError(PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_psFits_TABLE_EMPTY); return false; } // find all the columns needed psArray* columns = psArrayAlloc(((psMetadata*)table->data[0])->list->size); columns->n=0; // find the unique items in the array of metadata 'rows' for (int row=0; row < rows; row++) { psMetadata* rowMeta = table->data[row]; if (rowMeta != NULL) { psListIterator* iter = psListIteratorAlloc(rowMeta->list, PS_LIST_HEAD,true); while ( (item=psListGetAndIncrement(iter)) != NULL) { if (PS_META_IS_PRIMITIVE(item->type)) { bool found = false; for (int n=0; n < columns->n && ! found; n++) { if (strcmp(item->name, ((psMetadataItem*)(columns->data[n]))->name) == 0) { found = true; } } if (! found) { psArrayAdd(columns, columns->nalloc, item); } } } psFree(iter); } } if (columns->n == 0) { // no table columns found // XXX: Error? return false; } //create list of column names and types. psArray* columnNames = psArrayAlloc(columns->n); psArray* columnTypes = psArrayAlloc(columns->n); for (int n=0; n < columns->n; n++) { char* fitsType; columnNames->data[n] = psMemIncrRefCounter(((psMetadataItem*)columns->data[n])->name); if ( ! convertMetadataTypeToBinaryTForm(((psMetadataItem*)columns->data[n])->type, &fitsType)) { // XXX: error message return false; } columnTypes->data[n] = fitsType; } fits_create_tbl(fits->p_fd, BINARY_TBL, table->n, // number of rows in table columns->n, // number of columns in table (char**)columnNames->data, // names of the columns (char**)columnTypes->data, // format of the columns NULL, // physical unit of columns extname, // extension name &status); psFree(columnNames); psFree(columnTypes); // fill in the table elements with data for (int n = 0; n < columns->n; n++) { int row; item = columns->data[n]; if (PS_META_IS_PRIMITIVE(item->type)) { psVector* col = NULL; switch (item->type) { case PS_META_S32: col = psVectorAlloc(table->n, PS_TYPE_S32); for (row = 0; row < table->n; row++) { col->data.S32[row] = psMetadataLookupS32(NULL, table->data[row], item->name); } fits_write_col_int(fits->p_fd, n+1, // column number 1, // firstrow 1, // firstelem table->n, // nelements col->data.S32, &status); break; case PS_META_F32: col = psVectorAlloc(table->n, PS_TYPE_F32); for (row = 0; row < table->n; row++) { col->data.F32[row] = psMetadataLookupF32(NULL, table->data[row], item->name); } fits_write_col_flt(fits->p_fd, n+1, // column number 1, // firstrow 1, // firstelem table->n, // nelements col->data.F32, &status); break; case PS_META_F64: col = psVectorAlloc(table->n, PS_TYPE_F64); for (row = 0; row < table->n; row++) { col->data.F64[row] = psMetadataLookupF64(NULL, table->data[row], item->name); } fits_write_col_dbl(fits->p_fd, n+1, // column number 1, // firstrow 1, // firstelem table->n, // nelements col->data.F64, &status); break; case PS_META_BOOL: col = psVectorAlloc(table->n, PS_TYPE_BOOL); for (row = 0; row < table->n; row++) { col->data.S8[row] = psMetadataLookupBool(NULL, table->data[row], item->name); } fits_write_col_log(fits->p_fd, n+1, // column number 1, // firstrow 1, // firstelem table->n, // nelements col->data.S8, &status); break; default: // XXX: error message? break; } psFree(col); } else if (item->type == PS_META_STR) { psArray* col = psArrayAlloc(table->n); for (row = 0; row < table->n; row++) { col->data[row] = item->data.V; } fits_write_col_str(fits->p_fd, n, // column number 1, // firstrow 1, // firstelem table->n, // nelements (char**)col->data, &status); psFree(col); } } psFree(columns); return true; } bool psFitsUpdateTable(const psFits* fits, psMetadata* data, int row) { int status = 0; if (fits == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_NULL); return false; } if (data == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, PS_ERRORTEXT_psFits_IMAGE_NULL); return false; } // check to see if we even are positioned on a table HDU int hdutype; if ( fits_get_hdu_type(fits->p_fd,&hdutype, &status) != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED, fitsErr); return false; } if (hdutype != ASCII_TBL && hdutype != BINARY_TBL) { psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_NOT_TABLE_TYPE); return false; } psMetadataIterator* iter = psMetadataIteratorAlloc(data,PS_LIST_HEAD,NULL); psMetadataItem* item; while ( (item=psMetadataGetAndIncrement(iter)) != NULL) { if (PS_META_IS_PRIMITIVE(item->type)) { // operating on primitive data type, i.e., not a complex object int colnum = 0; if ( fits_get_colnum(fits->p_fd, CASESEN, item->name, &colnum, &status) == 0) { // cooresponding column found in table int dataType = 0; convertPsTypeToFits(item->type, NULL, NULL, &dataType); if (fits_write_col(fits->p_fd, dataType, colnum, row+1, 1, 1, &item->data,&status) != 0) { char fitsErr[MAX_STRING_LENGTH]; (void)fits_get_errstatus(status, fitsErr); psError(PS_ERR_IO, true, PS_ERRORTEXT_psFits_WRITE_FAILED, fits->filename, fitsErr); psFree(iter); return false; } } else { // the column was not found. psWarning("No column with the name '%s' exists in the table.", item->name); } } } psFree(iter); return true; }