Changeset 34053 for branches/eam_branches/ipp-20120601/psphot/src

Timestamp:

Jun 22, 2012, 3:33:32 PM (14 years ago)

Author:

eugene

Message:

read the time table in psTime.c; updates to the variance model tests; better trace options for output masks; convert to new PM_SOURCE_PHOTFIT version of weighting

Location:

branches/eam_branches/ipp-20120601/psphot/src

Files:

• psphot.h (modified) (2 diffs)
• psphotEfficiency.c (modified) (1 diff)
• psphotFitSourcesLinear.c (modified) (21 diffs)
• psphotMaskReadout.c (modified) (2 diffs)

branches/eam_branches/ipp-20120601/psphot/src/psphot.h

@@ -61 +61 @@
 bool            psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
 
+bool            psphotUpdateVariance (pmConfig *config, const pmFPAview *view, const char *filerule);
+bool            psphotUpdateVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
+
 bool            psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filerule);
 bool            psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
@@ -99 +102 @@
 
 bool            psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final);
-
-# if (HAVE_MODEL_VAR)
 bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode);
-# else
-bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final);
-# endif
 
 bool            psphotSourceSize (pmConfig *config, const pmFPAview *view, const char *filerule, bool getPSFsize);

branches/eam_branches/ipp-20120601/psphot/src/psphotEfficiency.c

@@ -420 +420 @@
 
     // psphotFitSourcesLinearReadout subtracts the model fits
-# if (HAVE_MODEL_VAR)
     if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true, PM_SOURCE_PHOTFIT_CONST)) {
-# else
-    if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) {
-# endif
         psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources.");
         psFree(fakeSources);

branches/eam_branches/ipp-20120601/psphot/src/psphotFitSourcesLinear.c

@@ -12 +12 @@
 static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal);
 
-# if (HAVE_MODEL_VAR)
 bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources);
 pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe);
 bool psphotFreeModelVariance (pmReadout *readout, psArray *sources);
-# endif
 
 // for now, let's store the detections on the readout->analysis for each readout
@@ -30 +28 @@
     assert (recipe);
 
-# if (HAVE_MODEL_VAR)
     pmSourceFitVarMode fitVarMode = psphotGetFitVarMode (recipe);
     if (!fitVarMode) {
@@ -40 +37 @@
     // do a single pass.
     pmSourceFitVarMode fitVarModePass1 = (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) ? PM_SOURCE_PHOTFIT_CONST : fitVarMode;
-# endif
 
     int num = psphotFileruleCount(config, filerule);
@@ -68 +64 @@
         psAssert (psf, "missing psf?");
 
-# if (HAVE_MODEL_VAR)
-        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1))
-# else
-        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final))
-# endif
-        {
+        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
             return false;
         }
 
-# if (HAVE_MODEL_VAR)
         // the MODEL_VAR weighting scheme requires knowledge of the model fluxes to generate the variance
         // after we have determined the initial set of fits, then we can generate the variance image and 
@@ -104 +94 @@
             }
         }
-# endif
 
         psphotVisualShowResidualImage (readout, (num > 0)); 
@@ -113 +102 @@
 }
 
-# if (HAVE_MODEL_VAR)
 // look up the fit variance mode from the recipe; older recipes do not have the value
 // 'LINEAR_FIT_VARIANCE_MODE'; in those cases, look for 'CONSTANT_PHOTOMETRIC_WEIGHTS' as a boolean and
@@ -136 +124 @@
         return PM_SOURCE_PHOTFIT_IMAGE_VAR;
     }
+    if (!strcasecmp(fitVarModeString, "SKY")   || !strcasecmp(fitVarModeString, "MODEL_SKY")) {
+        return PM_SOURCE_PHOTFIT_MODEL_SKY;
+    }
     if (!strcasecmp(fitVarModeString, "MODEL") || !strcasecmp(fitVarModeString, "MODEL_VAR")) {
         return PM_SOURCE_PHOTFIT_MODEL_VAR;
@@ -142 +133 @@
     return PM_SOURCE_PHOTFIT_NONE;
 }
