Changeset 34399 for branches/eam_branches/ipp-20120805/psphot/src/psphotKronIterate.c

Timestamp:

Sep 5, 2012, 4:04:34 PM (14 years ago)

Author:

eugene

Message:

merge changes from trunk

Location:

branches/eam_branches/ipp-20120805/psphot

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/psphotKronIterate.c (modified) (13 diffs)

branches/eam_branches/ipp-20120805/psphot/src/psphotKronIterate.c

@@ -5 +5 @@
 
 
-bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule)
+bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule, int pass)
 {
     bool status = true;
@@ -42 +42 @@
         // psAssert (psf, "missing psf?");
 
-        if (!psphotKronIterateReadout (config, recipe, view, filerule, readout, sources, psf, i)) {
+        if (!psphotKronIterateReadout (config, recipe, view, filerule, readout, sources, psf, i, pass)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
             return false;
@@ -54 +54 @@
 bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max);
 
-bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index) {
+#ifdef DUMP_KRS
+FILE *dumpFile = NULL;
+#endif
+
+bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index, int pass) {
 
     bool status = false;
+
+#ifdef DUMP_KRS
+    if (!dumpFile) {
+        dumpFile = fopen("kr.txt", "w");
+        psAssert (dumpFile, "failed to open kr.txt");
+    }
+    fprintf(dumpFile, "\n\n Input %d\n", index);
+#endif
 
     if (!sources->n) {
@@ -64 +76 @@
 
     psTimerStart ("psphot.kron");
-
 
     // determine the number of allowed threads
@@ -160 +171 @@
 
     // generate the window image: multiply the flux by this to downweight neighbors
+    // XXX: we don't need this image if we aren't going to apply the window
     psImage *kronWindow = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
     psImageInit (kronWindow, 1.0);
@@ -235 +247 @@
             PS_ARRAY_ADD_SCALAR(job->args, KRON_SMOOTH_SIGMA,  PS_TYPE_F32);
             PS_ARRAY_ADD_SCALAR(job->args, KRON_SB_MIN_DIVISOR,PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, pass,               PS_TYPE_S32);
 
 // set this to 0 to run without threading
@@ -297 +310 @@
     float KRON_SMOOTH_SIGMA         = PS_SCALAR_VALUE(job->args->data[11],F32);
     float KRON_SB_MIN_DIVISOR       = PS_SCALAR_VALUE(job->args->data[12],F32);
+    int pass                        = PS_SCALAR_VALUE(job->args->data[13],S32);
+#ifndef REVERT_ON_BAD_MEASUREMENT
+    (void) pass;
+#endif
 
     for (int j = 0; j < KRON_ITERATIONS; j++) {
@@ -308 +325 @@
             if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue;
             if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+            if (!isfinite(source->moments->Mrf)) {
+                // Once we save a bad Mrf measurement we give up on this source
+                // checking here allows us to avoid adding and subtracting the model
+                continue;
+            }
 
             // skip saturated stars modeled with a radial profile 
@@ -324 +346 @@
             }
 
-            // On first iteration set window radius to sky radius (if valid) on second iteration
-            // use a factor times the previous radial moment value up to a maximum value that
-            // depends on the surface brightness of the source
-            float maxWindow;
-            if (j == 0) {
-                maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
-            } else {
+            // On first iteration set window radius to sky radius (if valid). We also use this on subsequent
+            // iterations if we cannot find a better limit
+            float maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
+            if (j > 0) {
+                // on subsequent iterations we use a factor times the previous radial moment value
+                // limited to a maximum value that depends on the surface brightness of the source
                 if (KRON_SB_MIN_DIVISOR) {
-                    if (isfinite(source->moments->KronFlux) && (source->moments->KronFlux > 0)) {
-                        // Limit window radius based on surface brightness
+                    // Limit window radius based on surface brightness if we have a good measurement of kron flux
+                    if (isfinite(source->moments->KronFlux)  && (source->moments->KronFlux > 0)) {
                         float Rmax = sqrt(source->moments->KronFlux) / KRON_SB_MIN_DIVISOR;
 
-                        if (source->moments->Mrf > 0) {
+                        if (isfinite(source->moments->Mrf) && source->moments->Mrf > 0) {
                             maxWindow = PS_MIN(6.0*source->moments->Mrf, Rmax);
                         } else {
-                            maxWindow = Rmax;
+                            maxWindow = PS_MIN(Rmax, maxWindow);
                         }
-                    } else {
-                        maxWindow = RADIUS;
                     }
                 } else {
-                    // old code
+                    // old recipe, no surface brightness cut
                     maxWindow = isfinite(source->moments->Mrf) ? 6.0*source->moments->Mrf : RADIUS;
                 }
             }
             float windowRadius = PS_MAX(RADIUS, maxWindow);
+
+#ifdef REVERT_ON_BAD_MEASURMENT
+            // save previous measurements. We might revert back to them if this round fails
+            float MrfPrior = source->moments->Mrf;
+            float KronFluxPrior = source->moments->KronFlux;
+            float KronFluxErrPrior = source->moments->KronFluxErr;
+#endif
 
             // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
@@ -376 +402 @@
             }
 
+#ifdef REVERT_ON_BAD_MEASUREMENT
+            // on pass 2 if we get an invalid measurement on a pass 1 source where we had a good one previously
+            // in pass 1 keep that measurement
+            bool reverted = false;
+            if (pass > 1 && !isfinite(source->moments->Mrf) && (source->mode2 & PM_SOURCE_MODE2_PASS1_SRC)) {
+                source->moments->Mrf = MrfPrior;    // This is finite otherwise we wouldn't have gotten here
+                source->moments->KronFlux = KronFluxPrior;
+                source->moments->KronFluxErr = KronFluxErrPrior;
+                reverted = true;
+            }
+#endif
+#ifdef DUMP_KRS
+#ifndef REVERT_ON_BAD_MEASUREMENT
+            bool reverted = false;
+#endif
+            fprintf(dumpFile, "%7d %1d %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f %2d\n", source->id, reverted, source->moments->Mrf, MrfPrior, maxWindow, windowRadius, source->peak->xf, source->peak->yf, source->imageID);
+#endif
+
             // if we subtracted it above, re-subtract the object, leave local sky
             if (reSubtract) {
@@ -494 +538 @@
     }
 
+    float MrfTry = RF/RS;
+    if (RF <= 0. || RS <= 0 || !isfinite(MrfTry)) {
+        // We did not get a good measurement
+        source->moments->Mrf = NAN;
+        source->moments->KronFlux  = NAN;
+        source->moments->KronFluxErr  = NAN;
+        return false;
+    }
+
+    float Mrf = MAX(minKronRadius, MrfTry);
     // Saturate the 1st radial moment
-    float Mrf = MAX(minKronRadius, RF/RS);
     if (sqrt(source->peak->detValue) < 10.0) {
         Mrf = MIN (radius, Mrf);
@@ -544 +597 @@
     }
 
-    source->moments->Mrf = Mrf;
+    source->moments->Mrf         = Mrf;
     source->moments->KronFlux    = Sum;
     source->moments->KronFluxErr = sqrt(Var);
@@ -562 +615 @@
     psAssert(kronWindow, "need a window");
 
+    // XXX: If we are not applying the window then we don't need to check for valid Mrf here.
+    // We should give the this module a chance to measure a good value. 
+    // However experiments show that it hardly ever succeeds in getting a better value
     if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false;