IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 33094


Ignore:
Timestamp:
Jan 11, 2012, 2:23:07 PM (15 years ago)
Author:
eugene
Message:

generate a background model of the variance image for use in the model-based variance; finish the model-based variance analysis

Location:
branches/eam_branches/ipp-20111122/psphot/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20111122/psphot/src/psphotDefineFiles.c

    r25983 r33094  
    8585        output->save = true;
    8686    }
     87//    // optionally save the smoothed variance model (small FITS image)
     88//    if (psMetadataLookupBool(NULL, recipe, "SAVE.BACKMDL")) {
     89//        int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN");
     90//        int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN");
     91//        pmFPAfile *output = pmFPAfileDefineFromFile (config, input, DX, DY, "PSPHOT.VARMDL");
     92//        if (!output) {
     93//            psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.VARMDL");
     94//            return false;
     95//        }
     96//        output->save = true;
     97//    }
     98//    // optionally save the smoothed variance model's standard deviation (small FITS image)
     99//    if (psMetadataLookupBool(NULL, recipe, "SAVE.VARMDL.STDEV")) {
     100//        int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN");
     101//        int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN");
     102//        pmFPAfile *output = pmFPAfileDefineFromFile (config, input, DX, DY, "PSPHOT.VARMDL.STDEV");
     103//        if (!output) {
     104//            psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.VARMDL.STDEV");
     105//            return false;
     106//        }
     107//        output->save = true;
     108//    }
    87109    // optionally save the PSF Model
    88110    if (psMetadataLookupBool(NULL, recipe, "SAVE.PSF")) {

  • branches/eam_branches/ipp-20111122/psphot/src/psphotFitSourcesLinear.c

    r33086 r33094  
    1212static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal);
    1313
    14 bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, int index, psMetadata *recipe, pmReadout *readout, psArray *sources);
     14bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources);
    1515pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe);
    1616bool psphotFreeModelVariance (pmReadout *readout, psArray *sources);
    … …  
    7474        if (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) {
    7575            // generate the model variance image & source pointers
    76             if (!psphotGenerateModelVariance (config, view, i, recipe, readout, sources)) {
     76            if (!psphotGenerateModelVariance (config, view, file, i, recipe, readout, sources)) {
    7777                psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
    7878                return false;
    … …  
    554554}
    555555
    556 bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, int index, psMetadata *recipe, pmReadout *readout, psArray *sources) {
     556bool psphotModelBackgroundReadout(psImage *model,  // Model image
     557                                  psImage *modelStdev, // Model stdev image
     558                                  psMetadata *analysis, // Analysis metadata for outputs
     559                                  pmReadout *readout, // Readout for which to generate a background model
     560                                  psImageBinning *binning, // Binning parameters
     561                                  const pmConfig *config,// Configuration
     562                                  bool useVarianceImage
     563    );
     564
     565bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources) {
    557566
    558567    bool status = false;
    … …  
    573582    psImage *modelVar = psImageCopy (NULL, readout->variance, PS_TYPE_F32);
    574583
    575     // find the pmFPAfile for the readout (background model is saved on PSPHOT.BACKMDL regardless of 'filerule')
    576     pmFPAfile *backModelFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL", index); // File of interest
    577     assert (backModelFile);
    578 
    579     // find the background model readout from the correct location
    580     pmReadout *backModel = READOUT_OR_INTERNAL(view, backModelFile);
    581     psAssert (backModel, "this must exist");
    582 
    583584    // find the binning information
    584     psImageBinning *backBinning = psMetadataLookupPtr(&status, backModel->analysis, "PSPHOT.BACKGROUND.BINNING");
     585    psImageBinning *backBinning = psphotBackgroundBinning (modelVar, config);
    585586    assert (backBinning);
     587   
     588    psImage *varModel = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
     589    psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
     590
     591    if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) {
     592        psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
     593        psFree (varModel);
     594        psFree (varModelStdev);
     595        return false;
     596    }
    586597
    587598    // linear interpolation to full-scale
    588     if (!psImageUnbin (modelVar, backModel->image, backBinning)) {
    589         psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning");
    590         return false;
    591     }
     599    if (!psImageUnbin (modelVar, varModel, backBinning)) {
     600        psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning");
     601        psFree (varModel);
     602        psFree (varModelStdev);
     603        return false;
     604    }
     605
     606    psFree (varModel);
     607    psFree (varModelStdev);
    592608
    593609    // XXX for a test:
    … …  
    609625        psAssert (!source->modelVar, "programming error : modelVar should be NULL here");
    610626        psAssert (source->modelFlux, "programming error : modelFlux should not be NULL here");
    611         source->modelVar = psImageSubset(readout->image, region);
     627        psAssert (source->modelFlux->data.F32, "programming error : modelFlux should not be NULL here");
     628        source->modelVar = psImageSubset(modelVar, region);
    612629
    613630        // add the source model to the model variance image
    … …  
    615632    }
    616633
    617     // XXX add the readnoise to the image?  This is a bit problematic because I don't have the
    618     // per-cell information.  Alternatively, I could use the median of the variance image or
    619     // something equivalent to a background model to define the base level of the variance
    620 
    621634    // XXX for a test:
    622635    psphotSaveImage (NULL, modelVar, "model.var.fits");
    623     psphotSaveImage (NULL, readout->image, "image.var.fits");
     636    psphotSaveImage (NULL, readout->variance, "image.var.fits");
    624637
    625638    // we save the model variance for future reference
    … …  
    638651    assert (modelVar);
    639652
     653    psMetadataRemoveKey (readout->analysis, "PSPHOT.MODEL.VAR");
     654
    640655    // clear modelVar pointers for all of the source models
    641656    for (int i = 0; i < sources->n; i++) {
    … …  
    645660        psFree (source->modelVar);
    646661    }
    647     psFree (modelVar);
    648662
    649663    return true;

  • branches/eam_branches/ipp-20111122/psphot/src/psphotModelBackground.c

    r33086 r33094  
    3737// corresponding to the model.  Other information about the background model is saved on the
    3838// readout->analysis
    39 static bool psphotModelBackgroundReadout(psImage *model,  // Model image
    40                                   psImage *modelStdev, // Model stdev image
    41                                   psMetadata *analysis, // Analysis metadata for outputs
    42                                   pmReadout *readout, // Readout for which to generate a background model
    43                                   psImageBinning *binning, // Binning parameters
    44                                   const pmConfig *config // Configuration
     39bool psphotModelBackgroundReadout(psImage *model,  // Model image
     40                                  psImage *modelStdev, // Model stdev image
     41                                  psMetadata *analysis, // Analysis metadata for outputs
     42                                  pmReadout *readout, // Readout for which to generate a background model
     43                                  psImageBinning *binning, // Binning parameters
     44                                  const pmConfig *config,// Configuration
     45                                  bool useVarianceImage
    4546    )
    4647{
    … …  
    4950    bool status = true;
    5051
    51     psImage *image = readout->image, *mask = readout->mask; // Image and mask for readout
     52    psImage *image = useVarianceImage ? readout->variance : readout->image;
     53    psImage *mask = readout->mask; // Image and mask for readout
    5254
    5355    // select the appropriate recipe information
    … …  
    143145
    144146    // we save the binning structure for use in psphotMagnitudes
    145     psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning);
     147    if (!useVarianceImage) {
     148        psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning);
     149    }
    146150
    147151    psVector *dQ = psVectorAllocEmpty (100, PS_TYPE_F32);
    … …  
    330334    psVectorStats (statsDQ, dQ, NULL, NULL, 0);
    331335
    332     psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_MEAN",  PS_META_REPLACE, "sky mean", Value);
    333     psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev);
    334     psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_DQ",    PS_META_REPLACE, "sky quartile slope", statsDQ->sampleMedian);
    335     psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f, stdev %f, dQ %f", Value, ValueStdev, statsDQ->sampleMedian);
    336 
    337     // measure image and background stats and save for later output
    338     psStats *statsBck = psStatsAlloc (PS_STAT_SAMPLE_MEAN |
    339                                       PS_STAT_SAMPLE_STDEV |
    340                                       PS_STAT_MIN |
    341                                       PS_STAT_MAX);
    342     psImageStats (statsBck, model, NULL, 0);
    343     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MN",  PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
    344     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_SIG", PS_META_REPLACE, "sky model stdev",         statsBck->sampleStdev);
    345     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DEV", PS_META_REPLACE, "sky stdev",               ValueStdev);
    346     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DQ",  PS_META_REPLACE, "sky quartile slope",      statsDQ->sampleMedian);
    347     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MAX", PS_META_REPLACE, "sky model maximum value", statsBck->max);
    348     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MIN", PS_META_REPLACE, "sky model minimum value", statsBck->min);
    349     psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NX", PS_META_REPLACE, "sky model size (x)",      model->numCols);
    350     psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NY", PS_META_REPLACE, "sky model size (y)",      model->numRows);
    351     psLogMsg ("psphot", PS_LOG_INFO, "background sky : min %f mean %f max %f stdev %f",
    352               statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev);
     336    if (!useVarianceImage) {
     337        psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_MEAN",  PS_META_REPLACE, "sky mean", Value);
     338        psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev);
     339        psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_DQ",    PS_META_REPLACE, "sky quartile slope", statsDQ->sampleMedian);
     340        psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f, stdev %f, dQ %f", Value, ValueStdev, statsDQ->sampleMedian);
     341
     342        // measure image and background stats and save for later output
     343        psStats *statsBck = psStatsAlloc (PS_STAT_SAMPLE_MEAN |
     344                                          PS_STAT_SAMPLE_STDEV |
     345                                          PS_STAT_MIN |
     346                                          PS_STAT_MAX);
     347        psImageStats (statsBck, model, NULL, 0);
     348        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MN",  PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
     349        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_SIG", PS_META_REPLACE, "sky model stdev",         statsBck->sampleStdev);
     350        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DEV", PS_META_REPLACE, "sky stdev",               ValueStdev);
     351        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_DQ",  PS_META_REPLACE, "sky quartile slope",      statsDQ->sampleMedian);
     352        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MAX", PS_META_REPLACE, "sky model maximum value", statsBck->max);
     353        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MIN", PS_META_REPLACE, "sky model minimum value", statsBck->min);
     354        psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NX", PS_META_REPLACE, "sky model size (x)",      model->numCols);
     355        psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NY", PS_META_REPLACE, "sky model size (y)",      model->numRows);
     356        psLogMsg ("psphot", PS_LOG_INFO, "background sky : min %f mean %f max %f stdev %f",
     357                  statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev);
     358        psFree(statsBck);
     359    } else {
     360        psLogMsg ("psphot", PS_LOG_INFO, "variance data : mean %f, stdev %f, dQ %f", Value, ValueStdev, statsDQ->sampleMedian);
     361    }
    353362
    354363    psFree(statsDQ);
    … …  
    356365
    357366    psFree(stats);
    358     psFree(statsBck);
    359367    psFree(statsDefaults);
    360368    psFree(binning);
    … …  
    375383    psImage *modelStdev = psImageAlloc(binning->nXruff, binning->nYruff, PS_TYPE_F32); // Standard deviation
    376384
    377     if (!psphotModelBackgroundReadout(model, modelStdev, readout->analysis, readout, binning, config)) {
     385    if (!psphotModelBackgroundReadout(model, modelStdev, readout->analysis, readout, binning, config, false)) {
    378386        psFree(model);
    379387        psFree(modelStdev);
    … …  
    400408    pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning, index);
    401409
    402     if (!psphotModelBackgroundReadout(model->image, modelStdev->image, model->analysis, readout, binning, config)) {
     410    if (!psphotModelBackgroundReadout(model->image, modelStdev->image, model->analysis, readout, binning, config, false)) {
    403411        psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
    404412        return false;
Note: See TracChangeset for help on using the changeset viewer.