Changeset 34542

Timestamp:

Oct 24, 2012, 8:52:34 AM (14 years ago)

Author:

bills

Message:

In psphotStack add third pass for psphotKronIterate. Here measure kron flux
but not radius for matched sources. The radial momement for these is set to the
minimum of the values for the detected sources for the object

Location:

trunk/psphot/src

Files:

• psphot.h (modified) (1 diff)
• psphotKronIterate.c (modified) (12 diffs)
• psphotSourceMatch.c (modified) (2 diffs)
• psphotStackReadout.c (modified) (1 diff)

trunk/psphot/src/psphot.h

@@ -372 +372 @@
 bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, const char *filerule, int index);
 bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS);
-bool psphotDropBadMatchedSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
+bool psphotFilterMatchedSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
 
 bool psphotFitSourcesLinearStack (pmConfig *config, psArray *objects, bool final);

trunk/psphot/src/psphotKronIterate.c

@@ -2 +2 @@
 
 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 psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, bool applyWeight, psImage *smoothedPixels, bool measureRadius);
 
 
@@ -185 +185 @@
     }
 
+    if (pass == 3) {
+        KRON_SMOOTH = false;
+        KRON_ITERATIONS = 1;
+    }
+
     // We measure the Kron Radius on a smoothed copy of the readout image
     psImage *smoothedImage = NULL;
@@ -311 +316 @@
     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
+
+    bool measureRadius = true;
+    if (pass == 3) {
+        measureRadius = false;
+    }
 
     for (int j = 0; j < KRON_ITERATIONS; j++) {
@@ -320 +327 @@
             pmSource *source = sources->data[i];
             if (!source->peak) continue; // XXX how can we have a peak-less source?
+
+            if (pass == 3) {
+                // in pass 3 we only measure the flux for matched sources
+                if (!(source->mode2 & PM_SOURCE_MODE2_MATCHED)) {
+                    continue;
+                }
+            }
 
 # if (0)
@@ -331 +345 @@
 # endif
 
-            // 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;
+            if (measureRadius) {
+                // 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;
+            }
+
             if (!isfinite(source->moments->Mrf)) {
                 // Once we save a bad Mrf measurement we give up on this source
@@ -358 +375 @@
             // 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;
+            float maxWindow;
+            if (measureRadius) {
+                maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
+            } else {
+                maxWindow = source->moments->Mrf;
+            }
             if (j > 0) {
                 // on subsequent iterations we use a factor times the previous radial moment value
@@ -402 +424 @@
             // does not meet the requirements for measuring the Kron Radius or Magnitude.
             // Note: that it performs the cuts even if KRON_APPLY_WINDOW is false
-            if (psphotKronWindowSetSource (source, kronWindow, (j > 0), false, KRON_APPLY_WINDOW)) {
+            if (!measureRadius || 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);
+                psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal, 
+                    KRON_APPLY_WEIGHT, smoothedPixels, measureRadius);
+
                 psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
 
@@ -444 +468 @@
 
 bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal,
-    bool applyWeight, psImage *smoothedPixels) {
+    bool applyWeight, psImage *smoothedPixels, bool measureRadius) {
 
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -495 +519 @@
     psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
 
-    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
-
-        psF32 yDiff = row - yCM;
-        if (fabs(yDiff) > radius) continue;
-
-        // coordinate of mirror pixel
-        int yFlip = yCM - yDiff;
-        if (yFlip < 0) continue;
-        if (yFlip >= source->pixels->numRows) continue;
-
-        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
-            // check mask and value for this pixel
-            if (vMsk && (vMsk[row][col] & maskVal)) continue;
-            if (isnan(vPix[row][col])) continue;
-
-            psF32 xDiff = col - xCM;
-            if (fabs(xDiff) > radius) continue;
-
-            // coordinate of mirror pixel
-            int xFlip = xCM - xDiff;
-            if (xFlip < 0) continue;
-            if (xFlip >= source->pixels->numCols) continue;
-
-            // check mask and value for mirror pixel
-            if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
-            if (isnan(vPix[yFlip][xFlip])) continue;
-
-            // radius is just a function of (xDiff, yDiff)
-            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
-            if (r2 > R2) continue;
-
-            // flux * window
-            float z = r2 * rsigma2;
-            assert (z >= 0.0);
-
-            // weight by window image and wide Gaussian
-            float weight1  = vWin[row+Ywo][col+Xwo]*exp(-z);
-            float weight2  = vWin[yFlip+Ywo][xFlip+Xwo]*exp(-z);
-
-            float fDiff1 = vPix[row][col]*weight1;
-            float fDiff2 = vPix[yFlip][xFlip]*weight2;
-
-            float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
-
-            // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?)
-            psF32 rf = pDiff * sqrt(r2);
-            psF32 rs = 0.5 * (fDiff1 + fDiff2);
-
-            RF  += rf;
-            RS  += rs;
-        }
-    }
-
-    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
-    if (sqrt(source->peak->detValue) < 10.0) {
-        Mrf = MIN (radius, Mrf);
+    if (measureRadius) {
+        for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+
+            psF32 yDiff = row - yCM;
+            if (fabs(yDiff) > radius) continue;
+
+            // coordinate of mirror pixel
+            int yFlip = yCM - yDiff;
+            if (yFlip < 0) continue;
+            if (yFlip >= source->pixels->numRows) continue;
+
+            for (psS32 col = 0; col < source->pixels->numCols ; col++) {
+                // check mask and value for this pixel
+                if (vMsk && (vMsk[row][col] & maskVal)) continue;
+                if (isnan(vPix[row][col])) continue;
+
+                psF32 xDiff = col - xCM;
+                if (fabs(xDiff) > radius) continue;
+
+                // coordinate of mirror pixel
+                int xFlip = xCM - xDiff;
+                if (xFlip < 0) continue;
+                if (xFlip >= source->pixels->numCols) continue;
+
+                // check mask and value for mirror pixel
+                if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
+                if (isnan(vPix[yFlip][xFlip])) continue;
+
+                // radius is just a function of (xDiff, yDiff)
+                psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+                if (r2 > R2) continue;
+
+                // flux * window
+                float z = r2 * rsigma2;
+                assert (z >= 0.0);
+
+                // weight by window image and wide Gaussian
+                float weight1  = vWin[row+Ywo][col+Xwo]*exp(-z);
+                float weight2  = vWin[yFlip+Ywo][xFlip+Xwo]*exp(-z);
+
+                float fDiff1 = vPix[row][col]*weight1;
+                float fDiff2 = vPix[yFlip][xFlip]*weight2;
+
+                float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
+
+                // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?)
+                psF32 rf = pDiff * sqrt(r2);
+                psF32 rs = 0.5 * (fDiff1 + fDiff2);
+
+                RF  += rf;
+                RS  += rs;
+            }
+        }
+
+        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
+        if (sqrt(source->peak->detValue) < 10.0) {
+            Mrf = MIN (radius, Mrf);
+        }
+        source->moments->Mrf = Mrf;
     }
 
     // Calculate the Kron magnitude (make this block optional?)
-    float radKron  = 2.5*Mrf;
+    float radKron  = 2.5 * source->moments->Mrf;
     float radKron2 = radKron*radKron;
 
@@ -575 +602 @@
     // smoothed image above)
     vPix = source->pixels->data.F32;
-
 
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
@@ -607 +633 @@
     }
 
