Changeset 17396 for trunk/psphot/src/psphotSourceSize.c

Timestamp:

Apr 8, 2008, 8:36:06 AM (18 years ago)

Author:

eugene

Message:

merging from eam_branch_20080324 : psphot work on extended source fitting and related I/O functions

File:

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

trunk/psphot/src/psphotSourceSize.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 # include <gsl/gsl_sf_gamma.h>
+
+float psphotModelContour (psImage *image, psImage *weight, psImage *mask, pmModel *model, float Ro);
 
 // we need to call this function after sources have been fitted to the PSF model and
@@ -14 +16 @@
     psTimerStart ("psphot");
 
-    float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CRNSIGMA.LIMIT");
+    float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.NSIGMA.LIMIT");
+    assert (status);
+
+    float EXT_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.LIMIT");
     assert (status);
 
@@ -24 +29 @@
         if (isfinite(source->crNsigma)) continue;
 
-        source->crNsigma  = -1.0;
-        source->extNsigma = -1.0;
-
         // source must have been subtracted
-        source->crNsigma  = -3.0;
         if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue;
 
-        psF32 **resid = source->pixels->data.F32;
+        psF32 **resid  = source->pixels->data.F32;
         psF32 **weight = source->weight->data.F32;
-        psU8 **mask = source->maskObj->data.U8;
+        psU8 **mask    = source->maskObj->data.U8;
+
+        // 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);
+
+        // XXX prevent a source from being both CR and EXT?
+        if (source->extNsigma > EXT_NSIGMA_LIMIT) {
+          source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
+        }
 
         int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
         int yPeak = source->peak->yf - source->pixels->row0 + 0.5;
 
-        // skip sources which are too close to a boundary
-        source->crNsigma  = -4.0;
+        // XXX for now, skip sources which are too close to a boundary
+        // XXX raise a flag?
         if (xPeak < 1) continue;
         if (xPeak > source->pixels->numCols - 2) continue;
@@ -46 +55 @@
 
         // XXX for now, just skip any sources with masked pixels
-        source->crNsigma  = -5.0;
+        // XXX raise a flag?
         bool keep = true;
         for (int iy = -1; (iy <= +1) && keep; iy++) {
@@ -113 +122 @@
         }
         source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
-        source->extNsigma = (nEXT > 0) ? fabs(fEXT) / nEXT : 0.0;
         // NOTE: abs needed to make the Nsigma value positive 
 
@@ -122 +130 @@
 
         if (source->crNsigma > CR_NSIGMA_LIMIT) {
-          source->mode |= PM_SOURCE_MODE_CRLIMIT;
+          source->mode |= PM_SOURCE_MODE_CR_LIMIT;
         }
     }
@@ -169 +177 @@
  */
 
+
+// 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) {
+
+    psF32 *PAR = model->params->data.F32;
+
+    // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY;
+    // y^2 (1/SYY^2) + y (x SXY) + (x / SXX)^2 - Ro = 0;
+    // y = [-(x SXY) +/- sqrt ((x SXY)^2 - 4 (1/SYY^2) ((x/SXX)^2 - Ro))] * [SYY^2 / 2];
+    // y = [-B +/- sqrt (B^2 - 4 A C)] / [2 A];
+
+    // min/max value of x is where T -> 0
+    // solve this for x2:
+    float Q = Ro * PS_SQR(PAR[PM_PAR_SXX]) / (1.0 - PS_SQR(PAR[PM_PAR_SXX]*PAR[PM_PAR_SYY]*PAR[PM_PAR_SXY]) / 4.0);
+    if (Q < 0.0) return NAN; // ellipse is imaginary
+
+    int xMax = sqrt(Q);
+    int xMin = -1.0*xMax;
+
+    int nPts = 0;
+    float nSigma = 0.0;
+
+    for (int x = xMin; x <= xMax; x++) {
+        float A = PS_SQR (1.0 / PAR[PM_PAR_SYY]);
+        float B = x * PAR[PM_PAR_SXY];
+        float C = PS_SQR (x / PAR[PM_PAR_SXX]) - Ro;
+
+        float T = PS_SQR(B) - 4*A*C;
+        if (T < 0.0) continue;
+    
+        float yP = (-B + sqrt (T)) / (2.0 * A);
+        float yM = (-B - sqrt (T)) / (2.0 * A);
+
+        int xPix  = x  + PAR[PM_PAR_XPOS] - image->col0 + 0.5;
+        int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
+        int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
+
+        if (xPix < 0) continue;
+        if (xPix >= image->numCols) continue;
+
+        if ((yPixM >= 0) && (yPixM < image->numRows)) {
+            if (!mask || !mask->data.U8[yPixM][xPix]) {
+                float dSigma = image->data.F32[yPixM][xPix] / sqrt (weight->data.F32[yPixM][xPix]);
+                nSigma += dSigma;
+                nPts ++;
+            }
+        }
+        
+        if (yPixM == yPixP) continue;
+
+        if ((yPixP >= 0) && (yPixP < image->numRows)) {
+            if (!mask || !mask->data.U8[yPixP][xPix]) {
+                float dSigma = image->data.F32[yPixP][xPix] / sqrt (weight->data.F32[yPixP][xPix]);
+                nSigma += dSigma;
+                nPts ++;
+            }
+        }
+    }   
+    nSigma /= nPts;
+    return nSigma;
+}