Changeset 19915

Timestamp:

Oct 6, 2008, 3:16:44 AM (18 years ago)

Author:

eugene

Message:

apply the maskVal to mask tests; handle masked pixels in the outer 8 of top 9 pixels for CR; drop old EXT measurements (superceded by model contour); add visualizations

File:

• trunk/psphot/src/psphotSourceSize.c (modified) (8 diffs)

trunk/psphot/src/psphotSourceSize.c

@@ -2 +2 @@
 # include <gsl/gsl_sf_gamma.h>
 
-float psphotModelContour (psImage *image, psImage *weight, psImage *mask, pmModel *model, float Ro);
+float psphotModelContour (psImage *image, psImage *weight, psImage *mask, psMaskType maskVal, pmModel *model, float Ro);
 
 // we need to call this function after sources have been fitted to the PSF model and
@@ -51 +51 @@
 
         // check for extendedness: measure the delta flux significance at the 1 sigma contour
-        source->extNsigma = psphotModelContour (source->pixels, source->weight, source->maskObj, source->modelPSF, 1.0);
+        source->extNsigma = psphotModelContour (source->pixels, source->weight, source->maskObj, maskVal, source->modelPSF, 1.0);
 
         // XXX prevent a source from being both CR and EXT?
@@ -78 +78 @@
         if (!keep) continue;
 
+        if (mask[yPeak][xPeak] & maskVal) continue;
+
         // XXX need to deal with edge peaks... and mask
-        float cX = 2*resid[yPeak][xPeak] - resid[yPeak+0][xPeak-1] - resid[yPeak+0][xPeak+1];
-        float cY = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak+0] - resid[yPeak+1][xPeak+0];
-        float cL = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak-1] - resid[yPeak+1][xPeak+1];
-        float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1];
-
-        float dcX = 4*weight[yPeak][xPeak] + weight[yPeak+0][xPeak-1] + weight[yPeak+0][xPeak+1];
-        float dcY = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak+0] + weight[yPeak+1][xPeak+0];
-        float dcL = 4*weight[yPeak][xPeak] + weight[yPeak-1][xPeak-1] + weight[yPeak+1][xPeak+1];
-        float dcR = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak-1] + weight[yPeak-1][xPeak+1];
-
-        float nX = cX / sqrt(dcX);
-        float nY = cY / sqrt(dcY);
-        float nL = cL / sqrt(dcL);
-        float nR = cR / sqrt(dcR);
+        float nX = 0;
+        float nY = 0;
+        float nL = 0;
+        float nR = 0;
+
+        if (!(mask[yPeak+0][xPeak-1] & maskVal) && !(mask[yPeak+0][xPeak+1] & maskVal)) {
+            float cX = 2*resid[yPeak][xPeak] - resid[yPeak+0][xPeak-1] - resid[yPeak+0][xPeak+1];
+            float dcX = 4*weight[yPeak][xPeak] + weight[yPeak+0][xPeak-1] + weight[yPeak+0][xPeak+1];
+            nX = cX / sqrt(dcX);
+        }
+
+        if (!(mask[yPeak-1][xPeak+0] & maskVal) && !(mask[yPeak+1][xPeak+0] & maskVal)) {
+            float cY = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak+0] - resid[yPeak+1][xPeak+0];
+            float dcY = 4*weight[yPeak][xPeak] + weight[yPeak-1][xPeak+0] + weight[yPeak+1][xPeak+0];
+            nY = cY / sqrt(dcY);
+        }
+
+        if (!(mask[yPeak-1][xPeak-1] & maskVal) && !(mask[yPeak+1][xPeak+1] & maskVal)) {
+            float dcL = 4*weight[yPeak][xPeak] + weight[yPeak-1][xPeak-1] + weight[yPeak+1][xPeak+1];
+            float cL = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak-1] - resid[yPeak+1][xPeak+1];
+            nL = cL / sqrt(dcL);
+        }
+
+        if (!(mask[yPeak+1][xPeak-1] & maskVal) && !(mask[yPeak-1][xPeak+1] & maskVal)) {
+            float dcR = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak-1] + weight[yPeak-1][xPeak+1];
+            float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1];
+            nR = cR / sqrt(dcR);
+        }
 
         // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2)
@@ -100 +116 @@
         // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq);
 
-        // not strictly accurate: overcounts the chisq contribution from the center pixel (by factor of 4)
+        // not strictly accurate: overcounts the chisq contribution from the center pixel (by
+        // factor of 4); also biases a bit low if any pixels are masked
+        // XXX I am not sure I want to keep this value...
         source->psfChisq = PS_SQR (nX) + PS_SQR (nY) + PS_SQR (nX) + PS_SQR (nR);
 
         float fCR = 0.0;
-        float fEXT = 0.0;
         int nCR = 0;
-        int nEXT = 0;
         if (nX > 0.0) {
             fCR += nX;
             nCR ++;
-        } else {
-            fEXT += nX;
-            nEXT ++;
         }
         if (nY > 0.0) {
             fCR += nY;
             nCR ++;
-        } else {
-            fEXT += nY;
-            nEXT ++;
         }
         if (nL > 0.0) {
             fCR += nL;
             nCR ++;
-        } else {
-            fEXT += nL;
-            nEXT ++;
         }
         if (nR > 0.0) {
             fCR += nR;
             nCR ++;
-        } else {
-            fEXT += nR;
-            nEXT ++;
         }
         source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
-        // NOTE: abs needed to make the Nsigma value positive
-
-        if (!isfinite(source->crNsigma)) {
-          fprintf (stderr, ".");
-          source->crNsigma = -6.0;
-        }
+        assert (isfinite(source->crNsigma));
 
         // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
@@ -212 +211 @@
     }
 
+    // now that we have masked pixels associated with CRs, we can grow the mask
     if (grow > 0) {
         psImage *newMask = psImageConvolveMask(NULL, readout->mask, crMask, crMask, -grow, grow, -grow, grow);
@@ -225 +225 @@
               sources->n - first, psTimerMark ("psphot"));
 
+    psphotVisualPlotSourceSize (sources);
+    psphotVisualShowSourceSize (sources);
+
     return true;
 }
-
-
-// some possible classes:
-// all 4 curvatures are highly negative : CR
-// all 4 curvatures are highly positive : EXT
-
-// at least 2 are significantly negative, none are significantly positive : CR
-// at least 2 are significantly positive, none are significantly negative : EXT
-
-// any are significantly negative, some may be significantly positive : CR
-// any are significantly positive, none may be significantly positive : EXT
-
-/* Nn  Np  No
-   4   0   0   CR_1
-   3   1   0   CR_1
-   3   0   1   CR_1
-   2   2   0   CR_2
-   2   1   1   CR_2
-   2   0   2   CR_2
-   1   3   0   CR_3
-   1   2   1   CR_3
-   1   1   2   CR_3
-   1   0   3   CR_3
-   0   4   0   EXT
-   0   3   1   EXT
-   0   2   2   EXT
-   0   1   3   PSF
-   0   0   4   PSF
-*/
-
-/* Alternatively, write a f(CR) = Sum(nX,etc if >0) */
-
-
-/* I can write the formal probability that the 4 measurements are consistent with a PSF
- * based on how large the error values are (Nsigma -> erf -> P(Nsigma)).
- * I should examine this value as a function of flux and also examine the distribution of
- * the probabilities for the PSF sources.
- */
-
 
 // given the PSF ellipse parameters, navigate around the 1sigma contour, return the total
 // deviation in sigmas.  This is measure on the residual image - should we ignore negative
 // deviations?
-float psphotModelContour (psImage *image, psImage *weight, psImage *mask, pmModel *model, float Ro) {
+    float psphotModelContour (psImage *image, psImage *weight, psImage *mask, psMaskType maskVal, pmModel *model, float Ro) {
 
     psF32 *PAR = model->params->data.F32;
@@ -309 +273 @@
 
         if ((yPixM >= 0) && (yPixM < image->numRows)) {
-            if (!mask || !mask->data.U8[yPixM][xPix]) {
+            if (!mask || !(mask->data.U8[yPixM][xPix] & maskVal)) {
                 float dSigma = image->data.F32[yPixM][xPix] / sqrt (weight->data.F32[yPixM][xPix]);
                 nSigma += dSigma;
@@ -319 +283 @@
 
         if ((yPixP >= 0) && (yPixP < image->numRows)) {
-            if (!mask || !mask->data.U8[yPixP][xPix]) {
+            if (!mask || !(mask->data.U8[yPixP][xPix] & maskVal)) {
                 float dSigma = image->data.F32[yPixP][xPix] / sqrt (weight->data.F32[yPixP][xPix]);
                 nSigma += dSigma;