Changeset 34041 for branches/meh_branches/ppstack_test/psphot/src/psphotFitSourcesLinear.c

Timestamp:

Jun 19, 2012, 5:24:19 PM (14 years ago)

Author:

mhuber

Message:

merging latest r34040 trunk changes to branch

Location:

branches/meh_branches/ppstack_test

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (1 prop)
• psphot/src/psphotFitSourcesLinear.c (modified) (13 diffs)

branches/meh_branches/ppstack_test/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
 bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final)
@@ -24 +30 @@
     assert (recipe);
 
+# if (HAVE_MODEL_VAR)
+    pmSourceFitVarMode fitVarMode = psphotGetFitVarMode (recipe);
+    if (!fitVarMode) {
+        psError (PSPHOT_ERR_CONFIG, true, "failed to get LINEAR_FIT_VARIANCE_MODE");
+        return false;
+    }
+    // MODEL_VAR requires 2 passes -- in the first, we get the rough fluxes; in the second, we
+    // use the flux to define the model variance before fitting the objects.  Other modes only
+    // do a single pass.
+    pmSourceFitVarMode fitVarModePass1 = (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) ? PM_SOURCE_PHOTFIT_CONST : fitVarMode;
+# endif
+
     int num = psphotFileruleCount(config, filerule);
 
@@ -50 +68 @@
         psAssert (psf, "missing psf?");
 
-        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final)) {
+# if (HAVE_MODEL_VAR)
+        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1))
+# else
+        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final))
+# endif
+        {
             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 
+        // re-run the fit against that variance.
+        if (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) {
+            // generate the model variance image & source pointers
+            if (!psphotGenerateModelVariance (config, view, file, i, recipe, readout, sources)) {
+                psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
+                return false;
+            }
+
+            // replace all sources (use TMPF_SUBTRACTED as test)
+            psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe, false);
+
+            // rerun fit with correct fitVarMode
+            if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarMode)) {
+                psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
+                return false;
+            }
+
+            // free the model variance image & source pointers
+            if (!psphotFreeModelVariance (readout, sources)) {
+                psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
+                return false;
+            }
+        }
+# endif
 
         psphotVisualShowResidualImage (readout, (num > 0)); 
@@ -62 +113 @@
 }
 
-bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final) {
-
+# 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
+// set the value to either CONST or IMAGE_VAR
+pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe) {
+
+    bool status = false;
+
+    char *fitVarModeString = psMetadataLookupStr(&status, recipe, "LINEAR_FIT_VARIANCE_MODE");
+    if (!status) {
+        bool CONSTANT_PHOTOMETRIC_WEIGHTS = psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
+        if (!status) {
+            psAbort("You must provide a value for LINEAR_FIT_VARIANCE_MODE or CONSTANT_PHOTOMETRIC_WEIGHTS");
+        }
+        pmSourceFitVarMode fitVarMode = CONSTANT_PHOTOMETRIC_WEIGHTS ? PM_SOURCE_PHOTFIT_CONST : PM_SOURCE_PHOTFIT_IMAGE_VAR;
+        return fitVarMode;
+    } 
+    if (!strcasecmp(fitVarModeString, "CONSTANT") || !strcasecmp(fitVarModeString, "CONST")) {
+        return PM_SOURCE_PHOTFIT_CONST;
+    }
+    if (!strcasecmp(fitVarModeString, "IMAGE") || !strcasecmp(fitVarModeString, "IMAGE_VAR")) {
+        return PM_SOURCE_PHOTFIT_IMAGE_VAR;
+    }
+    if (!strcasecmp(fitVarModeString, "MODEL") || !strcasecmp(fitVarModeString, "MODEL_VAR")) {
+        return PM_SOURCE_PHOTFIT_MODEL_VAR;
+    }
+    psError (PSPHOT_ERR_CONFIG, false, "Invalid value for LINEAR_FIT_VARIANCE_MODE (%s)", fitVarModeString);
+    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 status;
     float x;
@@ -99 +185 @@
     if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
 
+# if (!HAVE_MODEL_VAR)
     bool CONSTANT_PHOTOMETRIC_WEIGHTS =
         psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
@@ -104 +191 @@
         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) {
@@ -231 +319 @@
     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: 
+
+    // 1) do 2 passes, where in the first case we use the CONST weighting, and in the second
+    // use the fitted model values to define the model
+
+    // 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)
     // SRCi is the current source of interest
@@ -238 +336 @@
 
         // 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) {
+            float var = pmSourceModelDotModel (SRCi, SRCi, PM_SOURCE_PHOTFIT_IMAGE_VAR, covarFactor, maskVal);
+            errors->data.F32[i] = 1.0 / sqrt(var);
+        } else {
+            errors->data.F32[i] = 1.0 / sqrt(f);
+        }
+# else
         // the formal error depends on the weighting scheme
         if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
@@ -248 +359 @@
             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);
 
@@ -257 +372 @@
         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;
@@ -283 +409 @@
 
             // 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);
         }
@@ -321 +451 @@
 
     // 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;
@@ -479 +613 @@
     return true;
 }
+
+# if (HAVE_MODEL_VAR)
+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
+                                  bool useVarianceImage
+    );
+
+bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources) {
+
+    bool status = false;
+    psRegion fullRegion = psRegionSet (0, 0, 0, 0);
+
+    // bit-masks to test for good/bad pixels
+    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
+    assert (maskVal);
+
+    // bit-mask to mark pixels not used in analysis
+    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
+    assert (markVal);
+
+    // maskVal is used to test for rejected pixels, and must include markVal
+    maskVal |= markVal;
+
+    // create a model variance image
+    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;
+    }
+
+    psFree (varModel);
+    psFree (varModelStdev);
+
+    // XXX for a test:
+    psphotSaveImage (NULL, modelVar, "model.bck.fits");
+
+    // insert all of the source models
+    for (int i = 0; i < sources->n; i++) {
+
+        // source of interest
+        pmSource *source = sources->data[i];
+
+        // skip sources which were not fitted already
+        if (!(source->mode & PM_SOURCE_MODE_LINEAR_FIT)) continue;
+
+        // pixel region appropriate for the source
+        psRegion region = psRegionForImage (source->pixels, fullRegion);
+
+        // define the source->modelVar pixels (view on modelVar image)
+        psAssert (!source->modelVar, "programming error : modelVar should be NULL here");
+        psAssert (source->modelFlux, "programming error : modelFlux should not be NULL here");
+        psAssert (source->modelFlux->data.F32, "programming error : modelFlux should not be NULL here");
+        source->modelVar = psImageSubset(modelVar, region);
+
+        // add the source model to the model variance image
+        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
+    psMetadataAddImage(readout->analysis, PS_LIST_TAIL, "PSPHOT.MODEL.VAR", PS_META_REPLACE, "model variance", modelVar);
+    psFree (modelVar);
+
+    return true;
+}
+
+bool psphotFreeModelVariance (pmReadout *readout, psArray *sources) {
+
+    bool status = false;
+
+    // find the binning information
+    psImage *modelVar = psMetadataLookupPtr(&status, readout->analysis, "PSPHOT.MODEL.VAR");
+    assert (modelVar);
+
+    psMetadataRemoveKey (readout->analysis, "PSPHOT.MODEL.VAR"); 
+
+    // clear modelVar pointers for all of the source models
+    for (int i = 0; i < sources->n; i++) {
+
+        // source of interest
+        pmSource *source = sources->data[i];
+        psFree (source->modelVar);
+    }
+
+    return true;
+}
+# endif