IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 24, 2007, 11:27:58 AM (19 years ago)
Author:
eugene
Message:

update from eam_branch_20070921

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects/pmPSF_IO.c

    r14936 r15000  
    66 *  @author EAM, IfA
    77 *
    8  *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2007-09-21 00:05:57 $
     8 *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2007-09-24 21:27:58 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    148148
    149149// for a pmPSF supplied on the analysis metadata, we write out
     150// if needed:
     151//   - a PHU blank header
    150152// - image header        : FITS Image NAXIS = 0
    151 // - psf table (+header) : FITS Table
    152 // - psf resid (+header) : FITS Image
    153 // if needed, we also write out a PHU blank header
     153// if (trendMode == MAP)
     154//   - psf resid (+header) : FITS Image
     155// else
     156//   - psf table (+header) : FITS Table
    154157bool pmPSFmodelWrite (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
    155158{
    … …  
    159162
    160163    if (!analysis) return false;
     164
     165    // select the current recipe
     166    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSPHOT");
     167    if (!recipe) {
     168        psError(PS_ERR_UNKNOWN, false, "missing recipe %s", "PSPHOT");
     169        return false;
     170    }
    161171
    162172    // write a PHU? (only if input image is MEF)
    … …  
    215225
    216226    // write out the IMAGE header segment
    217     // this header block is new, write it to disk
     227    // if this header block is new, write it to disk
    218228    if (hdu->header != file->header) {
    219229        // add EXTNAME, EXTHEAD, EXTTYPE to header
    … …  
    250260        psMetadataAddStr (header, PS_LIST_TAIL, "PSF_NAME", 0, "PSF model name", modelName);
    251261
    252         psMetadataAddBool (header, PS_LIST_TAIL, "POISSON",  0, "Use Poisson errors in fits?", psf->poissonErrors);
     262        psMetadataAddBool (header, PS_LIST_TAIL, "ERR_LMM",  0, "Use Poisson errors in fits?", psf->poissonErrorsPhotLMM);
     263        psMetadataAddBool (header, PS_LIST_TAIL, "ERR_LIN",  0, "Use Poisson errors in fits?", psf->poissonErrorsPhotLin);
     264        psMetadataAddBool (header, PS_LIST_TAIL, "ERR_PAR",  0, "Use Poisson errors in fits?", psf->poissonErrorsParams);
    253265
    254266        int nPar = pmModelClassParameterCount (psf->type)    ;
    255267        psMetadataAdd (header, PS_LIST_TAIL, "PSF_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
     268
     269        psMetadataAddS32 (header, PS_LIST_TAIL, "IMAXIS1", 0, "Image X Size", psf->fieldNx);
     270        psMetadataAddS32 (header, PS_LIST_TAIL, "IMAXIS2", 0, "Image Y Size", psf->fieldNy);
     271        psMetadataAddS32 (header, PS_LIST_TAIL, "IMREF1",  0, "Image X Ref",  psf->fieldXo);
     272        psMetadataAddS32 (header, PS_LIST_TAIL, "IMREF2",  0, "Image Y Ref",  psf->fieldYo);
     273
     274        // extract PSF Clump info
     275        pmPSFClump psfClump;
     276
     277        psfClump.X  = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.X");   assert (status);
     278        psfClump.Y  = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.Y");   assert (status);
     279        psfClump.dX = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DX");  assert (status);
     280        psfClump.dY = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DY");  assert (status);
     281
     282        psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLX", 0, "psf clump center", psfClump.X);
     283        psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLY", 0, "psf clump center", psfClump.Y);
     284        psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLDX", 0, "psf clump size", psfClump.dX);
     285        psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLDY", 0, "psf clump size", psfClump.dY);
    256286
    257287        // save the dimensions of each parameter
    258288        for (int i = 0; i < nPar; i++) {
    259289            char name[9];
    260             psPolynomial2D *poly = psf->params->data[i];
    261             if (poly == NULL) continue;
     290            int nX, nY;
     291
     292            pmTrend2D *trend = psf->params->data[i];
     293            if (trend == NULL) continue;
     294
     295            if (trend->mode == PM_TREND_MAP) {
     296              nX = trend->map->map->numCols;
     297              nY = trend->map->map->numRows;
     298            } else {
     299              nX = trend->poly->nX;
     300              nY = trend->poly->nY;
     301            }
    262302            snprintf (name, 9, "PAR%02d_NX", i);
    263             psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nX);
     303            psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", nX);
    264304            snprintf (name, 9, "PAR%02d_NY", i);
    265             psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nY);
     305            psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", nY);
     306            snprintf (name, 9, "PAR%02d_MD", i);
     307            char *modeName = pmTrend2DModeToString (trend->mode);
     308            psMetadataAddStr (header, PS_LIST_TAIL, name, 0, "", modeName);
     309            psFree (modeName);
    266310        }
    267311
    268312        // 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);
     313        psMetadataAddF32 (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid);
     314        psMetadataAddF32 (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid);
     315        psMetadataAddF32 (header, PS_LIST_TAIL, "CHISQ",    0, "chi-square for fit", psf->chisq);
     316        psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS",   0, "number of stars used to measure PSF", psf->nPSFstars);
    273317
    274318        // XXX can we drop this now?
    … …  
    278322        psArray *psfTable = psArrayAllocEmpty (100);
    279323        for (int i = 0; i < nPar; i++) {
    280             psPolynomial2D *poly = psf->params->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);
     324            pmTrend2D *trend = psf->params->data[i];
     325            if (trend == NULL) continue; // skip unset parameters (eg, XPOS)
     326
     327            if (trend->mode == PM_TREND_MAP) {
     328                // write the image components into a table: this is needed because they may each be a different size
     329                psImageMap *map = trend->map;
     330                for (int ix = 0; ix < map->map->numCols; ix++) {
     331                    for (int iy = 0; iy < map->map->numRows; iy++) {
     332                        psMetadata *row = psMetadataAlloc ();
     333                        psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
     334                        psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
     335                        psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
     336                        psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", map->map->data.F32[iy][ix]);
     337                        psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", map->error->data.F32[iy][ix]);
     338                        psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", 0); // no cells are masked
     339
     340                        psArrayAdd (psfTable, 100, row);
     341                        psFree (row);
     342                    }
     343                }
     344            } else {
     345                // write the polynomial components into a table
     346                psPolynomial2D *poly = trend->poly;
     347                for (int ix = 0; ix <= poly->nX; ix++) {
     348                    for (int iy = 0; iy <= poly->nY; iy++) {
     349                        psMetadata *row = psMetadataAlloc ();
     350                        psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
     351                        psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
     352                        psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
     353                        psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
     354                        psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
     355                        psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->mask[ix][iy]);
     356
     357                        psArrayAdd (psfTable, 100, row);
     358                        psFree (row);
     359                    }
    295360                }
    296361            }
    … …  
    356421    // XXX save the growth curve
    357422    // XXX save ApTrend (as image?)
     423    // XXX write the ApTrend with the same API as will be used for the PSF parameters above
    358424
    359425# if (0)
    … …  
    474540    psTrace ("psModules.objects", 5, "read psf model for %s\n", file->filename);
    475541
     542    // select the current recipe
     543    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSPHOT");
     544    if (!recipe) {
     545        psError(PS_ERR_UNKNOWN, false, "missing recipe %s", "PSPHOT");
     546        return false;
     547    }
     548
    476549    // Menu of EXTNAME rules
    477550    psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
    … …  
    506579    if (!header) psAbort("cannot read table header");
    507580
     581    pmPSFOptions *options = pmPSFOptionsAlloc();
     582
    508583    // load the PSF model parameters from the FITS table
    509584    char *modelName = psMetadataLookupStr (&status, header, "PSF_NAME");
    510     pmModelType type = pmModelClassGetType (modelName);
    511     if (type == -1) {
     585    options->type = pmModelClassGetType (modelName);
     586    if (options->type == -1) {
    512587        psError(PS_ERR_UNKNOWN, true, "invalid model name %s in psf file %s", modelName, file->filename);
    513588        return false;
    514589    }
    515590
    516     bool poissonErrors = psMetadataLookupBool (&status, header, "POISSON");
     591    // psf clump data
     592    pmPSFClump psfClump;
     593
     594    psfClump.X  = psMetadataLookupF32 (&status, header, "PSF_CLX" );  assert(status);
     595    psfClump.Y  = psMetadataLookupF32 (&status, header, "PSF_CLY" );  assert(status);
     596    psfClump.dX = psMetadataLookupF32 (&status, header, "PSF_CLDX");  assert(status);
     597    psfClump.dY = psMetadataLookupF32 (&status, header, "PSF_CLDY");  assert(status);
     598
     599    psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.X" , 0, "psf clump center", psfClump.X);
     600    psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.Y" , 0, "psf clump center", psfClump.Y);
     601    psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.DX", 0, "psf clump size",   psfClump.dX);
     602    psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.DY", 0, "psf clump size",   psfClump.dY);
     603
     604    options->poissonErrorsPhotLMM = psMetadataLookupBool (&status, header, "ERR_LMM");
     605    options->poissonErrorsPhotLin = psMetadataLookupBool (&status, header, "ERR_LIN");
     606    options->poissonErrorsParams  = psMetadataLookupBool (&status, header, "ERR_PAR");
     607
     608    options->psfFieldNx = psMetadataLookupS32 (&status, header, "IMAXIS1");
     609    options->psfFieldNy = psMetadataLookupS32 (&status, header, "IMAXIS2");
     610    options->psfFieldXo = psMetadataLookupS32 (&status, header, "IMREF1");
     611    options->psfFieldYo = psMetadataLookupS32 (&status, header, "IMREF2");
     612
     613    psImageBinning *binning = psImageBinningAlloc();
     614    binning->nXfine = options->psfFieldNx;
     615    binning->nYfine = options->psfFieldNy;
    517616
    518617    // we determine the PSF parameter polynomials from the MD-defined polynomials
    519     pmPSF *psf = pmPSFAlloc (type, poissonErrors, NULL);
     618    pmPSF *psf = pmPSFAlloc (options);
    520619
    521620    // check the number of expected parameters
    … …  
    524623        psAbort("mismatch model par count");
    525624
    526     // load the dimensions of each parameter
     625    // load the trend mode and dimensions of each parameter
    527626    for (int i = 0; i < nPar; i++) {
    528627        char name[9];
    529628        snprintf (name, 9, "PAR%02d_NX", i);
    530         int nXorder = psMetadataLookupS32 (&status, header, name);
     629        binning->nXruff = psMetadataLookupS32 (&status, header, name);
    531630        if (!status) continue;          // not all parameters are defined
     631
    532632        snprintf (name, 9, "PAR%02d_NY", i);
    533         int nYorder = psMetadataLookupS32 (&status, header, name);
     633        binning->nYruff = psMetadataLookupS32 (&status, header, name);
    534634        if (!status) {
    535             psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX defined for PAR %d, but no NY", i);
     635            psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX defined for PAR %d, but not NY", i);
    536636            return false;
    537637        }
    538         psf->params->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder);
    539     }
     638
     639        snprintf (name, 9, "PAR%02d_MD", i);
     640        char *modeName = psMetadataLookupStr (&status, header, name);
     641        if (!status) {
     642            psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX & NY defined for PAR %d, but not MD", i);
     643            return false;
     644        }
     645        pmTrend2DMode psfTrendMode = pmTrend2DModeFromString (modeName);
     646        if (psfTrendMode == PM_TREND_NONE) {
     647            psfTrendMode = PM_TREND_POLY_ORD;
     648        }
     649
     650        psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
     651        psImageBinningSetSkipByOffset (binning, options->psfFieldXo, options->psfFieldYo);
     652        psf->params->data[i] = pmTrend2DNoImageAlloc (psfTrendMode, binning, NULL);
     653    }
     654    psFree (binning);
    540655
    541656    // other required information describing the PSF
    … …  
    557672        int xPow = psMetadataLookupS32 (&status, row, "X_POWER");
    558673        int yPow = psMetadataLookupS32 (&status, row, "Y_POWER");
    559         // XXX sanity check here
    560 
    561         psPolynomial2D *poly = psf->params->data[iPar];
    562         if (poly == NULL) {
    563             psError(PS_ERR_UNKNOWN, true, "values for parameter %d, but missing NX", iPar);
     674
     675        pmTrend2D *trend = psf->params->data[iPar];
     676        if (trend == NULL) {
     677            psError(PS_ERR_UNKNOWN, true, "parameter %d not available", iPar);
    564678            return false;
    565679        }
    566680
    567         poly->coeff[xPow][yPow]    = psMetadataLookupF32 (&status, row, "VALUE");
    568         poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR");
    569         poly->mask[xPow][yPow]     = psMetadataLookupU8  (&status, row, "MASK");
     681        if (trend->mode == PM_TREND_MAP) {
     682            psImageMap *map = trend->map;
     683            assert (map);
     684            assert (map->map);
     685            assert (map->error);
     686            assert (xPow >= 0);
     687            assert (yPow >= 0);
     688            assert (xPow < map->map->numCols);
     689            assert (yPow < map->map->numRows);
     690            map->map->data.F32[yPow][xPow]    = psMetadataLookupF32 (&status, row, "VALUE");
     691            map->error->data.F32[yPow][xPow]  = psMetadataLookupF32 (&status, row, "ERROR");
     692        } else {
     693            psPolynomial2D *poly = trend->poly;
     694            assert (poly);
     695            assert (xPow >= 0);
     696            assert (yPow >= 0);
     697            assert (xPow <= poly->nX);
     698            assert (yPow <= poly->nY);
     699            poly->coeff[xPow][yPow]    = psMetadataLookupF32 (&status, row, "VALUE");
     700            poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR");
     701            poly->mask[xPow][yPow]     = psMetadataLookupU8  (&status, row, "MASK");
     702        }
    570703    }
    571704    psFree (header);
    … …  
    618751
    619752// create a psMetadata representation (human-readable) of a psf model
    620 psMetadata *pmPSFtoMetadata (psMetadata *metadata, pmPSF *psf)
    621 {
    622 
    623     if (metadata == NULL) {
    624         metadata = psMetadataAlloc ();
    625     }
    626 
    627     char *modelName = pmModelClassGetName (psf->type);
    628     psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName);
    629 
    630     int nPar = pmModelClassParameterCount (psf->type)    ;
    631     psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
    632 
    633     for (int i = 0; i < nPar; i++) {
    634         psPolynomial2D *poly = psf->params->data[i];
    635         if (poly == NULL)
    636             continue;
    637         psPolynomial2DtoMetadata (metadata, poly, "PSF_PAR%02d", i);
    638     }
    639 
    640     // XXX fix this
    641     psWarning ("APTREND is currently missing");
    642     // psPolynomial4DtoMetadata (metadata, psf->ApTrend, "APTREND");
    643 
    644     psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_AP_RESID", PS_DATA_F32, "aperture residual", psf->ApResid);
    645     psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_dAP_RESID", PS_DATA_F32, "aperture residual scatter", psf->dApResid);
    646     psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias);
    647 
    648     psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "chi-square for fit", psf->chisq);
    649     psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_NSTARS", PS_DATA_S32, "number of stars used to measure PSF", psf->nPSFstars);
    650     psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_POISSON_ERRORS", PS_DATA_BOOL, "Poisson errors for fits", psf->poissonErrors);
    651 
    652     return metadata;
    653 }
    654 
    655 // parse a psMetadata representation (human-readable) of a psf model
    656 pmPSF *pmPSFfromMetadata (psMetadata *metadata)
    657 {
    658 
    659     bool status;
    660     char keyword[80];
    661 
    662     char *modelName = psMetadataLookupPtr (&status, metadata, "PSF_MODEL_NAME");
    663     pmModelType type = pmModelClassGetType (modelName);
    664 
    665     bool poissonErrors = psMetadataLookupPtr (&status, metadata, "PSF_POISSON_ERRORS");
    666     if (!status)
    667         poissonErrors = true;
    668 
    669     // we determine the PSF parameter polynomials from the MD-defined polynomials
    670     pmPSF *psf = pmPSFAlloc (type, poissonErrors, NULL);
    671 
    672     int nPar = psMetadataLookupS32 (&status, metadata, "PSF_MODEL_NPAR");
    673     if (nPar != pmModelClassParameterCount (psf->type))
    674         psAbort("mismatch model par count");
    675 
    676     // un-fitted terms, not in the Metadata, are left NULL
    677     // XXX add a double-check of the expected number?
    678     for (int i = 0; i < nPar; i++) {
    679         sprintf (keyword, "PSF_PAR%02d", i);
    680         psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword);
    681         if (!status)
    682             continue;
    683         psPolynomial2D *poly = psPolynomial2DfromMetadata (folder);
    684         psFree (psf->params->data[i]);
    685         psf->params->data[i] = poly;
    686     }
    687 
    688     // load the APTREND data
    689     // XXX fix this to work with pmTrend2D
    690     psWarning ("APTREND is not being read");
    691     # if (0)
    692     sprintf (keyword, "APTREND");
    693     psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword);
    694     psPolynomial4D *poly = psPolynomial4DfromMetadata (folder);
    695     psFree (psf->ApTrend);
    696     psf->ApTrend = poly;
    697     # endif
    698 
    699     psf->ApResid = psMetadataLookupF32 (&status, metadata, "PSF_AP_RESID");
    700     psf->dApResid = psMetadataLookupF32 (&status, metadata, "PSF_dAP_RESID");
    701     psf->skyBias = psMetadataLookupF32 (&status, metadata, "PSF_SKY_BIAS");
    702 
    703     psf->chisq = psMetadataLookupF32 (&status, metadata, "PSF_CHISQ");
    704     psf->nPSFstars = psMetadataLookupS32 (&status, metadata, "PSF_NSTARS");
    705 
    706     psFree (metadata);
    707     return (psf);
    708 }
     753// XXX pmPSF to/from Metadata functions were defined for 1.22 and earlier, but were dropped
Note: See TracChangeset for help on using the changeset viewer.