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

Timestamp:

Apr 9, 2010, 4:41:16 PM (16 years ago)

Author:

eugene

Message:

working on psphotStack (nearly done)

File:

• branches/eam_branches/stackphot.20100406/psphot/src/psphotSourceMatch.c (modified) (1 diff)

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

@@ -14 +14 @@
         if (!psphotMatchSourcesReadout (objects, config, view, "PSPHOT.INPUT", i)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to merge sources for PSPHOT.INPUT entry %d", i);
-            return false;
+            psFree (objects);
+            return NULL;
         }
     }
+
+    // psphotMatchSourcesGenerate (objects, config, view, "PSPHOT.INPUT", i);
+
+    return objects;
+}
+
+bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, char *filename, int index) { 
+ 
+    bool status = false;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    int RADIUS = psMetadataLookupF32 (&status, recipe, "PSPHOT.STACK.MATCH.RADIUS");
+    psAssert (status, "programming error: must define PSPHOT.STACK.MATCH.RADIUS");
+
+    // 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?");
+    psAssert (detections->newSources, "new sources not defined?");
+    psAssert (!detections->allSources, "all sources already defined?");
+
+    // XXX TEST:
+    if (detections->newSources) {
+        psphotMatchSourcesToObjects(objects, detections->newSources, RADIUS);
+    }
+
     return true;
 }
 
-// we need a couple of functions to distinguish coincident sources: 
-// XXX these are similar (identical?) to the goals of pmSourceMatch.c 
- 
 # define NEXT1 { i++; continue; } 
 # define NEXT2 { j++; continue; } 
- 
-bool psphotSourceMerge (psArray *objects, pmConfig *config, const pmFPAview *view, char *filename, int index) { 
+bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS) { 
  
     float dx, dy; 
  
+    float RADIUS2 = RADIUS*RADIUS;
+
     // sort the source list by X 
-    pmSourceSortByX (sources1); 
-    pmSourceSortByX (sources2); 
+    sources = psArraySort (sources, pmSourceSortByX); 
+    objects = psArraySort (objects, pmPhotObjSortByX); 
  
+    psVector *found = psVectorAlloc(sources->n, PS_TYPE_U8);
+    psVectorInit (found, 0);
+
+    // match sources to existing objects
+
+    psLogMsg ("psphot", PS_LOG_DETAIL, "attempt to match sources (%ld vs %ld)", sources->n, objects->n);
+
     int i, j; 
-    for (i = j = 0; (i < sources1->n) && (j < sources2->n); ) { 
+    for (i = j = 0; (i < sources->n) && (j < objects->n); ) { 
  
-        pmSource *src1 = sources1->data[i]; 
-        pmSource *src2 = sources2->data[j]; 
+        pmSource  *src = sources->data[i]; 
+        pmPhotObj *obj = objects->data[j]; 
  
-        if (!src1) NEXT1; 
-        if (!src1->peak) NEXT1; 
-        if (!finite(src1->peak->xf)) NEXT1; 
-        if (!finite(src1->peak->yf)) NEXT1; 
+        if (!src) NEXT1; 
+        if (!src->peak) NEXT1; 
+        if (!finite(src->peak->xf)) NEXT1; 
+        if (!finite(src->peak->yf)) NEXT1; 
  
-        if (!src2) NEXT2; 
-        if (!src2->peak) NEXT2; 
-        if (!finite(src2->peak->xf)) NEXT2; 
-        if (!finite(src2->peak->yf)) NEXT2; 
+        if (!obj) NEXT2; 
+        if (!finite(obj->x)) NEXT2; 
+        if (!finite(obj->y)) NEXT2; 
  
-        dx = src1->peak->xf - src2->peak->xf; 
+        dx = src->peak->xf - obj->x; 
         if (dx < -1.02*RADIUS) NEXT1; 
         if (dx > +1.02*RADIUS) NEXT2; 
  
         // we are within match range, look for matches: 
-        for (int J = j; (dx > -1.02*radius) && (J < sources2->n); J++) { 
+        int Jmin = -1;
+        float Rmin = RADIUS2;
+        for (int J = j; (dx > -1.02*RADIUS) && (J < objects->n); J++) { 
  
-            dx = src1->peak->xf - src2->peak->xf; 
-            dy = src1->peak->yf - src2->peak->yf; 
+            obj = objects->data[J]; 
+            
+            dx = src->peak->xf - obj->x; 
+            dy = src->peak->yf - obj->y; 
  
-            dr = dx*dx + dy*dy; 
+            float dr = dx*dx + dy*dy; 
             if (dr > RADIUS2) continue; 
- 
-            // add to group? 
-        } 
+            if (dr < Rmin) {
+                Rmin = dr;
+                Jmin  = J;
+            }
+        }
+
+        // no match, try next source
+        if (Jmin == -1) {
+            i++;
+            continue;
+        }
+        obj = objects->data[Jmin]; 
+
+        // add to object
+        pmPhotObjAddSource (obj, src);
+        found->data.U8[i] = 1;
         i++; 
     } 
+
+    // add missed sources to new objects
+
+    for (i = 0; i < sources->n; i++) {
+
+        if (found->data.U8[i]) continue;
+
+        pmSource *src = sources->data[i]; 
+
+        pmPhotObj *obj = pmPhotObjAlloc();
+        pmPhotObjAddSource(obj, src);
+        psArrayAdd (objects, 100, obj);
+    }
+    psLogMsg ("psphot", PS_LOG_DETAIL, "matched sources (%ld vs %ld)", sources->n, objects->n);
+
+    return true;
 }