Changeset 27650

Timestamp:

Apr 11, 2010, 2:48:58 PM (16 years ago)

Author:

eugene

Message:

create sources for matched detections found on other images (psphotStack)

Location:

branches/eam_branches/stackphot.20100406/psphot/src

Files:

• . (modified) (1 prop)
• psphotMergeSources.c (modified) (3 diffs)
• psphotSourceMatch.c (modified) (3 diffs)
• psphotSourceStats.c (modified) (10 diffs)

branches/eam_branches/stackphot.20100406/psphot/src

@@ -21 +21 @@
 psphotForced
 psphotMakePSF
+psphotStack

@@ -21 +21 @@
 psphotForced
 psphotMakePSF
+psphotStack

branches/eam_branches/stackphot.20100406/psphot/src/psphotMergeSources.c

@@ -76 +76 @@
     pmDetections *extCMF = NULL;
     psArray *extSourcesTXT = NULL;
-
-    // find the currently selected readout
-    pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // File of interest
+    int index = 0;
+
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest
     psAssert (file, "missing file?");
 
@@ -118 +119 @@
                 psFree (source->modelPSF);
                 source->modelPSF = NULL;
+                source->imageID = index;
 
                 psArrayAdd (detections->newSources, 100, source);
@@ -140 +142 @@
                 psFree (source->modelPSF);
                 source->modelPSF = NULL;
+                source->imageID = index;
 
                 psArrayAdd (detections->newSources, 100, source);

branches/eam_branches/stackphot.20100406/psphot/src/psphotSourceMatch.c

@@ -1 +1 @@
 # include "psphotInternal.h" 
+
+bool psphotMatchSourcesGenerate (pmConfig *config, const pmFPAview *view, psArray *objects);
  
 psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view) 
@@ -19 +21 @@
     }
 
-    // psphotMatchSourcesGenerate (objects, config, view, "PSPHOT.INPUT", i);
+    psphotMatchSourcesGenerate (config, view, objects);
 
     return objects;
@@ -135 +137 @@
         pmPhotObjAddSource(obj, src);
         psArrayAdd (objects, 100, obj);
+        psFree (obj);
     }
     psLogMsg ("psphot", PS_LOG_DETAIL, "matched sources (%ld vs %ld)", sources->n, objects->n);
 
+    psFree (found);
     return true;
 } 
+
+bool psphotMatchSourcesGenerate (pmConfig *config, const pmFPAview *view, psArray *objects) {
+
+    bool status = false;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    // determine properties (sky, moments) of initial sources
+    float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
+    psAssert (status, "missing SKY_OUTER_RADIUS in recipe?");
+
+    int nImages = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // generate look-up arrays for detections and readouts
+    psArray *detArrays = psArrayAlloc(nImages);
+    psArray *readouts = psArrayAlloc(nImages);
+
+    for (int i = 0; i < nImages; i++) {
+
+        // find the currently selected readout
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // 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?");
+
+        detArrays->data[i] = psMemIncrRefCounter(detections);
+        readouts->data[i] = psMemIncrRefCounter(readout);
+    }
+
+    // vector to track if source for an image is found
+    psVector *found = psVectorAlloc(nImages, PS_TYPE_U8);
+
+    for (int i = 0; i < objects->n; i++) { 
+        pmPhotObj *obj = objects->data[i]; 
+
+        // mark the images for which sources have been found
+        psVectorInit (found, 0);
+        for (int j = 0; j < obj->sources->n; j++) {
+
+            pmSource *src = obj->sources->data[j]; 
+            int index = src->imageID;
+            psAssert (index >= 0, "invalid index");
+            psAssert (index < found->n, "invalid index");
+
+            found->data.U8[index] = 1;
+        }
+
+        // generate new sources for the image that are missing
+        for (int index = 0; index < found->n; index++) {
+            if (found->data.U8[index]) continue;
+
+            pmDetections *detections = detArrays->data[index];
+            pmReadout *readout = readouts->data[index];
+            int row0 = readout->image->row0;
+            int col0 = readout->image->col0;
+
+            // XXX the peak type is not really used in psphot
+            // PM_PEAK_LONE is certainly not true, but irrelevant
+            float peakFlux = readout->image->data.F32[(int)(obj->y-row0-0.5)][(int)(obj->x-col0-0.5)];
+            pmPeak *peak = pmPeakAlloc(obj->x, obj->y, peakFlux, PM_PEAK_LONE);
+            peak->flux = peakFlux;
+            peak->SN = 1.0;
+            peak->xf = obj->x;
+            peak->yf = obj->y;
+            peak->dx = NAN;
+            peak->dy = NAN;
+            
+            // XXX assign to a footprint?
+
+            // create a new source
+            pmSource *source = pmSourceAlloc();
+            source->imageID = index;
+
+            // add the peak
+            source->peak = psMemIncrRefCounter(peak);
+
+            // allocate space for moments
+            source->moments = pmMomentsAlloc();
+
+            // allocate image, weight, mask arrays for each peak (square of radius OUTER)
+            pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
+
+            peak->assigned = true;
+            pmPhotObjAddSource(obj, source);
+            psArrayAdd (detections->newSources, 100, source);
+            psFree (source);
+            psFree (peak);
+        }
+    }
+
+    // how many sources do we have now?
+    int nSources = 0;
+    for (int i = 0; i < objects->n; i++) { 
+        pmPhotObj *obj = objects->data[i]; 
+        nSources += obj->sources->n;
+    }
+    psLogMsg ("psphot", PS_LOG_DETAIL, "total of %d sources for %d images", nSources, nImages);
+
+
+    psFree (found);
+    psFree (detArrays);
+    psFree (readouts);
+    return true;
+}

