IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 13013


Ignore:
Timestamp:
Apr 24, 2007, 3:29:40 PM (19 years ago)
Author:
rhl
Message:

Save PSPHOT.BACKMDL.STDEV; change e.g. SKY_MEAN to SKY_MODEL_MEAN

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotImageMedian.c

    r12792 r13013  
    22static int npass = 0;
    33
    4 // generate the median in NxN boxes, clipping heavily
    5 // linear interpolation to generate full-scale model
    6 bool psphotImageMedian (pmConfig *config, pmFPAview *view)
    7 {
    8     bool status;
    9     pmFPA *inFPA;
    10     pmFPAfile *file;
    11     static char *defaultStatsName = "FITTED_MEAN";
    12 
     4// we have 4 possibilities: (INTERNAL or I/O file) and (exists or not)
     5// select model pixels (from output background model file, or create internal file)
     6static pmReadout *get_model_readout(const char *name, // name of internal/external file
     7                                    const pmConfig *config, // configuration information
     8                                    pmFPAview *view,
     9                                    pmFPA *inFPA,
     10                                    const psImageBinning *binning) {
    1311    pmReadout *model = NULL;
    14     pmReadout *readout = NULL;
    15     pmReadout *background = NULL;
    16     pmReadout *backSub = NULL;
    17 
    18     psTimerStart ("psphot");
    19 
    20     // select the appropriate recipe information
    21     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    22 
    23     // user supplied seed, if available
    24     unsigned long seed = psMetadataLookupS32 (&status, recipe, "IMSTATS_SEED");
    25     if (!status) {
    26         seed = 0;
    27     }
    28     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, seed);
    29 
    30     // subtract this amount extra from the sky
    31     float SKY_BIAS = psMetadataLookupF32 (&status, recipe, "SKY_BIAS");
    32     if (!status) {
    33         SKY_BIAS = 0;
    34     }
    35 
    36     // supply the sky background statistics options
    37     char *statsName = psMetadataLookupStr (&status, recipe, "SKY_STAT");
    38     if (statsName == NULL) {
    39         statsName = defaultStatsName;
    40     }
    41     psStatsOptions statsOption = psStatsOptionFromString (statsName);
    42     if (!(statsOption & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_QUARTILE | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_MEAN_V3))) {
    43         statsOption = PS_STAT_FITTED_MEAN;
    44     }
    45     psStats *stats = psStatsAlloc (statsOption);
    46 
    47     // set range for old-version of sky statistic
    48     if (statsOption & PS_STAT_ROBUST_QUARTILE) {
    49         stats->min = 0.25;
    50         stats->max = 0.75;
    51     }
    52 
    53     // set user-option for number of pixels per region
    54     stats->nSubsample = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
    55     if (!status) {
    56         stats->nSubsample = 1000;
    57     }
    58 
    59     // optionally set the binsize
    60     stats->binsize = psMetadataLookupF32 (&status, recipe, "SKY_HISTOGRAM_BINS");
    61     if (status) {
    62         stats->options |= PS_STAT_USE_BINSIZE;
    63     }
    64 
    65     // optionally set the sigma clipping
    66     stats->clipSigma = psMetadataLookupF32 (&status, recipe, "SKY_CLIP_SIGMA");
    67     if (!status) {
    68         if ((stats->options & PS_STAT_FITTED_MEAN) || (stats->options & PS_STAT_FITTED_MEAN_V2)) {
    69             stats->clipSigma = 1.0;
    70         } else {
    71             stats->clipSigma = 3.0;
    72         }
    73     }
    74 
    75     // stats is not initialized by psStats???  use this to save the input options
    76     psStats *statsDefaults = psStatsAlloc (statsOption);
    77     *statsDefaults = *stats;
    78 
    79     // find the currently selected readout
    80     file = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
    81     inFPA = file->fpa;
    82     readout = pmFPAviewThisReadout (view, inFPA);
    83 
    84     psImage *image = readout->image;
    85     psImage *mask  = readout->mask;
    86 
    87     // I have the fine image size, I know the binning factor, determine the ruff image size
    88     psImageBinning *binning = psImageBinningAlloc();
    89     binning->nXfine = image->numCols;
    90     binning->nYfine = image->numRows;
    91     binning->nXbin  = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN");
    92     binning->nYbin  = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN");
    93 
    94     psImageBinningSetRuffSize(binning, PS_IMAGE_BINNING_CENTER);
    95     psImageBinningSetSkip(binning, image);
    96     status = psMetadataAddPtr(recipe, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning);
    97     PS_ASSERT (status, false);
    98 
    99     // we save the binning structure for use in psphotMagnitudes
    100 
    101     // we have 4 possibilities: (INTERNAL or I/O file) and (exists or not)
    102     // select model pixels (from output background model file, or create internal file)
    103     file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKMDL");
     12
     13    bool status = true;
     14    pmFPAfile *file = psMetadataLookupPtr(&status, config->files, name);
    10415    if (file == NULL) {
    10516        // we are not using PSPHOT.BACKMDL as an I/O file: define an internal version
    106         model = pmFPAfileDefineInternal (config->files, "PSPHOT.BACKMDL", binning->nXruff, binning->nYruff, PS_TYPE_F32);
     17        model = pmFPAfileDefineInternal (config->files, name, binning->nXruff, binning->nYruff, PS_TYPE_F32);
    10718    } else {
    10819        if (file->mode == PM_FPA_MODE_INTERNAL) {
    … …  
    12334        }
    12435    }
     36
     37    return model;
     38}
     39
     40// generate the median in NxN boxes, clipping heavily
     41// linear interpolation to generate full-scale model
     42bool psphotImageMedian (pmConfig *config, pmFPAview *view)
     43{
     44    bool status = true;
     45    static char *defaultStatsName = "FITTED_MEAN";
     46
     47    pmReadout *readout = NULL;
     48    pmReadout *background = NULL;
     49    pmReadout *backSub = NULL;
     50
     51    psTimerStart ("psphot");
     52
     53    // select the appropriate recipe information
     54    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     55
     56    // user supplied seed, if available
     57    unsigned long seed = psMetadataLookupS32 (&status, recipe, "IMSTATS_SEED");
     58    if (!status) {
     59        seed = 0;
     60    }
     61    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, seed);
     62
     63    // subtract this amount extra from the sky
     64    float SKY_BIAS = psMetadataLookupF32 (&status, recipe, "SKY_BIAS");
     65    if (!status) {
     66        SKY_BIAS = 0;
     67    }
     68
     69    // supply the sky background statistics options
     70    char *statsName = psMetadataLookupStr (&status, recipe, "SKY_STAT");
     71    if (statsName == NULL) {
     72        statsName = defaultStatsName;
     73    }
     74    psStatsOptions statsOptionLocation = psStatsOptionFromString(statsName);
     75    if (!(statsOptionLocation & (PS_STAT_SAMPLE_MEAN |
     76                                 PS_STAT_SAMPLE_MEDIAN |
     77                                 PS_STAT_ROBUST_MEDIAN |
     78                                 PS_STAT_ROBUST_QUARTILE |
     79                                 PS_STAT_CLIPPED_MEAN |
     80                                 PS_STAT_FITTED_MEAN |
     81                                 PS_STAT_FITTED_MEAN_V2 |
     82                                 PS_STAT_FITTED_MEAN_V3))) {
     83        statsOptionLocation = PS_STAT_FITTED_MEAN;
     84    }
     85
     86    psStatsOptions statsOptionWidth = PS_STAT_NONE;
     87    if (statsOptionLocation & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN)) {
     88        statsOptionWidth = PS_STAT_SAMPLE_STDEV;
     89    } else if (statsOptionLocation & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_QUARTILE)) {
     90#if 1
     91        statsOptionWidth = PS_STAT_ROBUST_STDEV; // not set; => NaN
     92#else
     93        statsOptionWidth = PS_STAT_FITTED_STDEV;
     94#endif
     95    } else if (statsOptionLocation & PS_STAT_FITTED_MEAN) {
     96        statsOptionWidth = PS_STAT_FITTED_STDEV;
     97    } else if (statsOptionLocation & PS_STAT_CLIPPED_MEAN) {
     98        statsOptionWidth = PS_STAT_CLIPPED_STDEV;
     99    } else if (statsOptionLocation & PS_STAT_FITTED_MEAN_V2) {
     100        statsOptionWidth = PS_STAT_FITTED_STDEV_V2;
     101    } else if (statsOptionLocation & PS_STAT_FITTED_MEAN_V3) {
     102        statsOptionWidth = PS_STAT_FITTED_STDEV_V3;
     103    } else {
     104        psAbort("Unable to estimate variance of selected statsOptionLocations 0x%x", statsOptionLocation);
     105    }
     106
     107    const psStatsOptions statsOption = statsOptionLocation | statsOptionWidth;
     108    psStats *stats = psStatsAlloc (statsOption);
     109
     110    // set range for old-version of sky statistic
     111    if (statsOptionLocation & PS_STAT_ROBUST_QUARTILE) {
     112        stats->min = 0.25;
     113        stats->max = 0.75;
     114    }
     115
     116    // set user-option for number of pixels per region
     117    stats->nSubsample = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
     118    if (!status) {
     119        stats->nSubsample = 1000;
     120    }
     121
     122    // optionally set the binsize
     123    stats->binsize = psMetadataLookupF32 (&status, recipe, "SKY_HISTOGRAM_BINS");
     124    if (status) {
     125        stats->options |= PS_STAT_USE_BINSIZE;
     126    }
     127
     128    // optionally set the sigma clipping
     129    stats->clipSigma = psMetadataLookupF32 (&status, recipe, "SKY_CLIP_SIGMA");
     130    if (!status) {
     131        if ((stats->options & PS_STAT_FITTED_MEAN) || (stats->options & PS_STAT_FITTED_MEAN_V2)) {
     132            stats->clipSigma = 1.0;
     133        } else {
     134            stats->clipSigma = 3.0;
     135        }
     136    }
     137
     138    // stats is not initialized by psStats???  use this to save the input options
     139    psStats *statsDefaults = psStatsAlloc (statsOption);
     140    *statsDefaults = *stats;
     141
     142    // find the currently selected readout
     143    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
     144    pmFPA *inFPA = file->fpa;
     145    readout = pmFPAviewThisReadout (view, inFPA);
     146
     147    psImage *image = readout->image;
     148    psImage *mask  = readout->mask;
     149
     150    // I have the fine image size, I know the binning factor, determine the ruff image size
     151    psImageBinning *binning = psImageBinningAlloc();
     152    binning->nXfine = image->numCols;
     153    binning->nYfine = image->numRows;
     154    binning->nXbin  = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN");
     155    binning->nYbin  = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN");
     156
     157    psImageBinningSetRuffSize(binning, PS_IMAGE_BINNING_CENTER);
     158    psImageBinningSetSkip(binning, image);
     159    status = psMetadataAddPtr(recipe, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning);
     160    PS_ASSERT (status, false);
     161
     162    // we save the binning structure for use in psphotMagnitudes
     163    pmReadout *model = get_model_readout("PSPHOT.BACKMDL", config, view, inFPA, binning);
     164    pmReadout *modelStdev = get_model_readout("PSPHOT.BACKMDL.STDEV", config, view, inFPA, binning);
     165
    125166    psF32 **modelData = model->image->data.F32;
     167    psF32 **modelStdevData = modelStdev->image->data.F32;
    126168
    127169    // measure clipped median for subimages
    … …  
    157199                    modelData[iy][ix] = stats->robustMedian;
    158200                } else {       
    159                     modelData[iy][ix] = psStatsGetValue(stats, statsOption);
     201                    modelData[iy][ix] = psStatsGetValue(stats, statsOptionLocation);
    160202                }
     203                modelStdevData[iy][ix] = psStatsGetValue(stats, statsOptionWidth);
    161204            } else {
    162205                psStatsOptions currentOptions = stats->options;
    163                 stats->options = PS_STAT_ROBUST_MEDIAN;
     206                stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV;
    164207                if (!psImageBackground(stats, subset, submask, 0xff, rng)) {
    165208                    psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for "
    166209                               "(%dx%d, (row0,col0) = (%d,%d)",
    167210                               subset->numRows, subset->numCols, subset->row0, subset->col0);
    168                     modelData[iy][ix] = NAN;
     211                    modelData[iy][ix] = modelStdevData[iy][ix] = NAN;
    169212                } else {
    170213                    modelData[iy][ix] = psStatsGetValue (stats, PS_STAT_ROBUST_MEDIAN);
     214                    modelStdevData[iy][ix] = psStatsGetValue(stats, PS_STAT_ROBUST_STDEV);
    171215                }
    172216                // drop errors caused by psImageBackground failures
    … …  
    183227    // patch over bad regions (use average of 8 possible neighbor pixels)
    184228    // XXX consider testing all pixels against the 8 neighbors and replacing outliers...
    185     float Count = 0;
    186     float Value = 0;
     229    double Count = 0;                   // number of good pixels
     230    double Value = 0;                   // sum of good pixel's value
     231    double ValueStdev = 0;              // sum of good pixel's standard deviations
    187232    for (int iy = 0; iy < model->image->numRows; iy++) {
    188233        for (int ix = 0; ix < model->image->numCols; ix++) {
    189234            if (!isnan(modelData[iy][ix])) {
    190235                Value += modelData[iy][ix];
    191                 Count += 1;
     236                ValueStdev += modelStdevData[iy][ix];
     237                Count++;
    192238                continue;
    193239            }
    194             float value = 0;
    195             float count = 0;
     240
     241            double value = 0;
     242            double count = 0;
    196243            for (int jy = iy - 1; jy <= iy + 1; jy++) {
    197244                if (jy <   0) continue;
    … …  
    206253            }
    207254            if (count > 0) modelData[iy][ix] = value / count;
     255        }
     256    }
     257    assert (Count > 0);
     258    Value /= Count;
     259    ValueStdev /= Count;
    208260           
    209         }
    210     }
    211     Value = Value / Count;
    212 
    213261    // patch over remaining bad regions (use global average)
    214262    for (int iy = 0; iy < model->image->numRows; iy++) {
    … …  
    216264            if (!isnan(modelData[iy][ix])) continue;
    217265            modelData[iy][ix] = Value;
     266            modelStdevData[iy][ix] = ValueStdev;
    218267        }
    219268    }
    … …  
    221270    psLogMsg ("psphot", PS_LOG_MINUTIA, "build median image: %f sec\n", psTimerMark ("psphot"));
    222271
    223     // measure background stats and save for later output
    224     psStats *statsBck = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_MIN | PS_STAT_MAX);
     272    psMetadataAddF32(recipe, PS_LIST_TAIL, "SKY_MEAN", PS_META_REPLACE, "sky mean", Value);
     273    psMetadataAddF32(recipe, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev);
     274    psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f stdev %f", Value, ValueStdev);
     275   
     276    // measure image and background stats and save for later output
     277    psStats *statsBck = psStatsAlloc (PS_STAT_SAMPLE_MEAN |
     278                                      PS_STAT_SAMPLE_STDEV |
     279                                      PS_STAT_MIN |
     280                                      PS_STAT_MAX);
    225281    psImageStats (statsBck, model->image, NULL, 0);
    226     psMetadataAddF32 (recipe, PS_LIST_TAIL, "SKY_MEAN", PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
    227     psMetadataAddF32 (recipe, PS_LIST_TAIL, "SKYSTDEV", PS_META_REPLACE, "sky model stdev",         statsBck->sampleStdev);
    228     psMetadataAddF32 (recipe, PS_LIST_TAIL, "SKY_MAX",  PS_META_REPLACE, "sky model maximum value", statsBck->max);
    229     psMetadataAddF32 (recipe, PS_LIST_TAIL, "SKY_MIN",  PS_META_REPLACE, "sky model minimum value", statsBck->min);
    230     psMetadataAddS32 (recipe, PS_LIST_TAIL, "SKY_NX",   PS_META_REPLACE, "sky model size (x)",      model->image->numCols);
    231     psMetadataAddS32 (recipe, PS_LIST_TAIL, "SKY_NY",   PS_META_REPLACE, "sky model size (y)",      model->image->numRows);
     282    psMetadataAddF32 (recipe, PS_LIST_TAIL, "SKY_MODEL_MEAN",
     283                      PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
     284    psMetadataAddF32 (recipe, PS_LIST_TAIL, "SKY_MODEL_STDEV",
     285                      PS_META_REPLACE, "sky model stdev",         statsBck->sampleStdev);
     286    psMetadataAddF32 (recipe, PS_LIST_TAIL, "SKY_MODEL_MAX",
     287                      PS_META_REPLACE, "sky model maximum value", statsBck->max);
     288    psMetadataAddF32 (recipe, PS_LIST_TAIL, "SKY_MODEL_MIN",
     289                      PS_META_REPLACE, "sky model minimum value", statsBck->min);
     290    psMetadataAddS32 (recipe, PS_LIST_TAIL, "SKY_MODEL_NX",
     291                      PS_META_REPLACE, "sky model size (x)",      model->image->numCols);
     292    psMetadataAddS32 (recipe, PS_LIST_TAIL, "SKY_MODEL_NY",
     293                      PS_META_REPLACE, "sky model size (y)",      model->image->numRows);
    232294    psLogMsg ("psphot", PS_LOG_INFO, "background sky : min %f mean %f max %f stdev %f",
    233295              statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev);
Note: See TracChangeset for help on using the changeset viewer.