IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36342


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

rules for performing the extended source fits and petrosian analysis are defined in psphotChooseAnalysisOptions.c; only the model-specific limits are applied here

Location:
branches/eam_branches/ipp-20130904/psphot/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r36198 r36342  
    22void psphotRadialProfileShowSkips ();
    33
    4 // measure the elliptical radial profile and use this to measure the petrosian parameters for the sources
     4// measure the petrosian parameters for the sources
    55
    66// for now, let's store the detections on the readout->analysis for each readout
    … …  
    1616    psAssert (recipe, "missing recipe?");
    1717
    18     // perform full non-linear fits / extended source analysis?
    19     if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANALYSIS")) {
     18    // measure petrosians?
     19    if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN")) {
    2020        psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n");
    2121        return true;
    … …  
    4242int  Nskip5 = 0;
    4343int  Nskip6 = 0;
    44 int  Nskip7 = 0;
    45 int  Nskip8 = 0;
    46 int  Nskip9 = 0;
    47 int Nskip10 = 0;
    48 int Nskip11 = 0;
    49 int Nskip12 = 0;
    50 int Nskip13 = 0;
    51 int Nskip14 = 0;
    5244
    5345# define SKIP(VALUE) { VALUE++; continue; }
    … …  
    5749
    5850    bool status;
    59     int NfaintEXT = 0;
    60     int NfaintPSF = 0;
    61 
    6251    int Next = 0;
    6352    int Npetro = 0;
    … …  
    131120            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npetro
    132121            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nannuli
    133             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NfaintEXT
    134             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NfaintPSF
    135122
    136123// set this to 0 to run without threading
    … …  
    154141            scalar = job->args->data[7];
    155142            Nannuli += scalar->data.S32;
    156             scalar = job->args->data[8];
    157             NfaintEXT += scalar->data.S32;
    158             scalar = job->args->data[9];
    159             NfaintPSF += scalar->data.S32;
    160143            psFree(job);
    161144# endif
    … …  
    182165                scalar = job->args->data[7];
    183166                Nannuli += scalar->data.S32;
    184                 scalar = job->args->data[8];
    185                 NfaintEXT += scalar->data.S32;
    186                 scalar = job->args->data[9];
    187                 NfaintPSF += scalar->data.S32;
    188167            }
    189168            psFree(job);
    … …  
    196175    psLogMsg ("psphot", PS_LOG_INFO, "  %d petrosian\n", Npetro);
    197176    psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
    198     psLogMsg ("psphot", PS_LOG_INFO, "  skipped: %d EXT, %d PSF\n", NfaintEXT, NfaintPSF);
    199177
    200178    fprintf (stderr, "ext analysis skipped @ 1  : %d\n", Nskip1);
    … …  
    204182    fprintf (stderr, "ext analysis skipped @ 5  : %d\n", Nskip5);
    205183    fprintf (stderr, "ext analysis skipped @ 6  : %d\n", Nskip6);
    206     fprintf (stderr, "ext analysis skipped @ 7  : %d\n", Nskip7);
    207     fprintf (stderr, "ext analysis skipped @ 8  : %d\n", Nskip8);
    208     fprintf (stderr, "ext analysis skipped @ 9  : %d\n", Nskip9);
    209     fprintf (stderr, "ext analysis skipped @ 10 : %d\n", Nskip10);
    210     fprintf (stderr, "ext analysis skipped @ 11 : %d\n", Nskip11);
    211     fprintf (stderr, "ext analysis skipped @ 12 : %d\n", Nskip12);
    212     fprintf (stderr, "ext analysis skipped @ 13 : %d\n", Nskip13);
    213     fprintf (stderr, "ext analysis skipped @ 14 : %d\n", Nskip14);
    214184
    215185    psphotRadialProfileShowSkips ();
    … …  
    217187    psphotVisualShowResidualImage (readout, false);
    218188
    219     bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
    220     if (doPetrosian) {
    221         psphotVisualShowPetrosians (sources);
    222     }
     189    psphotVisualShowPetrosians (sources);
    223190
    224191    return true;
    … …  
    228195
    229196    bool status;
    230 
    231     int NfaintEXT = 0;
    232     int NfaintPSF = 0;
    233197
    234198    int Next = 0;
    … …  
    244208    float skynoise          = PS_SCALAR_VALUE(job->args->data[4],F32);
    245209
    246     // S/N limit to perform full non-linear fits
    247     float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");
    248 
    249     // which extended source analyses should we perform?
    250     bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
    251     bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
    252     bool doPetroStars   = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS");
    253 
    254210    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    255211    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    … …  
    263219        Nall ++;
    264220
    265         // if we have checked the source validity on the basis of the object set, then
    266         // we either skip these tests below or we skip the source completely
    267         if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_SKIP) SKIP (Nskip1);
    268         if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_KEEP) goto keepSource;
    269 
    270         // skip PSF-like and non-astronomical objects
    271         if (source->type == PM_SOURCE_TYPE_DEFECT) SKIP (Nskip2);
    272         if (source->type == PM_SOURCE_TYPE_SATURATED) SKIP (Nskip3);
    273         if (source->mode & PM_SOURCE_MODE_DEFECT) SKIP (Nskip4);
    274         if (source->mode & PM_SOURCE_MODE_SATSTAR) SKIP (Nskip5);
    275 
    276         // skip saturated stars modeled with a radial profile
    277         if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) SKIP (Nskip6);
    278 
    279         // optionally allow non-extended objects to get petrosians as well
    280         if (!doPetroStars) {
    281             if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) SKIP (Nskip7);
    282             // if (source->type == PM_SOURCE_TYPE_STAR) SKIP (Nskip8); -- XXX this might be set for input lists (instead of EXT_LIMIT)
    283         }
    284 
    285         // limit selection to some SN limit
    286         // assert (source->peak); // how can a source not have a peak?
    287         // limit selection to some SN limit
    288         bool skipSource = false;
    289         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
    290             skipSource = (source->moments->KronFlux < SN_LIM * source->moments->KronFluxErr);
    291             NfaintEXT ++;
    292         } else {
    293             skipSource = (sqrt(source->peak->detValue) < SN_LIM);
    294             NfaintPSF ++;
    295         }
    296         if (skipSource) SKIP (Nskip9);
    297 
    298         // limit selection by analysis region (this automatically apply
    299         if (source->peak->x < region->x0) SKIP (Nskip10);
    300         if (source->peak->y < region->y0) SKIP (Nskip11);
    301         if (source->peak->x > region->x1) SKIP (Nskip12);
    302         if (source->peak->y > region->y1) SKIP (Nskip13);
    303 
    304     keepSource:
     221        // rules for measuring petrosian parameters for specific objects are set in
     222        // psphotChooseAnalysisOptions.c
     223        if (!(source->tmpFlags & PM_SOURCE_TMPF_PETRO)) SKIP (Nskip1);
     224
     225        // limit selection by analysis region (XXX move this into psphotChooseAnalysisOption?)
     226        if (source->peak->x < region->x0) SKIP (Nskip2);
     227        if (source->peak->y < region->y0) SKIP (Nskip3);
     228        if (source->peak->x > region->x1) SKIP (Nskip4);
     229        if (source->peak->y > region->y1) SKIP (Nskip5);
    305230
    306231        // replace object in image
    … …  
    317242        pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius);
    318243
    319         // if we request any of these measurements, we require the radial profile
    320         if (doPetrosian || doAnnuli) {
    321             if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) {
    322                 // all measurements below require the radial profile; skip them all
    323                 // re-subtract the object, leave local sky
    324                 psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    325                 pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    326                 SKIP (Nskip14);
    327             }
    328             Nannuli ++;
    329             source->mode |= PM_SOURCE_MODE_RADIAL_FLUX;
    330         }
     244        // measure the radial profile
     245        if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) {
     246          // re-subtract the object, leave local sky
     247          psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
     248          pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     249          SKIP (Nskip6);
     250        }
     251
     252        Nannuli ++;
     253        source->mode |= PM_SOURCE_MODE_RADIAL_FLUX;
    331254
    332255        // Petrosian Mags
    333         if (doPetrosian) {
    334             if (!psphotPetrosian (source, recipe, skynoise, maskVal)) {
    335                 psTrace ("psphot", 5, "FAILED petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    336             } else {
    337                 psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    338                 Npetro ++;
    339                 source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
    340             }
     256        if (!psphotPetrosian (source, recipe, skynoise, maskVal)) {
     257          psTrace ("psphot", 5, "FAILED petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
     258        } else {
     259          psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
     260          Npetro ++;
     261          source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
    341262        }
    342263
    … …  
    362283    scalar->data.S32 = Nannuli;
    363284
    364     scalar = job->args->data[8];
    365     scalar->data.S32 = NfaintEXT;
    366 
    367     scalar = job->args->data[9];
    368     scalar->data.S32 = NfaintPSF;
    369 
    370285    return true;
    371286}

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

    r36333 r36342  
    6464    psLogMsg("psphot", PS_LOG_INFO, "extended source fits for image %d", index);
    6565    psphotVisualShowImage(readout);
    66 
    67     // psphotSaveImage (NULL, readout->image, "test.01.fits");
    6866
    6967    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    … …  
    302300            psArrayAdd(job->args, 1, cells->data[j]); // sources
    303301            psArrayAdd(job->args, 1, models);
     302
    304303            // Allocate a metadata iterator here because psMetadataIteratorAlloc/Free are not thread safe
    305304            psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL);
    … …  
    321320
    322321// XXX TEST
    323             if (true && !isInteractive) {
     322            if (!isInteractive) {
    324323                if (!psThreadJobAddPending(job)) {
    325324                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    … …  
    403402    psphotSersicModelClassCleanup();
    404403
    405     // psphotSaveImage (NULL, readout->image, "test.02.fits");
    406 
    407404    psphotVisualShowResidualImage (readout, false);
    408405
    … …  
    458455        pmSource *source = sources->data[i];
    459456
    460         // skip PSF-like and non-astronomical objects
    461         if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
    462         if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
    463         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
    464 
    465         // skip saturated stars modeled with a radial profile
    466         if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
    467 
    468         // XXX this should use peak?
     457        // rules for measuring petrosian parameters for specific objects are set in
     458        // psphotChooseAnalysisOptions.c
     459        if (!(source->tmpFlags & PM_SOURCE_TMPF_EXT_FIT)) continue;
     460
     461        // limit selection by analysis region (XXX move this into psphotChooseAnalysisOption?)
    469462        if (source->peak->x < region->x0) continue;
    470463        if (source->peak->y < region->y0) continue;
    471464        if (source->peak->x > region->x1) continue;
    472465        if (source->peak->y > region->y1) continue;
    473 
    474 
    475         // XXX for a test, just do the obvious trail
    476         // XXX if (source->peak->xf < 1100) continue;
    477         // XXX if (source->peak->xf > 1400) continue;
    478         // XXX if (source->peak->yf >  245) continue;
    479466
    480467        // replace object in image
    … …  
    487474        psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal); // NOTE : 6 allocs
    488475
    489         // XXX WATCH OUT HERE!!
    490         // fitRadius = 30;
    491 
    492476        // UPDATE : we have changed the moments calculation.  There is now an iteration within
    493477        // psphotKronMasked to determine moments appropriate for a larger object.  The values
    … …  
    502486# ifdef TEST_OBJECT
    503487        bool testObject = false;
    504         // testObject |= ((fabs(source->peak->xf -  179) < 5) && (fabs(source->peak->yf - 1138) < 5));
    505         // testObject |= ((fabs(source->peak->xf - 5047) < 5) && (fabs(source->peak->yf -  151) < 5));
    506         // testObject |= ((fabs(source->peak->xf - 3929) < 5) && (fabs(source->peak->yf - 4109) < 5));
    507         // testObject |= ((fabs(source->peak->xf -  915) < 5) && (fabs(source->peak->yf - 5998) < 5));
    508         // testObject |= ((fabs(source->peak->xf - 5406) < 5) && (fabs(source->peak->yf -  326) < 5));
     488        testObject |= ((fabs(source->peak->xf -  179) < 5) && (fabs(source->peak->yf - 1138) < 5));
    509489        if (testObject) {
    510490            fprintf (stderr, "test object @ %f, %f\n", source->peak->xf, source->peak->yf);
    … …  
    528508          assert (status);
    529509
    530           // limit selection to some SN limit
    531           bool skipSource = false;
    532           if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
    533               skipSource = (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr);
    534           } else {
    535               skipSource = (sqrt(source->peak->detValue) < FIT_SN_LIM);
    536           }
    537           if (skipSource) {
     510          // limit selection to some SN limit for specific models (this value only applies if > EXTENDED_SOURCE_SN_LIM)
     511          if (isfinite(FIT_SN_LIM)) {
     512            if (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr) {
    538513              Nfaint ++;
    539514              continue;
     515            }
    540516          }
    541517
    … …  
    612588              }
    613589          }
    614           // XXX really need to do this in a cleaner way:
    615           // XXX I'm not sure I can get to this statement with a null modelFit
    616           if (!modelFit) continue;
     590          psAssert (modelFit, "modelFit not set?");
    617591
    618592          // test for fit quality / result
Note: See TracChangeset for help on using the changeset viewer.