Changeset 27650 for branches/eam_branches/stackphot.20100406/psphot/src/psphotSourceMatch.c

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)
• psphotSourceMatch.c (modified) (3 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/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;
+}

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