Changeset 32612

Timestamp:

Nov 3, 2011, 7:41:16 PM (15 years ago)

Author:

eugene

Message:

fix up profile wings to use flux and slope as criteria for finish; save on source

File:

• branches/eam_branches/ipp-20110906/psphot/src/psphotRadialProfileWings.c (modified) (13 diffs)

branches/eam_branches/ipp-20110906/psphot/src/psphotRadialProfileWings.c

@@ -4 +4 @@
 // 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 +51 @@
 }
 
+// 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 +75 @@
 
     // 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 +109 @@
     }
 
+    // file = fopen ("radii.dat", "w");
+
     // threaded measurement of the source magnitudes
     // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
@@ -116 +130 @@
             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 +163 @@
     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 +174 @@
     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 +196 @@
 
         // 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 +206 @@
 }
 
+# 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 +235 @@
     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 +267 @@
         // \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 +287 @@
         // 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 +297 @@
     }
 
+    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?)