-    source->moments->Mrf         = Mrf;
     source->moments->KronFlux    = Sum;
     source->moments->KronFluxErr = sqrt(Var);
@@ -621 +646 @@
     if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
     if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
+    if (source->mode2 & PM_SOURCE_MODE2_MATCHED) return true;
+
     if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) return false;
     if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) return false;

trunk/psphot/src/psphotSourceMatch.c

@@ -559 +559 @@
 }
 
-bool psphotDropBadMatchedSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects) {
+bool psphotFilterMatchedSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects) {
 
     bool status = false;
 
-    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Drop Bad Matched Sources ---");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Filter Matched Sources ---");
 
     psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
@@ -615 +615 @@
     psFree(dropped);
 
+    // find the "best" Mrf from the detected sources. 
+    // Currently we use the smallest positive value
+    for (int i=0; i< objects->n; i++) {
+        pmPhotObj *obj = objects->data[i]; 
+
+        float minMrf = 1000.;
+        bool hasMatched = false;
+        for (int j = 0; j < obj->sources->n; j++) {
+            pmSource *source = obj->sources->data[j]; 
+            if (source->mode2 & PM_SOURCE_MODE2_MATCHED) {
+                hasMatched = true;
+                continue;
+            }
+            float Mrf = source->moments->Mrf;
+            if (isfinite(Mrf) && Mrf < minMrf && Mrf > 0) {
+                minMrf = Mrf;
+            }
+        }
+
+        if (!hasMatched || minMrf > 120.) {
+            continue;
+        }
+
+        // set Mrf for matched sources to the value found above
+        for (int j = 0; j < obj->sources->n; j++) {
+            pmSource *source = obj->sources->data[j]; 
+            if (source->mode2 & PM_SOURCE_MODE2_MATCHED) {
+                source->moments->Mrf = minMrf;
+            }
+        }
+    }
+
     return true;
 }

trunk/psphot/src/psphotStackReadout.c

@@ -369 +369 @@
     psMemDump("psfstats");
 
-    // drop matched sources without any useful measurements
-    psphotDropBadMatchedSources (config, view, STACK_SRC, objects);
+    // drop matched sources without any useful measurements and set kron radii for the ones
+    // we decide to keep
+    psphotFilterMatchedSources (config, view, STACK_SRC, objects);
+
+    // measure kron fluxes for the matched sources only
+    psphotKronIterate(config, view, STACK_SRC, 3);
 
     // measure elliptical apertures, petrosians (objects sorted by S/N)