Changeset 34399 for branches/eam_branches/ipp-20120805/psphot/src/psphotSourceMatch.c

Timestamp:

Sep 5, 2012, 4:04:34 PM (14 years ago)

Author:

eugene

Message:

merge changes from trunk

Location:

branches/eam_branches/ipp-20120805/psphot

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/psphotSourceMatch.c (modified) (13 diffs)

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

@@ -1 +1 @@
 # include "psphotInternal.h" 
 
-bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
+// Structure for storing the contents of the PSPHOT.STACK.MATCH.FILTERS list from recipe
+#define MATCH_INFO_FILTER_STR_LEN 16
+typedef struct {
+    int     inputNum;
+    int     order;
+    bool    matchAll;
+    char    filterID[MATCH_INFO_FILTER_STR_LEN];
+} psphotStackMatchInfo;
+
+static psphotStackMatchInfo * psphotStackGetMatchInfo (pmConfig *config, const pmFPAview *view, const char *filerule);
+
+// functions for sorting the match info array
+static int compareMatchInfoByInputNum(const void *a, const void *b);
+static int compareMatchInfoByOrder(const void *a, const void *b);
+
+bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects, psphotStackMatchInfo *matchInfo);
 bool psphotMatchSourcesSetIDs (psArray *objects);
 
@@ -14 +29 @@
     int num = psphotFileruleCount(config, filerule);
 
-    // loop over the available readouts
-    for (int i = 0; i < num; i++) {
+    psphotStackMatchInfo *matchInfo = psphotStackGetMatchInfo(config, view, filerule);
+    
+    // loop over the available readouts matching sources found to objects. The inputs are processed
+    // in the order specified by the recipe
+    for (int j = 0; j < num; j++) {
+        int i = matchInfo[j].inputNum;
         if (!psphotMatchSourcesReadout (objects, config, view, filerule, i)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to merge sources for %s entry %d", filerule, i);
@@ -23 +42 @@
     }
 
+    // Now re-order the matchInfo array by input number so we can find each inputs entry easily
+    qsort(matchInfo, num, sizeof(psphotStackMatchInfo), compareMatchInfoByInputNum);
+
     // create sources for images where an object has been detected in the other images
-    psphotMatchSourcesAddMissing (config, view, filerule, objects);
+    psphotMatchSourcesAddMissing (config, view, filerule, objects, matchInfo);
 
     // choose a consistent position; set common sequence values
     psphotMatchSourcesSetIDs (objects);
 
+    psFree(matchInfo);
+
     return objects;
 }
 
-bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, const char *filerule, int index) { 
+bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, const char *filerule, int index) {
  
     bool status = false;
@@ -61 +85 @@
 # define NEXT1 { i++; continue; } 
 # define NEXT2 { j++; continue; } 
-bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS) { 
+bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS) {
  
     float dx, dy; 
@@ -70 +94 @@
     sources = psArraySort (sources, pmSourceSortByX); 
     objects = psArraySort (objects, pmPhotObjSortByX); 
- 
+
     psVector *foundSrc = psVectorAlloc(sources->n, PS_TYPE_U8);
     psVectorInit (foundSrc, 0);
@@ -87 +111 @@
         pmPhotObj *obj = objects->data[j]; 
  
-        if (!src) NEXT1; 
+        if (!src) NEXT1;
         if (!src->peak) NEXT1; 
         if (!isfinite(src->peak->xf)) NEXT1; 
         if (!isfinite(src->peak->yf)) NEXT1; 
- 
+
         if (!obj) NEXT2; 
         if (!isfinite(obj->x)) NEXT2; 
@@ -141 +165 @@
     } 
 
-    // add missed sources to new objects
-
+    // create new objects for unmatched sources
     for (i = 0; i < sources->n; i++) {
 
-        if (foundSrc->data.U8[i]) continue;
+        if (foundSrc->data.U8[i]) continue;
 
         pmSource *src = sources->data[i]; 
 
-        pmPhotObj *obj = pmPhotObjAlloc();
-        pmPhotObjAddSource(obj, src);
-        psArrayAdd (objects, 100, obj);
-        psFree (obj);
+        pmPhotObj *obj = pmPhotObjAlloc();
+        pmPhotObjAddSource(obj, src);
+        psArrayAdd (objects, 100, obj);
+        psFree (obj);
     }
     psLogMsg ("psphot", PS_LOG_DETAIL, "matched sources (%ld vs %ld)", sources->n, objects->n);
@@ -161 +184 @@
 } 
 
-bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects) {
+bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects, psphotStackMatchInfo *matchInfo) {
 
     bool status = false;
@@ -169 +192 @@
     psAssert (recipe, "missing recipe?");
 
+    int minDetectionsForForced = psMetadataLookupS32 (&status, recipe, "PSPHOT.STACK.MIN.DETECT.FOR.FORCED");
+    if (!status) {
+        minDetectionsForForced = 2;
+    }
+
     // determine properties (sky, moments) of initial sources
     float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
@@ -205 +233 @@
         int nSpansMax = 0;
         int iSpansMax = -1;
+
+        bool matchAll = false;
 
         // mark the images for which sources have been found
@@ -219 +249 @@
                 iSpansMax = j;
             }
