Index: /branches/eam_branches/ipp-20120627/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphot.h (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphot.h (revision 34247) @@ -102,6 +102,6 @@ bool psphotMergeSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index); -bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final); -bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode); +bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final, bool skipNegativeFluxSources); +bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode, bool skipNegativeFluxSources); bool psphotSourceSize (pmConfig *config, const pmFPAview *view, const char *filerule, bool getPSFsize); @@ -265,4 +265,5 @@ bool psphotEllipticalProfile (pmSource *source, bool RAW_RADIUS); bool psphotEllipticalContour (pmSource *source); +bool psphotLimitRadialApertures(psMetadata *recipe, long nSources); // psphotVisual functions @@ -423,4 +424,6 @@ // bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index); +bool psphotStackSetInputsToSkip(pmConfig *config, const pmFPAview *view, const char *filerule, bool set) ; + bool psphotStackRenormaliseVariance(const pmConfig *config, pmReadout *readout); Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotChoosePSF.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotChoosePSF.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotChoosePSF.c (revision 34247) @@ -47,4 +47,9 @@ if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) { psLogMsg ("psphot", PS_LOG_DETAIL, "psf model supplied for input file %d", index); + return true; + } + + if (psMetadataLookupBool (&status, readout->analysis, "PSPHOT.SKIP.INPUT")) { + psLogMsg ("psphot", PS_LOG_DETAIL, "skipping choose PSF for input file %d", index); return true; } Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotCullPeaks.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotCullPeaks.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotCullPeaks.c (revision 34247) @@ -63,4 +63,15 @@ # if (PM_PEAKS_CULL_WITH_SMOOTHED_IMAGE) psLogMsg ("psphot", PS_LOG_INFO, "Culling peaks from footprints using the smoothed image"); + bool setNaNsToZero = psMetadataLookupF32 (&status, recipe, "FOOTPRINT_SET_NANS_TO_ZERO"); + if (setNaNsToZero) { + // Set NaN pixels to zero so that they do not considered in the footprint analysis + // XXX: Currently the caller does nothing further with the significance readout so + // modifiying the image does no harm + // XXX: Note: signifR->weight is not used by pmFootprintCullPeaks so we don't + // need to touch it + psImageClipNaN(signifR->image, 0); + } + + # else psLogMsg ("psphot", PS_LOG_INFO, "Culling peaks from footprints using the raw (unsmoothed) image"); @@ -73,5 +84,5 @@ if (fp->npix > 30000) { - fprintf (stderr, "big footprint: %f %f to %f %f (%d pix)\n", fp->bbox.x0, fp->bbox.y0, fp->bbox.x1, fp->bbox.y1, fp->npix); + psLogMsg("psphot.cull.footprints", PS_LOG_WARN, "big footprint: %f %f to %f %f (%d pix)\n", fp->bbox.x0, fp->bbox.y0, fp->bbox.x1, fp->bbox.y1, fp->npix); } psTimerStart ("psphot.cull.footprints"); @@ -97,7 +108,8 @@ # endif + float dtime = psTimerMark ("psphot.cull.footprints"); if (dtime > 1.0) { - fprintf (stderr, "slow cull for %d (%f sec)\n", i, dtime); + psLogMsg("psphot.cull.footprints", PS_LOG_WARN, "slow cull for %d (%f sec)\n", i, dtime); } } Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotEfficiency.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotEfficiency.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotEfficiency.c (revision 34247) @@ -420,5 +420,5 @@ // psphotFitSourcesLinearReadout subtracts the model fits - if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true, PM_SOURCE_PHOTFIT_CONST)) { + if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true, PM_SOURCE_PHOTFIT_CONST, false)) { psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources."); psFree(fakeSources); Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotExtendedSourceFits.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotExtendedSourceFits.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotExtendedSourceFits.c (revision 34247) @@ -175,4 +175,7 @@ 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); + psArrayAdd(job->args, 1, iter); psArrayAdd(job->args, 1, AnalysisRegion); // XXX make a pointer @@ -203,18 +206,19 @@ } psScalar *scalar = NULL; - scalar = job->args->data[7]; + scalar = job->args->data[8]; Next += scalar->data.S32; - scalar = job->args->data[8]; + scalar = job->args->data[9]; Nconvolve += scalar->data.S32; - scalar = job->args->data[9]; + scalar = job->args->data[10]; NconvolvePass += scalar->data.S32; - scalar = job->args->data[10]; + scalar = job->args->data[11]; Nplain += scalar->data.S32; - scalar = job->args->data[11]; + scalar = job->args->data[12]; NplainPass += scalar->data.S32; - scalar = job->args->data[12]; + scalar = job->args->data[13]; Nfaint += scalar->data.S32; - scalar = job->args->data[13]; + scalar = job->args->data[14]; Nfail += scalar->data.S32; + psFree(job->args->data[3]); // iterator allocated above psFree(job); # endif @@ -235,18 +239,19 @@ } else { psScalar *scalar = NULL; - scalar = job->args->data[7]; + scalar = job->args->data[8]; Next += scalar->data.S32; - scalar = job->args->data[8]; + scalar = job->args->data[9]; Nconvolve += scalar->data.S32; - scalar = job->args->data[9]; + scalar = job->args->data[10]; NconvolvePass += scalar->data.S32; - scalar = job->args->data[10]; + scalar = job->args->data[11]; Nplain += scalar->data.S32; - scalar = job->args->data[11]; + scalar = job->args->data[12]; NplainPass += scalar->data.S32; - scalar = job->args->data[12]; + scalar = job->args->data[13]; Nfaint += scalar->data.S32; - scalar = job->args->data[13]; + scalar = job->args->data[14]; Nfail += scalar->data.S32; + psFree(job->args->data[3]); // metadata iterator allocated above } psFree(job); @@ -285,8 +290,9 @@ psArray *sources = job->args->data[1]; psMetadata *models = job->args->data[2]; - psRegion *region = job->args->data[3]; - int psfSize = PS_SCALAR_VALUE(job->args->data[4],S32); - psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA); - psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA); + psMetadataIterator *iter = job->args->data[3]; + psRegion *region = job->args->data[4]; + int psfSize = PS_SCALAR_VALUE(job->args->data[5],S32); + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA); // Define source fitting parameters for extended source fits @@ -338,5 +344,6 @@ // loop here over the models chosen for each source (exclude by S/N) - psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL); + // Reset the iterator + psMetadataIteratorSet(iter, PS_LIST_HEAD); psMetadataItem *item = NULL; while ((item = psMetadataGetAndIncrement (iter)) != NULL) { @@ -428,5 +435,4 @@ psFree (modelFit); } - psFree (iter); // evaluate the relative quality of the models, choose one @@ -489,23 +495,23 @@ // change the value of a scalar on the array (wrap this and put it in psArray.h) - scalar = job->args->data[7]; + scalar = job->args->data[8]; scalar->data.S32 = Next; - scalar = job->args->data[8]; + scalar = job->args->data[9]; scalar->data.S32 = Nconvolve; - scalar = job->args->data[9]; + scalar = job->args->data[10]; scalar->data.S32 = NconvolvePass; - scalar = job->args->data[10]; + scalar = job->args->data[11]; scalar->data.S32 = Nplain; - scalar = job->args->data[11]; + scalar = job->args->data[12]; scalar->data.S32 = NplainPass; - scalar = job->args->data[12]; + scalar = job->args->data[13]; scalar->data.S32 = Nfaint; - scalar = job->args->data[13]; + scalar = job->args->data[14]; scalar->data.S32 = Nfail; Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotFitSourcesLinear.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotFitSourcesLinear.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotFitSourcesLinear.c (revision 34247) @@ -13,5 +13,5 @@ // for now, let's store the detections on the readout->analysis for each readout -bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final) +bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final, bool skipNegativeFluxSources) { bool status = true; @@ -51,4 +51,10 @@ psAssert (readout, "missing readout?"); + if (psMetadataLookupBool (&status, readout->analysis, "PSPHOT.SKIP.INPUT")) { + psLogMsg ("psphot", PS_LOG_DETAIL, "skipping fit sources for input file %d", i); + continue; + } + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); psAssert (detections, "missing detections?"); @@ -60,5 +66,5 @@ psAssert (psf, "missing psf?"); - if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1)) { + if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1, skipNegativeFluxSources)) { psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i); return false; @@ -79,5 +85,6 @@ // rerun fit with correct fitVarMode - if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarMode)) { + // XXX: does skipNegativeFlux work here? + if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarMode, skipNegativeFluxSources)) { psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i); return false; @@ -130,5 +137,5 @@ } -bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode) { +bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode, bool skipNegativeFluxSources) { bool status; float x; @@ -177,4 +184,7 @@ float MIN_VALID_FLUX = psMetadataLookupF32(&status, recipe, "PSF_FIT_MIN_VALID_FLUX"); if (!status) { + MIN_VALID_FLUX = 0.0; + } + if (skipNegativeFluxSources && MIN_VALID_FLUX < 0) { MIN_VALID_FLUX = 0.0; } @@ -252,9 +262,9 @@ if (modelSum < 0.5) continue; // skip sources with no model constraint (somewhat arbitrary limit) if (modelSum < 0.8) { - fprintf (stderr, "low-sig model @ %f, %f (%f masked sum, %f sum, %f peak)\n", + psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "low-sig model @ %f, %f (%f masked sum, %f sum, %f peak)\n", source->peak->xf, source->peak->yf, maskedSum, modelSum, source->peak->rawFlux); } if (maskedSum < 0.5) { - fprintf (stderr, "worrying model @ %f, %f (%f masked sum, %f sum, %f peak)\n", + psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "worrying model @ %f, %f (%f masked sum, %f sum, %f peak)\n", source->peak->xf, source->peak->yf, maskedSum, modelSum, source->peak->rawFlux); } Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotForcedReadout.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotForcedReadout.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotForcedReadout.c (revision 34247) @@ -64,5 +64,5 @@ // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) - psphotFitSourcesLinear (config, view, filerule, false); + psphotFitSourcesLinear (config, view, filerule, false, false); // identify CRs and extended sources Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotRadialApertures.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotRadialApertures.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotRadialApertures.c (revision 34247) @@ -71,4 +71,10 @@ pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); psAssert (readout, "missing readout?"); + + if (psMetadataLookupBool (&status, readout->analysis, "PSPHOT.SKIP.INPUT")) { + psLogMsg ("psphot", PS_LOG_DETAIL, "skipping radial aptertures for input file %d", index); + return true; + } + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotRadialBins.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotRadialBins.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotRadialBins.c (revision 34247) @@ -202,4 +202,69 @@ } +// If the number of sources is large, edit the recipe to remove radial bins beyond +// a value specified in the recipe +bool psphotLimitRadialApertures(psMetadata *recipe, long nSources) { + + bool status = false; + long sourceLimit = psMetadataLookupS32 (&status, recipe, "RADIAL.NUM.SOURCES.LIMIT"); + if (!status) { + sourceLimit = 50000; + } + if (nSources < sourceLimit) { + return true; + } + psF32 maxRadius = psMetadataLookupF32 (&status, recipe, "RADIAL.SOURCES.OVER.LIMIT.RADIUS"); + if (!status) { + maxRadius = 20.0; + } + psLogMsg ("psphot", PS_LOG_INFO, "Number of objects: %ld is greater than limit %ld. Limiting radial annular bins to %.3f pixels\n", + nSources, sourceLimit, maxRadius); + + psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); + psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + if (!radMin || !radMin->n) { + psError (PSPHOT_ERR_CONFIG, true, "error in definition of annular bins (radMin missing or empty)"); + return false; + } + if (!radMax || !radMax->n) { + psError (PSPHOT_ERR_CONFIG, true, "error in definition of annular bins (radMax missing or empty)"); + return false; + } + if (radMax->n != radMin->n) { + psError (PSPHOT_ERR_CONFIG, true, "length of radMin %ld and radMax %ld not equal)", radMin->n, radMax->n); + return false; + } + int i = 0; + psF32 lastRadius = 0; + for (; i < radMax->n; i++) { + if (radMax->data.F32[i] > maxRadius) { + break; + } + lastRadius = radMax->data.F32[i]; + } + if (i == radMax->n) { + psLogMsg ("psphot", PS_LOG_INFO, "Radius of all bins is within the limit. lastRadius: %.3f\n", lastRadius); + return true; + } + psLogMsg ("psphot", PS_LOG_INFO, " radius limit exceeded at bin %d of %ld. New max radius is %.3f\n", + i, radMax->n, lastRadius); + + psVector *radMinNew = psVectorRealloc(radMin, i); + if (radMinNew != radMin) { + psMetadataAddVector(recipe, PS_LIST_TAIL, "RADIAL.ANNULAR.BINS.LOWER", PS_META_REPLACE, "", radMinNew); + // XXX: I don't need this so I? + // psFree(radMinNew); + } + + psVector *radMaxNew = psVectorRealloc(radMax, i); + if (radMaxNew != radMax) { + psMetadataAddVector(recipe, PS_LIST_TAIL, "RADIAL.ANNULAR.BINS.UPPER", PS_META_REPLACE, "", radMaxNew); + // XXX: I don't need to free this do I? + // psFree(radMaxNew); + } + + return true; +} + // the area-weighted mean radius is given by: Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotReadout.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotReadout.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotReadout.c (revision 34247) @@ -185,5 +185,5 @@ // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) - psphotFitSourcesLinear (config, view, filerule, false); // pass 1 (detections->allSources) + psphotFitSourcesLinear (config, view, filerule, false, false); // pass 1 (detections->allSources) // measure the radial profiles to the sky @@ -212,5 +212,5 @@ // linear fit to include all sources (subtract again) // NOTE : apply to ALL sources (extended + psf) - psphotFitSourcesLinear (config, view, filerule, true); // pass 2 (detections->allSources) + psphotFitSourcesLinear (config, view, filerule, true, false); // pass 2 (detections->allSources) // if we only do one pass, skip to extended source analysis @@ -261,5 +261,5 @@ // NOTE: apply to ALL sources - psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources) + psphotFitSourcesLinear (config, view, filerule, true, false); // pass 3 (detections->allSources) } @@ -302,5 +302,5 @@ // NOTE: apply to ALL sources - psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources) + psphotFitSourcesLinear (config, view, filerule, true, false); // pass 3 (detections->allSources) } Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotReadoutForcedKnownSources.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotReadoutForcedKnownSources.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotReadoutForcedKnownSources.c (revision 34247) @@ -45,5 +45,5 @@ // linear PSF fit to source peaks - psphotFitSourcesLinear (config, view, filerule, false); + psphotFitSourcesLinear (config, view, filerule, false, false); // calculate source magnitudes Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotReadoutKnownSources.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotReadoutKnownSources.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotReadoutKnownSources.c (revision 34247) @@ -57,5 +57,5 @@ // linear PSF fit to source peaks - psphotFitSourcesLinear (config, view, filerule, false); + psphotFitSourcesLinear (config, view, filerule, false, false); // measure aperture photometry corrections Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotReadoutMinimal.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotReadoutMinimal.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotReadoutMinimal.c (revision 34247) @@ -81,5 +81,5 @@ // linear PSF fit to source peaks - psphotFitSourcesLinear (config, view, filerule, false); + psphotFitSourcesLinear (config, view, filerule, false, false); // measure the radial profiles to the sky Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotReplaceUnfit.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotReplaceUnfit.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotReplaceUnfit.c (revision 34247) @@ -55,4 +55,9 @@ pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); psAssert (readout, "missing readout?"); + + if (psMetadataLookupBool (&status, readout->analysis, "PSPHOT.SKIP.INPUT")) { + psLogMsg ("psphot", PS_LOG_DETAIL, "skipping replace all sources for input file %d", index); + return true; + } pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); @@ -304,4 +309,10 @@ psAssert (readout, "missing readout?"); + if (psMetadataLookupBool (&status, readout->analysis, "PSPHOT.SKIP.INPUT")) { + psLogMsg ("psphot", PS_LOG_DETAIL, "skipping reset models for input file %d", index); + return true; + } + + // XXX the sources have already been copied (merge into here?) pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotSetThreads.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotSetThreads.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotSetThreads.c (revision 34247) @@ -40,5 +40,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 14); + task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 15); task->function = &psphotExtendedSourceFits_Threaded; psThreadTaskAdd(task); Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotSourceSize.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotSourceSize.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotSourceSize.c (revision 34247) @@ -18,4 +18,5 @@ float sizeLimitCR; float magLimitCR; + int maxWindowCR; } psphotSourceSizeOptions; @@ -26,5 +27,5 @@ bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options); bool psphotSourceSelectCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options); -bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal); +bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal, int maxWindowCR); bool psphotMaskCosmicRayFootprintCheck (psArray *sources); int psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh); @@ -139,4 +140,9 @@ if (!status) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSPHOT.CRMASK.APPLY is not defined."); + return false; + } + options.maxWindowCR = psMetadataLookupS32 (&status, recipe, "PSPHOT.CR.MAX.WINDOW"); + if (!status) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSPHOT.CR.MAX.WINDOW is not defined."); return false; } @@ -647,5 +653,5 @@ psTrace("psphot", 6, "mask cosmic ray at %f, %f\n", source->peak->xf, source->peak->yf); if (options->applyCRmask) { - psphotMaskCosmicRay(readout, source, options->crMask); + psphotMaskCosmicRay(readout, source, options->crMask, options->maxWindowCR); } else { source->mode |= PM_SOURCE_MODE_CR_LIMIT; @@ -687,5 +693,5 @@ // does no repair or recovery of the CR pixels, it only masks them out. My test code can be // found at /data/ipp031.0/watersc1/psphot.20091209/algo_check.c -bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal) { +bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal, int maxWindowCR) { // Get the actual images and information about the peak. @@ -699,4 +705,16 @@ int ys = footprint->bbox.y0; int ye = footprint->bbox.y1 + 1; + + // We occasionally get very large footprints. When the footprint bounding box is large this function will + // do an incredible amount of work for no benefit. + // Limit the size of the area we examine. + if (xe - xs > maxWindowCR) { + xs = peak->x - maxWindowCR / 2; + xe = xs + maxWindowCR; + } + if (ye - ys > maxWindowCR) { + ys = peak->y - maxWindowCR / 2; + ye = ys + maxWindowCR; + } LIMIT_XRANGE(xs, mask); @@ -741,7 +759,10 @@ for (int i = 0; i < footprint->spans->n; i++) { pmSpan *sp = footprint->spans->data[i]; - for (int j = sp->x0; j <= sp->x1; j++) { - int y = sp->y - ys; - int x = j - xs; + int y = sp->y - ys; + if (y < 0 || y >= mymask->numRows) { + // can we break if y >= numRows? + continue; + } + for (int x = PS_MAX(sp->x0 - xs, 0); x <= PS_MIN(sp->x1 - xs, mymask->numCols-1); x++) { mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x01; } @@ -883,5 +904,6 @@ bool psphotMaskCosmicRayFootprintCheck (psArray *sources) { - +#ifdef CHECK_FOOTPRINTS + // This gets really expensive for complex images for (int i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; @@ -894,4 +916,5 @@ } } +#endif return true; } Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotStackMatchPSFs.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotStackMatchPSFs.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotStackMatchPSFs.c (revision 34247) @@ -87,9 +87,4 @@ } - // copy the header data from Src to Out - // pmHDU *hduSrc = pmHDUFromReadout(readoutSrc); - // psAssert (hduSrc, "input missing hdu?"); - - // set NAN pixels to 'SAT' psImageMaskType maskSat = pmConfigMaskGet("SAT", config); @@ -116,9 +111,14 @@ // save the output fwhm values in the readout->analysis. we may have / will have multiple output PSF sizes, // so we save this in a vector. if the vector is not yet defined, create it + // Skip this if the readout deconvolution fraction was over the limit. // NOTE: fwhmValues as defined here has 1 + nMatched PSF : 0 == unmatched psVector *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32); psVectorAppend(fwhmValues, NAN); // XXX this corresponds to the unmatched image set - for (int i = 0; i < options->targetSeeing->n; i++) { - psVectorAppend(fwhmValues, options->targetSeeing->data.F32[i]); + + bool overLimit = psMetadataLookupBool(NULL, readoutOut->analysis, "DECONV.OVERLIMIT"); + if (!overLimit) { + for (int i = 0; i < options->targetSeeing->n; i++) { + psVectorAppend(fwhmValues, options->targetSeeing->data.F32[i]); + } } psMetadataAddVector(readoutSrc->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues); Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotStackMatchPSFsNext.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotStackMatchPSFsNext.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotStackMatchPSFsNext.c (revision 34247) @@ -7,18 +7,23 @@ int nRadialEntries = 0; - // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, 0); // File of interest - psAssert (file, "missing file?"); + // find the numer of fwhmValues in the first non-skipped input + int num = psphotFileruleCount(config, filerule); + for (int i=0; ifiles, filerule, i); // File of interest + psAssert (file, "missing file?"); - pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); - psAssert (readout, "missing readout?"); - - bool status = false; - psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES"); - if (!fwhmValues) { - return 1; + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + bool status = false; + if (!psMetadataLookupBool (&status, readout->analysis, "PSPHOT.SKIP.INPUT")) { + psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES"); + if (fwhmValues) { + nRadialEntries = fwhmValues->n; + } else { + nRadialEntries = 1; + } + break; + } } - - nRadialEntries = fwhmValues->n; return nRadialEntries; } @@ -60,4 +65,9 @@ pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); psAssert (readout, "missing readout?"); + + if (psMetadataLookupBool (&status, readout->analysis, "PSPHOT.SKIP.INPUT")) { + psLogMsg("psphot", PS_LOG_INFO, "skipping smooth %d to next psf", index); + return true; + } psLogMsg("psphot", PS_LOG_INFO, "smooth %d to next psf", index); Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotStackMatchPSFsUtils.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotStackMatchPSFsUtils.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotStackMatchPSFsUtils.c (revision 34247) @@ -310,6 +310,14 @@ float deconv = psMetadataLookupF32(NULL, readoutOut->analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction if (deconv > deconvLimit) { +#if (1) + // XXX: don't reject the image set + psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image for matched psf analysis%d\n", deconv, deconvLimit, index); + psMetadataAddBool(readoutOut->analysis, PS_LIST_TAIL, "DECONV.OVERLIMIT", PS_META_REPLACE, "", true); +#else psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n", deconv, deconvLimit, index); goto escape; +#endif + } else { + psMetadataAddBool(readoutOut->analysis, PS_LIST_TAIL, "DECONV.OVERLIMIT", PS_META_REPLACE, "", false); } @@ -461,2 +469,29 @@ return optWidths; } + +// Set input to be skipped if the decovolution fraction was overlimit. Use for radial apertures +// This interface can be potentiall be extended for other uses +bool psphotStackSetInputsToSkip(pmConfig *config, const pmFPAview *view, const char *filerule, bool set) { + int num = psphotFileruleCount(config, filerule); + bool status; + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + for (int i = 0; i < num; i++) { + if (i == chisqNum) { + continue; + } + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + if (set) { + bool overLimit = psMetadataLookupBool(&status, readout->analysis, "DECONV.OVERLIMIT"); + psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSPHOT.SKIP.INPUT", PS_META_REPLACE, "Skip analysis", overLimit); + } else { + psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSPHOT.SKIP.INPUT", PS_META_REPLACE, "Skip analysis", false); + } + } + + return true; +} Index: /branches/eam_branches/ipp-20120627/psphot/src/psphotStackReadout.c =================================================================== --- /branches/eam_branches/ipp-20120627/psphot/src/psphotStackReadout.c (revision 34246) +++ /branches/eam_branches/ipp-20120627/psphot/src/psphotStackReadout.c (revision 34247) @@ -1,3 +1,5 @@ # include "psphotInternal.h" + +static bool psphotStackLoadWCS(pmConfig *config, const pmFPAview *view, const char *filerule); // we have 3 possible real filesets: @@ -67,4 +69,10 @@ char *STACK_SRC = useRaw ? STACK_RAW : STACK_CNV; char *STACK_DET = STACK_RAW; + + // load WCS + if (!psphotStackLoadWCS(config, view, STACK_SRC)) { + psError (PSPHOT_ERR_CONFIG, false, "trouble loading WCS for %s", STACK_SRC); + return false; + } // set the photcode for each image @@ -179,5 +187,5 @@ // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) - psphotFitSourcesLinear (config, view, STACK_SRC, false); + psphotFitSourcesLinear (config, view, STACK_SRC, false, false); psphotStackVisualFilerule(config, view, STACK_SRC); @@ -207,5 +215,5 @@ // NOTE : apply to ALL sources (extended + psf) // NOTE 2 : this function subtracts the models from the given filerule (SRC), not DET - psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 2 (detections->allSources) + psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 2 (detections->allSources) // NOTE: possibly re-measure background model here with objects subtracted / or masked @@ -278,5 +286,16 @@ } + // gcc doesn't like the label to refer to a declaration so we declare this here + bool splitLinearFit; + pass1finish: + + splitLinearFit = psMetadataLookupBool(NULL, recipe, "PSPHOT.STACK.SPLIT.LINEAR.FIT"); + if (splitLinearFit) { + psLogMsg ("psphot", 3, "splitting fit of detected and matched soures\n"); + // Fit the detected sources separately from matched that wea are about to create. + // NOTE: apply to ALL sources but only include sources with postitive flux in the fit + psphotFitSourcesLinear (config, view, STACK_SRC, true, true); // pass 3 (detections->allSources) + } // generate the objects (objects unify the sources from the different images) NOTE: could @@ -286,4 +305,8 @@ psMemDump("matchsources"); + // check the source density. If it too high change the number of radial bins + // in the recipe. + psphotLimitRadialApertures(recipe, objects->n); + // Construct an initial model for each object, set the radius to fitRadius, set circular // fit mask. NOTE: only applied to sources without guess models @@ -294,6 +317,12 @@ psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects); - // NOTE: apply to ALL sources - psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 3 (detections->allSources) + if (splitLinearFit) { + // NOTE: apply to Matched sources. Since the sources that we fit above are subtracted, they will + // not be included in this fit. + psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 4 (detections->allSources) + } else { + // Fit all sources together + psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 3 (detections->allSources) + } // measure the radial profiles to the sky (only measures new objects) @@ -344,9 +373,11 @@ } } - psphotRadialApertures (config, view, STACK_CNV, 0); // XXX entry 0 == unmatched? + psphotRadialApertures (config, view, STACK_CNV, 0); // entry 0 == unmatched psMemDump("extmeas"); + // mark any inputs that we want to skip the matched apertures for + psphotStackSetInputsToSkip(config, view, STACK_OUT, true); + int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT); - for (int entry = 1; entry < nRadialEntries; entry++) { // NOTE: entry 0 is the unmatched image set @@ -359,8 +390,7 @@ // this is necessary to get the right normalization for the new models - psphotFitSourcesLinear (config, view, STACK_OUT, false); + psphotFitSourcesLinear (config, view, STACK_OUT, false, false); // measure circular, radial apertures (objects sorted by S/N) - // entry 0 == unmatched? pass entry + 1? psphotRadialApertures (config, view, STACK_OUT, entry); @@ -373,4 +403,5 @@ } } + psphotStackSetInputsToSkip(config, view, STACK_OUT, false); // measure aperture photometry corrections @@ -413,4 +444,26 @@ // create the exported-metadata and free local data return psphotReadoutCleanup (config, view, STACK_SRC); +} + +static bool psphotStackLoadWCS(pmConfig *config, const pmFPAview *view, const char *filerule) { + + int num = psphotFileruleCount(config, filerule); + + for (int i=0; ifiles, filerule, i); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmHDU *hdu = pmHDUFromReadout(readout); + psAssert (hdu, "input missing hdu?"); + + if (!pmAstromReadWCS(file->fpa, readout->parent->parent, hdu->header, 1.0)) { + psError (PSPHOT_ERR_UNKNOWN, false, "failed to read WCS from header for input %d", i); + return false; + } + } + return true; }