Changeset 31337

Timestamp:

Apr 21, 2011, 12:17:34 PM (15 years ago)

Author:

eugene

Message:

set seq for sources loaded from text; set the peak flux values for sources load from text; add masked kron flux measurements

Location:

branches/eam_branches/ipp-20110404/psphot/src

Files:

• psphotKronMasked.c (modified) (5 diffs)
• psphotLoadSRCTEXT.c (modified) (1 diff)
• psphotMergeSources.c (modified) (1 diff)
• psphotReadout.c (modified) (3 diffs)

branches/eam_branches/ipp-20110404/psphot/src/psphotKronMasked.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+
+bool psphotKronMag (pmSource *source, float radius, psImageMaskType maskVal);
 
 bool psphotKronMasked (pmConfig *config, const pmFPAview *view, const char *filerule)
@@ -29 +31 @@
         psAssert (detections, "missing detections?");
 
-        psArray *sources = detections->newSources;
+        psArray *sources = detections->allSources;
         psAssert (sources, "missing sources?");
 
@@ -58 +60 @@
     if (!status) {
         nThreads = 0;
+    }
+
+    float RADIUS = psMetadataLookupF32 (&status, readout->analysis, "PSF_MOMENTS_RADIUS");
+    if (!status) {
+        RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     }
 
@@ -84 +91 @@
         if (!source->moments) continue;
 
+        // replace object in image
+        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+        }
+
+        psphotKronMag (source, RADIUS, maskVal);
+
+        // re-subtract the object, leave local sky
+        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+
+        continue;
+
+        // XXX skip this code
         // generate the pixel masks
         // int Xo = source->moments->Mx;
@@ -117 +137 @@
     return true;
 }
+
+bool psphotKronMag (pmSource *source, float radius, 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(radius, 0.0, false);
+
+    psF32 Sum = 0.0;
+    psF32 Var = 0.0;
+    psF32 R2 = PS_SQR(radius);
+
+    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
+    // image coordinates.  the source->pixels image has an offset relative to its parent of
+    // col0,row0: a pixel (x,y) in the primary image has coordinates of (x-col0, y-row0) in
+    // this subimage.  we subtract off the peak coordinates, adjusted to this subimage, to have
+    // minimal round-off error in the sums.  since these values are subtracted just to minimize
+    // the dynamic range and are added back below, the exact value does not matter. these are
+    // (int) so they can be used in the image index below.
+
+    // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates
+    // Xn  = SUM (x - xc)^n * (z - sky)
+
+    psF32 RF = 0.0;
+    psF32 RS = 0.0;
+
+# if (1)
+    float dX = source->moments->Mx - source->peak->xf;
+    float dY = source->moments->My - source->peak->yf;
+    float dR = hypot(dX, dY);
+    float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
+    float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
+# else
+    float Xo = source->moments->Mx;
+    float Yo = source->moments->My;
+# endif
+
+    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
+    // xCM, yCM from pixel coords to pixel index here.
+    psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
+    psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
+
+    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+
+        psF32 yDiff = row - yCM;
+        if (fabs(yDiff) > radius) continue;
+
+        psF32 *vPix = source->pixels->data.F32[row];
+
+        psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
+
+        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++) {
+            if (vMsk) {
+                if (*vMsk & maskVal) {
+                    vMsk++;
+                    continue;
+                }
+                vMsk++;
+            }
+            if (isnan(*vPix)) continue;
+
+            psF32 xDiff = col - xCM;
+            if (fabs(xDiff) > radius) continue;
+
+            // radius is just a function of (xDiff, yDiff)
+            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+            if (r2 > R2) continue;
+
+            psF32 pDiff = *vPix;
+
+            Sum += pDiff;
+
+            // Kron Flux uses the 1st radial moment (NOT Gaussian windowed)
+            psF32 rf = pDiff * sqrt(r2);
+            psF32 rs = pDiff;
+
+            RF  += rf;
+            RS  += rs;
+        }
+    }
+
+    // Saturate the 1st radial moment
+    float Mrf = MIN(radius, MAX(0.25*radius, RF/RS));
+
+    // Calculate the Kron magnitude (make this block optional?)
+    float radKron  = 2.5*Mrf;
+    float radKron2 = radKron*radKron;
+
+    int nKronPix = 0;
+    Sum = Var = 0.0;
+
+    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+
+        psF32 yDiff = row - yCM;
+        if (fabs(yDiff) > radKron) continue;
+
+        psF32 *vPix = source->pixels->data.F32[row];
+        psF32 *vWgt = source->variance->data.F32[row];
+
+        psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
+
+        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
+            if (vMsk) {
+                if (*vMsk & maskVal) {
+                    vMsk++;
+                    continue;
+                }
+                vMsk++;
+            }
+            if (isnan(*vPix)) continue;
+
+            psF32 xDiff = col - xCM;
+            if (fabs(xDiff) > radKron) continue;
+
+            // radKron is just a function of (xDiff, yDiff)
+            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+            if (r2 > radKron2) continue;
+
+            psF32 pDiff = *vPix;
+            psF32 wDiff = *vWgt;
+
+            Sum += pDiff;
+            Var += wDiff;
+            nKronPix ++;
+        }
+    }
+    // XXX for a test, save the old values here:
+    source->moments->KronCore    = source->moments->KronFlux;
+    source->moments->KronCoreErr = source->moments->KronFluxErr;
+
+    source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
+    source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
+
+    return true;
+}

branches/eam_branches/ipp-20110404/psphot/src/psphotLoadSRCTEXT.c

@@ -52 +52 @@
 
             // NOTE: most of these values are irrelevant for loaded source positions
-            source->seq       = 0;
+            source->seq       = source->id;
             PAR[PM_PAR_XPOS]  = X;
             PAR[PM_PAR_YPOS]  = Y;

branches/eam_branches/ipp-20110404/psphot/src/psphotMergeSources.c

@@ -137 +137 @@
                 pmSource *source = extSourcesTXT->data[i];
                 source->mode |= PM_SOURCE_MODE_EXTERNAL;
+
+                // the supplied peak flux needs to be re-normalized
+                source->peak->rawFlux = 1.0;
+                source->peak->smoothFlux = 1.0;
+                source->peak->detValue = 1.0;
 
                 // drop the loaded source modelPSF

branches/eam_branches/ipp-20110404/psphot/src/psphotReadout.c

@@ -146 +146 @@
     }
 
-    // psphotKronMasked(config, view, filerule);
-
     // find blended neighbors of very saturated stars (detections->newSources)
     // if (!psphotDeblendSatstars (config, view, filerule)) {
@@ -201 +199 @@
     psphotFitSourcesLinear (config, view, filerule, false); // pass 1 (detections->allSources)
     psphotDumpChisqs (config, view, filerule);
+
+    // XXX re-measure the kron mags with models subtracted
+    psphotKronMasked(config, view, filerule);
 
     // identify CRs and extended sources (only unmeasured sources are measured)
@@ -311 +312 @@
 pass1finish:
 
+    // XXX re-measure the kron mags with models subtracted
+    psphotKronMasked(config, view, filerule);
+
     // measure source size for the remaining sources
     // NOTE: applies only to NEW (unmeasured) sources