+            // If this detection was on an input that has the matchAll flag set we create matched sources
+            // irespective of the number of inputs in which the object was found.
+            if (matchInfo[index].matchAll) {
+                matchAll = true;
+            }
 
             found->data.U8[index] = 1;
         }
+
+        // skip adding matched sources for this object if the number of detections for less than
+        // the supplied mininum unless one of the detections was in a band for which the matchAll flag is set
+        if (!matchAll && obj->sources->n < minDetectionsForForced) {
+            continue;
+        }
 
         // we make a copy of the largest footprint; this will be used for all new sources associated with this object
@@ -299 +340 @@
         pmPhotObj *obj = objects->data[i]; 
         nSources += obj->sources->n;
-        psAssert (obj->sources->n == nImages, "failed to match sources?");
+        if (minDetectionsForForced <= 1) {
+            psAssert (obj->sources->n == nImages, "failed to match sources?");
+        }
     }
     psLogMsg ("psphot", PS_LOG_DETAIL, "total of %d sources for %d images", nSources, nImages);
@@ -375 +418 @@
     return copy;
 }
+
+// qsort function to sort match info by order
+static int compareMatchInfoByOrder(const void *a, const void *b) {
+    psphotStackMatchInfo *infoA = (psphotStackMatchInfo *) a;
+    psphotStackMatchInfo *infoB = (psphotStackMatchInfo *) b;
+
+    int     result;
+    if (infoA->order < infoB->order) {
+        result = -1;
+    } else if (infoA->order == infoB->order) {
+        result = 0;
+    } else {
+        result = 1;
+    }
+    return result;
+}
+
+// qsort function to sort match info by input number
+static int compareMatchInfoByInputNum(const void *a, const void *b) {
+    psphotStackMatchInfo *infoA = (psphotStackMatchInfo *) a;
+    psphotStackMatchInfo *infoB = (psphotStackMatchInfo *) b;
+
+    int     result;
+    if (infoA->inputNum < infoB->inputNum) {
+        result = -1;
+    } else if (infoA->inputNum == infoB->inputNum) {
+        result = 0;
+    } else {
+        result = 1;
+    }
+    return result;
+}
+
+// psphotStackGetMatchInfo
+// process the recipe PSPHOT.STACK.MATCH.FILTERS which controls the order that the inputs
+// are processed for source matching and the rules for creating matched sources for sources that
+// are only detected in a single band.
+static psphotStackMatchInfo *psphotStackGetMatchInfo (pmConfig *config, const pmFPAview *view, const char *filerule) {
+
+    bool status = false;
+    psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+
+    psMetadata *filterInfo = psMetadataLookupPtr (&status, recipe, "PSPHOT.STACK.MATCH.FILTERS");
+    if (!status || !filterInfo) {   
+        psLogMsg ("psphot", PS_LOG_WARN, "PSPHOT.MATCH.STACK.FILTERS not found in the recipe. Will process inputs in order supplied.");
+        filterInfo = NULL;
+    }
+
+    int num = psphotFileruleCount(config, filerule);
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
+    // Find the filter for each of the inputs in fpa concepts
+    psArray *inputFilters = psArrayAlloc(num);
+    for (int i = 0 ; i < num; i++) {
+        if (i != chisqNum) {
+            pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i);
+            psAssert (file, "missing file?");
+
+            psString filterid = psMetadataLookupStr(&status, file->fpa->concepts, "FPA.FILTERID");
+            psAssert (filterid, "missing FPA.FILTERID?");
+            psArraySet(inputFilters, i, filterid);
+        } else {
+            // The chisq image inherits its filterid from the first input so we shouldn't use that.
+            psString tmp  = psStringCopy("chisq");
+            psArraySet(inputFilters, i, tmp);
+            psFree(tmp);
+        }
+    }
+    // Allocate the match info array
+    psphotStackMatchInfo *matchInfo = psAlloc(num *sizeof(psphotStackMatchInfo));
+    memset(matchInfo, 0, num*sizeof(psphotStackMatchInfo));
+    int highest_order = -1;
+    if (filterInfo) {
+        // PSPHOT.STACK.MATCH.FILTERS is a metadata which contains a list of metadata objects each
+        // entry pertaining to a filter.
+        // The objects contained in each filter's metadata are strings.
+      
+        // Loop over the entries in the recipe and find the entry for each our our input readouts.
+        // XXX: Will this work if more than one input with a given filter is supplied?
+        // I think so.
+        psMetadataIterator *iter = psMetadataIteratorAlloc(filterInfo, PS_LIST_HEAD, NULL);
+        psMetadataItem *item = NULL;
+        while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
+            if (item->type != PS_DATA_METADATA) {
+                psAbort ("Invalid type for PSPHOT.STACK.MATCH.FILTER: %s, not a metadata folder", item->name);
+            }
+            psString thisFilter = psMetadataLookupStr (&status, item->data.md, "FILTER.ID");
+            psAssert(thisFilter, "missing FILTER.ID");
+
+            psString orderStr = psMetadataLookupStr (&status, item->data.md, "ORDER");
+            psAssert(orderStr, "missing ORDER");
+            psS32 order = atoi(orderStr);
+
+            psString matchAllStr = psMetadataLookupStr (&status, item->data.md, "MATCH.ALL");
+            psAssert(matchAllStr, "missing MATCH.ALL");
+            bool matchAll = (strcasecmp("true", matchAllStr) == 0);
+
+            // look for this filter in the inputs
+            for (int i = 0; i < num; i++) {
+                if (!strcmp((char *)inputFilters->data[i], thisFilter)) {
+                    // We have an input for this one
+                    matchInfo[i].inputNum = i;
+                    strncpy(matchInfo[i].filterID, thisFilter, MATCH_INFO_FILTER_STR_LEN - 1 );
+                    matchInfo[i].order = order;
+                    matchInfo[i].matchAll = matchAll;
+                    psLogMsg ("psphot", PS_LOG_DETAIL, "input: %d %s match order: %d match all: %d\n",
+                                                               i, thisFilter, order, matchAll);
+                    if (order > highest_order) {
+                        highest_order = order;
+                    }
+                    break;
+                }
+            }
+        }
+        psFree(iter);
+    }
+
+    // Make sure we have a matchInfo for all of the inputs. Fill in an entry for any that were not found.
+    for (int i = 0; i < num; i++) {
+        if (!*matchInfo[i].filterID) {
+            matchInfo[i].inputNum = i;
+            strncpy(matchInfo[i].filterID, inputFilters->data[i], MATCH_INFO_FILTER_STR_LEN - 1);
+            matchInfo[i].order = ++highest_order;
+            matchInfo[i].matchAll = false;
+            psLogMsg ("psphot", PS_LOG_WARN, "Entry in PSPHOT.MATCH.STACK.FILTERS not found for input: %d filter: %s using order %d",
+                i, (char *) inputFilters->data[i], highest_order);
+        }
+    }
+    psFree(inputFilters);
+
+    if (filterInfo) {
+        // Sort the array by ORDER.
+        qsort(matchInfo, num, sizeof(psphotStackMatchInfo), compareMatchInfoByOrder);
+    } else {
+        // no need to sort we just built the list in the order that the inputs were supplied
+    }
+
+    return matchInfo;
+}
+
+bool psphotDropBadMatchedSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects) {
+
+    bool status = false;
+
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Drop Bad Matched Sources ---");
+
+    psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+
+    // select the appropriate recipe information
+    bool keepBadMatches = psMetadataLookupBool (&status, recipe, "PSPHOT.STACK.KEEP.BAD.MATCHES");
+    if (!status) {
+        keepBadMatches = true;
+    }
+
+    if (keepBadMatches) {
+        psLogMsg ("psphot", PS_LOG_INFO, "keeping bad matches");
+        return true;
+    }
+
+    int numImages = psphotFileruleCount(config, filerule);
+    psVector *dropped = psVectorAlloc(numImages, PS_TYPE_U32);
+    psVectorInit(dropped, 0);
+
+    int nDropped = 0;
+    for (int i = 0; i < objects->n; i++) { 
+        pmPhotObj *obj = objects->data[i]; 
+
+        // traverse the array from the end so that sources don't move until after we've processed them
+        for (int j = obj->sources->n - 1; j >= 0; j--) {
+
+            pmSource *source = obj->sources->data[j]; 
+            // This applies only to matched sources
+            if (!(source->mode2 & PM_SOURCE_MODE2_MATCHED)) continue;
+
+            if (isfinite(source->apFlux)) continue;
+
+            psTrace ("psphot", 7, "Dropping matched source from image %d at (%d, %d) no valid flux",
+                source->imageID, source->peak->x, source->peak->y);
+
+            psAssert(source->imageID >= 0 && source->imageID < numImages, "bad imageID %d", source->imageID);
+
+            dropped->data.U32[source->imageID]++;
+            
+            nDropped++;
+
+            psArrayRemoveIndex(obj->sources, j);
+        }
+    }
+
+    psLogMsg ("psphot", PS_LOG_DETAIL, "Dropped %d matched sources with no valid flux", nDropped);
+    psLogMsg ("psphot", PS_LOG_DETAIL, "       Input  Num Dropped");
+    for (int i=0; i<numImages; i++) {
+        psLogMsg ("psphot", PS_LOG_DETAIL, "    %8d     %8d", i, dropped->data.U32[i]);
+    }
+    psFree(dropped);
+
+    return true;
+}
+
+