Index: /branches/eam_branches/ipp-20120405/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphot.h (revision 33952) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphot.h (revision 33953) @@ -189,5 +189,5 @@ // used by psphotFindDetections pmReadout *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal); -psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax); +psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax, int *totalPeaks); bool psphotFindFootprints (pmDetections *detections, pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal); psErrorCode psphotCullPeaks(const pmReadout *readout, const pmReadout *signifRO, const psMetadata *recipe, psArray *footprints); @@ -472,5 +472,5 @@ bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule); -bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf); +bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index); bool psphotKronIterate_Threaded (psThreadJob *job); Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotFindDetections.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotFindDetections.c (revision 33952) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotFindDetections.c (revision 33953) @@ -48,4 +48,5 @@ // Use the new pmFootprints approach? const bool useFootprints = psMetadataLookupBool(NULL, recipe, "USE_FOOTPRINTS"); + const bool footprintUseUnsubtracted = psMetadataLookupBool(NULL, recipe, "FOOTPRINT_USE_UNSUBTRACTED"); pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); @@ -64,4 +65,5 @@ } + bool replaceSourcesForFootprints = false; if (firstPass) { pass = 1; @@ -70,4 +72,5 @@ } else { pass = 2; + replaceSourcesForFootprints = footprintUseUnsubtracted; NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2"); PS_ASSERT (status, NULL); NMAX = 0; // unlimited number of peaks in final pass: allow a limit (PEAKS_NMAX_2) ? @@ -96,5 +99,6 @@ // detect the peaks in the significance image - detections->peaks = psphotFindPeaks (significance, readout, recipe, threshold, NMAX); + int totalPeaks = 0; + detections->peaks = psphotFindPeaks (significance, readout, recipe, threshold, NMAX, &totalPeaks); psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PEAK_THRESHOLD", PS_META_REPLACE, "Peak Detection Threshold", threshold); if (!detections->peaks) { @@ -105,8 +109,26 @@ return false; } + // hard limit on number of peaks we will accept. (To avoid memory overload in psphotStack) + int maxPeaks = psMetadataLookupS32 (&status, recipe, "PEAKS_NMAX_TOTAL"); PS_ASSERT (status, NULL); + if (maxPeaks && (totalPeaks > maxPeaks)) { + psFree (detections); + psError (PSPHOT_ERR_DATA, true, "Too many peaks %d found PEAKS_NMAX_TOTAL: %d", totalPeaks, maxPeaks); + return false; + } // optionally merge peaks into footprints if (useFootprints) { + if (replaceSourcesForFootprints) { + psphotAddOrSubNoiseReadout(config, view, filerule, index, recipe, false); + psphotReplaceAllSourcesReadout (config, view, filerule, index, recipe, false); + psFree (significance); + significance = psphotSignificanceImage (readout, recipe, maskVal); + } + psphotFindFootprints (detections, significance, readout, recipe, threshold, pass, maskVal); + + if (replaceSourcesForFootprints) { + psphotRemoveAllSourcesReadout (config, view, filerule, index, recipe, false); + } } Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotFindPeaks.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotFindPeaks.c (revision 33952) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotFindPeaks.c (revision 33953) @@ -4,5 +4,5 @@ // image must be constructed to represent (S/N)^2. If nMax is non-zero, only return a maximum // of nMax peaks -psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax) { +psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax, int *totalPeaks) { bool status = false; @@ -18,4 +18,8 @@ psError(PSPHOT_ERR_DATA, false, "no peaks found in this image"); return NULL; + } + // return the total number of peaks found before the nMax limit is applied + if (totalPeaks) { + *totalPeaks = peaks->n; } Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotKronIterate.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotKronIterate.c (revision 33952) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotKronIterate.c (revision 33953) @@ -1,9 +1,7 @@ # include "psphotInternal.h" -# ifndef ROUND -# define ROUND(X) ((int) ((X) + 0.5*SIGN(X))) -# endif - -bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert); -bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal); + +bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool oldWindow); +bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, bool applyWeight, psImage *smoothedPixels); + bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule) @@ -11,6 +9,4 @@ bool status = true; - // return true; - fprintf (stdout, "\n"); psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Kron Iterate ---"); @@ -46,5 +42,5 @@ // psAssert (psf, "missing psf?"); - if (!psphotKronIterateReadout (config, recipe, view, readout, sources, psf)) { + if (!psphotKronIterateReadout (config, recipe, view, filerule, readout, sources, psf, i)) { psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i); return false; @@ -58,5 +54,5 @@ bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max); -bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf) { +bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index) { bool status = false; @@ -69,4 +65,5 @@ psTimerStart ("psphot.kron"); + // determine the number of allowed threads int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads @@ -87,5 +84,27 @@ int KRON_ITERATIONS = psMetadataLookupS32 (&status, recipe, "KRON_ITERATIONS"); if (!status) { - KRON_ITERATIONS = 2; + KRON_ITERATIONS = 1; + } + + bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT"); + if (!status) { + KRON_APPLY_WEIGHT = true; + } + + bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW"); + if (!status) { + KRON_APPLY_WINDOW = false; + } + bool KRON_SMOOTH = psMetadataLookupBool (&status, recipe, "KRON_SMOOTH"); + if (!status) { + KRON_SMOOTH = false; + } + float KRON_SMOOTH_SIGMA = psMetadataLookupF32 (&status, recipe, "KRON_SMOOTH_SIGMA"); + if (!status) { + KRON_SMOOTH_SIGMA = 1.7; + } + float KRON_SMOOTH_NSIGMA = psMetadataLookupS32 (&status, recipe, "KRON_SMOOTH_NSIGMA"); + if (!status) { + KRON_SMOOTH_NSIGMA = 2; } @@ -119,5 +138,38 @@ // set a window function for each source based on the moments // (this skips really bad sources (no peak, no moments, DEFECT) - psphotKronWindowSetSource (source, kronWindow, false, true); + psphotKronWindowSetSource (source, kronWindow, false, true, KRON_APPLY_WINDOW); + } + + // We measure the Kron Radius on a smoothed copy of the readout image + psImage *smoothedImage = NULL; + if (KRON_SMOOTH) { + // Build the smoothed source image + // Replace the subtracted sources + psphotReplaceAllSourcesReadout(config, view, filerule, index, recipe, false); + // Copy the image and smooth + psTimerStart ("psphot.kron.smooth"); + smoothedImage = psImageCopy(NULL, readout->image, PS_TYPE_F32); + psImageSmooth(smoothedImage, KRON_SMOOTH_SIGMA, KRON_SMOOTH_NSIGMA); + psLogMsg ("psphot.kron", PS_LOG_INFO, "smoothed image %f sec\n", psTimerMark ("psphot.kron.smooth")); + + // remove the sources + psphotRemoveAllSourcesReadout( config, view, filerule, index, recipe, false ); + // Now subtract smooth versions of the sources from the smoothed image + psTimerStart ("psphot.kron.smooth.sources"); + for (int i=0; i< sources->n; i++) { + pmSource *source = sources->data[i]; + // If source has been subtracted from the readout image subtract a "smoothed" version from the smoothedImage + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + // cache copy of smoothedPixels in the source + // tmpPtr is for use by a single "module" and must be null otherwise + psAssert(source->tmpPtr == NULL, "source->tmpPtr is not null!"); + + psImage *smoothedPixels = psImageSubset(smoothedImage, source->region); + source->tmpPtr = (psPtr) smoothedPixels; + pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, false, maskVal, 0, 0); + } + } + psLogMsg ("psphot.kron", PS_LOG_INFO, "removed %ld smoothed sources %f sec\n", sources->n, psTimerMark ("psphot.kron.smooth.sources")); + } @@ -146,7 +198,11 @@ PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS, PS_TYPE_F32); PS_ARRAY_ADD_SCALAR(job->args, KRON_ITERATIONS, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WEIGHT, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WINDOW, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, KRON_SMOOTH_SIGMA, PS_TYPE_F32); + psArrayAdd(job->args, 1, smoothedImage); // set this to 0 to run without threading -# if (0) +# if (1) if (!psThreadJobAddPending(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); @@ -156,4 +212,5 @@ if (!psphotKronIterate_Threaded(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + // psFree(AnalysisRegion); return false; } @@ -179,4 +236,5 @@ psFree (cellGroups); psFree (kronWindow); + psFree (smoothedImage); psLogMsg ("psphot.kron", PS_LOG_WARN, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n); @@ -194,6 +252,10 @@ float MIN_KRON_RADIUS = PS_SCALAR_VALUE(job->args->data[6],F32); int KRON_ITERATIONS = PS_SCALAR_VALUE(job->args->data[7],S32); - - // XXX TEST : set iteration to 1 + bool KRON_APPLY_WEIGHT = PS_SCALAR_VALUE(job->args->data[8],S32); + bool KRON_APPLY_WINDOW = PS_SCALAR_VALUE(job->args->data[9],S32); + float KRON_SMOOTH_SIGMA = PS_SCALAR_VALUE(job->args->data[10],F32); + psImage *smoothedImage = job->args->data[11]; + + // psImage *smoothedPixels = NULL; for (int j = 0; j < KRON_ITERATIONS; j++) { for (int i = 0; i < sources->n; i++) { @@ -202,11 +264,21 @@ if (!source->peak) continue; // XXX how can we have a peak-less source? - // allocate space for moments + // check status of this source's moments if (!source->moments) continue; + if (!source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED) continue; + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; // replace object in image bool reSubtract = false; + psImage *smoothedPixels = NULL; if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { - pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + smoothedPixels = (psImage *) source->tmpPtr; + if (smoothedPixels) { + // psFree(source->tmpPtr); + // smoothedPixels = psImageSubset(smoothedImage, source->region); + pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, true, maskVal, 0, 0); + // source->tmpPtr = smoothedPixels; + } reSubtract = true; } @@ -223,47 +295,60 @@ // XXX float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), maxWindow); - // Sextractor apparently takes something like the skyRadius window, measures the moments inside that window, - // then uses an elliptical outline based on the 2nd moments, but 6x the size - // XXX TEST : use a window based on the radial profile numbers: max is skyRadius, min is RADIUS // if we lack the skyRadius (eg MATCHED sources), go to the default value - float maxWindow, windowRadius; - - if (j == 0) { - maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS; - windowRadius = PS_MAX(RADIUS, maxWindow); - } else { - maxWindow = isfinite(source->moments->Mrf) ? 6.0*source->moments->Mrf : RADIUS; - windowRadius = PS_MAX(RADIUS, maxWindow); - } - -#ifdef notdef - fprintf(stderr, "Redefining pixels for source: %d %4d %4d new radius: %f\n", - i, source->peak->x, source->peak->y, windowRadius+2); -#endif + float maxWindow; + if (j == 0) { + maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS; + } else { + maxWindow = isfinite(source->moments->Mrf) ? 6.0*source->moments->Mrf : RADIUS; + } + float windowRadius = PS_MAX(RADIUS, maxWindow); + // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS - pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2); + bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2); psAssert (source->pixels, "WTF?"); + if (extend && smoothedPixels) { + psFree(source->tmpPtr); + smoothedPixels = psImageSubset(smoothedImage, source->region); + psAssert (smoothedPixels, "WTF?"); + source->tmpPtr = (psPtr) smoothedPixels ; + } + // clear the window function for this source based on the moments - psphotKronWindowSetSource (source, kronWindow, (j > 0), false); - - // this function populates moments->Mrf,KronFlux,KronFluxErr - psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal); - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); - - // set a window function for each source based on the moments - psphotKronWindowSetSource (source, kronWindow, true, true); + // Note this function also applies cuts on the source and returns false if it + // does not meet the requirements for measuring the Kron Radius or Magnitude. + // Not that it performs that function even if KRON_APPLY_WINDOW is false + if (psphotKronWindowSetSource (source, kronWindow, (j > 0), false, KRON_APPLY_WINDOW)) { + + // this function populates moments->Mrf,KronFlux,KronFluxErr + psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal, KRON_APPLY_WEIGHT, smoothedPixels); + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + + // set a window function for each source based on the moments + psphotKronWindowSetSource (source, kronWindow, true, true, KRON_APPLY_WINDOW); + } // if we subtracted it above, re-subtract the object, leave local sky if (reSubtract) { - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + if (smoothedPixels) { + pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, false, maskVal, 0, 0); + if (j + 1 == KRON_ITERATIONS) { + // We're done with the smoothedPixels + psFree(source->tmpPtr); + source->tmpPtr = NULL; + smoothedPixels = NULL; + } + } } } } + // psFree(smoothedPixels); return true; } -bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal) { +bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, + bool applyWeight, psImage *smoothedPixels) { PS_ASSERT_PTR_NON_NULL(source, false); @@ -273,5 +358,5 @@ psF32 R2 = PS_SQR(radius); - float rsigma2 = 0.5 / PS_SQR(radius/2.0); + float rsigma2 = applyWeight ? 0.5 / R2 : 0; // a note about coordinates: coordinates of objects throughout psphot refer to the primary @@ -286,4 +371,5 @@ // Xn = SUM (x - xc)^n * (z - sky) + psF32 RF = 0.0; psF32 RS = 0.0; @@ -304,5 +390,10 @@ int Ywo = source->pixels->row0; - psF32 **vPix = source->pixels->data.F32; + psF32 **vPix; + if (smoothedPixels) { + vPix = smoothedPixels->data.F32; + } else { + vPix = source->pixels->data.F32; + } psF32 **vWin = kronWindow->data.F32; psF32 **vWgt = source->variance->data.F32; @@ -377,4 +468,9 @@ float Var = 0.0; float Win = 0.0; + + // set vPix to the source pixels (it may have been set to the + // smoothed image above) + vPix = source->pixels->data.F32; + for (psS32 row = 0; row < source->pixels->numRows ; row++) { @@ -415,5 +511,5 @@ } -bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert) { +bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool applyWindow) { if (!source) return false; @@ -422,5 +518,13 @@ if (source->type == PM_SOURCE_TYPE_DEFECT) return false; if (source->type == PM_SOURCE_TYPE_SATURATED) return false; + if (!source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED) return false; + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) return false; psAssert(kronWindow, "need a window"); + + if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false; + + if (!applyWindow) { + return true; + } // we have a source with moments Mx, My, Mxx, Mxy, Myy. we just need to define a Gaussian that has @@ -433,19 +537,29 @@ float Yo = source->moments->My; - float Mxx = source->moments->Mxx; - float Mxy = source->moments->Mxy; - float Myy = source->moments->Myy; - - float Mmajor = 0.5*(Mxx + Myy) + 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy)); - float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy)); - - // Mxx, Mxy, Myy define the elliptical shape, but Mrf defines the width - float scale = PS_SQR(0.5 * source->moments->Mrf) / Mmajor; - - float Sxx = scale * Mmajor * Mminor / Myy; // sigma_x^2 - float Sxy = Mxy / (scale * Mmajor * Mminor); - float Syy = scale * Mmajor * Mminor / Mxx; // sigma_y^2 - - float Smajor = sqrt(Mmajor); + psEllipseMoments moments; + moments.x2 = source->moments->Mxx; + moments.y2 = source->moments->Myy; + moments.xy = source->moments->Mxy; + + psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0); + if (! (isfinite(axes.major) && isfinite(axes.minor) && isfinite(axes.theta)) ) { + // Shall we log a proper warning? This happens often with matched sources (forced photometry) + // fprintf(dump, "invalid axes found id: %4d major: %f minor: %f theta: %f\n", source->id, axes.major, axes.minor, axes.theta); + return false; + } + + // Why this factor of 0.5 ? + float scale = 0.5 * source->moments->Mrf / axes.major; + axes.major *= scale; + axes.minor *= scale; + + psEllipseShape shape = psEllipseAxesToShape(axes); + if (! (isfinite(shape.sx) && isfinite(shape.sy) && isfinite(shape.sxy)) ) { + // Shall we log a proper warning? This happens often with matched sources (forced photometry) + // fprintf(dump, "invalid shape found id: %d sx: %f sy %f sxy: %f\n", source->id, shape.sx, shape.sy, shape.sxy); + return false; + } + + float Smajor = axes.major; int minX = PS_MIN(PS_MAX(Xo - 5*Smajor, 0), Nx - 1); @@ -454,7 +568,9 @@ int maxY = PS_MIN(PS_MAX(Yo + 5*Smajor, 0), Ny - 1); - float rMxx = 0.5 / Sxx; - float rMyy = 0.5 / Syy; - + float rMxx = 0.5 / PS_SQR(shape.sx); + float rMyy = 0.5 / PS_SQR(shape.sy); + float Sxy = -1. * shape.sxy; // factor of -1 is included to match the previous window function + // implementation. XXX: Is this correct? + for (int iy = minY; iy < maxY; iy++) { for (int ix = minX; ix < maxX; ix++) { @@ -467,7 +583,8 @@ float f = insert ? 1.001 - exp(-z) : 1.0 / (1.001 - exp(-z)); - kronWindow->data.F32[iy][ix] *= f; + kronWindow->data.F32[iy][ix] *= f; } } + return true; } Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotMergeSources.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotMergeSources.c (revision 33952) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotMergeSources.c (revision 33953) @@ -744,6 +744,13 @@ psAssert (readoutOut, "missing readout?"); + pmDetections *detectionsOutOld = psMetadataLookupPtr (&status, readoutOut->analysis, "PSPHOT.DETECTIONS"); + psArray *oldFootprints = detectionsOutOld ? detectionsOutOld->footprints : NULL; + // replace any existing DETECTION container on readoutOut->analysis with the new one pmDetections *detectionsOut = pmDetectionsAlloc(); + if (oldFootprints) { + // ... but hang on to any existing footprints so that they can be merged with new footprints in pass 2 + detectionsOut->footprints = psMemIncrRefCounter(oldFootprints); + } detectionsOut->allSources = psArrayAllocEmpty (100); if (!psMetadataAddPtr (readoutOut->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", detectionsOut)) { Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotRadialProfileWings.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotRadialProfileWings.c (revision 33952) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotRadialProfileWings.c (revision 33953) @@ -55,4 +55,5 @@ static float MIN_RADIUS = NAN; static float SKY_STDEV = NAN; +static float SKY_SLOPE_MIN = NAN; // static FILE *file = NULL; @@ -81,5 +82,5 @@ MIN_RADIUS = psMetadataLookupF32 (&status, readout->analysis, "PSF_MOMENTS_RADIUS"); - if (!status) { + if (!status || (MIN_RADIUS > MAX_RADIUS)) { MIN_RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); } @@ -89,4 +90,10 @@ if (!status) { SKY_STDEV = 1.0; // a crude default value (why would this not exist?) + } + + // SKY_SLOPE_MIN is the sigma of the sky model (ie, smoothed on large scales) + SKY_SLOPE_MIN = psMetadataLookupF32 (&status, recipe, "SKY_SLOPE_MIN"); + if (!status) { + SKY_SLOPE_MIN = 3.0; } @@ -302,10 +309,10 @@ limitSlope = slope; } - if (!limit && isfinite(slope) && (fabs(slope) < 3.0)) { + if (!limit && isfinite(slope) && (fabs(slope) < SKY_SLOPE_MIN)) { // SB no longer changing. limit = true; // linearly interpolate to the radius at which we hit the sky, using the last flux and the limiting slope if (isfinite(lastFlux)) { - float interpolatedRadius = lastRadius + lastFlux / 3.0; + float interpolatedRadius = lastRadius + lastFlux / SKY_SLOPE_MIN; if (interpolatedRadius < MAX_RADIUS) { limitRadius = interpolatedRadius; Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotReadout.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotReadout.c (revision 33952) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotReadout.c (revision 33953) @@ -185,5 +185,9 @@ // remove noise for subtracted objects (ie, return to normal noise level) // NOTE: this needs to operate only on the OLD sources - psphotSubNoise (config, view, filerule); // pass 1 (detections->allSources) + bool footprintUseUnsubtracted = psMetadataLookupBool(NULL, recipe, "FOOTPRINT_USE_UNSUBTRACTED"); + // Note: if footprintUseUnsubtracted is true the noise was already subtracted in psphotFindDetections() + if (!footprintUseUnsubtracted) { + psphotSubNoise (config, view, filerule); // pass 1 (detections->allSources) + } // define new sources based on only the new peaks & measure moments Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotReadoutMinimal.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotReadoutMinimal.c (revision 33952) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotReadoutMinimal.c (revision 33953) @@ -64,5 +64,8 @@ // Construct an initial model for each object, set the radius to fitRadius, set circular // fit mask. NOTE: only applied to sources without guess models - psphotGuessModels (config, view, filerule); + if (!psphotGuessModels (config, view, filerule)) { + psLogMsg ("psphot", 3, "failure to Guess Model"); + return psphotReadoutCleanup (config, view, filerule); + } // linear PSF fit to source peaks Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotSetThreads.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotSetThreads.c (revision 33952) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotSetThreads.c (revision 33953) @@ -30,5 +30,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 8); + task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 12); task->function = &psphotKronIterate_Threaded; psThreadTaskAdd(task); Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotSourceStats.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotSourceStats.c (revision 33952) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotSourceStats.c (revision 33953) @@ -605,4 +605,7 @@ } psAssert (isfinite(Sigma), "did we miss a case?"); + if (!isfinite(minKronRadius)) { + minKronRadius = Sigma; + } // choose a grid scale that is a fixed fraction of the psf sigma^2 Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotStackImageLoop.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotStackImageLoop.c (revision 33952) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotStackImageLoop.c (revision 33953) @@ -42,4 +42,5 @@ // XXX for now, just load the full set of images up front except for EXPNUM which we defer pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.RAW"); + pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.CNV"); if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot."); Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotStackMatchPSFsUtils.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotStackMatchPSFsUtils.c (revision 33952) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotStackMatchPSFsUtils.c (revision 33953) @@ -105,9 +105,12 @@ psAssert(psphotRecipe, "Need PSPHOT recipe"); +#if (0) + // DON'T OVERWRITE PSPHOT recipe value psString maskBadStr = psMetadataLookupStr(NULL, ppStackRecipe, "MASK.BAD");// Name of bits to mask for bad psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) - psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad); + psMetadataAddImageMask(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad); +#endif psImage *binned = psphotModelBackgroundReadoutNoFile(ro, config); // Binned background model Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotStackObjects.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotStackObjects.c (revision 33952) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotStackObjects.c (revision 33953) @@ -119,6 +119,4 @@ } if (!skipSourcePetro) keepObjectPetro = true; - - keepObjectPetro = true; } Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotStackParseCamera.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotStackParseCamera.c (revision 33952) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotStackParseCamera.c (revision 33953) @@ -129,4 +129,9 @@ } + psS64 stack_id = psMetadataLookupS64(&status, input, "STACK_ID"); + if (!status) { + psError(PS_ERR_UNKNOWN, false, "Unable to find STACK_ID from sources %d", i); + return false; + } // generate an pmFPAimage for the output convolved image // XXX output of these files should be optional @@ -139,5 +144,5 @@ } outputImage->save = true; - outputImage->fileID = i; // this is used to generate output names + outputImage->fileID = stack_id; // this is used to generate output names pmFPAfile *outputMask = pmFPAfileDefineOutput(config, outputImage->fpa, "PSPHOT.STACK.OUTPUT.MASK"); @@ -151,5 +156,5 @@ } outputMask->save = true; - outputMask->fileID = i; // this is used to generate output names + outputMask->fileID = stack_id; // this is used to generate output names pmFPAfile *outputVariance = pmFPAfileDefineOutput(config, outputImage->fpa, "PSPHOT.STACK.OUTPUT.VARIANCE"); @@ -163,5 +168,5 @@ } outputVariance->save = true; - outputVariance->fileID = i; // this is used to generate output names + outputVariance->fileID = stack_id; // this is used to generate output names // the output sources are carried on the outputImage->fpa structures @@ -172,5 +177,5 @@ } outsources->save = true; - outsources->fileID = i; // this is used to generate output names + outsources->fileID = stack_id; // this is used to generate output names } } Index: /branches/eam_branches/ipp-20120405/psphot/src/psphotStackReadout.c =================================================================== --- /branches/eam_branches/ipp-20120405/psphot/src/psphotStackReadout.c (revision 33952) +++ /branches/eam_branches/ipp-20120405/psphot/src/psphotStackReadout.c (revision 33953) @@ -89,6 +89,16 @@ return psphotReadoutCleanup (config, view, STACK_SRC); } - if (!psphotSubtractBackground (config, view, STACK_SRC)) { - return psphotReadoutCleanup (config, view, STACK_SRC); + if (strcmp(STACK_SRC, STACK_DET)) { +#define MODEL_BACKGROUND_SRC 1 +#ifdef MODEL_BACKGROUND_SRC + // work around the fact that the background levels on the convolved + // and unconvolved stacks can be different + if (!psphotModelBackground (config, view, STACK_SRC)) { + return psphotReadoutCleanup (config, view, STACK_SRC); + } +#endif + if (!psphotSubtractBackground (config, view, STACK_SRC)) { + return psphotReadoutCleanup (config, view, STACK_SRC); + } } if (!strcasecmp (breakPt, "BACKMDL")) { @@ -190,4 +200,5 @@ psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources) + // if we only do one pass, skip to extended source analysis if (!strcasecmp (breakPt, "PASS1")) goto pass1finish; @@ -211,5 +222,5 @@ // subtract all sources from DET (this will subtract using the psf model for SRC, which // will somewhat oversubtract the sources -- this is OK - psphotRemoveAllSources (config, view, STACK_DET, false); // ignore subtraction state for sources + psphotRemoveAllSources (config, view, STACK_DET, false); // do not ignore subtraction state for sources } @@ -223,5 +234,9 @@ // remove noise for subtracted objects (ie, return to normal noise level) // NOTE: this needs to operate only on the OLD sources - psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources) + // NOTE: if fooprintsUseUnsubtracted, the noise has already been removed by psphotFindDetections + bool footprintsUseUnsubtracted = psMetadataLookupBool(NULL, recipe, "FOOTPRINT_USE_UNSUBTRACTED"); + if (!footprintsUseUnsubtracted) { + psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources) + } // if DET and SRC are different images, copy the detections from DET to SRC @@ -303,6 +318,8 @@ // create source children for the OUT filerule (for radial aperture photometry) - psArray *objectsRadial = psphotSourceChildrenByObject (config, view, STACK_OUT, objects); - if (!objectsRadial) { + // (These are not just for radial aperture photomoetry. If they aren't defined we + // get no sources output + psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects); + if (!objectsOut) { psFree(objects); psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis"); @@ -310,33 +327,37 @@ } - // measure circular, radial apertures (objects sorted by S/N) - // this forces photometry on the undetected sources from other images - psphotRadialApertures (config, view, STACK_SRC, 0); // XXX entry 0 == unmatched? - psMemDump("extmeas"); - - int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT); - - for (int entry = 1; entry < nRadialEntries; entry++) { - // NOTE: entry 0 is the unmatched image set - - // re-measure the PSF for the smoothed image (using entries in 'allSources') - psphotChoosePSF (config, view, STACK_OUT, false); - - // this is necessary to update the models based on the new PSF - psphotResetModels (config, view, STACK_OUT); - - // this is necessary to get the right normalization for the new models - psphotFitSourcesLinear (config, view, STACK_OUT, false); - - // measure circular, radial apertures (objects sorted by S/N) - // entry 0 == unmatched? pass entry + 1? - psphotRadialApertures (config, view, STACK_OUT, entry); - - // replace the flux in the image so it is returned to its original state - psphotReplaceAllSources (config, view, STACK_OUT, false); - - // smooth to the next FWHM, or set 'smoothAgain' to false if no more - psphotStackMatchPSFsNext(config, view, STACK_OUT, entry); - psMemDump("matched"); + + bool radial_apertures = psMetadataLookupBool(NULL, recipe, "RADIAL_APERTURES"); + if (radial_apertures) { + // measure circular, radial apertures (objects sorted by S/N) + // this forces photometry on the undetected sources from other images + psphotRadialApertures (config, view, STACK_SRC, 0); // XXX entry 0 == unmatched? + psMemDump("extmeas"); + + int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT); + + for (int entry = 1; entry < nRadialEntries; entry++) { + // NOTE: entry 0 is the unmatched image set + + // re-measure the PSF for the smoothed image (using entries in 'allSources') + psphotChoosePSF (config, view, STACK_OUT, false); + + // this is necessary to update the models based on the new PSF + psphotResetModels (config, view, STACK_OUT); + + // this is necessary to get the right normalization for the new models + psphotFitSourcesLinear (config, view, STACK_OUT, false); + + // measure circular, radial apertures (objects sorted by S/N) + // entry 0 == unmatched? pass entry + 1? + psphotRadialApertures (config, view, STACK_OUT, entry); + + // replace the flux in the image so it is returned to its original state + psphotReplaceAllSources (config, view, STACK_OUT, false); + + // smooth to the next FWHM, or set 'smoothAgain' to false if no more + psphotStackMatchPSFsNext(config, view, STACK_OUT, entry); + psMemDump("matched"); + } } @@ -344,4 +365,5 @@ if (!psphotApResid (config, view, STACK_SRC)) { psFree (objects); + psFree (objectsOut); psLogMsg ("psphot", 3, "failed on psphotApResid"); return psphotReadoutCleanup (config, view, STACK_SRC); @@ -375,5 +397,5 @@ psFree (objects); - psFree (objectsRadial); + psFree (objectsOut); // create the exported-metadata and free local data