Index: trunk/psphot/src/psphotSourceStats.c =================================================================== --- trunk/psphot/src/psphotSourceStats.c (revision 29936) +++ trunk/psphot/src/psphotSourceStats.c (revision 31154) @@ -122,4 +122,8 @@ // allocate space for moments source->moments = pmMomentsAlloc(); + + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) { + fprintf (stderr, "moment failure\n"); + } // allocate image, weight, mask arrays for each peak (square of radius OUTER) @@ -361,5 +365,4 @@ bool status = false; - float BIG_RADIUS; psScalar *scalar = NULL; @@ -373,4 +376,9 @@ 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); + + // if no valid pixels, massive swing or very large Mrf, likely saturated source, + // try a much larger box + float BIG_RADIUS = 3.0*RADIUS; + float BIG_SIGMA = 3.0*SIGMA; // maskVal is used to test for rejected pixels, and must include markVal @@ -390,5 +398,5 @@ // skip faint sources for moments measurement - if (source->peak->SN < MIN_SN) { + if (sqrt(source->peak->detValue) < MIN_SN) { source->mode |= PM_SOURCE_MODE_BELOW_MOMENTS_SN; Nfaint++; @@ -416,39 +424,40 @@ } - // measure basic source moments (no S/N clipping on input pixels) - status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, maskVal); - // XXX moments / aperture test: - if (0) { - // clear the mask bit and set the circular mask pixels - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, 4.0*source->moments->Mrf, "OR", markVal); - - float apMag = NAN; - pmSourcePhotometryAper (&apMag, NULL, NULL, NULL, source->pixels, source->variance, source->maskObj, maskVal); - fprintf (stderr, "apMag: %f, kronMag: %f\n", apMag, -2.5*log10(source->moments->KronFlux)); - - // clear the mask bit - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); - } - if (status) { - Nmoments ++; - continue; - } - - // if no valid pixels, or massive swing, likely saturated source, - // try a much larger box - BIG_RADIUS = PS_MIN (INNER, 3*RADIUS); - psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y); - status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, 0.0, maskVal); - if (status) { - source->mode |= PM_SOURCE_MODE_BIG_RADIUS; - Nmoments ++; - continue; - } - - source->mode |= PM_SOURCE_MODE_MOMENTS_FAILURE; - Nfail ++; - psErrorClear(); - continue; + // XXX how can this be set: this is raised in pmSourceRoughClass, called after this function? + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + fprintf (stderr, "satstar: %f,%f\n", source->peak->xf, source->peak->yf); + } + + if (!(source->peak->type == PM_PEAK_SUSPECT_SATURATION)) { + // measure basic source moments (no S/N clipping on input pixels) + status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, maskVal); + } else { + // For saturated stars, choose a much larger box NOTE this is slightly sleazy, but + // only slightly: pmSourceRedefinePixels uses the readout to pass the pointers to + // the parent image data. I guess the API could be simplified: we could recover + // this from the source in the function + + pmReadout tmpReadout; + tmpReadout.image = (psImage *)source->pixels->parent; + tmpReadout.mask = (psImage *)source->maskView->parent; + tmpReadout.variance = (psImage *)source->variance->parent; + + // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS + pmSourceRedefinePixels (source, &tmpReadout, source->peak->x, source->peak->y, BIG_RADIUS + 2); + + psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y); + status = pmSourceMoments (source, BIG_RADIUS, BIG_SIGMA, 0.0, maskVal); + source->mode |= PM_SOURCE_MODE_BIG_RADIUS; + } + if (!status) { + source->mode |= PM_SOURCE_MODE_MOMENTS_FAILURE; + Nfail ++; + psErrorClear(); + continue; + } + Nmoments ++; + + // re-try big sources or not?? + // if (source->moments->Mrf < 2.0*SIGMA) } @@ -473,6 +482,9 @@ float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); psAssert (status, "missing MOMENTS_SN_MIN"); - float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM"); psF32 MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX"); + + // when we set the window, we are not attempting to measure spatial variations; we can use a somewhat higher S/N limit + // since we are using all sources (true?) + float PSF_SN_LIM = 2.0*psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM"); // XXX move this to a config file? @@ -480,7 +492,8 @@ float sigma[NSIGMA] = {1.0, 2.0, 3.0, 4.5, 6.0, 9.0, 12.0, 18.0}; float Sout[NSIGMA]; - - // this sorts by peak->SN - sources = psArraySort (sources, pmSourceSortBySN); + int Nout[NSIGMA]; // number of stars found in clump : use this to control the number of regions measured by psphotRoughClass + + // this sorts by peak->rawFlux + sources = psArraySort (sources, pmSourceSortByFlux); // loop over radii: @@ -495,5 +508,5 @@ // skip faint sources for moments measurement - if (source->peak->SN < MIN_SN) { + if (sqrt(source->peak->detValue) < MIN_SN) { continue; } @@ -510,5 +523,5 @@ // determine the PSF parameters from the source moment values pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_AR_MAX); - psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X, Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]); + psLogMsg ("psphot", 3, "radius %.1f, nStars: %d of %d in clump, nSigma: %5.2f, X, Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nTotal, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]); #if 0 @@ -531,7 +544,9 @@ Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i]; + Nout[i] = psfClump.nStars; } // we are looking for sigma for which Sout = 0.65 (or some other value) + int Nstars = 0; float Sigma = NAN; float minS = Sout[0]; @@ -549,4 +564,5 @@ if ((Sout[i] < 0.65) && (Sout[i+1] < 0.65)) continue; Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]); + Nstars = 0.5*(Nout[i] + Nout[i+1]); } psAssert (isfinite(Sigma), "did we miss a case?"); @@ -558,4 +574,5 @@ psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma); psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma); + psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_CLUMP_NSTARS", PS_META_REPLACE, "number of stars in clump", Nstars); psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma);