Index: /branches/eam_branches/ipp-20110906/psphot/src/psphotRadialProfileWings.c =================================================================== --- /branches/eam_branches/ipp-20110906/psphot/src/psphotRadialProfileWings.c (revision 32611) +++ /branches/eam_branches/ipp-20110906/psphot/src/psphotRadialProfileWings.c (revision 32612) @@ -4,5 +4,5 @@ // reach sky + X sigma -bool psphotRadialProfileWingsSource (pmSource *source, pmReadout *readout, float minRadius, float maxRadius, psImageMaskType maskVal); +bool psphotRadialProfileWingsSource (pmSource *source, pmReadout *readout, psImageMaskType maskVal); bool psphotRadialProfileFluxAtRadius (psVector *flux, psVector *fluxVar, pmSource *source, pmReadout *readout, float Radius, float dRadius, psImageMaskType maskVal); @@ -51,4 +51,10 @@ } +// these are set before we fork off to threads and used by all threads as constant values +static float MAX_RADIUS = NAN; +static float MIN_RADIUS = NAN; +static float SKY_STDEV = NAN; +// static FILE *file = NULL; + bool psphotRadialProfileWingsReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources) { @@ -69,12 +75,18 @@ // XXX is this a good recipe value to use for MAX RADIUS?? - float MAX_RADIUS = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_RADIUS"); + MAX_RADIUS = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_RADIUS"); if (!status) { MAX_RADIUS = 50.0; } - float MIN_RADIUS = 0.25*psMetadataLookupF32 (&status, readout->analysis, "PSF_MOMENTS_RADIUS"); + MIN_RADIUS = 0.25*psMetadataLookupF32 (&status, readout->analysis, "PSF_MOMENTS_RADIUS"); if (!status) { MIN_RADIUS = 0.25*psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); + } + + // SKY_STDEV is the sigma of the sky model (ie, smoothed on large scales) + SKY_STDEV = psMetadataLookupF32 (&status, readout->analysis, "MSKY_SIG"); + if (!status) { + SKY_STDEV = 1.0; // a crude default value (why would this not exist?) } @@ -97,4 +109,6 @@ } + // file = fopen ("radii.dat", "w"); + // threaded measurement of the source magnitudes // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) @@ -116,6 +130,4 @@ psArrayAdd(job->args, 1, cells->data[j]); // sources PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); - PS_ARRAY_ADD_SCALAR(job->args, MIN_RADIUS, PS_TYPE_F32); - PS_ARRAY_ADD_SCALAR(job->args, MAX_RADIUS, PS_TYPE_F32); // set this to 0 to run without threading @@ -151,4 +163,6 @@ psFree (cellGroups); + // fclose (file); + psLogMsg ("psphot.wings", PS_LOG_WARN, "measure radial profile wings : %f sec for %ld objects\n", psTimerMark ("psphot.wings"), sources->n); return true; @@ -160,6 +174,4 @@ psArray *sources = job->args->data[1]; psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[2],PS_TYPE_IMAGE_MASK_DATA); - float MIN_RADIUS = PS_SCALAR_VALUE(job->args->data[3],F32); - float MAX_RADIUS = PS_SCALAR_VALUE(job->args->data[4],F32); for (int i = 0; i < sources->n; i++) { @@ -184,5 +196,5 @@ // this function populates moments->Mrf,KronFlux,KronFluxErr - psphotRadialProfileWingsSource (source, readout, MIN_RADIUS, MAX_RADIUS, maskVal); + psphotRadialProfileWingsSource (source, readout, maskVal); // if we subtracted it above, re-subtract the object, leave local sky @@ -194,27 +206,28 @@ } +# define TEST_X 3158 +# define TEST_Y 3096 + // XXX use integer radius values? the rings assume integer values, right? or do they? -bool psphotRadialProfileWingsSource (pmSource *source, pmReadout *readout, float minRadius, float maxRadius, psImageMaskType maskVal) { +bool psphotRadialProfileWingsSource (pmSource *source, pmReadout *readout, psImageMaskType maskVal) { PS_ASSERT_PTR_NON_NULL(source, false); PS_ASSERT_PTR_NON_NULL(source->peak, false); PS_ASSERT_PTR_NON_NULL(source->pixels, false); - PS_ASSERT_FLOAT_LARGER_THAN(minRadius, 0.0, false); - PS_ASSERT_FLOAT_LARGER_THAN(maxRadius, 0.0, false); // psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; - // radii will be minRadius to maxRadius in NN log steps: + // radii will be MIN_RADIUS to MAX_RADIUS in NN log steps: float NSTEP = 10.0; float MIN_DR = 2; float NSIGMA = 1.0; - float SKY = 0.0; - float alpha = pow ((maxRadius / minRadius), 1.0/NSTEP) - 1.0; - float dRmax = maxRadius * alpha / (1.0 + alpha); // approximate size of last annulus, to get a rough size for vector allocation + float THRESHOLD = 2.0*SKY_STDEV; + float alpha = pow ((MAX_RADIUS / MIN_RADIUS), 1.0/NSTEP) - 1.0; + float dRmax = MAX_RADIUS * alpha / (1.0 + alpha); // approximate size of last annulus, to get a rough size for vector allocation int iter = 0; - psVector *flux = psVectorAllocEmpty(7*maxRadius*dRmax, PS_TYPE_F32); - psVector *fluxVar = psVectorAllocEmpty(7*maxRadius*dRmax, PS_TYPE_F32); + psVector *flux = psVectorAllocEmpty(7*MAX_RADIUS*dRmax, PS_TYPE_F32); + psVector *fluxVar = psVectorAllocEmpty(7*MAX_RADIUS*dRmax, PS_TYPE_F32); // should I just use sample median here? @@ -222,8 +235,19 @@ psStats *varStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); - for (float radius = minRadius; radius < maxRadius; iter ++) { + float lastFlux = NAN; + float lastRadius = NAN; + float lastSlope = NAN; + bool limit = false; + float limitRadius = NAN; + float limitFlux = NAN; + float limitSlope = NAN; + + // note: radius is the inner radius of the annulus; outer radius = inner radius + dR + for (float radius = MIN_RADIUS; !limit && (radius < MAX_RADIUS); iter ++) { float dR = (int)(radius * alpha); if (dR < MIN_DR) dR = MIN_DR; + float outerRadius = radius + dR; + float meanRadius = (2.0/3.0) * (outerRadius*outerRadius*outerRadius - radius*radius*radius) / (PS_SQR(outerRadius) - PS_SQR(radius)); // extract a vector of the pixel values (signal, variance) at this radius + dR @@ -243,10 +267,19 @@ // \mean(variance) * Npts and cancelling the Npts term inside and out of the sqrt() - fprintf (stderr, "%f %f : %f : %f %f\n", source->peak->xf, source->peak->yf, radius, meanFlux, meanFluxError); - - if (meanFlux - NSIGMA * meanFluxError < SKY) { - fprintf (stderr, "hit the sky at %f\n", radius); - } - + float slope = NAN; + if (isfinite(lastFlux)) { + slope = (meanFlux - lastFlux) / (meanRadius - lastRadius); + } + + // fprintf (stderr, "%f %f : %f : %f %f : %f\n", source->peak->xf, source->peak->yf, radius, meanFlux, meanFluxError, slope); + + if (!limit) { + limit |= (meanFlux - NSIGMA * meanFluxError < THRESHOLD); // dropped to sky level + limit |= isfinite(slope) && (fabs(slope) < 3.0); // SB no longer changing. + limitRadius = meanRadius; + limitFlux = meanFlux; + limitSlope = slope; + } + // completion criteria: // 1) flux - NSIGMA * dflux <= sky @@ -254,4 +287,8 @@ // 3) flux flat? + lastFlux = meanFlux; + lastRadius = meanRadius; + lastSlope = slope; + // reset the flux & fluxVar vector length to zero for re-use above: flux->n = 0; @@ -260,8 +297,19 @@ } + if (!limit) { + limitRadius = lastRadius; + limitFlux = lastFlux; + limitSlope = lastSlope; + } + // fprintf (file, "%f %f : %f %f : %f\n", source->peak->xf, source->peak->yf, limitRadius, limitFlux, limitSlope); + psFree (flux); psFree (fluxVar); psFree (fluxStats); psFree (varStats); + + source->skyRadius = limitRadius; + source->skyFlux = limitFlux; + source->skySlope = limitSlope; // save the max radius (and anything else?)