IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36343


Ignore:
Timestamp:
Dec 4, 2013, 2:00:56 PM (13 years ago)
Author:
eugene
Message:

rules for performing the extended source fits and petrosian analysis are defined in psphotChooseAnalysisOptions.c

Location:
branches/eam_branches/ipp-20130904/psphot/src
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20130904/psphot/src/Makefile.am

    r36271 r36343  
    205205        psphotSourceFits.c             \
    206206        psphotRadiusChecks.c           \
     207        psphotChooseAnalysisOptions.c  \
    207208        psphotOutput.c                 \
    208209        psphotFakeSources.c            \

  • branches/eam_branches/ipp-20130904/psphot/src/psphot.h

    r36301 r36343  
    125125bool            psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool add);
    126126bool            psphotAddOrSubNoise_Threaded (psThreadJob *job);
     127
     128bool            psphotChooseAnalysisOptions (pmConfig *config, const pmFPAview *view, const char *filerule);
     129bool            psphotChooseAnalysisOptionsReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     130bool            psphotChooseAnalysisOptionsByObject(pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
    127131
    128132bool            psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule);

  • branches/eam_branches/ipp-20130904/psphot/src/psphotChooseAnalysisOptions.c

    r36336 r36343  
    33// choose which sources will be processed for the petrosian radii and/or galaxy/trail model fits
    44// a source for which we want galaxy models or petrosians gets
     5
     6// currently we also skip the following items (regardless of recipe values):
     7// TYPE_DEFECT, TYPE_SATURATED, MODE_DEFECT, MODE_SATSTAR, SATSTAR_PROFILE
    58
    69// RULES:
    … …  
    1114// * select by Star/Galaxy?
    1215// ** PSPHOT.EXT.NSIGMA.LIMIT (sets EXT_LIMIT bit -- star or galaxy)
    13 // ** PSPHOT.FIT.EXT_LIMIT (bool : ignore or not EXT_LIMIT)
     16// -- PSPHOT.FIT.EXT_LIMIT (bool : ignore or not EXT_LIMIT)
    1417
    1518// * select by S/N?
    … …  
    1720
    1821// * select by Flux?
     22// * needs to depend on the filter
     23// * define a metadata block
     24// -- EXT.ANALYSIS.MAG.LIMITS (metadata block per filter)
     25
    1926// * select by galactic latitude (global override)?
     27// * if (|b| > LIM) do not limit the number of sources by density
     28// -- GLAT_MAX
     29
    2030// * de-select by density (global override)?
    21 
    22 // Fit all sources?
    23 // * PSPHOT.EXT.FIT.ALL.SOURCES
    24 // * PSPHOT.EXT.FIT.ALL.THRESH (density limit)
    25 
    26 // galactic lat limit?
    27 //
    28 
    29 // density limit?
    30 // magnitude limit?
    31 // S/N limit?
    32 
    33 // DIFF MODE :
    34 // *
     31// ??
    3532
    3633// for which sources do we want to run the extended source analysis (petrosian and/or galaxy fits)?
    … …  
    5552    for (int i = 0; i < num; i++) {
    5653        if (i == chisqNum) continue; // skip chisq image
    57         if (!psphotSourceSizeReadout (config, view, filerule, i, recipe)) {
     54        if (!psphotChooseAnalysisOptionsReadout (config, view, filerule, i, recipe)) {
    5855            psError (PSPHOT_ERR_CONFIG, false, "failed on source size analysis for %s entry %d", filerule, i);
    5956            return false;
    … …  
    6259    return true;
    6360}
     61
     62bool GetGalacticCoords (psSphere *ptGal, psSphere *ptSky, psSphereRot *toGal, pmChip *chip, float xPos, float yPos);
    6463
    6564// this function use an internal flag to mark sources which have already been measured
    … …  
    7776    psAssert (readout, "missing readout?");
    7877
     78    pmChip *chip = pmFPAviewThisChip(view, file->fpa);
     79    psAssert (chip, "missing chip?");
     80
    7981    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    8082    psAssert (detections, "missing detections?");
    … …  
    8890    }
    8991
     92    // we do the petrosian analysis for the same sources as the extended source fits IFF
     93    // this recipe value is turned on (otherwise we only make a selection based on the mag limits)
     94    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
     95
    9096    // Option to enable fitting of all objects with extended model.
    91     options.extFitAll = psMetadataLookupBool(&status, recipe, "PSPHOT.EXT.FIT.ALL.SOURCES");
     97    bool extFitAll = psMetadataLookupBool(&status, recipe, "PSPHOT.EXT.FIT.ALL.SOURCES");
    9298    assert (status);
    9399
    94100    // Fitting everything is fine, but if the source density is high, we probably shouldn't.
    95     options.extFitAllThresh = psMetadataLookupF32(&status, recipe, "PSPHOT.EXT.FIT.ALL.THRESH");
     101    float extFitAllThresh = psMetadataLookupF32(&status, recipe, "PSPHOT.EXT.FIT.ALL.THRESH");
    96102    assert (status);
    97103   
    98104    // Determine if this readout is above the threshold to ext fit all sources
    99     options->extFitAllReadout = false;
    100     if (options->extFitAll) {
     105    if (extFitAll) {
    101106      float maskFrac = psMetadataLookupF32(&status,readout->analysis,"READOUT.MASK.FRAC");
    102       if (status) {
    103         maskFrac = 0.0;
     107      if (status) maskFrac = 0.0;
     108      if (sources->n * (1.0 - maskFrac) > extFitAllThresh) {
     109        extFitAll = false;
    104110      }
    105       if (sources->n * (1.0 - maskFrac) > options->extFitAllThresh) {
    106         options->extFitAllReadout = false;
    107       }
    108       else {
    109         options->extFitAllReadout = true;
    110       }
    111     }
    112    
     111    }
     112
     113    // use EXT_LIMIT bit to select objects?
     114    bool useEXT_LIMIT = psMetadataLookupBool(&status, recipe, "EXT.NSIGMA.LIMIT.USE");
     115    assert (status);
     116
     117    float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");
     118    assert (status);
     119
     120    // use GAL_LIMIT to select / skip objects?
     121    psSphereRot *toGal = NULL;
     122    float GAL_LIMIT = 0.0;
     123    bool useGAL_LIMIT = psMetadataLookupBool(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.USE");
     124    assert (status);
     125    if (useGAL_LIMIT) {
     126        GAL_LIMIT = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT");
     127        toGal = psSphereRotICRSToGalactic();
     128    }
     129
     130    // use MAG.LIMITS
     131    psMetadata *magLimits = psMetadataLookupPtr (&status, recipe, "EXT.ANALYSIS.MAG.LIMITS");
     132    if (!status || !magLimits) {   
     133        psLogMsg ("psphot", PS_LOG_WARN, "EXT.ANALYSIS.MAG.LIMITS not found in the recipe, will use other criteria.");
     134        magLimits = NULL;
     135    }
     136
     137    float extFitFluxLim = NAN;
     138
     139    if (magLimits) {
     140        float extFitMagLimDefault = NAN;
     141        float extFitMagLim = NAN;
     142
     143        // match to the given filter
     144        psString filterID = psMetadataLookupStr(&status, file->fpa->concepts, "FPA.FILTERID");
     145        psAssert (filterID, "missing FPA.FILTERID?");
     146
     147        psMetadataIterator *iter = psMetadataIteratorAlloc(magLimits, PS_LIST_HEAD, NULL);
     148        psMetadataItem *item = NULL;
     149        while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
     150            if (item->type != PS_DATA_METADATA) {
     151                psAbort ("Invalid type for EXT.ANALYSIS.MAG.LIMITS: %s, not a metadata folder", item->name);
     152            }
     153            psString thisFilter = psMetadataLookupStr (&status, item->data.md, "FILTER.ID");
     154            psAssert(thisFilter, "missing FILTER.ID");
     155
     156            // find a matching filter or default to 'any'
     157            if (!strcasecmp (thisFilter, "any")) {
     158                psString extFitMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT");
     159                psAssert(extFitMagLimStr, "missing MAG.LIMIT");
     160                extFitMagLimDefault = atof (extFitMagLimStr);
     161            }
     162
     163            // find a matching filter or default to 'any'
     164            if (!strcasecmp (thisFilter, filterID)) {
     165                psString extFitMagLimStr = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT");
     166                psAssert(extFitMagLimStr, "missing MAG.LIMIT");
     167                extFitMagLim = atof (extFitMagLimStr);
     168                break;
     169            }
     170        }
     171        if (!isfinite (extFitMagLim)) extFitMagLim = extFitMagLimDefault;
     172
     173        // now I need to convert the mag limits into instrumental flux limits
     174        // I need to get a zero point and exposure time for this image
     175       
     176        // select the exposure time
     177        float exptime = psMetadataLookupF32(&status, file->fpa->concepts, "CELL.EXPOSURE");
     178        if (!status) {
     179            exptime = psMetadataLookupF32(&status, file->fpa->concepts, "FPA.EXPOSURE");
     180            psAssert (status, "missing CELL.EXPOSURE and FPA.EXPOSURE?");
     181        }
     182
     183        // select the exposure time
     184        float zeropt = psMetadataLookupF32(&status, file->fpa->concepts, "FPA.ZP");
     185        psAssert (status, "missing FPA.ZP?");
     186
     187        extFitFluxLim = exptime * pow (10.0, 0.4*(zeropt - extFitMagLim));
     188    }
     189
     190    pmSourceTmpF clearBits = ~(PM_SOURCE_TMPF_EXT_FIT | PM_SOURCE_TMPF_PETRO);
    113191    for (psS32 i = 0 ; i < sources->n ; i++) {
    114192
    115193        pmSource *source = (pmSource *) sources->data[i];
     194
     195        // clear the 2 relevant bits
     196        source->tmpFlags &= clearBits;
     197
     198        // skip PSF-like and non-astronomical objects
     199        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
     200        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
     201        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
     202        if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
    116203
    117204        // skip saturated stars modeled with a radial profile
    118205        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
    119206
    120         // If the recipe requests we do extended source fits to everything, set
    121         // the EXT_LIMIT flag
    122         if (options->extFitAllReadout) {
    123           psTrace("psphotSourceClassRegion.EXTALT",10,"In extFitAll: %d %d\n",options->extFitAll,options->extFitAllReadout);
    124           source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
    125         }
    126 
    127         // saturated star (too many saturated pixels or peak above saturation limit).  These
    128         // may also be saturated galaxies, or just large saturated regions.
    129         if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    130             source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    131             continue;
    132         }
    133 
    134         // any sources missing a large fraction should just be treated as PSFs
    135         if ((source->pixWeightNotBad < 0.9) || (source->pixWeightNotPoor < 0.9)) {
    136             source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    137             continue;
    138         }
    139 
    140         bool isCR = (source->moments->SN > 7.0) && (Mminor < options->sizeLimitCR);
    141         if (isCR) {
    142             psTrace("psphotSourceClassRegion.CR",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g CR\n",
    143                     source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid);
    144             source->mode |= PM_SOURCE_MODE_DEFECT;
    145             source->tmpFlags |= PM_SOURCE_TMPF_SIZE_CR_CANDIDATE;
    146             source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    147             Ncr ++;
    148             continue;
    149         }
    150 
    151         source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
    152 
    153         // Alternate extended source limit calculation
    154         if (options->altDiffExt) {
    155           // ratio of major to minor axes
    156           // MRV = Major / Minor
    157           // MRV = (0.5 * (Mxx + Myy) + 0.5 * sqrt( (Mxx + Myy)^2 + 4 Mxy^2)) /
    158           //       (0.5 * (Mxx + Myy) - 0.5 * sqrt( (Mxx + Myy)^2 + 4 Mxy^2))
    159           // MRV = (2 * 0.5 * (Mxx + Myy) - Minor) / Minor
    160           float momentRatioVeres = (Mxx + Myy - Mminor) / Mminor;
    161           bool  isAltEXT = (momentRatioVeres > options->altDiffExtThresh);
    162           if (isAltEXT) {
    163             psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g ALTEXT\t%g %g\n",
    164                     source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid,
    165                     momentRatioVeres,options->altDiffExtThresh);
    166             source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
    167             source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    168             Next ++;
    169             continue;
    170           }
    171         }
    172        
    173         // Everything else should just be treated as a PSF
     207        // Do the fits if the recipe requests we do extended source fits to everything
     208        if (extFitAll) {
     209          source->tmpFlags |= PM_SOURCE_TMPF_EXT_FIT;
     210          if (doPetrosian) source->tmpFlags |= PM_SOURCE_TMPF_PETRO;
     211          continue;
     212        }
     213
     214        if (useEXT_LIMIT) {
     215            if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) {
     216                // not extended so skip
     217                continue;
     218            }
     219        }
     220
     221        if (useGAL_LIMIT) {
     222            psSphere ptGal, ptSky;
     223            GetGalacticCoords (&ptGal, &ptSky, toGal, chip, source->peak->xf, source->peak->yf);
     224            if (fabs(ptGal.d) < GAL_LIMIT) continue;
     225            // include an exception for low density skycells below the limit?
     226        }
     227
     228        // for petro and extFit, we will either use the mag limits or the S/N
     229        if (isfinite(extFitFluxLim)) {
     230            if (source->moments->KronFlux > extFitFluxLim) {
     231                source->tmpFlags |= PM_SOURCE_TMPF_EXT_FIT;
     232                if (doPetrosian) source->tmpFlags |= PM_SOURCE_TMPF_PETRO;
     233            }
     234        } else {
     235            if (source->moments->KronFlux > SN_LIM * source->moments->KronFluxErr) {
     236                source->tmpFlags |= PM_SOURCE_TMPF_EXT_FIT;
     237                if (doPetrosian) source->tmpFlags |= PM_SOURCE_TMPF_PETRO;
     238            }
     239        }
    174240    }
    175241
    … …  
    178244    return true;
    179245}
     246
     247// this function use an internal flag to mark sources which have already been measured
     248bool psphotChooseAnalysisOptionsByObject(pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects)
     249{
     250    bool status;
     251
     252    psTimerStart ("psphot.options");
     253
     254    fprintf (stdout, "\n");
     255    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Choose Analysis Options (By Object) ---");
     256
     257    // select the appropriate recipe information
     258    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     259    psAssert (recipe, "missing recipe?");
     260
     261    // S/N lim to perform radial aperture analysis
     262    float SN_LIM_RADIAL = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
     263
     264    // we do the petrosian analysis for the same sources as the extended source fits IFF
     265    // this recipe value is turned on (otherwise we only make a selection based on the mag limits)
     266    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
     267
     268    // Option to enable fitting of all objects with extended model.
     269    bool extFitAll = psMetadataLookupBool(&status, recipe, "PSPHOT.EXT.FIT.ALL.SOURCES");
     270    assert (status);
     271
     272    // Fitting everything is fine, but if the source density is high, we probably shouldn't.
     273    float extFitAllThresh = psMetadataLookupF32(&status, recipe, "PSPHOT.EXT.FIT.ALL.THRESH");
     274    assert (status);
     275   
     276    // use EXT_LIMIT bit to select objects?
     277    bool useEXT_LIMIT = psMetadataLookupBool(&status, recipe, "EXT.NSIGMA.LIMIT.USE");
     278    assert (status);
     279
     280    float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");
     281    assert (status);
     282
     283    // use GAL_LIMIT to select / skip objects?
     284    psSphereRot *toGal = NULL;
     285    float GAL_LIMIT = 0.0;
     286    bool useGAL_LIMIT = psMetadataLookupBool(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT.USE");
     287    assert (status);
     288    if (useGAL_LIMIT) {
     289        GAL_LIMIT = psMetadataLookupF32(&status, recipe, "EXT.FIT.MIN.GAL.LIMIT");
     290        toGal = psSphereRotICRSToGalactic();
     291    }
     292
     293    // use MAG.LIMITS
     294    psMetadata *magLimits = psMetadataLookupPtr (&status, recipe, "EXT.ANALYSIS.MAG.LIMITS");
     295    if (!status || !magLimits) {   
     296        psLogMsg ("psphot", PS_LOG_WARN, "EXT.ANALYSIS.MAG.LIMITS not found in the recipe, will use other criteria.");
     297        magLimits = NULL;
     298    }
     299
     300    int num = psphotFileruleCount(config, filerule);
     301
     302    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
     303    if (!status) chisqNum = -1;
     304
     305    // Find the filter for each of the inputs in fpa concepts
     306    psArray *inputFilters = psArrayAlloc(num);
     307    psArray *chips = psArrayAlloc(num);
     308    psVector *zeropt = psVectorAlloc (inputFilters->n, PS_TYPE_F32);
     309    psVector *exptime = psVectorAlloc (inputFilters->n, PS_TYPE_F32);
     310
     311    // we will not use the chisq image to set the fitting limits
     312    if (chisqNum >= 0) {
     313      inputFilters->data[chisqNum] = psStringCopy("chisq");
     314      chips->data[chisqNum] = NULL;
     315      zeropt->data.F32[chisqNum] = NAN;
     316      exptime->data.F32[chisqNum] = NAN;
     317    }
     318
     319    // get the needed metadata for the non-chisq images
     320    for (int i = 0 ; i < num; i++) {
     321      if (i == chisqNum) continue;
     322
     323      pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i);
     324      psAssert (file, "missing file?");
     325
     326      pmChip *chip = pmFPAviewThisChip(view, file->fpa);
     327      psAssert (chip, "missing chip?");
     328      chips->data[i] = psMemIncrRefCounter (chip);
     329
     330      pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     331      psAssert (readout, "missing readout?");
     332
     333      pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     334      psAssert (detections, "missing detections?");
     335
     336      psArray *sources = detections->allSources;
     337      psAssert (sources, "missing sources?");
     338
     339      float maskFrac = psMetadataLookupF32(&status,readout->analysis,"READOUT.MASK.FRAC");
     340      if (!status) maskFrac = 0.0;
     341      if (sources->n * (1.0 - maskFrac) > extFitAllThresh) {
     342        extFitAll = false;
     343      }
     344
     345      // select the filterID for this image
     346      char *filterID = psMetadataLookupStr(&status, file->fpa->concepts, "FPA.FILTERID");
     347      psAssert (status, "missing FPA.FILTERID?");
     348      inputFilters->data[i] = psStringCopy (filterID);
     349
     350      // select the exposure time for this image
     351      exptime->data.F32[i] = psMetadataLookupF32(&status, file->fpa->concepts, "CELL.EXPOSURE");
     352      if (!status) {
     353        exptime->data.F32[i] = psMetadataLookupF32(&status, file->fpa->concepts, "FPA.EXPOSURE");
     354        psAssert (status, "missing CELL.EXPOSURE and FPA.EXPOSURE?");
     355      }
     356
     357      // select the zero point for this image
     358      zeropt->data.F32[i] = psMetadataLookupF32(&status, file->fpa->concepts, "FPA.ZP");
     359      psAssert (status, "missing FPA.ZP?");
     360    }
     361
     362    // find extFitFluxLim->data.F32[i] for i == image number
     363    psVector *extFitFluxLim = NULL;
     364    if (magLimits) {
     365        extFitFluxLim = psVectorAlloc (inputFilters->n, PS_TYPE_F32);
     366        psVectorInit (extFitFluxLim, NAN);
     367
     368        float extFitMagLimDefault = NAN;
     369
     370        // match mag limits (flux limits) to the filters
     371        psMetadataIterator *iter = psMetadataIteratorAlloc(magLimits, PS_LIST_HEAD, NULL);
     372        psMetadataItem *item = NULL;
     373        while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
     374            if (item->type != PS_DATA_METADATA) {
     375                psAbort ("Invalid type for EXT.ANALYSIS.MAG.LIMITS: %s, not a metadata folder", item->name);
     376            }
     377            psString thisFilter = psMetadataLookupStr (&status, item->data.md, "FILTER.ID");
     378            psAssert(thisFilter, "missing FILTER.ID");
     379
     380            // save the default value to assign to unset filters
     381            if (!strcasecmp (thisFilter, "any")) {
     382                psString magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT");
     383                psAssert(magString, "missing MAG.LIMIT");
     384                extFitMagLimDefault = atof (magString);
     385                continue;
     386            }
     387
     388            // not every entry in the metadata block needs to match to an image in our list
     389            for (int i = 0; i < num; i++) {
     390                if (i == chisqNum) continue;
     391                if (!strcasecmp (thisFilter, inputFilters->data[i])) {
     392                    psString magString = psMetadataLookupStr (&status, item->data.md, "MAG.LIMIT");
     393                    psAssert(magString, "missing MAG.LIMIT");
     394                    float magvalue = atof (magString);
     395                    extFitFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - magvalue));
     396                    break;
     397                }
     398            }
     399        }
     400
     401        for (int i = 0; i < num; i++) {
     402            if (i == chisqNum) continue;
     403            if (isfinite(extFitFluxLim->data.F32[i])) continue;
     404            extFitFluxLim->data.F32[i] = exptime->data.F32[i] * pow (10.0, 0.4*(zeropt->data.F32[i] - extFitMagLimDefault));
     405        }
     406    }
     407
     408    pmSourceTmpF clearBits = ~(PM_SOURCE_TMPF_EXT_FIT | PM_SOURCE_TMPF_PETRO);
     409
     410    for (int i = 0; i < objects->n; i++) {
     411        pmPhotObj *object = objects->data[i];
     412        if (!object) continue;
     413        if (!object->sources) continue;
     414
     415        // we check each source for an object and keep the object if any source is valid
     416
     417        bool doObjectRadial = false;
     418        bool doObjectExtFit = false;
     419        for (int j = 0; !doObjectExtFit && !doObjectRadial && (j < object->sources->n); j++) {
     420
     421            pmSource *source = object->sources->data[j];
     422            if (!source) continue;
     423            if (!source->peak) continue;
     424            if (source->imageID < 0) continue; // skip sources which come from other images?
     425            if (source->imageID >= num) continue; // skip sources which come from other images?
     426
     427            // skip PSF-like and non-astronomical objects
     428            if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
     429            if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
     430            if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
     431            if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
     432       
     433            // skip saturated stars modeled with a radial profile
     434            if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
     435       
     436            // XXX should I fit all even if one of the detections matches the above?
     437
     438            // check on radial aperture analysis (fewer options for now)
     439            if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
     440                bool doSourceRadial = (source->moments->KronFlux > SN_LIM_RADIAL * source->moments->KronFluxErr);
     441                doObjectRadial = doObjectRadial | doSourceRadial;
     442            } else {
     443                bool doSourceRadial = (sqrt(source->peak->detValue) > SN_LIM_RADIAL);
     444                doObjectRadial = doObjectRadial | doSourceRadial;
     445            }
     446
     447            // Do the fits if the recipe requests we do extended source fits to everything
     448            if (extFitAll) {
     449                doObjectExtFit = true;
     450                continue;
     451            }
     452
     453            if (useEXT_LIMIT) {
     454                if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) {
     455                    continue; // not extended so skip
     456                }
     457            }
     458
     459            int imageID = source->imageID;
     460
     461            if (useGAL_LIMIT) {
     462                psSphere ptGal, ptSky;
     463                GetGalacticCoords (&ptGal, &ptSky, toGal, chips->data[imageID], source->peak->xf, source->peak->yf);
     464                if (fabs(ptGal.d) < GAL_LIMIT) continue;
     465                // include an exception for low density skycells below the limit?
     466            }
     467
     468            float fluxLim = extFitFluxLim ? extFitFluxLim->data.F32[imageID] : NAN;
     469            // for petro and extFit, we will either use the mag limits or the S/N
     470            if (isfinite(fluxLim)) {
     471                if (source->moments->KronFlux > extFitFluxLim->data.F32[imageID]) {
     472                    doObjectExtFit = true;
     473                }
     474            } else {
     475                if (source->moments->KronFlux > SN_LIM * source->moments->KronFluxErr) {
     476                    doObjectExtFit = true;
     477                }
     478            }
     479        }
     480
     481        for (int j = 0; j < object->sources->n; j++) {
     482               
     483            pmSource *source = object->sources->data[j];
     484            if (!source) continue;
     485            if (!source->peak) continue;
     486               
     487            // clear the 2 relevant bits
     488            source->tmpFlags &= clearBits;
     489               
     490            // skip PSF-like and non-astronomical objects??
     491            if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
     492            if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
     493            if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
     494            if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
     495               
     496            // skip saturated stars modeled with a radial profile
     497            if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
     498               
     499            // Do the fits if the recipe requests we do extended source fits to everything
     500            if (doObjectExtFit) {
     501                source->tmpFlags |= PM_SOURCE_TMPF_EXT_FIT;
     502                if (doPetrosian) source->tmpFlags |= PM_SOURCE_TMPF_PETRO;
     503            }
     504
     505            // Do the fits if the recipe requests we do extended source fits to everything
     506            if (doObjectRadial) {
     507                source->tmpFlags |=  PM_SOURCE_TMPF_RADIAL_KEEP;
     508                source->tmpFlags &= ~PM_SOURCE_TMPF_RADIAL_SKIP;
     509            } else {
     510                source->tmpFlags |=  PM_SOURCE_TMPF_RADIAL_SKIP;
     511                source->tmpFlags &= ~PM_SOURCE_TMPF_RADIAL_KEEP;
     512            }
     513        }
     514    }
     515
     516    psLogMsg ("psphot.options", PS_LOG_WARN, "choose analysis options for %ld objects: %f sec\n", objects->n, psTimerMark ("psphot.options"));
     517
     518    return true;
     519}
     520
     521bool GetGalacticCoords (psSphere *ptGal, psSphere *ptSky, psSphereRot *toGal, pmChip *chip, float xPos, float yPos) {
     522
     523    pmFPA *fpa = chip->parent;
     524
     525    if (!chip->toFPA) goto escape;
     526    if (!fpa->toTPA) goto escape;
     527    if (!fpa->toSky) goto escape;
     528
     529    psPlane ptCH, ptFP, ptTP;
     530
     531    // calculate the astrometry for the coordinate of interest
     532    ptCH.x = xPos;
     533    ptCH.y = yPos;
     534    psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
     535    psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
     536    psDeproject (ptSky, &ptTP, fpa->toSky);
     537    psSphereRotApply (ptGal, toGal, ptSky);
     538
     539    return true;
     540
     541escape:
     542    // no astrometry calibration, give up
     543    ptSky->r = NAN;
     544    ptSky->d = NAN;
     545
     546    ptGal->r = NAN;
     547    ptGal->d = NAN;
     548
     549    return false;
     550}
Note: See TracChangeset for help on using the changeset viewer.