IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 34338


Ignore:
Timestamp:
Aug 23, 2012, 9:40:34 AM (14 years ago)
Author:
bills
Message:

Changes to the creation of "matched" sources for forced photometry in psphotStack to
reduce the number of false sources.
New recipe values PSPHOT.STACK.MIN.DETECT.FOR.FORCED and PSPHOT.STACK.MATCH.FILTERS are added..
The latter is an array (actually metadata) of metadatas an entry for each filter.
Contained in each entry are two values ORDER and MATCH.ALL.
When matching objects the inputs are processed in the defined order. (The intent is that
more sensitive bands are to be processed first.)
If the number of detections for each object is less than PSPHOT.STACK.MIN.DETECT.FOR.FORCED
(set to 2 here) matched sources are not created unless one of the detections was in a
bad for which MATCH.ALL is set. This is currently only true for y band.
This reduces the number of suprious sources.
Note that the entries defined in PSPHOT.STACK.MIN.DETECT.FOR.FORCED are the PS1 bands.
If psphotStack is supplied inputs whose FILTER.ID is not in the defined set, then will be processed
fine, being processed in order supplied following any inputs which do have entries.

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/ippconfig/recipes/psphot.config

    r34307 r34338  
    366366PSPHOT.STACK.SPLIT.LINEAR.FIT       BOOL  T    # require positive flux for pass 3 linear fit then fit matched sources separately
    367367
     368# psphotStack source matching parameters
     369PSHOT.STACK.MIN.DETECT.FOR.FORCED   S32   2   # number of bands that an object must be detected in before creating
     370                                              # sources in undetected bands
     371                                             
     372# PSPHOT.STACK.MATCH.FILTER
     373# list of metadata to control how source matching is performed for each band
     374# FILTER.ID is used to match to the FPA's concept value FPA.FILTERID
     375# ORDER controls the order in which the bands are processed. Lower numbers are processed first.
     376#       The intent is to process in order of sensitivity
     377# MATCH.ALL if true causes forced photometry sources for objects detected in that band but don't meet
     378#       the minimum number of detections cut
     379# TYPE is used in the recipe bookkeeping and has no semantic (but must be unique)
     380# Note: if an input has a filter not in this list it will be processed after all of the others with MATCH.ALL = FALSE
     381PSPHOT.STACK.MATCH.FILTERS   METADATA
     382    TYPE        PSPHOT_STACK.MATCH.FILTER FILTER.ID ORDER MATCH.ALL
     383    gband       PSPHOT_STACK.MATCH.FILTER   g        4       FALSE
     384    rband       PSPHOT_STACK.MATCH.FILTER   r        2       FALSE
     385    iband       PSPHOT_STACK.MATCH.FILTER   i        1       FALSE
     386    zband       PSPHOT_STACK.MATCH.FILTER   z        3       FALSE
     387    yband       PSPHOT_STACK.MATCH.FILTER   y        5       TRUE
     388    wband       PSPHOT_STACK.MATCH.FILTER   w        0       FALSE
     389END
     390
    368391RADIAL_APERTURES                    BOOL  F    # calculate flux in circular radial apertures?
    369392RADIAL_APERTURES_SN_LIM             F32   0.0  # S/N limit for radial aperture calculation

  • trunk/psphot/src/psphotSourceMatch.c

    r33587 r34338  
    11# include "psphotInternal.h"
    22
    3 bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
     3// Structure for storing the contents of the PSPHOT.STACK.MATCH.FILTERS list from recipe
     4#define MATCH_INFO_FILTER_STR_LEN 16
     5typedef struct {
     6    int     inputNum;
     7    int     order;
     8    bool    matchAll;
     9    char    filterID[MATCH_INFO_FILTER_STR_LEN];
     10} psphotStackMatchInfo;
     11
     12static psphotStackMatchInfo * psphotStackGetMatchInfo (pmConfig *config, const pmFPAview *view, const char *filerule);
     13
     14// functions for sorting the match info array
     15static int compareMatchInfoByInputNum(const void *a, const void *b);
     16static int compareMatchInfoByOrder(const void *a, const void *b);
     17
     18bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects, psphotStackMatchInfo *matchInfo);
    419bool psphotMatchSourcesSetIDs (psArray *objects);
    520
    … …  
    1429    int num = psphotFileruleCount(config, filerule);
    1530
    16     // loop over the available readouts
    17     for (int i = 0; i < num; i++) {
     31    psphotStackMatchInfo *matchInfo = psphotStackGetMatchInfo(config, view, filerule);
     32   
     33    // loop over the available readouts matching sources found to objects. The inputs are processed
     34    // in the order specified by the recipe
     35    for (int j = 0; j < num; j++) {
     36        int i = matchInfo[j].inputNum;
    1837        if (!psphotMatchSourcesReadout (objects, config, view, filerule, i)) {
    1938            psError (PSPHOT_ERR_CONFIG, false, "failed to merge sources for %s entry %d", filerule, i);
    … …  
    2342    }
    2443
     44    // Now re-order the matchInfo array by input number so we can find each inputs entry easily
     45    qsort(matchInfo, num, sizeof(psphotStackMatchInfo), compareMatchInfoByInputNum);
     46
    2547    // create sources for images where an object has been detected in the other images
    26     psphotMatchSourcesAddMissing (config, view, filerule, objects);
     48    psphotMatchSourcesAddMissing (config, view, filerule, objects, matchInfo);
    2749
    2850    // choose a consistent position; set common sequence values
    2951    psphotMatchSourcesSetIDs (objects);
    3052
     53    psFree(matchInfo);
     54
    3155    return objects;
    3256}
    3357
    34 bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, const char *filerule, int index) { 
     58bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, const char *filerule, int index) {
    3559 
    3660    bool status = false;
    … …  
    6185# define NEXT1 { i++; continue; }
    6286# define NEXT2 { j++; continue; }
    63 bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS) { 
     87bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS) {
    6488 
    6589    float dx, dy;
    … …  
    7094    sources = psArraySort (sources, pmSourceSortByX);
    7195    objects = psArraySort (objects, pmPhotObjSortByX);
    72  
     96
    7397    psVector *foundSrc = psVectorAlloc(sources->n, PS_TYPE_U8);
    7498    psVectorInit (foundSrc, 0);
    … …  
    87111        pmPhotObj *obj = objects->data[j];
    88112 
    89         if (!src) NEXT1; 
     113        if (!src) NEXT1;
    90114        if (!src->peak) NEXT1;
    91115        if (!isfinite(src->peak->xf)) NEXT1;
    92116        if (!isfinite(src->peak->yf)) NEXT1;
    93  
     117
    94118        if (!obj) NEXT2;
    95119        if (!isfinite(obj->x)) NEXT2;
    … …  
    141165    }
    142166
    143     // add missed sources to new objects
    144 
     167    // create new objects for unmatched sources
    145168    for (i = 0; i < sources->n; i++) {
    146169
    147         if (foundSrc->data.U8[i]) continue;
     170        if (foundSrc->data.U8[i]) continue;
    148171
    149172        pmSource *src = sources->data[i];
    150173
    151         pmPhotObj *obj = pmPhotObjAlloc();
    152         pmPhotObjAddSource(obj, src);
    153         psArrayAdd (objects, 100, obj);
    154         psFree (obj);
     174        pmPhotObj *obj = pmPhotObjAlloc();
     175        pmPhotObjAddSource(obj, src);
     176        psArrayAdd (objects, 100, obj);
     177        psFree (obj);
    155178    }
    156179    psLogMsg ("psphot", PS_LOG_DETAIL, "matched sources (%ld vs %ld)", sources->n, objects->n);
    … …  
    161184}
    162185
    163 bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects) {
     186bool psphotMatchSourcesAddMissing (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects, psphotStackMatchInfo *matchInfo) {
    164187
    165188    bool status = false;
    … …  
    169192    psAssert (recipe, "missing recipe?");
    170193
     194    int minDetectionsForForced = psMetadataLookupS32 (&status, recipe, "PSPHOT.STACK.MIN.DETECT.FOR.FORCED");
     195    if (!status) {
     196        minDetectionsForForced = 2;
     197    }
     198
    171199    // determine properties (sky, moments) of initial sources
    172200    float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
    … …  
    205233        int nSpansMax = 0;
    206234        int iSpansMax = -1;
     235
     236        bool matchAll = false;
    207237
    208238        // mark the images for which sources have been found
    … …  
    219249                iSpansMax = j;
    220250            }
     251            // If this detection was on an input that has the matchAll flag set we create matched sources
     252            // irespective of the number of inputs in which the object was found.
     253            if (matchInfo[index].matchAll) {
     254                matchAll = true;
     255            }
    221256
    222257            found->data.U8[index] = 1;
    223258        }
     259
     260        // skip adding matched sources for this object if the number of detections for less than
     261        // the supplied mininum unless one of the detections was in a band for which the matchAll flag is set
     262        if (!matchAll && obj->sources->n < minDetectionsForForced) {
     263            continue;
     264        }
    224265
    225266        // we make a copy of the largest footprint; this will be used for all new sources associated with this object
    … …  
    299340        pmPhotObj *obj = objects->data[i];
    300341        nSources += obj->sources->n;
    301         psAssert (obj->sources->n == nImages, "failed to match sources?");
     342        if (minDetectionsForForced <= 1) {
     343            psAssert (obj->sources->n == nImages, "failed to match sources?");
     344        }
    302345    }
    303346    psLogMsg ("psphot", PS_LOG_DETAIL, "total of %d sources for %d images", nSources, nImages);
    … …  
    375418    return copy;
    376419}
     420
     421// qsort function to sort match info by order
     422static int compareMatchInfoByOrder(const void *a, const void *b) {
     423    psphotStackMatchInfo *infoA = (psphotStackMatchInfo *) a;
     424    psphotStackMatchInfo *infoB = (psphotStackMatchInfo *) b;
     425
     426    int     result;
     427    if (infoA->order < infoB->order) {
     428        result = -1;
     429    } else if (infoA->order == infoB->order) {
     430        result = 0;
     431    } else {
     432        result = 1;
     433    }
     434    return result;
     435}
     436
     437// qsort function to sort match info by input number
     438static int compareMatchInfoByInputNum(const void *a, const void *b) {
     439    psphotStackMatchInfo *infoA = (psphotStackMatchInfo *) a;
     440    psphotStackMatchInfo *infoB = (psphotStackMatchInfo *) b;
     441
     442    int     result;
     443    if (infoA->inputNum < infoB->inputNum) {
     444        result = -1;
     445    } else if (infoA->inputNum == infoB->inputNum) {
     446        result = 0;
     447    } else {
     448        result = 1;
     449    }
     450    return result;
     451}
     452
     453// psphotStackGetMatchInfo
     454// process the recipe PSPHOT.STACK.MATCH.FILTERS which controls the order that the inputs
     455// are processed for source matching and the rules for creating matched sources for sources that
     456// are only detected in a single band.
     457static psphotStackMatchInfo *psphotStackGetMatchInfo (pmConfig *config, const pmFPAview *view, const char *filerule) {
     458
     459    bool status = false;
     460    psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     461
     462    psMetadata *filterInfo = psMetadataLookupPtr (&status, recipe, "PSPHOT.STACK.MATCH.FILTERS");
     463    if (!status || !filterInfo) {   
     464        psLogMsg ("psphot", PS_LOG_WARN, "PSPHOT.MATCH.STACK.FILTERS not found in the recipe. Will process inputs in order supplied.");
     465        filterInfo = NULL;
     466    }
     467
     468    int num = psphotFileruleCount(config, filerule);
     469    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
     470    if (!status) chisqNum = -1;
     471
     472    // Find the filter for each of the inputs in fpa concepts
     473    psArray *inputFilters = psArrayAlloc(num);
     474    for (int i = 0 ; i < num; i++) {
     475        if (i != chisqNum) {
     476            pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i);
     477            psAssert (file, "missing file?");
     478
     479            psString filterid = psMetadataLookupStr(&status, file->fpa->concepts, "FPA.FILTERID");
     480            psAssert (filterid, "missing FPA.FILTERID?");
     481            psArraySet(inputFilters, i, filterid);
     482        } else {
     483            // The chisq image inherits its filterid from the first input so we shouldn't use that.
     484            psString tmp  = psStringCopy("chisq");
     485            psArraySet(inputFilters, i, tmp);
     486            psFree(tmp);
     487        }
     488    }
     489    // Allocate the match info array
     490    psphotStackMatchInfo *matchInfo = psAlloc(num *sizeof(psphotStackMatchInfo));
     491    memset(matchInfo, 0, num*sizeof(psphotStackMatchInfo));
     492    int highest_order = -1;
     493    if (filterInfo) {
     494        // PSPHOT.STACK.MATCH.FILTERS is a metadata which contains a list of metadata objects each
     495        // entry pertaining to a filter.
     496        // The objects contained in each filter's metadata are strings.
     497     
     498        // Loop over the entries in the recipe and find the entry for each our our input readouts.
     499        // XXX: Will this work if more than one input with a given filter is supplied?
     500        // I think so.
     501        psMetadataIterator *iter = psMetadataIteratorAlloc(filterInfo, PS_LIST_HEAD, NULL);
     502        psMetadataItem *item = NULL;
     503        while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
     504            if (item->type != PS_DATA_METADATA) {
     505                psAbort ("Invalid type for PSPHOT.STACK.MATCH.FILTER: %s, not a metadata folder", item->name);
     506            }
     507            psString thisFilter = psMetadataLookupStr (&status, item->data.md, "FILTER.ID");
     508            psAssert(thisFilter, "missing FILTER.ID");
     509
     510            psString orderStr = psMetadataLookupStr (&status, item->data.md, "ORDER");
     511            psAssert(orderStr, "missing ORDER");
     512            psS32 order = atoi(orderStr);
     513
     514            psString matchAllStr = psMetadataLookupStr (&status, item->data.md, "MATCH.ALL");
     515            psAssert(matchAllStr, "missing MATCH.ALL");
     516            bool matchAll = (strcasecmp("true", matchAllStr) == 0);
     517
     518            // look for this filter in the inputs
     519            for (int i = 0; i < num; i++) {
     520                if (!strcmp((char *)inputFilters->data[i], thisFilter)) {
     521                    // We have an input for this one
     522                    matchInfo[i].inputNum = i;
     523                    strncpy(matchInfo[i].filterID, thisFilter, MATCH_INFO_FILTER_STR_LEN - 1 );
     524                    matchInfo[i].order = order;
     525                    matchInfo[i].matchAll = matchAll;
     526                    psLogMsg ("psphot", PS_LOG_DETAIL, "input: %d %s match order: %d match all: %d\n",
     527                                                               i, thisFilter, order, matchAll);
     528                    if (order > highest_order) {
     529                        highest_order = order;
     530                    }
     531                    break;
     532                }
     533            }
     534        }
     535        psFree(iter);
     536    }
     537
     538    // Make sure we have a matchInfo for all of the inputs. Fill in an entry for any that were not found.
     539    for (int i = 0; i < num; i++) {
     540        if (!*matchInfo[i].filterID) {
     541            matchInfo[i].inputNum = i;
     542            strncpy(matchInfo[i].filterID, inputFilters->data[i], MATCH_INFO_FILTER_STR_LEN - 1);
     543            matchInfo[i].order = ++highest_order;
     544            matchInfo[i].matchAll = false;
     545            psLogMsg ("psphot", PS_LOG_WARN, "Entry in PSPHOT.MATCH.STACK.FILTERS not found for input: %d filter: %s using order %d",
     546                i, (char *) inputFilters->data[i], highest_order);
     547        }
     548    }
     549    psFree(inputFilters);
     550
     551    if (filterInfo) {
     552        // Sort the array by ORDER.
     553        qsort(matchInfo, num, sizeof(psphotStackMatchInfo), compareMatchInfoByOrder);
     554    } else {
     555        // no need to sort we just built the list in the order that the inputs were supplied
     556    }
     557
     558    return matchInfo;
     559}
Note: See TracChangeset for help on using the changeset viewer.