Changeset 32236

Timestamp:

Aug 31, 2011, 3:23:36 AM (15 years ago)

Author:

eugene

Message:

cleanup test code, allow for usage with new or allSources, only subtract source if we added it in, do not symmetrize the mag, just the radius

File:

• branches/eam_branches/ipp-20110710/psphot/src/psphotKronIterate.c (modified) (10 diffs)

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

@@ -37 +37 @@
         psAssert (detections, "missing detections?");
 
-        psArray *sources = detections->allSources;
+        psArray *sources = detections->newSources ? detections->newSources : detections->allSources;
         psAssert (sources, "missing sources?");
 
@@ -127 +127 @@
     // psphotVisualRangeImage (kapa, kronWindow, "kronwin", 0, 0.0, 1.0);
 
-    psphotSaveImage (NULL, kronWindow, "kron.window.v0.fits");
+    // psphotSaveImage (NULL, kronWindow, "kron.window.v0.fits");
 
     for (int j = 0; j < 5; j++) {
@@ -139 +139 @@
 
             // replace object in image
+            bool reSubtract = false;
             if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
                 pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+                reSubtract = true;
             }
 
@@ -151 +153 @@
             // clear the window function for this source based on the moments
             psphotKronWindowSetSource (source, kronWindow, (j > 0), false);
-            // psphotKronWindowSetSource (source, kronWindow, false, false);
             // psphotVisualRangeImage (kapa, kronWindow, "kronwin", 1, 0.0, 1.0);
-
-            // 165, 539;
-            if ((fabs(source->peak->xf - 165) < 3) && (fabs(source->peak->yf - 539) < 3)) {
-                fprintf (stderr, "test obj\n");
-            }
 
             // this function populates moments->Mrf,KronFlux,KronFluxErr
@@ -165 +161 @@
             // set a window function for each source based on the moments
             psphotKronWindowSetSource (source, kronWindow, true, true);
-            // psphotKronWindowSetSource (source, kronWindow, false, true);
-
-            // test source fluxes
-            pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
-            float kmag = -2.5*log10(source->moments->KronFlux);
-# define TEST_X1 167
-# define TEST_Y1 299
-# define TEST_X2 180
-# define TEST_Y2 300
-            if ((fabs(source->peak->xf - TEST_X1) < 3) && (fabs(source->peak->yf - TEST_Y1) < 3)) {
-                fprintf (stderr, "R1: %f vs %f  (%f) (%f)\n", source->moments->KronRadiusPSF, source->moments->Mrf, kmag, windowRadius);
+
+            // if we subtracted it above, re-subtract the object, leave local sky
+            if (reSubtract) {
+                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
             }
-            if ((fabs(source->peak->xf - TEST_X2) < 3) && (fabs(source->peak->yf - TEST_Y2) < 3)) {
-                fprintf (stderr, "R2: %f vs %f  (%f) (%f)\n", source->moments->KronRadiusPSF, source->moments->Mrf, kmag, windowRadius);
-            }
-
-            // re-subtract the object, leave local sky
-            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
         }
-        char name[64];
-        sprintf (name, "kron.window.v%d.fits", j+1);
-        psphotSaveImage (NULL, kronWindow, name);
+        // char name[64];
+        // sprintf (name, "kron.window.v%d.fits", j+1);
+        // psphotSaveImage (NULL, kronWindow, name);
     }
     psFree (kronWindow);
@@ -275 +258 @@
             assert (z >= 0.0);
 
-            // float weight  = 1.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 weight1  = vWin[row+Ywo][col+Xwo];
-            // float weight2  = vWin[yFlip+Ywo][xFlip+Xwo];
 
             float fDiff1 = vPix[row][col]*weight1;
@@ -315 +296 @@
         if (fabs(yDiff) > radKron) 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
@@ -328 +304 @@
             if (fabs(xDiff) > radKron) 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;
-
             // radKron is just a function of (xDiff, yDiff)
             psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
             if (r2 > radKron2) continue;
 
-            // float z = r2 * rsigma2;
-            // assert (z >= 0.0);
-
-            // float weight  = 1.0;
-            // float weight1  = vWin[row+Ywo][col+Xwo]*exp(-z);
-            // float weight2  = vWin[yFlip+Ywo][xFlip+Xwo]*exp(-z);
             float weight1  = vWin[row+Ywo][col+Xwo];
-            float weight2  = vWin[yFlip+Ywo][xFlip+Xwo];
-
             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));
+
+            float pDiff = fDiff1;
             psF32 wDiff = vWgt[row][col] * weight1;
 
@@ -365 +323 @@
     source->moments->Mrf = Mrf;
     source->moments->KronFlux    = Sum;
-    source->moments->KronFluxErr = Sum * (nKronPix / Win);
-    // source->moments->KronFluxErr = sqrt(Var / Win);
+    source->moments->KronFluxErr = sqrt(Var);
 
     return true;
@@ -393 +350 @@
     float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
 
-    // float kratio = source->moments->KronFinner / source->moments->KronFlux;
-
+    // Mxx, Mxy, Myy define the elliptical shape, but Mrf defines the width 
     float scale = PS_SQR(0.5 * source->moments->Mrf) / Mmajor;
-    // float scale = useKronRadius ? 2.0 * source->moments->Mrf / Mmajor : 2.0;
-    // float scale = (kratio > 0.4) ? 9.0 * source->moments->Mrf / Mmajor : 3.0 * source->moments->Mrf / Mmajor;
 
     float Sxx = scale * Mmajor * Mminor / Myy; // sigma_x^2