Changeset 32628

Timestamp:

Nov 8, 2011, 2:41:42 PM (15 years ago)

Author:

eugene

Message:

match sources before final linear fit (to enable forced photometry of PSF mags); only a single detection per image allowed for an object; distinguish conditions for Radial Aperture analysis and Petrosian analysis; no Radial Aperture analysis for chisq image

Location:

branches/eam_branches/ipp-20110906/psphot/src

Files:

• psphotExtendedSourceAnalysis.c (modified) (1 diff)
• psphotRadialApertures.c (modified) (3 diffs)
• psphotSourceMatch.c (modified) (7 diffs)
• psphotStackObjects.c (modified) (5 diffs)
• psphotStackReadout.c (modified) (5 diffs)

branches/eam_branches/ipp-20110906/psphot/src/psphotExtendedSourceAnalysis.c

@@ -209 +209 @@
         // if we have checked the source validity on the basis of the object set, then 
         // we either skip these tests below or we skip the source completely
-        if (source->tmpFlags & PM_SOURCE_TMPF_STACK_SKIP) continue;
-        if (source->tmpFlags & PM_SOURCE_TMPF_STACK_KEEP) goto keepSource;
+        if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_SKIP) continue;
+        if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_KEEP) goto keepSource;
 
         // skip PSF-like and non-astronomical objects

branches/eam_branches/ipp-20110906/psphot/src/psphotRadialApertures.c

@@ -26 +26 @@
     int num = psphotFileruleCount(config, filerule);
 
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
+
         if (!psphotRadialAperturesReadout (config, view, filerule, i, recipe, entry)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for %s entry %d", filerule, i);
@@ -218 +224 @@
         pmSource *source = sources->data[i];
 
+        // if we have checked the source validity on the basis of the object set, then 
+        // we either skip these tests below or we skip the source completely
+        if (source->tmpFlags & PM_SOURCE_TMPF_RADIAL_SKIP) continue;
+        if (source->tmpFlags & PM_SOURCE_TMPF_RADIAL_KEEP) goto keepSource;
+
         // skip PSF-like and non-astronomical objects
         if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
@@ -235 +246 @@
         if (source->peak->x > region->x1) continue;
         if (source->peak->y > region->y1) continue;
+
+    keepSource:
 
         // allocate pmSourceExtendedParameters, if not already defined

branches/eam_branches/ipp-20110906/psphot/src/psphotSourceMatch.c

@@ -67 +67 @@
     objects = psArraySort (objects, pmPhotObjSortByX); 
  
-    psVector *found = psVectorAlloc(sources->n, PS_TYPE_U8);
-    psVectorInit (found, 0);
+    psVector *foundSrc = psVectorAlloc(sources->n, PS_TYPE_U8);
+    psVectorInit (foundSrc, 0);
+
+    psVector *foundObj = psVectorAlloc(sources->n, PS_TYPE_U8);
+    psVectorInit (foundObj, 0);
 
     // match sources to existing objects
@@ -93 +96 @@
         if (dx > +1.02*RADIUS) NEXT2; 
  
+        /* this block will match a given detection to the closest object within range of that detection.
+           XXX note that this matches ALL detections within range of the single object to that same object 
+           this is bad, but I cannot just go in linear order (ie, mark off each object as they are
+           used).  I should make a list of all Nobj * Ndet pairs in range and choose the matches
+           based on their separations.  UGH
+        */
+    
         // we are within match range, look for matches: 
         int Jmin = -1;
@@ -98 +108 @@
         for (int J = j; (dx > -1.02*RADIUS) && (J < objects->n); J++) { 
  
+            // skip objects that are already assigned:
+            if (foundObj->data.U8[J]) continue;
             obj = objects->data[J]; 
             
@@ -117 +129 @@
         }
         obj = objects->data[Jmin]; 
+        foundObj->data.U8[Jmin] = 1;
 
         // add to object
         pmPhotObjAddSource (obj, src);
-        found->data.U8[i] = 1;
+        foundSrc->data.U8[i] = 1;
         i++; 
     } 
@@ -128 +141 @@
     for (i = 0; i < sources->n; i++) {
 
-        if (found->data.U8[i]) continue;
+        if (foundSrc->data.U8[i]) continue;
 
         pmSource *src = sources->data[i]; 
@@ -139 +152 @@
     psLogMsg ("psphot", PS_LOG_DETAIL, "matched sources (%ld vs %ld)", sources->n, objects->n);
 
-    psFree (found);
+    psFree (foundSrc);
+    psFree (foundObj);
     return true;
 } 
@@ -270 +284 @@
         pmPhotObj *obj = objects->data[i]; 
         nSources += obj->sources->n;
+        psAssert (obj->sources->n == nImages, "failed to match sources?");
     }
     psLogMsg ("psphot", PS_LOG_DETAIL, "total of %d sources for %d images", nSources, nImages);

branches/eam_branches/ipp-20110906/psphot/src/psphotStackObjects.c

@@ -70 +70 @@
 
     // S/N limit to perform full non-linear fits
-    float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");
+    float SN_LIM_PETRO  = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");
+    float SN_LIM_RADIAL = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
 
     bool doPetroStars   = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS");
@@ -81 +82 @@
         // we check each source for an object and keep the object if any source is valid
 
-        bool keepObject = false;
+        bool keepObjectRadial = false;
+        bool keepObjectPetro = false;
         for (int j = 0; j < object->sources->n; j++) {
 
@@ -94 +96 @@
             if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
             
-        // optionally allow non-extended objects to get petrosians as well
-        if (!doPetroStars) {
-            if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
-            if (source->type == PM_SOURCE_TYPE_STAR) continue;
-        }
-
-            // limit selection to some SN limit
-            // assert (source->peak); // how can a source not have a peak?
-            // limit selection to some SN limit
-            bool skipSource = false;
-            if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
-                skipSource = (source->moments->KronFlux < SN_LIM * source->moments->KronFluxErr);
-            } else {
-                skipSource = (sqrt(source->peak->detValue) < SN_LIM);
-            }
-            if (skipSource) continue;
-
             // limit selection by analysis region (this automatically apply
             if (source->peak->x < AnalysisRegion.x0) continue;
