Index: branches/eam_branches/ipp-20110906/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- branches/eam_branches/ipp-20110906/psphot/src/psphotExtendedSourceAnalysis.c (revision 32619) +++ branches/eam_branches/ipp-20110906/psphot/src/psphotExtendedSourceAnalysis.c (revision 32628) @@ -209,6 +209,6 @@ // 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_STACK_SKIP) continue; - if (source->tmpFlags & PM_SOURCE_TMPF_STACK_KEEP) goto keepSource; + if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_SKIP) continue; + if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_KEEP) goto keepSource; // skip PSF-like and non-astronomical objects Index: branches/eam_branches/ipp-20110906/psphot/src/psphotRadialApertures.c =================================================================== --- branches/eam_branches/ipp-20110906/psphot/src/psphotRadialApertures.c (revision 32619) +++ branches/eam_branches/ipp-20110906/psphot/src/psphotRadialApertures.c (revision 32628) @@ -26,6 +26,12 @@ int num = psphotFileruleCount(config, filerule); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image + if (!psphotRadialAperturesReadout (config, view, filerule, i, recipe, entry)) { psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for %s entry %d", filerule, i); @@ -218,4 +224,9 @@ pmSource *source = sources->data[i]; + // 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_RADIAL_SKIP) continue; + if (source->tmpFlags & PM_SOURCE_TMPF_RADIAL_KEEP) goto keepSource; + // skip PSF-like and non-astronomical objects if (source->type == PM_SOURCE_TYPE_DEFECT) continue; @@ -235,4 +246,6 @@ if (source->peak->x > region->x1) continue; if (source->peak->y > region->y1) continue; + + keepSource: // allocate pmSourceExtendedParameters, if not already defined Index: branches/eam_branches/ipp-20110906/psphot/src/psphotSourceMatch.c =================================================================== --- branches/eam_branches/ipp-20110906/psphot/src/psphotSourceMatch.c (revision 32619) +++ branches/eam_branches/ipp-20110906/psphot/src/psphotSourceMatch.c (revision 32628) @@ -67,6 +67,9 @@ objects = psArraySort (objects, pmPhotObjSortByX); - psVector *found = psVectorAlloc(sources->n, PS_TYPE_U8); - psVectorInit (found, 0); + psVector *foundSrc = psVectorAlloc(sources->n, PS_TYPE_U8); + psVectorInit (foundSrc, 0); + + psVector *foundObj = psVectorAlloc(sources->n, PS_TYPE_U8); + psVectorInit (foundObj, 0); // match sources to existing objects @@ -93,4 +96,11 @@ if (dx > +1.02*RADIUS) NEXT2; + /* this block will match a given detection to the closest object within range of that detection. + XXX note that this matches ALL detections within range of the single object to that same object + this is bad, but I cannot just go in linear order (ie, mark off each object as they are + used). I should make a list of all Nobj * Ndet pairs in range and choose the matches + based on their separations. UGH + */ + // we are within match range, look for matches: int Jmin = -1; @@ -98,4 +108,6 @@ for (int J = j; (dx > -1.02*RADIUS) && (J < objects->n); J++) { + // skip objects that are already assigned: + if (foundObj->data.U8[J]) continue; obj = objects->data[J]; @@ -117,8 +129,9 @@ } obj = objects->data[Jmin]; + foundObj->data.U8[Jmin] = 1; // add to object pmPhotObjAddSource (obj, src); - found->data.U8[i] = 1; + foundSrc->data.U8[i] = 1; i++; } @@ -128,5 +141,5 @@ for (i = 0; i < sources->n; i++) { - if (found->data.U8[i]) continue; + if (foundSrc->data.U8[i]) continue; pmSource *src = sources->data[i]; @@ -139,5 +152,6 @@ psLogMsg ("psphot", PS_LOG_DETAIL, "matched sources (%ld vs %ld)", sources->n, objects->n); - psFree (found); + psFree (foundSrc); + psFree (foundObj); return true; } @@ -270,4 +284,5 @@ pmPhotObj *obj = objects->data[i]; nSources += obj->sources->n; + psAssert (obj->sources->n == nImages, "failed to match sources?"); } psLogMsg ("psphot", PS_LOG_DETAIL, "total of %d sources for %d images", nSources, nImages); Index: branches/eam_branches/ipp-20110906/psphot/src/psphotStackObjects.c =================================================================== --- branches/eam_branches/ipp-20110906/psphot/src/psphotStackObjects.c (revision 32619) +++ branches/eam_branches/ipp-20110906/psphot/src/psphotStackObjects.c (revision 32628) @@ -70,5 +70,6 @@ // S/N limit to perform full non-linear fits - float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM"); + float SN_LIM_PETRO = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM"); + float SN_LIM_RADIAL = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM"); bool doPetroStars = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS"); @@ -81,5 +82,6 @@ // we check each source for an object and keep the object if any source is valid - bool keepObject = false; + bool keepObjectRadial = false; + bool keepObjectPetro = false; for (int j = 0; j < object->sources->n; j++) { @@ -94,21 +96,4 @@ if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; - // optionally allow non-extended objects to get petrosians as well - if (!doPetroStars) { - if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue; - if (source->type == PM_SOURCE_TYPE_STAR) continue; - } - - // 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); - } else { - skipSource = (sqrt(source->peak->detValue) < SN_LIM); - } - if (skipSource) continue; - // limit selection by analysis region (this automatically apply if (source->peak->x < AnalysisRegion.x0) continue; @@ -117,5 +102,23 @@ if (source->peak->y > AnalysisRegion.y1) continue; - keepObject = true; + // SN limit tests for RADIAL APERTURES: + bool skipSourceRadial = false; + if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { + skipSourceRadial = (source->moments->KronFlux < SN_LIM_RADIAL * source->moments->KronFluxErr); + } else { + skipSourceRadial = (sqrt(source->peak->detValue) < SN_LIM_RADIAL); + } + if (!skipSourceRadial) keepObjectRadial = true; + + // SN limit tests for PETRO + bool skipSourcePetro = false; + if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { + skipSourcePetro = (source->moments->KronFlux < SN_LIM_PETRO * source->moments->KronFluxErr); + } else { + skipSourcePetro = doPetroStars ? (sqrt(source->peak->detValue) < SN_LIM_PETRO) : true; + } + if (!skipSourcePetro) keepObjectPetro = true; + + keepObjectPetro = true; } @@ -128,10 +131,18 @@ // avoid the single-detection tests - if (keepObject) { - source->tmpFlags |= PM_SOURCE_TMPF_STACK_KEEP; - source->tmpFlags &= ~PM_SOURCE_TMPF_STACK_SKIP; + if (keepObjectPetro) { + source->tmpFlags |= PM_SOURCE_TMPF_PETRO_KEEP; + source->tmpFlags &= ~PM_SOURCE_TMPF_PETRO_SKIP; } else { - source->tmpFlags |= PM_SOURCE_TMPF_STACK_SKIP; - source->tmpFlags &= ~PM_SOURCE_TMPF_STACK_KEEP; + source->tmpFlags |= PM_SOURCE_TMPF_PETRO_SKIP; + source->tmpFlags &= ~PM_SOURCE_TMPF_PETRO_KEEP; + } + + if (keepObjectRadial) { + source->tmpFlags |= PM_SOURCE_TMPF_RADIAL_KEEP; + source->tmpFlags &= ~PM_SOURCE_TMPF_RADIAL_SKIP; + } else { + source->tmpFlags |= PM_SOURCE_TMPF_RADIAL_SKIP; + source->tmpFlags &= ~PM_SOURCE_TMPF_RADIAL_KEEP; } } Index: branches/eam_branches/ipp-20110906/psphot/src/psphotStackReadout.c =================================================================== --- branches/eam_branches/ipp-20110906/psphot/src/psphotStackReadout.c (revision 32619) +++ branches/eam_branches/ipp-20110906/psphot/src/psphotStackReadout.c (revision 32628) @@ -43,4 +43,6 @@ bool psphotStackReadout (pmConfig *config, const pmFPAview *view) { + + psArray *objects = NULL; // used below after 'pass1finish' label // measure the total elapsed time in psphotReadout. dtime is the elapsed time used jointly @@ -187,10 +189,10 @@ psphotReplaceAllSources (config, view, STACK_SRC); // pass 1 (detections->allSources) + // if we only do one pass, skip to extended source analysis + if (!strcasecmp (breakPt, "PASS1")) goto pass1finish; + // linear fit to include all sources (subtract again) // NOTE : apply to ALL sources (extended + psf) psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 2 (detections->allSources) - - // if we only do one pass, skip to extended source analysis - if (!strcasecmp (breakPt, "PASS1")) goto pass1finish; // NOTE: possibly re-measure background model here with objects subtracted / or masked @@ -241,10 +243,20 @@ // XXX check on free of sources... psphotMergeSources (config, view, STACK_SRC); // (detections->newSources + detections->allSources -> detections->allSources) - - // NOTE: apply to ALL sources - psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 3 (detections->allSources) } pass1finish: + + // generate the objects (objects unify the sources from the different images) NOTE: could + // this just match the detections for the chisq image, and not bother measuring the source + // stats in that case...? + objects = psphotMatchSources (config, view, STACK_SRC); + psMemDump("matchsources"); + + psphotStackObjectsUnifyPosition (objects); + + psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects); + + // NOTE: apply to ALL sources + psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 3 (detections->allSources) // measure the radial profiles to the sky (only measures new objects) @@ -260,17 +272,4 @@ psMemDump("psfstats"); - - // XXX *** NOTE : this should happen above before the last psphotFitSourcesLinear in order - // to generate psf forced photometry for the undetected sources - - // generate the objects (objects unify the sources from the different images) - // XXX this could just match the detections for the chisq image, and not bother measuring the - // source stats in that case... - psArray *objects = psphotMatchSources (config, view, STACK_SRC); - psMemDump("matchsources"); - - psphotStackObjectsUnifyPosition (objects); - - psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects); // measure elliptical apertures, petrosians (objects sorted by S/N) @@ -330,5 +329,7 @@ psphotMagnitudes(config, view, STACK_SRC); - if (0 && !psphotEfficiency(config, view, STACK_DET)) { + // XXX NOTE: this function wants to have the PSF of the image, but we (so far) only measure the + // PSF of the SRC image. can we fake it by generating the PSF for DET as well (up above)? + if (false && !psphotEfficiency(config, view, STACK_DET)) { psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources"); psErrorClear();