-# endif
-
-# if (HAVE_MODEL_VAR)
-bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode)
-# else
-bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final)
-# endif
-{
+
+bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode) {
     bool status;
     float x;
@@ -185 +170 @@
     if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
 
-# if (!HAVE_MODEL_VAR)
-    bool CONSTANT_PHOTOMETRIC_WEIGHTS =
-        psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
-    if (!status) {
-        psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
-    }
-# endif
     int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER");
     if (!status) {
@@ -319 +297 @@
     psSparseBorder *border = psSparseBorderAlloc (sparse, nBorder);
 
-# if (HAVE_MODEL_VAR)
     // if fitVarMode is MODEL_VAR, then we need to generate the model image variance
     // XXX we have two possibilities here: 
@@ -327 +304 @@
 
     // 2) do a single pass, and use the model guess to define the model variance (but do I trust the Model Guess?)
-# endif
 
     // fill out the sparse matrix elements and border elements (B)
@@ -336 +312 @@
 
         // diagonal elements of the sparse matrix (auto-cross-product)
-# if (HAVE_MODEL_VAR)
         f = pmSourceModelDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
-# else
-        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
-# endif
         psSparseMatrixElement (sparse, i, i, f);
 
-# if (HAVE_MODEL_VAR)
-        // if we have used CONSTANT errors, then we need to calculate the value of the parameter error
-        if (fitVarMode != PM_SOURCE_PHOTFIT_IMAGE_VAR) {
+        // if we have used CONSTANT errors, then we need to re-calculate the value of the parameter error
+        if (fitVarMode == PM_SOURCE_PHOTFIT_CONST) {
             float var = pmSourceModelDotModel (SRCi, SRCi, PM_SOURCE_PHOTFIT_IMAGE_VAR, covarFactor, maskVal);
             errors->data.F32[i] = 1.0 / sqrt(var);
@@ -351 +322 @@
             errors->data.F32[i] = 1.0 / sqrt(f);
         }
-# else
-        // the formal error depends on the weighting scheme
-        if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
-            float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor, maskVal);
-            errors->data.F32[i] = 1.0 / sqrt(var);
-        } else {
-            errors->data.F32[i] = 1.0 / sqrt(f);
-        }
-# endif
 
         // find the image x model value
-# if (HAVE_MODEL_VAR)
         f = pmSourceDataDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
-# else
-        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
-# endif
         psSparseVectorElement (sparse, i, f);
 
@@ -372 +330 @@
         switch (SKY_FIT_ORDER) {
           case 1:
-# if (HAVE_MODEL_VAR)
             f = pmSourceModelWeight (SRCi, 1, fitVarMode, covarFactor, maskVal);
             psSparseBorderElementB (border, i, 1, f);
             f = pmSourceModelWeight (SRCi, 2, fitVarMode, covarFactor, maskVal);
             psSparseBorderElementB (border, i, 2, f);
-# else
-            f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
-            psSparseBorderElementB (border, i, 1, f);
-            f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
-            psSparseBorderElementB (border, i, 2, f);
-# endif
 
           case 0:
-# if (HAVE_MODEL_VAR)
             f = pmSourceModelWeight (SRCi, 0, fitVarMode, covarFactor, maskVal);
-# else
-            f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
-# endif
             psSparseBorderElementB (border, i, 0, f);
             break;
@@ -409 +356 @@
 
             // got an overlap; calculate cross-product and add to output array
-# if (HAVE_MODEL_VAR)
             f = pmSourceModelDotModel (SRCi, SRCj, fitVarMode, covarFactor, maskVal);
-# else
-            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
-# endif
             psSparseMatrixElement (sparse, j, i, f);
         }
@@ -451 +394 @@
 
     // set the sky, sky_x, sky_y components of border matrix
-# if (HAVE_MODEL_VAR)
     SetBorderMatrixElements (border, readout, fitSources, (fitVarMode == PM_SOURCE_PHOTFIT_CONST), SKY_FIT_ORDER, markVal);
-# else
-    SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER, markVal);
-# endif
 
     psSparseConstraint constraint;
@@ -614 +553 @@
 }
 
-# if (HAVE_MODEL_VAR)
 bool psphotModelBackgroundReadout(psImage *model,  // Model image
                                   psImage *modelStdev, // Model stdev image
@@ -650 +588 @@
     psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
 
+    // generate an image of the mean variance image in DN
     if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
@@ -668 +607 @@
     psFree (varModelStdev);
 
-    // XXX for a test:
-    psphotSaveImage (NULL, modelVar, "model.bck.fits");
+    float gain = 1.0;  // accept 1.0 as a default since it is not critical to the analysis
+    pmCell *cell = readout->parent; // The parent cell
+    if (cell) {
+      gain = psMetadataLookupF32(&status, cell->concepts, "CELL.GAIN"); // Cell gain
+      if (!status) {
+        gain = 1.0;           // set note above
+      }
+    }
+    if (gain > 2.0) { /* warn? */ }
 
     // insert all of the source models
@@ -690 +636 @@
 
         // add the source model to the model variance image
+        // XXX note that this should be added with gain applied
+        // var_DN = flux_DN / gain [e/DN]
+        // to do this requires an API upgrade...
         pmSourceAdd (source, PM_MODEL_OP_MODELVAR, maskVal);
     }
-
-    // XXX for a test:
-    psphotSaveImage (NULL, modelVar, "model.var.fits");
-    psphotSaveImage (NULL, readout->variance, "image.var.fits");
 
     // we save the model variance for future reference
@@ -724 +669 @@
     return true;
 }
-# endif

branches/eam_branches/ipp-20120601/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) {
+
+    bool status;
+
+    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;
+}