Changeset 27838 for branches/tap_branches/psphot/src/psphotModelBackground.c

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (1 prop)
• psphot/src/psphotModelBackground.c (modified) (11 diffs)

branches/tap_branches/psphot/src

@@ -19 +19 @@
 psphotMomentsStudy
 psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

@@ -19 +19 @@
 psphotMomentsStudy
 psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

branches/tap_branches/psphot/src/psphotModelBackground.c

@@ -5 +5 @@
 // Determine the binning characteristics for the background model
 // Sets the ruff image size and skip based on recipe values for the binning
-static psImageBinning *backgroundBinning(const psImage *image, // Image for which to generate a bg model
+psImageBinning *psphotBackgroundBinning(const psImage *image, // Image for which to generate a bg model
                                          const pmConfig *config // Configuration
                                          )
@@ -28 +28 @@
 }
 
+
 // Generate the background model
-static bool backgroundModel(psImage *model,  // Model image
-                            psImage *modelStdev, // Model stdev image
-                            psMetadata *analysis, // Analysis metadata for outputs
-                            pmReadout *readout, // Readout for which to generate a background model
-                            psImageBinning *binning, // Binning parameters
-                            const pmConfig *config // Configuration
-                            )
+// generate the median in NxN boxes, clipping heavily
+// linear interpolation to generate full-scale model
+//
+// NOTE that the 'analysis' metedata pass in here is used to store the binning information.
+// This may be the analysis for this readout, but it may be the analysis for the pmFPAfile
+// corresponding to the model.  Other information about the background model is saved on the
+// readout->analysis
+static bool psphotModelBackgroundReadout(psImage *model,  // Model image
+                                  psImage *modelStdev, // Model stdev image
+                                  psMetadata *analysis, // Analysis metadata for outputs
+                                  pmReadout *readout, // Readout for which to generate a background model
+                                  psImageBinning *binning, // Binning parameters
+                                  const pmConfig *config // Configuration
+    )
 {
     psTimerStart ("psphot.background");
@@ -137 +145 @@
 
     // we save the binning structure for use in psphotMagnitudes
-    psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING",
-                     PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning);
+    psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning);
 
     psF32 **modelData = model->data.F32;
@@ -144 +151 @@
 
     // XXXX we can thread this here by running blocks in parallel
+
+    int nFailures = 0;
+    psImageBackgroundInit();
 
     // measure clipped median for subimages
@@ -211 +221 @@
                 stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV;
                 if (!psImageBackground(stats, &sample, subset, submask, maskVal, rng)) {
-                    psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for "
-                               "(%dx%d, (row0,col0) = (%d,%d)",
-                               subset->numRows, subset->numCols, subset->row0, subset->col0);
+                    if ((nFailures < 3) || (nFailures % 100 == 0)) {
+                        psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for "
+                                  "(%dx%d, (row0,col0) = (%d,%d)",
+                                  subset->numRows, subset->numCols, subset->row0, subset->col0);
+                    }
+                    nFailures ++;
                     modelData[iy][ix] = modelStdevData[iy][ix] = NAN;
                 } else {
@@ -232 +245 @@
             psFree (submask);
         }
+    }
+
+    if (nFailures) {
+        psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background for %d of %d subimages", nFailures, (model->numRows*model->numCols));
     }
 
@@ -293 +310 @@
     psLogMsg ("psphot", PS_LOG_INFO, "built background image: %f sec\n", psTimerMark ("psphot.background"));
 
-    psMetadataAddF32(recipe, PS_LIST_TAIL, "SKY_MEAN", PS_META_REPLACE, "sky mean", Value);
-    psMetadataAddF32(recipe, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev);
+    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_MEAN", PS_META_REPLACE, "sky mean", Value);
+    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev);
     psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f stdev %f", Value, ValueStdev);
 
@@ -303 +320 @@
                                       PS_STAT_MAX);
     psImageStats (statsBck, model, NULL, 0);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MN",
