Changeset 21108

Timestamp:

Jan 12, 2009, 4:14:01 PM (17 years ago)

Author:

Paul Price

Message:

Adding API for background model calculation that doesn't require definition of pmFPAfiles.

Location:

trunk/psphot/src

Files:

• psphot.h (modified) (2 diffs)
• psphotDiagnosticPlots.c (modified) (1 diff)
• psphotModelBackground.c (modified) (10 diffs)

trunk/psphot/src/psphot.h

@@ -213 +213 @@
 bool psphotKapaClose ();
 bool psphotImageBackgroundCellHistogram (psVector *values, float mean, float sigma, int ix, int iy);
-bool psphotDiagnosticPlots (pmConfig *config, char *name, ...);
+bool psphotDiagnosticPlots (const pmConfig *config, const char *name, ...);
 
 bool psphotMakeFluxScale (psImage *image, psMetadata *recipe, pmPSF *psf);
@@ -222 +222 @@
 int psphotSupplementStars (psArray *stars, psArray *sources, psImageBinning *binning, int ix, int iy);
 
+
+// Create a background model for a readout
+// Essentially identical to psphotModelBackground, but with a more convenient API for outsiders.
+psImage *psphotBackgroundModel(pmReadout *ro, // Readout for which to generate background model
+                               const pmConfig *config // Configuration
+    );
+
+
 #endif

trunk/psphot/src/psphotDiagnosticPlots.c

@@ -135 +135 @@
 }
 
