Index: trunk/ippconfig/recipes/psphot.config =================================================================== --- trunk/ippconfig/recipes/psphot.config (revision 34337) +++ trunk/ippconfig/recipes/psphot.config (revision 34338) @@ -366,4 +366,27 @@ PSPHOT.STACK.SPLIT.LINEAR.FIT BOOL T # require positive flux for pass 3 linear fit then fit matched sources separately +# psphotStack source matching parameters +PSHOT.STACK.MIN.DETECT.FOR.FORCED S32 2 # number of bands that an object must be detected in before creating + # sources in undetected bands + +# PSPHOT.STACK.MATCH.FILTER +# list of metadata to control how source matching is performed for each band +# FILTER.ID is used to match to the FPA's concept value FPA.FILTERID +# ORDER controls the order in which the bands are processed. Lower numbers are processed first. +# The intent is to process in order of sensitivity +# MATCH.ALL if true causes forced photometry sources for objects detected in that band but don't meet +# the minimum number of detections cut +# TYPE is used in the recipe bookkeeping and has no semantic (but must be unique) +# Note: if an input has a filter not in this list it will be processed after all of the others with MATCH.ALL = FALSE +PSPHOT.STACK.MATCH.FILTERS METADATA + TYPE PSPHOT_STACK.MATCH.FILTER FILTER.ID ORDER MATCH.ALL + gband PSPHOT_STACK.MATCH.FILTER g 4 FALSE + rband PSPHOT_STACK.MATCH.FILTER r 2 FALSE + iband PSPHOT_STACK.MATCH.FILTER i 1 FALSE + zband PSPHOT_STACK.MATCH.FILTER z 3 FALSE + yband PSPHOT_STACK.MATCH.FILTER y 5 TRUE + wband PSPHOT_STACK.MATCH.FILTER w 0 FALSE +END + RADIAL_APERTURES BOOL F # calculate flux in circular radial apertures? RADIAL_APERTURES_SN_LIM F32 0.0 # S/N limit for radial aperture calculation Index: trunk/psphot/src/psphotSourceMatch.c =================================================================== --- trunk/psphot/src/psphotSourceMatch.c (revision 34337) +++ trunk/psphot/src/psphotSourceMatch.c (revision 34338) @@ -1,5 +1,20 @@ # 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,6 +29,10 @@ 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,14 +42,19 @@ } + // 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,5 +85,5 @@ # 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,5 +94,5 @@ sources = psArraySort (sources, pmSourceSortByX); objects = psArraySort (objects, pmPhotObjSortByX); - + psVector *foundSrc = psVectorAlloc(sources->n, PS_TYPE_U8); psVectorInit (foundSrc, 0); @@ -87,9 +111,9 @@ 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,16 +165,15 @@ } - // 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,5 +184,5 @@ } -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,4 +192,9 @@ 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,4 +233,6 @@ int nSpansMax = 0; int iSpansMax = -1; + + bool matchAll = false; // mark the images for which sources have been found @@ -219,7 +249,18 @@ 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,5 +340,7 @@ 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,2 +418,142 @@ 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; +}