Changeset 10399

Timestamp:

Dec 1, 2006, 1:18:58 PM (20 years ago)

Author:

eugene

Message:

added SKY_HISTOGRAM_BINS and moved arguments to psImageBackground to a psStats argument

File:

• trunk/psphot/src/psphotImageMedian.c (modified) (9 diffs)

trunk/psphot/src/psphotImageMedian.c

@@ -9 +9 @@
     pmFPAfile *file;
     psRegion region;
-    int MAX_SAMPLE_PIXELS;
     static char *defaultStatsName = "FITTED_MEAN";
 
@@ -22 +21 @@
     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
 
-    MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
-    if (!status) {
-        MAX_SAMPLE_PIXELS = 1000;
-    }
+    // user supplied seed, if available
     unsigned long seed = psMetadataLookupS32 (&status, recipe, "IMSTATS_SEED");
     if (!status) {
@@ -38 +34 @@
     }
 
-    // supply the sky background statistic name:
+    // supply the sky background statistics options
     char *statsName = psMetadataLookupStr (&status, recipe, "SKY_STAT");
     if (statsName == NULL) {
@@ -46 +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))) {
         statsOption = PS_STAT_FITTED_MEAN;
+    }
+    psStats *stats = psStatsAlloc (statsOption);
+
+    // set range for old-version of sky statistic
+    if (statsOption & PS_STAT_ROBUST_QUARTILE) {
+        stats->min = 0.25;
+        stats->max = 0.75;
+    }
+
+    // set user-option for number of pixels per region
+    stats->nSubsample = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
+    if (!status) {
+        stats->nSubsample = 1000;
+    }
+
+    // optionally set the binsize
+    stats->binsize = psMetadataLookupF32 (&status, recipe, "SKY_HISTOGRAM_BINS");
+    if (status) {
+        stats->options |= PS_STAT_USE_BINSIZE;
     }
 
@@ -116 +131 @@
 
             // XXX the value of the upper and lower cuts probably should be studied...
-            psStats *stats = psImageBackground(subset, submask, 0xff, 0.25, 0.75, MAX_SAMPLE_PIXELS,
-                                               statsOption, rng);
-            if (stats != NULL) {
+            if (psImageBackground(stats, subset, submask, 0xff, rng)) {
                 if (stats->options & PS_STAT_ROBUST_QUARTILE) {
                     modelData[iy][ix] = stats->robustMedian;
@@ -125 +138 @@
                 }
             } else {
-                psStatsOptions fallback_option = PS_STAT_ROBUST_MEDIAN;
-                stats = psImageBackground(subset, submask, 0xff, 0.25, 0.75, MAX_SAMPLE_PIXELS,
-                                          fallback_option, rng);
-                if (stats == NULL) {
+                psStatsOptions currentOptions = stats->options;
+                stats->options |= PS_STAT_ROBUST_MEDIAN;
+                if (!psImageBackground(stats, subset, submask, 0xff, rng)) {
                     psAbort(PS_FILE_LINE, "Failed to estimate background using ROBUST_MEDIAN for "
                             "(%dx%d, (row0,col0) = (%d,%d)",
@@ -134 +146 @@
                 }
                 modelData[iy][ix] = psStatsGetValue (stats, PS_STAT_ROBUST_MEDIAN);
+                stats->options = currentOptions;
             }
             modelData[iy][ix] += SKY_BIAS;
 
-            psFree (stats);
             psFree (subset);
             psFree (submask);
@@ -150 +162 @@
 
     // measure background stats and save for later output
-    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_MIN | PS_STAT_MAX);
-    stats = psImageStats (stats, model->image, NULL, 0);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_MEAN", PS_DATA_F32 | PS_META_REPLACE, "sky model mean",          stats->sampleMean);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_STDEV", PS_DATA_F32 | PS_META_REPLACE, "sky model stdev",        stats->sampleStdev);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_MAX",  PS_DATA_F32 | PS_META_REPLACE, "sky model maximum value", stats->max);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_MIN",  PS_DATA_F32 | PS_META_REPLACE, "sky model minimum value", stats->min);
+    psStats *statsBck = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_MIN | PS_STAT_MAX);
+    psImageStats (statsBck, model->image, NULL, 0);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_MEAN", PS_DATA_F32 | PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_STDEV", PS_DATA_F32 | PS_META_REPLACE, "sky model stdev",        statsBck->sampleStdev);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_MAX",  PS_DATA_F32 | PS_META_REPLACE, "sky model maximum value", statsBck->max);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_MIN",  PS_DATA_F32 | PS_META_REPLACE, "sky model minimum value", statsBck->min);
     psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_NX",   PS_DATA_S32 | PS_META_REPLACE, "sky model size (x)",      nx);
     psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_NY",   PS_DATA_S32 | PS_META_REPLACE, "sky model size (y)",      ny);
-    psLogMsg ("psphot", 3, "background sky : min %f mean %f max %f stdev %f", stats->min, stats->sampleMean, stats->max, stats->sampleStdev);
-    psFree (stats);
+    psLogMsg ("psphot", 3, "background sky : min %f mean %f max %f stdev %f", statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev);
+    psFree (statsBck);
 
     // select background pixels, from output background file, or create
@@ -216 +228 @@
 
     psLogMsg ("psphot", 3, "subtracted background model: %f sec\n", psTimerMark ("psphot"));
+    psFree(stats);
     psFree(rng);