branches/eam_branches/stackphot.20100406/psphot/src/psphotSourceStats.c

@@ -2 +2 @@
 
 // convert detections to sources and measure their basic properties (moments, local sky, sky
-// variance) Note: this function only generates sources for the new peaks (peak->assigned)
+// variance) Note: this function only generates sources for the new peaks (peak->assigned).
+// The new sources are added to any existing sources on detections->newSources.  The sources 
+// on detections->allSources are ignored.
 bool psphotSourceStats (pmConfig *config, const pmFPAview *view, bool setWindow)
 {
@@ -40 +42 @@
     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     psAssert (detections, "missing detections?");
-    psAssert (!detections->newSources, "new sources already defined?");
 
     // XXX TEST:
@@ -84 +85 @@
 
     // generate the array of sources, define the associated pixel
-    sources = psArrayAllocEmpty (peaks->n);
+    if (!detections->newSources) {
+        detections->newSources = psArrayAllocEmpty (peaks->n);
+    }
+    sources = detections->newSources;
 
     // if there are no peaks, we save the empty source array and return
     if (!peaks->n) {
-        // save the new sources on the detection structure:
-        detections->newSources = sources;
         return true;
     }
@@ -100 +102 @@
         // create a new source
         pmSource *source = pmSourceAlloc();
+        source->imageID = index;
 
         // add the peak
@@ -119 +122 @@
         psLogMsg ("psphot", PS_LOG_INFO, "break point PEAKS, skipping MOMENTS\n");
         psphotVisualShowMoments (sources);
-        detections->newSources = sources;
         return true;
     }
@@ -126 +128 @@
         if (!psphotSetMomentsWindow(recipe, readout->analysis, sources)) {
             psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
-            psFree(sources);
+            psFree(detections->newSources);
             return false;
         }
@@ -178 +180 @@
                 psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS");
                 psFree (job);
-                psFree(sources);
+                psFree(detections->newSources);
                 return false;
             }
@@ -187 +189 @@
         if (!psThreadPoolWait (false)) {
             psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS");
-            psFree(sources);
+            psFree(detections->newSources);
             return false;
         }
@@ -215 +217 @@
     psphotVisualShowMoments (sources);
 
-    // save the new sources on the detection structure:
-    detections->newSources = sources;
-
-
     if (detections->allSources) {
         psphotMaskCosmicRayFootprintCheck(detections->allSources);
@@ -377 +375 @@
         pmSource *source = sources->data[i];
 
+        if (source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED) continue;
+        source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED;
+
         // skip faint sources for moments measurement
         if (source->peak->SN < MIN_SN) {

@@ -21 +21 @@
 psphotForced
 psphotMakePSF
+psphotStack