Changeset 8412 for trunk/psLib/src/fits/psFitsImage.c
- Timestamp:
- Aug 17, 2006, 12:15:18 PM (20 years ago)
- File:
-
- 1 edited
-
trunk/psLib/src/fits/psFitsImage.c (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/fits/psFitsImage.c
r8232 r8412 7 7 * @author Robert DeSonia, MHPCC 8 8 * 9 * @version $Revision: 1.1 2$ $Name: not supported by cvs2svn $10 * @date $Date: 2006-08- 08 23:32:23$9 * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2006-08-17 22:15:17 $ 11 11 * 12 12 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 14 14 15 15 #include <unistd.h> 16 16 #include <assert.h> 17 #include <string.h> 18 19 #include "psAssert.h" 17 20 #include "psFits.h" 18 21 #include "psFitsImage.h" 19 #include "string.h"20 22 #include "psError.h" 21 23 … … 31 33 #define MAX_STRING_LENGTH 256 // maximum length string for FITS routines 32 34 33 psImage* psFitsReadImage(const psFits* fits, // the psFits object 35 // Information required to read a FITS file 36 typedef struct 37 { 38 int nAxis; // Number of axes 39 int bitPix; // Bits per pixel 40 long nAxes[3]; // Length of each axis 41 long firstPixel[3]; // lower-left corner of image subset 42 long lastPixel[3]; // upper-right corner of image subset 43 long increment[3]; // increment for image subset 44 int fitsDatatype; // cfitsio data type 45 int psDatatype; // psLib data type 46 } 47 p_psFitsReadInfo; 48 49 static p_psFitsReadInfo *p_psFitsReadInfoAlloc(const psFits *fits, // The FITS file handle 50 psRegion region, // Region to read 51 int z // z-plane to read in cube 52 ) 53 { 54 PS_ASSERT_PTR_NON_NULL(fits, NULL); 55 PS_ASSERT_INT_NONNEGATIVE(z, NULL); 56 57 p_psFitsReadInfo *info = psAlloc(sizeof(p_psFitsReadInfo)); 58 memset(info, 0, sizeof(p_psFitsReadInfo)); 59 60 int status = 0; // CFITSIO status 61 char fitsErr[80]; // CFITSIO error message string 62 63 // check to see if we even are positioned on an image HDU 64 int hduType; // Type of HDU 65 if (fits_get_hdu_type(fits->fd, &hduType, &status) != 0) { 66 (void)fits_get_errstatus(status, fitsErr); 67 psError(PS_ERR_IO, true, 68 _("Could not determine the HDU type. CFITSIO Error: %s"), 69 fitsErr); 70 goto bad; 71 } 72 if (hduType != IMAGE_HDU) { 73 psError(PS_ERR_IO, true, 74 _("Current FITS HDU type must be an image.")); 75 goto bad; 76 } 77 78 // Get the data type 'bitPix' from the FITS image 79 if (fits_get_img_equivtype(fits->fd, &info->bitPix, &status) != 0) { 80 fits_get_errstatus(status, fitsErr); 81 psError(PS_ERR_IO, true, 82 _("Could not determine image data type. CFITSIO Error: %s"), 83 fitsErr); 84 goto bad; 85 } 86 87 /* Get the dimensions 'nAxis' from the FITS image */ 88 if (fits_get_img_dim(fits->fd, &info->nAxis, &status) != 0) { 89 (void)fits_get_errstatus(status, fitsErr); 90 psError(PS_ERR_IO, true, 91 _("Could not determine image dimensions. CFITSIO Error: %s"), 92 fitsErr); 93 goto bad; 94 } 95 96 /* Validate the number of axis */ 97 if ((info->nAxis < 2) || (info->nAxis > 3)) { 98 psError(PS_ERR_IO, true, 99 _("Image number of dimensions, %d, is not valid. Only two or three dimensions supported for FITS I/O."), 100 info->nAxis); 101 goto bad; 102 } 103 104 /* Get the Image size from the FITS file */ 105 if (fits_get_img_size(fits->fd, info->nAxis, info->nAxes, &status) != 0) { 106 (void)fits_get_errstatus(status, fitsErr); 107 psError(PS_ERR_IO, true, 108 _("Could not determine image size. CFITSIO Error: %s"), 109 fitsErr); 110 goto bad; 111 } 112 113 info->firstPixel[0] = region.x0 + 1; 114 info->firstPixel[1] = region.y0 + 1; 115 info->firstPixel[2] = z + 1; 116 117 if (region.x1 > 0) { 118 info->lastPixel[0] = region.x1; 119 } else { 120 info->lastPixel[0] = info->nAxes[0] + region.x1; // n.b., region.x1 < 0 121 } 122 if (region.y1 > 0) { 123 info->lastPixel[1] = region.y1; 124 } else { 125 info->lastPixel[1] = info->nAxes[1] + region.y1; // n.b., region.y1 < 0 126 } 127 info->lastPixel[2] = z + 1; 128 129 info->increment[0] = 1; 130 info->increment[1] = 1; 131 info->increment[2] = 1; 132 133 switch (info->bitPix) { 134 case BYTE_IMG: 135 info->psDatatype = PS_TYPE_U8; 136 info->fitsDatatype = TBYTE; 137 break; 138 case SBYTE_IMG: 139 info->psDatatype = PS_TYPE_S8; 140 info->fitsDatatype = TSBYTE; 141 break; 142 case USHORT_IMG: 143 info->psDatatype = PS_TYPE_U16; 144 info->fitsDatatype = TUSHORT; 145 break; 146 case SHORT_IMG: 147 info->psDatatype = PS_TYPE_S16; 148 info->fitsDatatype = TSHORT; 149 break; 150 case ULONG_IMG: 151 info->psDatatype = PS_TYPE_U32; 152 info->fitsDatatype = TUINT; 153 break; 154 case LONG_IMG: 155 info->psDatatype = PS_TYPE_S32; 156 info->fitsDatatype = TINT; 157 break; 158 case LONGLONG_IMG: 159 info->psDatatype = PS_TYPE_S64; 160 info->fitsDatatype = TLONGLONG; 161 break; 162 case FLOAT_IMG: 163 info->psDatatype = PS_TYPE_F32; 164 info->fitsDatatype = TFLOAT; 165 break; 166 case DOUBLE_IMG: 167 info->psDatatype = PS_TYPE_F64; 168 info->fitsDatatype = TDOUBLE; 169 break; 170 default: 171 psError(PS_ERR_IO, true, 172 _("FITS image type, BITPIX=%d, is not supported."), 173 info->bitPix); 174 goto bad; 175 } 176 177 return info; 178 179 // Common error cleanup 180 bad: 181 psFree(info); 182 return NULL; 183 } 184 185 // Read into an extant image of just the right size 186 static bool fitsReadImage(psImage *output, // Output image 187 const psFits *fits, // FITS file handle 188 p_psFitsReadInfo *info // Info on how to read 189 ) 190 { 191 // n.b., this assumes contiguous image buffer 192 assert(output); 193 assert(fits); 194 assert(info); 195 assert(output->numCols == info->lastPixel[0] - info->firstPixel[0] + 1); 196 assert(output->numRows == info->lastPixel[1] - info->firstPixel[1] + 1); // Right size 197 assert(!output->parent); // No parents means the buffer is contiguous 198 199 int anynull = 0; // Are there any NULLs in the data? 200 int status = 0; // cfitsio status 201 if (fits_read_subset(fits->fd, info->fitsDatatype, info->firstPixel, info->lastPixel, 202 info->increment, NULL, output->data.V[0], &anynull, &status) != 0) { 203 char fitsErr[80]; // CFITSIO error message string 204 (void)fits_get_errstatus(status, fitsErr); 205 psError(PS_ERR_IO, true, 206 _("Reading FITS file failed. CFITSIO Error: %s"), 207 fitsErr); 208 return false; 209 } 210 211 return true; 212 } 213 214 psImage* psFitsReadImage(const psFits *fits, // the psFits object 34 215 psRegion region, // the region in the FITS image to read 35 int z) // the z-plane in the FITS image cube to read 36 { 37 psS32 status = 0; /* CFITSIO file vars */ 38 psS32 nAxis = 0; 39 psS32 anynull = 0; 40 psS32 bitPix = 0; /* Pixel type */ 41 long nAxes[3]; 42 long firstPixel[3]; /* lower-left corner of image subset */ 43 long lastPixel[3]; /* upper-right corner of image subset */ 44 long increment[3]; /* increment for image subset */ 45 char fitsErr[80] = ""; /* CFITSIO error message string */ 46 psS32 fitsDatatype = 0; 47 psS32 datatype = 0; 48 49 if (fits == NULL) { 50 psError(PS_ERR_BAD_PARAMETER_NULL, true, 51 _("The input psFits object can not NULL.")); 216 int z // the z-plane in the FITS image cube to read 217 ) 218 { 219 p_psFitsReadInfo *info = p_psFitsReadInfoAlloc(fits, region, z); 220 221 psImage *output = psImageAlloc(info->lastPixel[0] - info->firstPixel[0] + 1, 222 info->lastPixel[1] - info->firstPixel[1] + 1, 223 info->psDatatype); 224 225 if (!fitsReadImage(output, fits, info)) { 226 psFree(info); 227 psFree(output); 52 228 return NULL; 53 229 } 54 230 55 // check to see if we even are positioned on an image HDU 56 int hdutype; 57 if ( fits_get_hdu_type(fits->fd,&hdutype, &status) != 0) { 58 char fitsErr[MAX_STRING_LENGTH]; 59 (void)fits_get_errstatus(status, fitsErr); 60 psError(PS_ERR_IO, true, 61 _("Could not determine the HDU type. CFITSIO Error: %s"), 62 fitsErr); 231 psFree(info); 232 return output; 233 } 234 235 psImage* psFitsReadImageBuffer(psImage *output, // Output image buffer 236 const psFits *fits, // the psFits object 237 psRegion region, // the region in the FITS image to read 238 int z // the z-plane in the FITS image cube to read 239 ) 240 { 241 if (output && output->parent) { 242 psError(PS_ERR_IO, true, "Unable to read into a buffer for a child image.\n"); 63 243 return NULL; 64 244 } 65 if (hdutype != IMAGE_HDU) { 66 psError(PS_ERR_IO, true, 67 _("Current FITS HDU type must be an image.")); 245 246 p_psFitsReadInfo *info = p_psFitsReadInfoAlloc(fits, region, z); 247 248 output = psImageRecycle(output, info->lastPixel[0] - info->firstPixel[0] + 1, 249 info->lastPixel[1] - info->firstPixel[1] + 1, 250 info->psDatatype); 251 252 if (!fitsReadImage(output, fits, info)) { 253 psFree(info); 254 psFree(output); 68 255 return NULL; 69 256 } 70 257 71 /* Get the data type 'bitPix' from the FITS image */ 72 if (fits_get_img_equivtype(fits->fd, &bitPix, &status) != 0) { 73 fits_get_errstatus(status, fitsErr); 74 psError(PS_ERR_IO, true, 75 _("Could not determine image data type. CFITSIO Error: %s"), 76 fitsErr); 77 return NULL; 78 } 79 80 /* Get the dimensions 'nAxis' from the FITS image */ 81 if (fits_get_img_dim(fits->fd, &nAxis, &status) != 0) { 82 (void)fits_get_errstatus(status, fitsErr); 83 psError(PS_ERR_IO, true, 84 _("Could not determine image dimensions. CFITSIO Error: %s"), 85 fitsErr); 86 return NULL; 87 } 88 89 /* Validate the number of axis */ 90 if ((nAxis < 2) || (nAxis > 3)) { 91 psError(PS_ERR_IO, true, 92 _("Image number of dimensions, %d, is not valid. Only two or three dimensions supported for FITS I/O."), 93 nAxis); 94 return NULL; 95 } 96 97 /* Get the Image size from the FITS file */ 98 if (fits_get_img_size(fits->fd, nAxis, nAxes, &status) != 0) { 99 (void)fits_get_errstatus(status, fitsErr); 100 psError(PS_ERR_IO, true, 101 _("Could not determine image size. CFITSIO Error: %s"), 102 fitsErr); 103 return NULL; 104 } 105 106 firstPixel[0] = region.x0 + 1; 107 firstPixel[1] = region.y0 + 1; 108 firstPixel[2] = z + 1; 109 110 if (region.x1 > 0) { 111 lastPixel[0] = region.x1; 112 } else { 113 lastPixel[0] = nAxes[0] + region.x1; // n.b., region.x1 < 0 114 } 115 if (region.y1 > 0) { 116 lastPixel[1] = region.y1; 117 } else { 118 lastPixel[1] = nAxes[1] + region.y1; // n.b., region.y1 < 0 119 } 120 lastPixel[2] = z + 1; 121 122 increment[0] = 1; 123 increment[1] = 1; 124 increment[2] = 1; 125 126 switch (bitPix) { 127 case BYTE_IMG: 128 datatype = PS_TYPE_U8; 129 fitsDatatype = TBYTE; 130 break; 131 case SBYTE_IMG: 132 datatype = PS_TYPE_S8; 133 fitsDatatype = TSBYTE; 134 break; 135 case USHORT_IMG: 136 datatype = PS_TYPE_U16; 137 fitsDatatype = TUSHORT; 138 break; 139 case SHORT_IMG: 140 datatype = PS_TYPE_S16; 141 fitsDatatype = TSHORT; 142 break; 143 case ULONG_IMG: 144 datatype = PS_TYPE_U32; 145 fitsDatatype = TUINT; 146 break; 147 case LONG_IMG: 148 datatype = PS_TYPE_S32; 149 fitsDatatype = TINT; 150 break; 151 case LONGLONG_IMG: 152 datatype = PS_TYPE_S64; 153 fitsDatatype = TLONGLONG; 154 break; 155 case FLOAT_IMG: 156 datatype = PS_TYPE_F32; 157 fitsDatatype = TFLOAT; 158 break; 159 case DOUBLE_IMG: 160 datatype = PS_TYPE_F64; 161 fitsDatatype = TDOUBLE; 162 break; 163 default: 164 psError(PS_ERR_IO, true, 165 _("FITS image type, BITPIX=%d, is not supported."), 166 bitPix); 167 return NULL; 168 } 169 170 psImage *output = psImageAlloc(lastPixel[0]-firstPixel[0]+1, 171 lastPixel[1]-firstPixel[1]+1, 172 datatype); 173 174 // n.b., this assumes contiguous image buffer 175 if (fits_read_subset(fits->fd, fitsDatatype, firstPixel, lastPixel, increment, 176 NULL, output->data.V[0], &anynull, &status) != 0) { 177 psFree(output); 178 (void)fits_get_errstatus(status, fitsErr); 179 psError(PS_ERR_IO, true, 180 _("Reading FITS file failed. CFITSIO Error: %s"), 181 fitsErr); 182 return NULL; 183 } 184 258 psFree(info); 185 259 return output; 186 187 260 } 188 261 … … 329 402 330 403 // check to see if we are positioned on an image HDU 331 int hdu type;332 if ( fits_get_hdu_type(fits->fd, &hdutype, &status) != 0) {404 int hduType; 405 if ( fits_get_hdu_type(fits->fd, &hduType, &status) != 0) { 333 406 char fitsErr[MAX_STRING_LENGTH]; 334 407 (void)fits_get_errstatus(status, fitsErr); … … 338 411 return NULL; 339 412 } 340 if (hdu type != IMAGE_HDU) {413 if (hduType != IMAGE_HDU) { 341 414 psError(PS_ERR_IO, true, 342 415 _("Current FITS HDU type must be an image."));
Note:
See TracChangeset
for help on using the changeset viewer.