@@ -117 +102 @@
             if (source->peak->y > AnalysisRegion.y1) continue;
             
-            keepObject = true;
+            // SN limit tests for RADIAL APERTURES:
+            bool skipSourceRadial = false;
+            if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
+                skipSourceRadial = (source->moments->KronFlux < SN_LIM_RADIAL * source->moments->KronFluxErr);
+            } else {
+                skipSourceRadial = (sqrt(source->peak->detValue) < SN_LIM_RADIAL);
+            }
+            if (!skipSourceRadial) keepObjectRadial = true;
+
+            // SN limit tests for PETRO
+            bool skipSourcePetro = false;
+            if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
+                skipSourcePetro = (source->moments->KronFlux < SN_LIM_PETRO * source->moments->KronFluxErr);
+            } else {
+                skipSourcePetro = doPetroStars ? (sqrt(source->peak->detValue) < SN_LIM_PETRO) : true;
+            }
+            if (!skipSourcePetro) keepObjectPetro = true;
+
+            keepObjectPetro = true;
         }
 
@@ -128 +131 @@
             // avoid the single-detection tests
 
-            if (keepObject) {
-                source->tmpFlags |=  PM_SOURCE_TMPF_STACK_KEEP;
-                source->tmpFlags &= ~PM_SOURCE_TMPF_STACK_SKIP;
+            if (keepObjectPetro) {
+                source->tmpFlags |=  PM_SOURCE_TMPF_PETRO_KEEP;
+                source->tmpFlags &= ~PM_SOURCE_TMPF_PETRO_SKIP;
             } else {
-                source->tmpFlags |=  PM_SOURCE_TMPF_STACK_SKIP;
-                source->tmpFlags &= ~PM_SOURCE_TMPF_STACK_KEEP;
+                source->tmpFlags |=  PM_SOURCE_TMPF_PETRO_SKIP;
+                source->tmpFlags &= ~PM_SOURCE_TMPF_PETRO_KEEP;
+            }       
+
+            if (keepObjectRadial) {
+                source->tmpFlags |=  PM_SOURCE_TMPF_RADIAL_KEEP;
+                source->tmpFlags &= ~PM_SOURCE_TMPF_RADIAL_SKIP;
+            } else {
+                source->tmpFlags |=  PM_SOURCE_TMPF_RADIAL_SKIP;
+                source->tmpFlags &= ~PM_SOURCE_TMPF_RADIAL_KEEP;
             }       
         }

branches/eam_branches/ipp-20110906/psphot/src/psphotStackReadout.c

@@ -43 +43 @@
 
 bool psphotStackReadout (pmConfig *config, const pmFPAview *view) {
+
+    psArray *objects = NULL; // used below after 'pass1finish' label
 
     // measure the total elapsed time in psphotReadout.  dtime is the elapsed time used jointly
@@ -187 +189 @@
     psphotReplaceAllSources (config, view, STACK_SRC); // pass 1 (detections->allSources)
 
+    // if we only do one pass, skip to extended source analysis
+    if (!strcasecmp (breakPt, "PASS1")) goto pass1finish;
+
     // linear fit to include all sources (subtract again)
     // NOTE : apply to ALL sources (extended + psf)
     psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 2 (detections->allSources)
-
-    // if we only do one pass, skip to extended source analysis
-    if (!strcasecmp (breakPt, "PASS1")) goto pass1finish;
 
     // NOTE: possibly re-measure background model here with objects subtracted / or masked
@@ -241 +243 @@
         // XXX check on free of sources...
         psphotMergeSources (config, view, STACK_SRC); // (detections->newSources + detections->allSources -> detections->allSources)
-
-        // NOTE: apply to ALL sources
-        psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 3 (detections->allSources)
     }
 
 pass1finish:
+
+    // generate the objects (objects unify the sources from the different images) NOTE: could
+    // this just match the detections for the chisq image, and not bother measuring the source
+    // stats in that case...?
+    objects = psphotMatchSources (config, view, STACK_SRC);
+    psMemDump("matchsources");
+
+    psphotStackObjectsUnifyPosition (objects);
+
+    psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
+
+    // NOTE: apply to ALL sources
+    psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 3 (detections->allSources)
 
     // measure the radial profiles to the sky (only measures new objects)
@@ -260 +272 @@
 
     psMemDump("psfstats");
-
-    // XXX *** NOTE : this should happen above before the last psphotFitSourcesLinear in order
-    // to generate psf forced photometry for the undetected sources
-
-    // generate the objects (objects unify the sources from the different images)
-    // XXX this could just match the detections for the chisq image, and not bother measuring the
-    // source stats in that case...
-    psArray *objects = psphotMatchSources (config, view, STACK_SRC);
-    psMemDump("matchsources");
-
-    psphotStackObjectsUnifyPosition (objects);
-
-    psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
 
     // measure elliptical apertures, petrosians (objects sorted by S/N)
@@ -330 +329 @@
     psphotMagnitudes(config, view, STACK_SRC);
 
-    if (0 && !psphotEfficiency(config, view, STACK_DET)) {
+    // XXX NOTE: this function wants to have the PSF of the image, but we (so far) only measure the 
+    // PSF of the SRC image.  can we fake it by generating the PSF for DET as well (up above)? 
+    if (false && !psphotEfficiency(config, view, STACK_DET)) {
         psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
         psErrorClear();