-bool psphotDiagnosticPlots (pmConfig *config, char *name, ...) {
+bool psphotDiagnosticPlots (const pmConfig *config, const char *name, ...) {
 
     va_list argPtr;

trunk/psphot/src/psphotModelBackground.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 static int npass = 0;
-
-// 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)
+static char *defaultStatsName = "FITTED_MEAN";
+
+// 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
+                                         const pmConfig *config // Configuration
+                                         )
 {
     bool status = true;
-    static char *defaultStatsName = "FITTED_MEAN";
-
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    assert (recipe);
+
+    // I have the fine image size, I know the binning factor, determine the ruff image size
+    psImageBinning *binning = psImageBinningAlloc();
+    binning->nXfine = image->numCols;
+    binning->nYfine = image->numRows;
+    binning->nXbin  = psMetadataLookupS32(&status, recipe, "BACKGROUND.XBIN");
+    binning->nYbin  = psMetadataLookupS32(&status, recipe, "BACKGROUND.YBIN");
+
+    psImageBinningSetRuffSize(binning, PS_IMAGE_BINNING_CENTER);
+    psImageBinningSetSkip(binning, image);
+
+    return binning;
+}
+
+// Generate the background model
+static bool backgroundModel(psImage *model,  // Model image
+                            psImage *modelStdev, // Model stdev image
+                            psImage *image, // Image for which to generate a background model
+                            psImage *mask, // Mask for image
+                            psImageBinning *binning, // Binning parameters
+                            const pmConfig *config // Configuration
+                            )
+{
     psTimerStart ("psphot.background");
 
-    // find the currently selected readout
-    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, filename);
-    pmFPA *inFPA = file->fpa;
-    pmReadout *readout = pmFPAviewThisReadout (view, inFPA);
-    psImage *image = readout->image;
-    psImage *mask  = readout->mask;
+    bool status = true;
 
     // select the appropriate recipe information
@@ -116 +139 @@
     *statsDefaults = *stats;
 
-    // I have the fine image size, I know the binning factor, determine the ruff image size
-    psImageBinning *binning = psImageBinningAlloc();
-    binning->nXfine = image->numCols;
-    binning->nYfine = image->numRows;
-    binning->nXbin  = psMetadataLookupS32 (&status, recipe, "BACKGROUND.XBIN");
-    binning->nYbin  = psMetadataLookupS32 (&status, recipe, "BACKGROUND.YBIN");
-
-    psImageBinningSetRuffSize(binning, PS_IMAGE_BINNING_CENTER);
-    psImageBinningSetSkip(binning, image);
+    // we save the binning structure for use in psphotMagnitudes
     status = psMetadataAddPtr(recipe, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning);
     PS_ASSERT (status, false);
 
-    // we save the binning structure for use in psphotMagnitudes
-    pmReadout *model      = pmFPAGenerateReadout (config, view, "PSPHOT.BACKMDL",       inFPA, binning);
-    pmReadout *modelStdev = pmFPAGenerateReadout (config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning);
-
-    psF32 **modelData = model->image->data.F32;
-    psF32 **modelStdevData = modelStdev->image->data.F32;
+
+    psF32 **modelData = model->data.F32;
+    psF32 **modelStdevData = modelStdev->data.F32;
 
     // XXXX we can thread this here by running blocks in parallel
@@ -140 +152 @@
     psRegion ruffRegion = {0,0,0,0};
     psRegion fineRegion = {0,0,0,0};
-    for (int iy = 0; iy < model->image->numRows; iy++) {
-        for (int ix = 0; ix < model->image->numCols; ix++) {
+    for (int iy = 0; iy < model->numRows; iy++) {
+        for (int ix = 0; ix < model->numCols; ix++) {
 
             // convert the ruff grid cell to the equivalent fine grid cell
@@ -229 +241 @@
         char name[256];
         sprintf (name, "backraw.%02d.fits", npass);
-        psphotSaveImage (NULL, model->image, name);
+        psphotSaveImage (NULL, model, name);
     }
 
@@ -237 +249 @@
     double Value = 0;                   // sum of good pixel's value
     double ValueStdev = 0;              // sum of good pixel's standard deviations
-    for (int iy = 0; iy < model->image->numRows; iy++) {
-        for (int ix = 0; ix < model->image->numCols; ix++) {
+    for (int iy = 0; iy < model->numRows; iy++) {
+        for (int ix = 0; ix < model->numCols; ix++) {
             if (!isnan(modelData[iy][ix])) {
                 Value += modelData[iy][ix];
@@ -250 +262 @@
             for (int jy = iy - 1; jy <= iy + 1; jy++) {
                 if (jy <   0) continue;
-                if (jy >= model->image->numRows) continue;
+                if (jy >= model->numRows) continue;
                 for (int jx = ix - 1; jx <= ix + 1; jx++) {
                     if (!jx && !jy) continue;
                     if (jx   <   0) continue;
-                    if (jx   >= model->image->numCols) continue;
+                    if (jx   >= model->numCols) continue;
                     value += modelData[jy][jx];
                     count += 1.0;
@@ -275 +287 @@
 
     // patch over remaining bad regions (use global average)
-    for (int iy = 0; iy < model->image->numRows; iy++) {
-        for (int ix = 0; ix < model->image->numCols; ix++) {
+    for (int iy = 0; iy < model->numRows; iy++) {
+        for (int ix = 0; ix < model->numCols; ix++) {
             if (!isnan(modelData[iy][ix])) continue;
             modelData[iy][ix] = Value;
@@ -294 +306 @@
                                       PS_STAT_MIN |
                                       PS_STAT_MAX);
-    psImageStats (statsBck, model->image, NULL, 0);
+    psImageStats (statsBck, model, NULL, 0);
     psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MN",
                       PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
@@ -304 +316 @@
                       PS_META_REPLACE, "sky model minimum value", statsBck->min);
     psMetadataAddS32 (recipe, PS_LIST_TAIL, "MSKY_NX",
-                      PS_META_REPLACE, "sky model size (x)",      model->image->numCols);
+                      PS_META_REPLACE, "sky model size (x)",      model->numCols);
     psMetadataAddS32 (recipe, PS_LIST_TAIL, "MSKY_NY",
-                      PS_META_REPLACE, "sky model size (y)",      model->image->numRows);
+                      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);
@@ -316 +328 @@
     psFree(rng);
 
+    return model;
+}
+
+
+
+psImage *psphotBackgroundModel(pmReadout *ro, const pmConfig *config)
+{
+    PM_ASSERT_READOUT_NON_NULL(ro, NULL);
+    PM_ASSERT_READOUT_IMAGE(ro, NULL);
+    PS_ASSERT_PTR_NON_NULL(config, NULL);
+
+    psImageBinning *binning = backgroundBinning(ro->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->image, ro->mask, binning, config)) {
+        psFree(model);
+        psFree(modelStdev);
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
+        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)
+{
+    bool status = true;
+
+    // find the currently selected readout
+    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, filename);
+    pmFPA *inFPA = file->fpa;
+    pmReadout *readout = pmFPAviewThisReadout (view, inFPA);
+    psImage *image = readout->image;
+    psImage *mask  = readout->mask;
+
+    psImageBinning *binning = backgroundBinning(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, image, mask, binning, config)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
+        return false;
+    }
+
     npass ++;
     return true;