Changeset 33086 for branches/eam_branches/ipp-20111122/psphot/src/psphotFitSourcesLinear.c

Timestamp:

Jan 11, 2012, 6:57:13 AM (15 years ago)

Author:

eugene

Message:

finish basic coding for model variance weighting

File:

• branches/eam_branches/ipp-20111122/psphot/src/psphotFitSourcesLinear.c (modified) (13 diffs)

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

@@ -12 +12 @@
 static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal);
 
+bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, int index, psMetadata *recipe, pmReadout *readout, psArray *sources);
+pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe);
+bool psphotFreeModelVariance (pmReadout *readout, psArray *sources);
+
 // 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)
@@ -29 +33 @@
         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;
 
@@ -66 +73 @@
         // re-run the fit against that variance.
         if (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) {
-            // generate the model variance image
-            if (!psphotGenerateModelVariance (recipe, readout, sources, psf, final, fitVarMode)) {
+            // generate the model variance image & source pointers
+            if (!psphotGenerateModelVariance (config, view, 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
@@ -77 +87 @@
                 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;
+            }
         }
 
@@ -86 +102 @@
 }
 
+// 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) {
 
@@ -96 +115 @@
             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;
+        pmSourceFitVarMode fitVarMode = CONSTANT_PHOTOMETRIC_WEIGHTS ? PM_SOURCE_PHOTFIT_CONST : PM_SOURCE_PHOTFIT_IMAGE_VAR;
         return fitVarMode;
     } 
@@ -109 +128 @@
     }
     psError (PSPHOT_ERR_CONFIG, false, "Invalid value for LINEAR_FIT_VARIANCE_MODE (%s)", fitVarModeString);
-    return PM_SOURCE_PHOTFIT_MODEL_NONE;
+    return PM_SOURCE_PHOTFIT_NONE;
 }
 
@@ -182 +201 @@
         // psAssert (source->peak, "source without peak??");
         // psAssert (source->peak->footprint, "peak without footprint??");
-
-        // XXX TEST
-        if (source->mode2 & PM_SOURCE_MODE2_MATCHED) {
-            fprintf (stderr, "forced photometry on supplied source\n");
-        }
 
         // turn this bit off and turn it on again if we pass this test
@@ -381 +395 @@
 
     // set the sky, sky_x, sky_y components of border matrix
-    SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER, markVal);
+    SetBorderMatrixElements (border, readout, fitSources, (fitVarMode == PM_SOURCE_PHOTFIT_CONST), SKY_FIT_ORDER, markVal);
 
     psSparseConstraint constraint;
@@ -540 +554 @@
 }
 
-bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
+bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, 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 pmFPAfile for the readout (background model is saved on PSPHOT.BACKMDL regardless of 'filename')
-    pmFPAfile *backModellFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL", index); // File of interest
+    // find the pmFPAfile for the readout (background model is saved on PSPHOT.BACKMDL regardless of 'filerule')
+    pmFPAfile *backModelFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL", index); // File of interest
     assert (backModelFile);
 
@@ -554 +582 @@
 
     // find the binning information
-    psImageBinning *backBinning = psMetadataLookupPtr(&status, model->analysis, "PSPHOT.BACKGROUND.BINNING");
+    psImageBinning *backBinning = psMetadataLookupPtr(&status, backModel->analysis, "PSPHOT.BACKGROUND.BINNING");
     assert (backBinning);
 
@@ -572 +600 @@
         pmSource *source = sources->data[i];
 
-        // XXX which sources need to be injected?
-        // skip non-astronomical objects (very likely defects)
-        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
-        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
-      
-        // do not include CRs in the full ensemble fit
-        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
-        
-        // do not include MOMENTS_FAILURES in the fit
-        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
-
-        // XXX need to identify the region appropriate for the source
+        // 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");
         source->modelVar = psImageSubset(readout->image, region);
 
         // add the source model to the model variance image
-        pmSourceAdd (source, PM_MODEL_OP_XXX, maskVal);
+        pmSourceAdd (source, PM_MODEL_OP_MODELVAR, maskVal);
     }
 
@@ -599 +621 @@
     // XXX for a test:
     psphotSaveImage (NULL, modelVar, "model.var.fits");
+    psphotSaveImage (NULL, readout->image, "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 (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
+bool psphotFreeModelVariance (pmReadout *readout, psArray *sources) {
+
+    bool status = false;
 
     // find the binning information
-    psImage *modelVar = psMetadataLookupPtr(&status, model->analysis, "PSPHOT.MODEL.VAR");
+    psImage *modelVar = psMetadataLookupPtr(&status, readout->analysis, "PSPHOT.MODEL.VAR");
     assert (modelVar);