-                      PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_SIG",
-                      PS_META_REPLACE, "sky model stdev",         statsBck->sampleStdev);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MAX",
-                      PS_META_REPLACE, "sky model maximum value", statsBck->max);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MIN",
-                      PS_META_REPLACE, "sky model minimum value", statsBck->min);
-    psMetadataAddS32 (recipe, PS_LIST_TAIL, "MSKY_NX",
-                      PS_META_REPLACE, "sky model size (x)",      model->numCols);
-    psMetadataAddS32 (recipe, PS_LIST_TAIL, "MSKY_NY",
-                      PS_META_REPLACE, "sky model size (y)",      model->numRows);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MN", PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_SIG", PS_META_REPLACE, "sky model stdev",         statsBck->sampleStdev);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MAX", PS_META_REPLACE, "sky model maximum value", statsBck->max);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MIN", PS_META_REPLACE, "sky model minimum value", statsBck->min);
+    psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NX", PS_META_REPLACE, "sky model size (x)",      model->numCols);
+    psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NY", PS_META_REPLACE, "sky model size (y)",      model->numRows);
     psLogMsg ("psphot", PS_LOG_INFO, "background sky : min %f mean %f max %f stdev %f",
               statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev);
@@ -327 +338 @@
 }
 
-
-
-psImage *psphotBackgroundModel(pmReadout *ro, const pmConfig *config)
+// generate a background model for a single readout. do not save an associated pmFPAfile for possible output
+psImage *psphotModelBackgroundReadoutNoFile(pmReadout *readout, const pmConfig *config)
 {
-    PM_ASSERT_READOUT_NON_NULL(ro, NULL);
-    PM_ASSERT_READOUT_IMAGE(ro, NULL);
+    PM_ASSERT_READOUT_NON_NULL(readout, NULL);
+    PM_ASSERT_READOUT_IMAGE(readout, NULL);
     PS_ASSERT_PTR_NON_NULL(config, NULL);
 
-    psImageBinning *binning = backgroundBinning(ro->image, config); // Image binning parameters
+    psImageBinning *binning = psphotBackgroundBinning(readout->image, config); // Image binning parameters
     psImage *model = psImageAlloc(binning->nXruff, binning->nYruff, PS_TYPE_F32); // Background model
     psImage *modelStdev = psImageAlloc(binning->nXruff, binning->nYruff, PS_TYPE_F32); // Standard deviation
 
-    if (!backgroundModel(model, modelStdev, ro->analysis, ro, binning, config)) {
+    if (!psphotModelBackgroundReadout(model, modelStdev, readout->analysis, readout, binning, config)) {
         psFree(model);
         psFree(modelStdev);
@@ -345 +355 @@
         return NULL;
     }
-
     psFree(modelStdev);
-
     return model;
 }
 
-
-// generate the median in NxN boxes, clipping heavily
-// linear interpolation to generate full-scale model
-bool psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filename)
+// generate a background model for readout number index; save an associated pmFPAfile for possible output
+bool psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filename, int index)
 {
-    bool status = true;
-
     // find the currently selected readout
-    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, filename);
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
+    psAssert (file, "missing file?");
+
     pmFPA *inFPA = file->fpa;
     pmReadout *readout = pmFPAviewThisReadout(view, inFPA);
-
-    psImageBinning *binning = backgroundBinning(readout->image, config); // Image binning parameters
-    pmReadout *model = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL", inFPA, binning);
-    pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning);
-
-    if (!backgroundModel(model->image, modelStdev->image, model->analysis, readout, binning, config)) {
+    psAssert (readout, "missing readout?");
+
+    psImageBinning *binning = psphotBackgroundBinning(readout->image, config); // Image binning parameters
+    pmReadout *model = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL", inFPA, binning, index);
+    pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning, index);
+
+    if (!psphotModelBackgroundReadout(model->image, modelStdev->image, model->analysis, readout, binning, config)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
         return false;
@@ -375 +382 @@
     return true;
 }
+
+// XXX supply filename or keep PSPHOT.INPUT fixed?
+bool psphotModelBackground (pmConfig *config, const pmFPAview *view)
+{
+    bool status = false;
+
+    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+        if (!psphotModelBackgroundReadoutFileIndex(config, view, "PSPHOT.INPUT", i)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to model background for PSPHOT.INPUT entry %d", i);
+            return false;
+        }
+    }
+    return true;
+}

@@ -19 +19 @@
 psphotMomentsStudy
 psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack