Changeset 32205 for branches/eam_branches/ipp-20110710/psphot/src/psphotKronMasked.c

Timestamp:

Aug 29, 2011, 2:59:09 AM (15 years ago)

Author:

eugene

Message:

find the footprint saddles; iterative search for kron (in masked); add iterative search for kron (downweighting); add footprint based mask

File:

• branches/eam_branches/ipp-20110710/psphot/src/psphotKronMasked.c (modified) (13 diffs)

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

@@ -1 +1 @@
 # include "psphotInternal.h"
-
+# ifndef ROUND
+# define ROUND(X) ((int) ((X) + 0.5*SIGN(X)))
+# endif
+
+bool psphotKronMask (pmSource *source, psImage *kronMask, psImageMaskType markVal);
 bool psphotKronMag (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal);
 
@@ -45 +49 @@
 }
 
+int psphotKapaChannel (int channel);
+bool psphotVisualShowMask (int kapaFD, psImage *inImage, const char *name, int channel);
+bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max);
+
 bool psphotKronMaskedReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf) {
 
@@ -72 +80 @@
     }
 
+    float EXT_FIT_MAX_RADIUS = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_RADIUS");
+
     // bit-masks to test for good/bad pixels
     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
@@ -83 +93 @@
     maskVal |= markVal;
 
+    // psphotSaveImage (NULL, readout->mask, "kron.unmasked.fits");
+
+    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_PSFONLY;
+
+    // XXX tmp visualization
+    // int kapa = psphotKapaChannel (1);
+
     // generate the mask image: increment counter for every source overlapping the pixel
-    psImage *kronMask = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_U8);
+    psImage *kronMask = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_S32);
     psImageInit (kronMask, 0);
-    int Nx = kronMask->numCols;
-    int Ny = kronMask->numRows;
     for (int i = 0; i < sources->n; i++) {
 
@@ -96 +111 @@
         if (!source->moments) continue;
 
+        // replace object in image
+        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+        }
+
+        // save the PSF-based values
+        // source->moments->KronFluxPSF    = source->moments->KronFlux;
+        // source->moments->KronFluxPSFErr = source->moments->KronFluxErr;
+        // source->moments->KronRadiusPSF  = source->moments->Mrf;
+
+        // iterate to the window radius
         float windowRadius = RADIUS;
-        if (source->moments->KronFlux / source->moments->KronFluxErr > 10.0) {
-            windowRadius = PS_MAX (RADIUS, 10*source->moments->Mrf);
-        }
-
-        // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
-        pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
-
-        pmMoments psfMoments = *source->moments;
-
-        // replace object in image
-        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
-            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-        }
-
-        psphotKronMag (source, windowRadius, MIN_KRON_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;
-        // int Yo = source->moments->My;
-        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;
-        
-        int Kr = 2.5*source->moments->Mrf;
-        int Kr2 = Kr*Kr;
-        
-        for (int iy = Yo - Kr; iy < Yo + Kr + 1; iy++) {
-            if (iy < 0) continue;
-            if (iy >= Ny) continue;
-            for (int ix = Xo - Kr; ix < Xo + Kr + 1; ix++) {
-                if (ix < 0) continue;
-                if (ix >= Nx) continue;
-
-                if (PS_SQR(ix - Xo) + PS_SQR(iy - Yo) > Kr2) continue;
-                if (kronMask->data.U8[iy][ix] < 0xff) {
-                    kronMask->data.U8[iy][ix] ++;
-                }
-            }
-        }
-    }           
-
+        for (int j = 0; j < 4; j++) {
+            // XXX use some S/N criterion to limit?
+            // if (source->moments->KronFlux / source->moments->KronFluxErr > 10.0) {
+            // windowRadius = PS_MAX (RADIUS, 10*source->moments->Mrf);
+            // }
+
+            // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
+            pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
+
+            // mask the pixels not contained by the footprint
+            // if (psphotMaskFootprint (readout, source, markVal)) {
+            //  // psphotVisualShowMask (kapa, source->maskObj, "mask", 2);
+            //  // psphotVisualRangeImage (kapa, source->pixels, "image", 1, -50, 300);
+            //  // fprintf (stderr, "masked\n");
+            // }            
+
+            // mask the pixels associated with near neighbors
+            if (psphotKronMask (source, kronMask, markVal)) {
+                // psphotVisualShowMask (kapa, source->maskObj, "mask", 2);
+                // psphotVisualRangeImage (kapa, source->pixels, "image", 1, -50, 300);
+                // fprintf (stderr, "masked\n");
+            }       
+            // this function populates moments->Mrf,KronFlux,KronFluxErr
+            // XXX what about KronInner, KronOuter, etc?
+            psphotKronMag (source, windowRadius, MIN_KRON_RADIUS, maskVal);
+            windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), EXT_FIT_MAX_RADIUS);
+            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+        }
+        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
+        float kmag = -2.5*log10(source->moments->KronFlux);
+        if (source->psfMag - kmag > 0.25) {
+            // psphotVisualShowMask (kapa, source->maskObj, "mask", 1);
+            // psphotVisualRangeImage (kapa, source->pixels, "image", 0, -50, 300);
+            // fprintf (stderr, "psf: %f, kron: %f, dmag: %f\n", source->psfMag, kmag, source->psfMag - kmag);
+            // fprintf (stderr, "continue\n");
+        }
+        // re-subtract the object, leave local sky
+        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+    }
+
+    // psphotSaveImage (NULL, readout->mask, "kron.masked.fits");
     // psphotSaveImage (NULL, kronMask, "kronmask.fits");
+
     psLogMsg ("psphot.kron", PS_LOG_DETAIL, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
 
@@ -153 +171 @@
     return true;
 }
+
+# define WEIGHTED 0
 
 bool psphotKronMag (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal) {
@@ -161 +181 @@
     PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
 
-    psF32 Sum = 0.0;
-    psF32 Var = 0.0;
     psF32 R2 = PS_SQR(radius);
+
+# if (WEIGHTED)
+    float rsigma2 = 0.5 / PS_SQR(radius / 4.0); // use a Gaussian window with sigma = R_window / 2
+# endif
 
     // a note about coordinates: coordinates of objects throughout psphot refer to the primary
@@ -179 +201 @@
     psF32 RS = 0.0;
 
-# if (1)
+    // the peak position is less accurate but less subject to extreme deviations
     float dX = source->moments->Mx - source->peak->xf;
     float dY = source->moments->My - source->peak->yf;
@@ -185 +207 @@
     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
@@ -221 +239 @@
             if (r2 > R2) continue;
 
+# if (WEIGHTED)
+            float z = r2 * rsigma2;
+            assert (z >= 0.0);
+            float weight  = exp(-z);
+            psF32 pDiff = *vPix * weight;
+# else
             psF32 pDiff = *vPix;
-
-            Sum += pDiff;
-
-            // Kron Flux uses the 1st radial moment (NOT Gaussian windowed)
+# endif
+
+            // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?)
             psF32 rf = pDiff * sqrt(r2);
             psF32 rs = pDiff;
@@ -245 +268 @@
 
     int nKronPix = 0;
-    Sum = Var = 0.0;
+    float Sum = 0.0;
+    float Var = 0.0;
 
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
@@ -274 +298 @@
             if (r2 > radKron2) continue;
 
+# if (WEIGHTED)
+            float z = r2 * rsigma2;
+            assert (z >= 0.0);
+            float weight  = exp(-z);
+            psF32 pDiff = *vPix * weight;
+            psF32 wDiff = *vWgt * weight;
+# else
             psF32 pDiff = *vPix;
             psF32 wDiff = *vWgt;
+# endif
 
             Sum += pDiff;
@@ -283 +315 @@
     }
 
-    source->moments->Mrh = Mrf;
-
-    // XXX for a test, save the old values here:
-    source->moments->KronFouter    = 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;
+    source->moments->Mrf = Mrf;
+    source->moments->KronFlux    = Sum;
+    source->moments->KronFluxErr = sqrt(Var);
+    // source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
+    // source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
 
     return true;
 }
+
+enum {Y_P, Y_M, X_P, X_M};
+
+bool psphotKronMask (pmSource *source, psImage *kronMask, psImageMaskType markVal) {
+
+    int X, Y, e, dXs, dYs;
+
+    if (!source->peak) return false;
+    if (!source->peak->saddlePoints) return false;
+    if (!source->peak->saddlePoints->n) return false;
+
+    // return true;
+
+    int xMax = source->maskObj->numCols;
+    int yMax = source->maskObj->numRows;
+    // int xOff = source->maskObj->col0;
+    // int yOff = source->maskObj->row0;
+
+    int Xo = ROUND(source->peak->xf) - source->maskObj->col0;
+    int Yo = ROUND(source->peak->yf) - source->maskObj->row0;
+
+    for (int i = 0; i < source->peak->saddlePoints->n; i++) {
+
+        psVector *saddlePoint = source->peak->saddlePoints->data[i];
+        int Xs = ROUND(saddlePoint->data.S32[0]) - source->maskObj->col0;
+        int Ys = ROUND(saddlePoint->data.S32[1]) - source->maskObj->row0;
+        
+        // fprintf (stderr, "mask %d,%d @ %d,%d (%d,%d)\n", Xo, Yo, Xs, Ys, saddlePoint->data.S32[0], saddlePoint->data.S32[1]);
+
+        // We want to mask the pixels between the edge of the image and the line perpendicular
+        // to the connecting line.  Call dX,dY = (Xs-Xo,Ys-Yo).  There are 4 cases: 
+        // +y : dY > 0, fabs(dY) > fabs(dX)
+        // -y : dY < 0, fabs(dY) > fabs(dX)
+        // +x : dX > 0, fabs(dX) > fabs(dY)
+        // -x : dX < 0, fabs(dX) > fabs(dY)
+
+        int dX = Xs - Xo;
+        int dY = Ys - Yo;
+
+        // +y : dY > 0, fabs(dY) > fabs(dX)
+        // -y : dY < 0, fabs(dY) > fabs(dX)
+        if (fabs(dY) > fabs(dX)) {
+            // points to right of saddle
+            Y = Ys;
+            e = 0;
+            dXs = (dY > 0) ? +dY : -dY;  // this forces dXs > 0
+            dYs = (dY > 0) ? -dX : +dX;
+            for (X = Xs; X < xMax; X++) {
+                if (dY > 0) {
+                    for (int Ym = Y; Ym < yMax; Ym++) {
+                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
+                        // kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
+                    }
+                } else {
+                    for (int Ym = Y; Ym >= 0; Ym--) {
+                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
+                        // kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
+                    }
+                }
+                e += dYs;
+                int e2 = 2 * e;
+                if (e2 > dXs) {
+                    Y ++;
+                    e -= dXs;
+                }
+                if (e2 < -dXs) {
+                    Y --;
+                    e += dXs;
+                }
+                if (Y >= yMax) break;
+                if (Y < 0) break;
+            }
+            // points to left of saddle
+            Y = Ys;
+            e = 0;
+            for (X = Xs; X >= 0; X--) {
+                if (dY > 0) {
+                    for (int Ym = Y; Ym < yMax; Ym++) {
+                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
+                        // kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
+                    }
+                } else {
+                    for (int Ym = Y; Ym >= 0; Ym--) {
+                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
+                        // kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
+                    }
+                }
+                e -= dYs;
+                int e2 = 2 * e;
+                if (e2 > dXs) {
+                    Y ++;
+                    e -= dXs;
+                }
+                if (e2 < -dXs) {
+                    Y --;
+                    e += dXs;
+                }
+                if (Y >= yMax) break;
+                if (Y < 0) break;
+            }
+        } else {
+            // points to right of saddle
+            X = Xs;
+            e = 0;
+            dXs = (dX > 0) ? -dY : +dY;
+            dYs = (dX > 0) ? +dX : -dX;  // this forces dYs > 0
+            for (Y = Ys; Y < yMax; Y++) {
+                if (dX > 0) {
+                    for (int Xm = X; Xm < xMax; Xm++) {
+                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
+                        // kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
+                    }
+                } else {
+                    for (int Xm = X; Xm >= 0; Xm--) {
+                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
+                        // kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
+                    }
+                }
+                e += dXs;
+                int e2 = 2 * e;
+                if (e2 > dYs) {
+                    X ++;
+                    e -= dYs;
+                }
+                if (e2 < -dYs) {
+                    X --;
+                    e += dYs;
+                }
+                if (X >= yMax) break;
+                if (X < 0) break;
+            }
+            // points to left of saddle
+            X = Xs;
+            e = 0;
+            for (Y = Ys; Y >= 0; Y--) {
+                if (dX > 0) {
+                    for (int Xm = X; Xm < xMax; Xm++) {
+                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
+                        // kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
+                    }
+                } else {
+                    for (int Xm = X; Xm >= 0; Xm--) {
+                        source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
+                        // kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
+                    }
+                }
+                e -= dXs;
+                int e2 = 2 * e;
+                if (e2 > dYs) {
+                    X ++;
+                    e -= dYs;
+                }
+                if (e2 < -dYs) {
+                    X --;
+                    e += dYs;
+                }
+                if (X >= yMax) break;
+                if (X < 0) break;
+            }
+        }
+    }
+    return true;
+}