Changeset 27840 for branches/simtest_nebulous_branches/psphot/src/psphotExtendedSourceAnalysis.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (1 prop)
• psphot/src/psphotExtendedSourceAnalysis.c (modified) (8 diffs)

branches/simtest_nebulous_branches/psphot/src

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

branches/simtest_nebulous_branches/psphot/src/psphotExtendedSourceAnalysis.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
+// for now, let's store the detections on the readout->analysis for each readout
+bool psphotExtendedSourceAnalysis (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?");
+
+    // perform full non-linear fits / extended source analysis?
+    if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANALYSIS")) {
+        psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n");
+        return true;
+    }
+
+    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 (!psphotExtendedSourceAnalysisReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for PSPHOT.INPUT entry %d", i);
+            return false;
+        }
+    }
+    return true;
+}
+
 // aperture-like measurements for extended sources
-bool psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe) {
+bool psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
 
     bool status;
@@ -11 +39 @@
     int Nkron = 0;
 
+    psTimerStart ("psphot.extended");
+
+    // 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 source size");
+        return true;
+    }
+
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
     assert (maskVal);
 
-    // perform full non-linear fits / extended source analysis?
-    if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANALYSIS")) {
-        psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n");
-        return true;
-    }
+    // XXX require petrosian analysis for non-linear fits? 
 
-    // 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");
+    // XXX temporary user-supplied systematic sky noise measurement (derive from background model)
+    float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE");
 
     // S/N limit to perform full non-linear fits
@@ -35 +77 @@
     bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
 
+# if (0)
+    // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN
+    // XXX use this to set skynoise
+    pmReadout *backModel = psphotSelectBackground (config, view);
+    pmReadout *backStdev = psphotSelectBackgroundStdev (config, view);
+# endif
+
     // source analysis is done in S/N order (brightest first)
     sources = psArraySort (sources, pmSourceSortBySN);
 
-    // XXX some init functions for the extended source recipe options?
+    // 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");
 
     // choose the sources of interest
@@ -49 +101 @@
         if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
+        if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
+        if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
 
         // limit selection to some SN limit
@@ -66 +121 @@
         Next ++;
 
+        // force source image to be a bit larger...
+        float radius = source->peak->xf - source->pixels->col0;
+        radius = PS_MAX (radius, source->peak->yf - source->pixels->row0);
+        radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0);
+        radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0);
+        pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius);
+
         // if we request any of these measurements, we require the radial profile
         if (doPetrosian || doIsophotal || doAnnuli || doKron) {
-            if (!psphotRadialProfile (source, recipe, maskVal)) {
+            if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) {
                 // all measurements below require the radial profile; skip them all
                 // re-subtract the object, leave local sky
                 psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
                 pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-                source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
                 continue;
             }
@@ -79 +140 @@
         }
 
+        // Petrosian Mags
+        if (doPetrosian) {
+            if (!psphotPetrosian (source, recipe, skynoise, maskVal)) {
+                psTrace ("psphot", 5, "FAILED petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
+            } else {
+                psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
+                Npetro ++;
+                source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
+            }
+        }
+
+# if (0)
         // Isophotal Mags
         if (doIsophotal) {
@@ -89 +162 @@
             }
         }
-
-        // Petrosian Mags
-        if (doPetrosian) {
-            if (!psphotPetrosian (source, recipe, maskVal)) {
-                psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
-            } else {
-                psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
-                Npetro ++;
-                source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
-            }
-        }
-
         // Kron Mags
         if (doKron) {
@@ -111 +172 @@
             }
         }
-
-        // Radial Annuli
-        if (doAnnuli) {
-            if (!psphotAnnuli (source, recipe, maskVal)) {
-                psError(PSPHOT_ERR_UNKNOWN, false, "failure in Annuli analysis");
-                return false;
-            } 
-            psTrace ("psphot", 5, "measured annuli for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
-            Nannuli ++;
-            source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
-        }
+# endif
 
         // re-subtract the object, leave local sky
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+
+        if (source->extpars) {
+            psFree(source->extpars->radFlux);
+            psFree(source->extpars->ellipticalFlux);
+        }
     }
 
-    psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot"), Next);
+    psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d petrosian\n", Npetro);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d isophotal\n", Nisophot);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d kron\n", Nkron);
+
+    psphotVisualShowResidualImage (readout);
+
+    if (doPetrosian) {
+        psphotVisualShowPetrosians (sources);
+    }
+
     return true;
 }

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack