Index: /branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- /branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceAnalysis.c (revision 36341) +++ /branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceAnalysis.c (revision 36342) @@ -2,5 +2,5 @@ void psphotRadialProfileShowSkips (); -// measure the elliptical radial profile and use this to measure the petrosian parameters for the sources +// measure the petrosian parameters for the sources // for now, let's store the detections on the readout->analysis for each readout @@ -16,6 +16,6 @@ psAssert (recipe, "missing recipe?"); - // perform full non-linear fits / extended source analysis? - if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANALYSIS")) { + // measure petrosians? + if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN")) { psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n"); return true; @@ -42,12 +42,4 @@ int Nskip5 = 0; int Nskip6 = 0; -int Nskip7 = 0; -int Nskip8 = 0; -int Nskip9 = 0; -int Nskip10 = 0; -int Nskip11 = 0; -int Nskip12 = 0; -int Nskip13 = 0; -int Nskip14 = 0; # define SKIP(VALUE) { VALUE++; continue; } @@ -57,7 +49,4 @@ bool status; - int NfaintEXT = 0; - int NfaintPSF = 0; - int Next = 0; int Npetro = 0; @@ -131,6 +120,4 @@ PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npetro PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nannuli - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NfaintEXT - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NfaintPSF // set this to 0 to run without threading @@ -154,8 +141,4 @@ scalar = job->args->data[7]; Nannuli += scalar->data.S32; - scalar = job->args->data[8]; - NfaintEXT += scalar->data.S32; - scalar = job->args->data[9]; - NfaintPSF += scalar->data.S32; psFree(job); # endif @@ -182,8 +165,4 @@ scalar = job->args->data[7]; Nannuli += scalar->data.S32; - scalar = job->args->data[8]; - NfaintEXT += scalar->data.S32; - scalar = job->args->data[9]; - NfaintPSF += scalar->data.S32; } psFree(job); @@ -196,5 +175,4 @@ psLogMsg ("psphot", PS_LOG_INFO, " %d petrosian\n", Npetro); psLogMsg ("psphot", PS_LOG_INFO, " %d annuli\n", Nannuli); - psLogMsg ("psphot", PS_LOG_INFO, " skipped: %d EXT, %d PSF\n", NfaintEXT, NfaintPSF); fprintf (stderr, "ext analysis skipped @ 1 : %d\n", Nskip1); @@ -204,12 +182,4 @@ fprintf (stderr, "ext analysis skipped @ 5 : %d\n", Nskip5); fprintf (stderr, "ext analysis skipped @ 6 : %d\n", Nskip6); - fprintf (stderr, "ext analysis skipped @ 7 : %d\n", Nskip7); - fprintf (stderr, "ext analysis skipped @ 8 : %d\n", Nskip8); - fprintf (stderr, "ext analysis skipped @ 9 : %d\n", Nskip9); - fprintf (stderr, "ext analysis skipped @ 10 : %d\n", Nskip10); - fprintf (stderr, "ext analysis skipped @ 11 : %d\n", Nskip11); - fprintf (stderr, "ext analysis skipped @ 12 : %d\n", Nskip12); - fprintf (stderr, "ext analysis skipped @ 13 : %d\n", Nskip13); - fprintf (stderr, "ext analysis skipped @ 14 : %d\n", Nskip14); psphotRadialProfileShowSkips (); @@ -217,8 +187,5 @@ psphotVisualShowResidualImage (readout, false); - bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - if (doPetrosian) { - psphotVisualShowPetrosians (sources); - } + psphotVisualShowPetrosians (sources); return true; @@ -228,7 +195,4 @@ bool status; - - int NfaintEXT = 0; - int NfaintPSF = 0; int Next = 0; @@ -244,12 +208,4 @@ float skynoise = PS_SCALAR_VALUE(job->args->data[4],F32); - // S/N limit to perform full non-linear fits - float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM"); - - // which extended source analyses should we perform? - bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - bool doPetroStars = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS"); - // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels @@ -263,44 +219,13 @@ Nall ++; - // if we have checked the source validity on the basis of the object set, then - // we either skip these tests below or we skip the source completely - if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_SKIP) SKIP (Nskip1); - if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_KEEP) goto keepSource; - - // skip PSF-like and non-astronomical objects - if (source->type == PM_SOURCE_TYPE_DEFECT) SKIP (Nskip2); - if (source->type == PM_SOURCE_TYPE_SATURATED) SKIP (Nskip3); - if (source->mode & PM_SOURCE_MODE_DEFECT) SKIP (Nskip4); - if (source->mode & PM_SOURCE_MODE_SATSTAR) SKIP (Nskip5); - - // skip saturated stars modeled with a radial profile - if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) SKIP (Nskip6); - - // optionally allow non-extended objects to get petrosians as well - if (!doPetroStars) { - if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) SKIP (Nskip7); - // if (source->type == PM_SOURCE_TYPE_STAR) SKIP (Nskip8); -- XXX this might be set for input lists (instead of EXT_LIMIT) - } - - // limit selection to some SN limit - // assert (source->peak); // how can a source not have a peak? - // limit selection to some SN limit - bool skipSource = false; - if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { - skipSource = (source->moments->KronFlux < SN_LIM * source->moments->KronFluxErr); - NfaintEXT ++; - } else { - skipSource = (sqrt(source->peak->detValue) < SN_LIM); - NfaintPSF ++; - } - if (skipSource) SKIP (Nskip9); - - // limit selection by analysis region (this automatically apply - if (source->peak->x < region->x0) SKIP (Nskip10); - if (source->peak->y < region->y0) SKIP (Nskip11); - if (source->peak->x > region->x1) SKIP (Nskip12); - if (source->peak->y > region->y1) SKIP (Nskip13); - - keepSource: + // rules for measuring petrosian parameters for specific objects are set in + // psphotChooseAnalysisOptions.c + if (!(source->tmpFlags & PM_SOURCE_TMPF_PETRO)) SKIP (Nskip1); + + // limit selection by analysis region (XXX move this into psphotChooseAnalysisOption?) + if (source->peak->x < region->x0) SKIP (Nskip2); + if (source->peak->y < region->y0) SKIP (Nskip3); + if (source->peak->x > region->x1) SKIP (Nskip4); + if (source->peak->y > region->y1) SKIP (Nskip5); // replace object in image @@ -317,26 +242,22 @@ pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius); - // if we request any of these measurements, we require the radial profile - if (doPetrosian || doAnnuli) { - if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) { - // all measurements below require the radial profile; skip them all - // re-subtract the object, leave local sky - psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - SKIP (Nskip14); - } - Nannuli ++; - source->mode |= PM_SOURCE_MODE_RADIAL_FLUX; - } + // measure the radial profile + if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) { + // re-subtract the object, leave local sky + psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + SKIP (Nskip6); + } + + Nannuli ++; + source->mode |= PM_SOURCE_MODE_RADIAL_FLUX; // Petrosian Mags - if (doPetrosian) { - if (!psphotPetrosian (source, recipe, skynoise, maskVal)) { - psTrace ("psphot", 5, "FAILED petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); - } else { - psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); - Npetro ++; - source->mode |= PM_SOURCE_MODE_EXTENDED_STATS; - } + if (!psphotPetrosian (source, recipe, skynoise, maskVal)) { + psTrace ("psphot", 5, "FAILED petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); + } else { + psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); + Npetro ++; + source->mode |= PM_SOURCE_MODE_EXTENDED_STATS; } @@ -362,10 +283,4 @@ scalar->data.S32 = Nannuli; - scalar = job->args->data[8]; - scalar->data.S32 = NfaintEXT; - - scalar = job->args->data[9]; - scalar->data.S32 = NfaintPSF; - return true; } Index: /branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceFits.c =================================================================== --- /branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceFits.c (revision 36341) +++ /branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceFits.c (revision 36342) @@ -64,6 +64,4 @@ psLogMsg("psphot", PS_LOG_INFO, "extended source fits for image %d", index); psphotVisualShowImage(readout); - - // psphotSaveImage (NULL, readout->image, "test.01.fits"); pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); @@ -302,4 +300,5 @@ psArrayAdd(job->args, 1, cells->data[j]); // sources psArrayAdd(job->args, 1, models); + // Allocate a metadata iterator here because psMetadataIteratorAlloc/Free are not thread safe psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL); @@ -321,5 +320,5 @@ // XXX TEST - if (true && !isInteractive) { + if (!isInteractive) { if (!psThreadJobAddPending(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); @@ -403,6 +402,4 @@ psphotSersicModelClassCleanup(); - // psphotSaveImage (NULL, readout->image, "test.02.fits"); - psphotVisualShowResidualImage (readout, false); @@ -458,23 +455,13 @@ pmSource *source = sources->data[i]; - // skip PSF-like and non-astronomical objects - if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue; - if (source->type == PM_SOURCE_TYPE_DEFECT) continue; - if (source->type == PM_SOURCE_TYPE_SATURATED) continue; - - // skip saturated stars modeled with a radial profile - if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; - - // XXX this should use peak? + // rules for measuring petrosian parameters for specific objects are set in + // psphotChooseAnalysisOptions.c + if (!(source->tmpFlags & PM_SOURCE_TMPF_EXT_FIT)) continue; + + // limit selection by analysis region (XXX move this into psphotChooseAnalysisOption?) if (source->peak->x < region->x0) continue; if (source->peak->y < region->y0) continue; if (source->peak->x > region->x1) continue; if (source->peak->y > region->y1) continue; - - - // XXX for a test, just do the obvious trail - // XXX if (source->peak->xf < 1100) continue; - // XXX if (source->peak->xf > 1400) continue; - // XXX if (source->peak->yf > 245) continue; // replace object in image @@ -487,7 +474,4 @@ psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal); // NOTE : 6 allocs - // XXX WATCH OUT HERE!! - // fitRadius = 30; - // UPDATE : we have changed the moments calculation. There is now an iteration within // psphotKronMasked to determine moments appropriate for a larger object. The values @@ -502,9 +486,5 @@ # ifdef TEST_OBJECT bool testObject = false; - // testObject |= ((fabs(source->peak->xf - 179) < 5) && (fabs(source->peak->yf - 1138) < 5)); - // testObject |= ((fabs(source->peak->xf - 5047) < 5) && (fabs(source->peak->yf - 151) < 5)); - // testObject |= ((fabs(source->peak->xf - 3929) < 5) && (fabs(source->peak->yf - 4109) < 5)); - // testObject |= ((fabs(source->peak->xf - 915) < 5) && (fabs(source->peak->yf - 5998) < 5)); - // testObject |= ((fabs(source->peak->xf - 5406) < 5) && (fabs(source->peak->yf - 326) < 5)); + testObject |= ((fabs(source->peak->xf - 179) < 5) && (fabs(source->peak->yf - 1138) < 5)); if (testObject) { fprintf (stderr, "test object @ %f, %f\n", source->peak->xf, source->peak->yf); @@ -528,14 +508,10 @@ assert (status); - // limit selection to some SN limit - bool skipSource = false; - if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { - skipSource = (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr); - } else { - skipSource = (sqrt(source->peak->detValue) < FIT_SN_LIM); - } - if (skipSource) { + // limit selection to some SN limit for specific models (this value only applies if > EXTENDED_SOURCE_SN_LIM) + if (isfinite(FIT_SN_LIM)) { + if (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr) { Nfaint ++; continue; + } } @@ -612,7 +588,5 @@ } } - // XXX really need to do this in a cleaner way: - // XXX I'm not sure I can get to this statement with a null modelFit - if (!modelFit) continue; + psAssert (modelFit, "modelFit not set?"); // test for fit quality / result