Changeset 20144

Timestamp:

Oct 14, 2008, 10:14:51 AM (18 years ago)

Author:

Paul Price

Message:

Cleaning up, adding comments.

File:

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

trunk/psphot/src/psphotSourceSize.c

@@ -2 +2 @@
 # include <gsl/gsl_sf_gamma.h>
 
-float psphotModelContour (psImage *image, psImage *weight, psImage *mask, psMaskType maskVal, pmModel *model, float Ro);
+static float psphotModelContour(const psImage *image, const psImage *weight, const psImage *mask,
+                                psMaskType maskVal, const pmModel *model, float Ro);
 
 // we need to call this function after sources have been fitted to the PSF model and
@@ -10 +11 @@
 // deviation from the psf model at the r = FWHM/2 position
 
-bool psphotSourceSize (pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, long first)
+bool psphotSourceSize(pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, long first)
 {
-
     bool status;
 
@@ -41 +41 @@
 
         // skip source if it was already measured
-        if (isfinite(source->crNsigma)) continue;
+        if (isfinite(source->crNsigma)) {
+            psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since already measured\n");
+            continue;
+        }
 
         // source must have been subtracted
-        if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue;
+        if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) {
+            psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since source is not subtracted\n");
+            continue;
+        }
 
         psF32 **resid  = source->pixels->data.F32;
@@ -51 +57 @@
 
         // check for extendedness: measure the delta flux significance at the 1 sigma contour
-        source->extNsigma = psphotModelContour (source->pixels, source->weight, source->maskObj, maskVal, 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?
         if (source->extNsigma > EXT_NSIGMA_LIMIT) {
-          source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
-        }
-
+            source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
+        }
+
+        // Integer position of peak
         int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
         int yPeak = source->peak->yf - source->pixels->row0 + 0.5;
@@ -63 +71 @@
         // 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;
-        if (yPeak < 1) continue;
-        if (yPeak > source->pixels->numRows - 2) continue;
+        if (xPeak < 1 || xPeak > source->pixels->numCols - 2 ||
+            yPeak < 1 || yPeak > source->pixels->numRows - 2) {
+            psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
+            continue;
+        }
 
         // XXX for now, just skip any sources with masked pixels
@@ -73 +82 @@
         for (int iy = -1; (iy <= +1) && keep; iy++) {
             for (int ix = -1; (ix <= +1) && keep; ix++) {
-                if (mask[yPeak+iy][xPeak+ix]) { keep &= false; }
-            }
-        }
-        if (!keep) continue;
-
-        if (mask[yPeak][xPeak] & maskVal) continue;
-
-        // XXX need to deal with edge peaks... and mask
-        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);
-        }
+                if (mask[yPeak+iy][xPeak+ix] & maskVal) {
+                    keep = false;
+                }
+            }
+        }
+        if (!keep) {
+            psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n");
+            continue;
+        }
+
+        // Compare the central pixel with those on either side, for the four possible lines through it.
+
+        // Across the middle: y = 0
+        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];
+        float nX = cX / sqrtf(dcX);
+
+        // Up the centre: x = 0
+        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];
+        float nY = cY / sqrtf(dcY);
+
+        // Diagonal: x = y
+        float cL = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak-1]  - resid[yPeak+1][xPeak+1];
+        float dcL = 4*weight[yPeak][xPeak] + weight[yPeak-1][xPeak-1] + weight[yPeak+1][xPeak+1];
+        float nL = cL / sqrtf(dcL);
+
+        // Diagonal: x = - y
+        float cR = 2*resid[yPeak][xPeak]   - resid[yPeak+1][xPeak-1]  - resid[yPeak-1][xPeak+1];
+        float dcR = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak-1] + weight[yPeak-1][xPeak+1];
+        float nR = cR / sqrtf(dcR);
 
         // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2)
@@ -118 +122 @@
         // 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);
+        // 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;
@@ -156 +160 @@
             psImage *weight = source->weight;
 
+            // XXX This should be a recipe variable
             # define SN_LIMIT 5.0
 
@@ -232 +237 @@
 
 // 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
+// deviation in sigmas.  This is measured on the residual image - should we ignore negative
 // deviations?
-    float psphotModelContour (psImage *image, psImage *weight, psImage *mask, psMaskType maskVal, pmModel *model, float Ro) {
-
-    psF32 *PAR = model->params->data.F32;
-
-    // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY;
+static float psphotModelContour(const psImage *image, const psImage *weight, const psImage *mask,
+                                psMaskType maskVal, const pmModel *model, float Ro)
+{
+    psF32 *PAR = model->params->data.F32; // Model parameters
+    float sxx = PAR[PM_PAR_SXX], sxy = PAR[PM_PAR_SXY], syy = PAR[PM_PAR_SYY]; // Ellipse parameters
+
+    // We treat the contour as an ellipse:
+    // 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;
-
+    // This is a quadratic, Ay^2 + By + C with A = 1/SYY^2, B = x*SXY, C = (x / SXX)^2 - Ro
+    // The solution is y = [-B +/- sqrt (B^2 - 4 A C)] / [2 A], so:
+    // y = [-(x SXY) +/- sqrt ((x SXY)^2 - 4 (1/SYY^2) ((x/SXX)^2 - Ro))] * [SYY^2 / 2]
+
+    // min/max value of x is where B^2 - 4AC = 0; solve this for x
+    float Q = Ro * PS_SQR(sxx) / (1.0 - PS_SQR(sxx * syy * sxy) / 4.0);
+    if (Q < 0.0) {
+        // ellipse is imaginary
+        return NAN;
+    }
+
+    int radius = sqrtf(Q) + 0.5;        // Radius of ellipse
+    int nPts = 0;                       // Number of points in ellipse
+    float nSigma = 0.0;                 //
+
+    for (int x = -radius; x <= radius; x++) {
+        // Polynomial coefficients
+        float A = PS_SQR (1.0 / syy);
+        float B = x * sxy;
+        float C = PS_SQR (x / sxx) - Ro;
         float T = PS_SQR(B) - 4*A*C;
-        if (T < 0.0) continue;
-
+        if (T < 0.0) {
+            continue;
+        }
+
+        // y position in source frame
         float yP = (-B + sqrt (T)) / (2.0 * A);
         float yM = (-B - sqrt (T)) / (2.0 * A);
 
+        // Get the closest pixel positions (image frame)
         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] & maskVal)) {
-                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] & maskVal)) {
-                float dSigma = image->data.F32[yPixP][xPix] / sqrt (weight->data.F32[yPixP][xPix]);
-                nSigma += dSigma;
-                nPts ++;
-            }
+        if (xPix < 0 || xPix >= image->numCols) {
+            continue;
+        }
+
+        if (yPixM >= 0 && yPixM < image->numRows &&
+            !(mask && (mask->data.PS_TYPE_MASK_DATA[yPixM][xPix] & maskVal))) {
+            float dSigma = image->data.F32[yPixM][xPix] / sqrtf(weight->data.F32[yPixM][xPix]);
+            nSigma += dSigma;
+            nPts++;
+        }
+
+        if (yPixM == yPixP) {
+            continue;
+        }
+
+        if (yPixP >= 0 && yPixP < image->numRows &&
+            !(mask && (mask->data.PS_TYPE_MASK_DATA[yPixP][xPix] & maskVal))) {
+            float dSigma = image->data.F32[yPixP][xPix] / sqrtf(weight->data.F32[yPixP][xPix]);
+            nSigma += dSigma;
+            nPts++;
         }
     }