Changeset 13526 for trunk/psModules/src/objects/pmPSF_IO.c
- Timestamp:
- May 25, 2007, 4:52:27 PM (19 years ago)
- File:
-
- 1 edited
-
trunk/psModules/src/objects/pmPSF_IO.c (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/objects/pmPSF_IO.c
r13414 r13526 6 6 * @author EAM, IfA 7 7 * 8 * @version $Revision: 1.1 8$ $Name: not supported by cvs2svn $9 * @date $Date: 2007-05- 18 00:17:12$8 * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2007-05-26 02:52:27 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 20 20 /* INCLUDE FILES */ 21 21 /*****************************************************************************/ 22 23 #include <string.h> 24 #include <strings.h> /* for strn?casecmp */ 22 25 23 26 #include <pslib.h> … … 212 215 } 213 216 214 // for each Readout (ie, analysed image), we write out: header + table with PSF model parameters, 215 // and header + image for the PSF residual images 217 // for each Readout (ie, analysed image), we write out 218 // - image header : FITS Image NAXIS = 0 219 // - psf table (+header) : FITS Table 220 // - psf resid (+header) : FITS Image 221 // if needed, we also write out a PHU blank header 216 222 bool pmReadoutWritePSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) 217 223 { 218 224 bool status; 225 pmHDU *hdu; 226 pmHDU *phu; 227 char *headName, *tableName, *residName; 228 229 // write a PHU? (only if input image is MEF) 230 // write a header? (only if this is the first readout for cell) 231 // note that the file->header is set to track the last hdu->header written 232 // write the data? (always?) 233 234 // get the current header 235 hdu = pmFPAviewThisHDU (view, file->fpa); 236 237 // if file does not yet have a PHU, attempt to write one to disk 238 // we only need a PHU if chips->n > 1 and file->fileLevel == FPA 239 // otherwise, the chip header fills the PHU location 240 // XXX this code could be placed in a 'pmPSF_WritePHU' function and called 241 // from pmFPAfileIO.c. 242 if ((file->fileLevel == PM_FPA_LEVEL_FPA) && (file->fpa->chips->n > 1) && !file->phu) { 243 244 // find the FPA phu 245 phu = pmFPAviewThisPHU (view, file->fpa); 246 247 // if there is no PHU, this is a single header+image (extension-less) file 248 // if there is a PHU, write it out as a 'blank' 249 psMetadata *outhead = psMetadataAlloc(); 250 if (phu) { 251 psMetadataCopy (outhead, phu->header); 252 } else { 253 pmConfigConformHeader (outhead, file->format); 254 255 psMetadata *fileData = psMetadataLookupMetadata(NULL, file->format, "FILE"); // File information 256 const char *fpaNameHdr = psMetadataLookupStr(NULL, fileData, "FPA.NAME"); 257 if (fpaNameHdr && strlen(fpaNameHdr) > 0) { 258 const char *fpaName = psMetadataLookupStr(NULL, file->fpa->concepts, "FPA.NAME"); 259 psMetadataAddStr(outhead, PS_LIST_TAIL, fpaNameHdr, PS_META_REPLACE, "FPA name", fpaName); 260 } 261 } 262 263 psMetadataAddBool (outhead, PS_LIST_TAIL, "EXTEND", PS_META_REPLACE, "this file has extensions", true); 264 psFitsWriteBlank (file->fits, outhead, ""); 265 file->phu = true; 266 psTrace ("pmFPAfile", 5, "wrote phu %s (type: %d)\n", file->filename, file->type); 267 psFree (outhead); 268 } 269 270 // define the EXTNAME values used for image header, table data, and residual image segments 271 { 272 // lookup the EXTNAME values used for table data and image header segments 273 char *rule = NULL; 274 275 // Menu of EXTNAME rules 276 psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES"); 277 if (!menu) { 278 psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config"); 279 return false; 280 } 281 282 // EXTNAME for image header 283 rule = psMetadataLookupStr(&status, menu, "PSF.HEAD"); 284 if (!rule) { 285 psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.HEAD in EXTNAME.RULES in camera.config"); 286 return false; 287 } 288 headName = pmFPAfileNameFromRule (rule, file, view); 289 290 // EXTNAME for table data 291 rule = psMetadataLookupStr(&status, menu, "PSF.TABLE"); 292 if (!rule) { 293 psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config"); 294 psFree (headName); 295 return false; 296 } 297 tableName = pmFPAfileNameFromRule (rule, file, view); 298 299 // EXTNAME for resid data 300 rule = psMetadataLookupStr(&status, menu, "PSF.RESID"); 301 if (!rule) { 302 psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config"); 303 psFree (headName); 304 psFree (tableName); 305 return false; 306 } 307 residName = pmFPAfileNameFromRule (rule, file, view); 308 } 309 310 // write out the IMAGE header segment 311 // this header block is new, write it to disk 312 if (hdu->header != file->header) { 313 // add EXTNAME, EXTHEAD, EXTTYPE to header 314 psMetadataAddStr (hdu->header, PS_LIST_TAIL, "EXTTABLE", PS_META_REPLACE, "name of table extension", tableName); 315 psMetadataAddStr (hdu->header, PS_LIST_TAIL, "EXTRESID", PS_META_REPLACE, "name of resid extension", residName); 316 psMetadataAddStr (hdu->header, PS_LIST_TAIL, "EXTTYPE", PS_META_REPLACE, "extension type", "IMAGE"); 317 if (!file->phu) { 318 // this hdu->header acts as the PHU: set EXTEND to be true 319 psMetadataAddBool (hdu->header, PS_LIST_TAIL, "EXTEND", PS_META_REPLACE, "this file has extensions", true); 320 } 321 322 psFitsWriteBlank (file->fits, hdu->header, headName); 323 psTrace ("pmFPAfile", 5, "wrote ext head %s (type: %d)\n", file->filename, file->type); 324 file->header = hdu->header; 325 psFree (headName); 326 } 219 327 220 328 // select the psf of interest 221 329 pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); 222 330 if (!psf) { 223 psError(PS_ERR_UNKNOWN, true, "missing PSF for this readout"); 224 return false; 225 } 226 227 // lookup the EXTNAME values used for table data and residual image segments 228 char *rule = NULL; 229 230 // Menu of EXTNAME rules 231 psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES"); 232 if (!menu) { 233 psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config"); 234 return false; 235 } 236 // EXTNAME for table data 237 rule = psMetadataLookupStr(&status, menu, "PSF.TABLE"); 238 if (!rule) { 239 psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config"); 240 return false; 241 } 242 char *tableName = pmFPAfileNameFromRule (rule, file, view); 331 psError(PS_ERR_UNKNOWN, true, "missing PSF for this readout"); 332 psFree (tableName); 333 psFree (residName); 334 return false; 335 } 243 336 244 337 // write the PSF model parameters in a FITS table 245 246 // we need to write a header for the table, 247 psMetadata *header = psMetadataAlloc(); 248 249 char *modelName = pmModelGetType (psf->type); 250 psMetadataAddStr (header, PS_LIST_TAIL, "PSF_NAME", 0, "PSF model name", modelName); 251 252 psMetadataAddBool (header, PS_LIST_TAIL, "POISSON", 0, "Use Poisson errors in fits?", psf->poissonErrors); 253 254 int nPar = pmModelParameterCount (psf->type) ; 255 psMetadataAdd (header, PS_LIST_TAIL, "PSF_NPAR", PS_DATA_S32, "PSF model parameter count", nPar); 256 257 // save the dimensions of each parameter 258 for (int i = 0; i < nPar; i++) { 259 char name[9]; 260 psPolynomial2D *poly = psf->params_NEW->data[i]; 261 if (poly == NULL) continue; 262 snprintf (name, 9, "PAR%02d_NX", i); 263 psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nX); 264 snprintf (name, 9, "PAR%02d_NY", i); 265 psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nY); 266 } 267 268 // other required information describing the PSF 269 psMetadataAddF32 (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid); 270 psMetadataAddF32 (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid); 271 psMetadataAddF32 (header, PS_LIST_TAIL, "CHISQ", 0, "chi-square for fit", psf->chisq); 272 psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS", 0, "number of stars used to measure PSF", psf->nPSFstars); 273 274 // XXX can we drop this now? 275 psMetadataAddF32 (header, PS_LIST_TAIL, "SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias); 276 277 // build a FITS table of the PSF parameters 278 psArray *psfTable = psArrayAllocEmpty (100); 279 for (int i = 0; i < nPar; i++) { 280 psPolynomial2D *poly = psf->params_NEW->data[i]; 281 if (poly == NULL) continue; // skip unset parameters (eg, XPOS) 282 for (int ix = 0; ix < poly->nX; ix++) { 283 for (int iy = 0; iy < poly->nY; iy++) { 284 285 psMetadata *row = psMetadataAlloc (); 286 psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i); 287 psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER", 0, "", ix); 288 psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER", 0, "", iy); 289 psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE", 0, "", poly->coeff[ix][iy]); 290 psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR", 0, "", poly->coeffErr[ix][iy]); 291 psMetadataAddU8 (row, PS_LIST_TAIL, "MASK", 0, "", poly->mask[ix][iy]); 292 293 psArrayAdd (psfTable, 100, row); 294 psFree (row); 295 } 296 } 297 } 298 299 // write an empty FITS segment if we have no PSF information 300 if (psfTable->n == 0) { 301 // XXX this is probably an error (if we have a PSF, how do we have no data?) 302 psFitsWriteBlank (file->fits, header, tableName); 303 } else { 304 psTrace ("pmFPAfile", 5, "writing psf data %s\n", tableName); 305 if (!psFitsWriteTable (file->fits, header, psfTable, tableName)) { 306 psError(PS_ERR_IO, false, "writing psf table data %s\n", tableName); 307 psFree (tableName); 308 psFree (psfTable); 309 psFree (header); 310 return false; 311 } 312 } 313 psFree (tableName); 314 psFree (psfTable); 315 psFree (header); 316 317 // EXTNAME for residual images 318 rule = psMetadataLookupStr(&status, menu, "PSF.RESID"); 319 if (!rule) { 320 psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config"); 321 return false; 322 } 323 char *imageName = pmFPAfileNameFromRule (rule, file, view); 338 { 339 // we need to write a header for the table, 340 psMetadata *header = psMetadataAlloc(); 341 342 char *modelName = pmModelGetType (psf->type); 343 psMetadataAddStr (header, PS_LIST_TAIL, "PSF_NAME", 0, "PSF model name", modelName); 344 345 psMetadataAddBool (header, PS_LIST_TAIL, "POISSON", 0, "Use Poisson errors in fits?", psf->poissonErrors); 346 347 int nPar = pmModelParameterCount (psf->type) ; 348 psMetadataAdd (header, PS_LIST_TAIL, "PSF_NPAR", PS_DATA_S32, "PSF model parameter count", nPar); 349 350 // save the dimensions of each parameter 351 for (int i = 0; i < nPar; i++) { 352 char name[9]; 353 psPolynomial2D *poly = psf->params_NEW->data[i]; 354 if (poly == NULL) continue; 355 snprintf (name, 9, "PAR%02d_NX", i); 356 psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nX); 357 snprintf (name, 9, "PAR%02d_NY", i); 358 psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nY); 359 } 360 361 // other required information describing the PSF 362 psMetadataAddF32 (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid); 363 psMetadataAddF32 (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid); 364 psMetadataAddF32 (header, PS_LIST_TAIL, "CHISQ", 0, "chi-square for fit", psf->chisq); 365 psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS", 0, "number of stars used to measure PSF", psf->nPSFstars); 366 367 // XXX can we drop this now? 368 psMetadataAddF32 (header, PS_LIST_TAIL, "SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias); 369 370 // build a FITS table of the PSF parameters 371 psArray *psfTable = psArrayAllocEmpty (100); 372 for (int i = 0; i < nPar; i++) { 373 psPolynomial2D *poly = psf->params_NEW->data[i]; 374 if (poly == NULL) continue; // skip unset parameters (eg, XPOS) 375 for (int ix = 0; ix <= poly->nX; ix++) { 376 for (int iy = 0; iy <= poly->nY; iy++) { 377 378 psMetadata *row = psMetadataAlloc (); 379 psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i); 380 psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER", 0, "", ix); 381 psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER", 0, "", iy); 382 psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE", 0, "", poly->coeff[ix][iy]); 383 psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR", 0, "", poly->coeffErr[ix][iy]); 384 psMetadataAddU8 (row, PS_LIST_TAIL, "MASK", 0, "", poly->mask[ix][iy]); 385 386 psArrayAdd (psfTable, 100, row); 387 psFree (row); 388 } 389 } 390 } 391 392 // write an empty FITS segment if we have no PSF information 393 if (psfTable->n == 0) { 394 // XXX this is probably an error (if we have a PSF, how do we have no data?) 395 psFitsWriteBlank (file->fits, header, tableName); 396 } else { 397 psTrace ("pmFPAfile", 5, "writing psf data %s\n", tableName); 398 if (!psFitsWriteTable (file->fits, header, psfTable, tableName)) { 399 psError(PS_ERR_IO, false, "writing psf table data %s\n", tableName); 400 psFree (tableName); 401 psFree (residName); 402 psFree (psfTable); 403 psFree (header); 404 return false; 405 } 406 } 407 psFree (tableName); 408 psFree (psfTable); 409 psFree (header); 410 } 324 411 325 412 // write the residual images (3D) 326 header = psMetadataAlloc (); 327 if (psf->residuals == NULL) { 328 // set some header keywords to make it clear there are no residuals? 329 psFitsWriteBlank (file->fits, header, imageName); 330 psFree (imageName); 331 psFree (header); 332 return true; 333 } 334 335 psMetadataAddS32 (header, PS_LIST_TAIL, "XBIN", 0, "", psf->residuals->xBin); 336 psMetadataAddS32 (header, PS_LIST_TAIL, "YBIN", 0, "", psf->residuals->yBin); 337 psMetadataAddS32 (header, PS_LIST_TAIL, "XCENTER", 0, "", psf->residuals->xCenter); 338 psMetadataAddS32 (header, PS_LIST_TAIL, "YCENTER", 0, "", psf->residuals->yCenter); 339 340 // write the residuals as three planes of the image 341 // this call creates an extension with NAXIS3 = 3 342 if (psf->residuals->Rx) { 343 // this call creates an extension with NAXIS3 = 3 344 psArray *images = psArrayAllocEmpty (3); 345 psArrayAdd (images, 1, psf->residuals->Ro); 346 psArrayAdd (images, 1, psf->residuals->Rx); 347 psArrayAdd (images, 1, psf->residuals->Ry); 348 349 psFitsWriteImageCube (file->fits, header, images, imageName); 350 psFree (images); 351 } else { 352 // this call creates an extension with NAXIS3 = 1 353 psFitsWriteImage(file->fits, header, psf->residuals->Ro, 0, imageName); 354 } 355 356 psFree (imageName); 357 psFree (header); 413 { 414 psMetadata *header = psMetadataAlloc (); 415 if (psf->residuals == NULL) { 416 // set some header keywords to make it clear there are no residuals? 417 psFitsWriteBlank (file->fits, header, residName); 418 psFree (residName); 419 psFree (header); 420 return true; 421 } 422 423 psMetadataAddS32 (header, PS_LIST_TAIL, "XBIN", 0, "", psf->residuals->xBin); 424 psMetadataAddS32 (header, PS_LIST_TAIL, "YBIN", 0, "", psf->residuals->yBin); 425 psMetadataAddS32 (header, PS_LIST_TAIL, "XCENTER", 0, "", psf->residuals->xCenter); 426 psMetadataAddS32 (header, PS_LIST_TAIL, "YCENTER", 0, "", psf->residuals->yCenter); 427 428 // write the residuals as three planes of the image 429 // this call creates an extension with NAXIS3 = 3 430 if (psf->residuals->Rx) { 431 // this call creates an extension with NAXIS3 = 3 432 psArray *images = psArrayAllocEmpty (3); 433 psArrayAdd (images, 1, psf->residuals->Ro); 434 psArrayAdd (images, 1, psf->residuals->Rx); 435 psArrayAdd (images, 1, psf->residuals->Ry); 436 437 psFitsWriteImageCube (file->fits, header, images, residName); 438 psFree (images); 439 } else { 440 // this call creates an extension with NAXIS3 = 1 441 psFitsWriteImage(file->fits, header, psf->residuals->Ro, 0, residName); 442 } 443 psFree (residName); 444 psFree (header); 445 } 446 358 447 return true; 359 448 … … 458 547 psMetadata *header = NULL; 459 548 549 psTrace ("psModules.objects", 5, "read psf model for %s\n", file->filename); 550 460 551 // Menu of EXTNAME rules 461 552 psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES"); … … 493 584 char *modelName = psMetadataLookupStr (&status, header, "PSF_NAME"); 494 585 pmModelType type = pmModelSetType (modelName); 586 if (type == -1) { 587 psError(PS_ERR_UNKNOWN, true, "invalid model name %s in psf file %s", modelName, file->filename); 588 return false; 589 } 495 590 496 591 bool poissonErrors = psMetadataLookupBool (&status, header, "POISSON"); … … 509 604 snprintf (name, 9, "PAR%02d_NX", i); 510 605 int nXorder = psMetadataLookupS32 (&status, header, name); 606 if (!status) continue; // not all parameters are defined 511 607 snprintf (name, 9, "PAR%02d_NY", i); 512 608 int nYorder = psMetadataLookupS32 (&status, header, name); 609 if (!status) { 610 psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX defined for PAR %d, but no NY", i); 611 return false; 612 } 513 613 psf->params_NEW->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder); 514 614 } … … 527 627 528 628 // fill in the matching psf->params entries 529 for (int i = 0; i >table->n; i++) {629 for (int i = 0; i < table->n; i++) { 530 630 psMetadata *row = table->data[i]; 531 631 int iPar = psMetadataLookupS32 (&status, row, "MODEL_TERM"); … … 535 635 536 636 psPolynomial2D *poly = psf->params_NEW->data[iPar]; 637 if (poly == NULL) { 638 psError(PS_ERR_UNKNOWN, true, "values for parameter %d, but missing NX", iPar); 639 return false; 640 } 537 641 538 642 poly->coeff[xPow][yPow] = psMetadataLookupF32 (&status, row, "VALUE"); … … 583 687 psFree (tableName); 584 688 psFree (imageName); 689 psFree (header); 585 690 586 691 return true;
Note:
See TracChangeset
for help on using the changeset viewer.
