Changeset 26894 for trunk/psphot/src/psphotBlendFit.c

Timestamp:

Feb 10, 2010, 7:36:29 PM (16 years ago)

Author:

eugene

Message:

updates from eam_branches/20091201 (substantially changes to the psphotReadout APIs to support future stack photometry; improvements to the CR masking code)

File:

• trunk/psphot/src/psphotBlendFit.c (modified) (5 diffs)

trunk/psphot/src/psphotBlendFit.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
+// for now, let's store the detections on the readout->analysis for each readout
+bool psphotBlendFit (pmConfig *config, const pmFPAview *view)
+{
+    bool status = true;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+        if (!psphotBlendFitReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (non-linear) for PSPHOT.INPUT entry %d", i);
+            return false;
+        }
+    }
+    return true;
+}
+
 // XXX I don't like this name
-bool psphotBlendFit (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) {
+bool psphotBlendFitReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
 
     int Nfit = 0;
@@ -12 +34 @@
     psTimerStart ("psphot.fit.nonlinear");
 
-    // select the appropriate recipe information
-    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
-    assert (recipe);
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
+    psAssert (file, "missing file?");
+
+    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+    psAssert (readout, "missing readout?");
+
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    psArray *sources = detections->allSources;
+    psAssert (sources, "missing sources?");
+
+    if (!sources->n) {
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend fit");
+        return true;
+    }
+
+    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
+    psAssert (psf, "missing psf?");
 
     // determine the number of allowed threads
@@ -40 +79 @@
     psphotInitLimitsPSF (recipe, readout);
     psphotInitLimitsEXT (recipe);
-    psphotInitRadiusPSF (recipe, psf->type);
+    psphotInitRadiusPSF (recipe, readout->analysis, psf->type);
 
     // starts the timer, sets up the array of fitSets
@@ -47 +86 @@
     // source analysis is done in S/N order (brightest first)
     sources = psArraySort (sources, pmSourceSortBySN);
+    if (!sources->n) {
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend");
+        return true;
+    }
 
     // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
@@ -272 +315 @@
 }
 
-# if (0)
-bool psphotBlendFit_Unthreaded (int *nfit, int *npsf, int *next, int *nfail, pmReadout *readout, psMetadata *recipe, psArray *sources, pmPSF *psf, psArray *newSources) {
-
-    bool status = false;
-    int Nfit = 0;
-    int Npsf = 0;
-    int Next = 0;
-    int Nfail = 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;
-
-    // S/N limit to perform full non-linear fits
-    float FIT_SN_LIM = psMetadataLookupF32 (&status, recipe, "FULL_FIT_SN_LIM");
-
-    // option to limit analysis to a specific region
-    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
-    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
-    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
-
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-
-        // skip non-astronomical objects (very likely defects)
-        if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
-        if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) continue;
-        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
-        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
-
-        // skip DBL second sources (ie, added by psphotFitBlob)
-        if (source->mode &  PM_SOURCE_MODE_PAIR) continue;
-
-        // limit selection to some SN limit
-        if (source->peak->SN < FIT_SN_LIM) continue;
-
-        // exclude sources outside optional analysis region
-        if (source->peak->xf < AnalysisRegion.x0) continue;
-        if (source->peak->yf < AnalysisRegion.y0) continue;
-        if (source->peak->xf > AnalysisRegion.x1) continue;
-        if (source->peak->yf > AnalysisRegion.y1) continue;
-
-        // if model is NULL, we don't have a starting guess
-        if (source->modelPSF == NULL) continue;
-
-        // skip sources which are insignificant flux?
-        // XXX this is somewhat ad-hoc
-        if (source->modelPSF->params->data.F32[1] < 0.1) {
-            psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n",
-                     source->modelPSF->params->data.F32[1],
-                     source->modelPSF->params->data.F32[2],
-                     source->modelPSF->params->data.F32[3]);
-            continue;
-        }
-
-        // replace object in image
-        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
-            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-        }
-        Nfit ++;
-
-        // try fitting PSFs or extended sources depending on source->mode
-        // these functions subtract the resulting fitted source
-        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
-            if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) {
-                source->type = PM_SOURCE_TYPE_EXTENDED;
-                psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
-                Next ++;
-                continue;
-            }
-        } else {
-            if (psphotFitBlend (readout, source, psf, maskVal, markVal)) {
-                source->type = PM_SOURCE_TYPE_STAR;
-                psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
-                Npsf ++;
-                continue;
-            }
-        }
-
-        psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf);
-        Nfail ++;
-
-        // re-subtract the object, leave local sky
-        pmSourceCacheModel (source, maskVal);
-        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-    }
-
-    // change the value of a scalar on the array (wrap this and put it in psArray.h)
-    *nfit  = Nfit;
-    *npsf  = Npsf;
-    *next  = Next;
-    *nfail = Nfail;
-
-    // moments are modified by the fit; re-display
-    psphotVisualPlotMoments (recipe, sources);
-    psphotVisualShowResidualImage (readout);
-
-    return true;
-}
-# endif