Changeset 34086 for trunk/psphot/src/psphotMaskReadout.c

Timestamp:

Jun 26, 2012, 11:33:10 AM (14 years ago)

Author:

eugene

Message:

re-enable MODEL_VAR option for linear photometry fit

Location:

trunk/psphot

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/psphotMaskReadout.c (modified) (2 diffs)

trunk/psphot/src/psphotMaskReadout.c

@@ -94 +94 @@
 
     // test output of files at this stage
-    if (psTraceGetLevel("psphot") >= 5) {
+    if (psTraceGetLevel("psphot.imsave") >= 5) {
         psphotSaveImage (NULL, readout->image,  "image.fits");
         psphotSaveImage (NULL, readout->mask,   "mask.fits");
@@ -105 +105 @@
     return true;
 }
+
+// XXX this function and support below was created to test the theory that the faint-end
+// bias results from the Poisson variation of the background pixels.  This is NOT the
+// case.  Using the code below maintains the faint-end bias.
+bool psphotUpdateVariance (pmConfig *config, const pmFPAview *view, const char *filerule) {
+
+    bool status = false;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    int num = psphotFileruleCount(config, filerule);
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+
+        // Generate the mask and weight images, including the user-defined analysis region of interest
+        if (!psphotUpdateVarianceReadout (config, view, filerule, i, recipe)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to generate mask for %s entry %d", filerule, i);
+            return false;
+        }
+    }
+    return true;
+}
+
+// determine the mean variance image (equivalent to the background model, but for the variance image)
+// set the variance image to the MAX(input, mean)
+bool psphotUpdateVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
+
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
+    psAssert (file, "missing file?");
+
+    // find the currently selected readout
+    pmReadout  *readout = pmFPAviewThisReadout (view, file->fpa);
+    psAssert (readout, "missing readout?");
+
+    pmSourceFitVarMode varMode = psphotGetFitVarMode (recipe);
+    if (varMode == PM_SOURCE_PHOTFIT_NONE) {
+      psError (PSPHOT_ERR_CONFIG, false, "need valid LINEAR_FIT_VARIANCE_MODE");
+      return false;
+    }
+
+    // make this an option via the recipe
+    if (varMode != PM_SOURCE_PHOTFIT_MODEL_SKY) return true;
+
+    // create a model variance image (full-scale image to take result of psImageUnbin below)
+    psImage *modelVar = psImageCopy (NULL, readout->variance, PS_TYPE_F32);
+
+    // find the binning information
+    psImageBinning *backBinning = psphotBackgroundBinning (modelVar, config);
+    assert (backBinning);
+    
+    psImage *varModel = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
+    psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
+
+    if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
+        psFree (varModel);
+        psFree (varModelStdev);
+        return false;
+    }
+
+    // linear interpolation to full-scale
+    if (!psImageUnbin (modelVar, varModel, backBinning)) {
+        psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning");
+        psFree (varModel);
+        psFree (varModelStdev);
+        return false;
+    }
+
+    // XXX save these?
+    psFree (varModel);
+    psFree (varModelStdev);
+
+    psImage *im = readout->image;
+    psImage *wt = readout->variance;
+    for (int j = 0; j < im->numRows; j++) {
+      for (int i = 0; i < im->numCols; i++) {
+        if (!isfinite(im->data.F32[j][i])) continue;
+        if (!isfinite(wt->data.F32[j][i])) continue;
+        // XXX for a test, make variance constant wt->data.F32[j][i] = PS_MAX(wt->data.F32[j][i], modelVar->data.F32[j][i]);
+        wt->data.F32[j][i] = modelVar->data.F32[j][i];
+      }
+    }
+
+    // test output of files at this stage
+    if (psTraceGetLevel("psphot.imsave") >= 5) {
+        psphotSaveImage (NULL, readout->image,  "image.varsky.fits");
+        psphotSaveImage (NULL, readout->mask,   "mask.varsky.fits");
+        psphotSaveImage (NULL, readout->variance, "variance.varsky.fits");
+    }
+
+    psFree (modelVar);
+
+    return true;
+}