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

Timestamp:

Feb 5, 2009, 5:03:33 PM (17 years ago)

Author:

Paul Price

Message:

Merging pap_branch_20090128. Resolved a small number of conflicts. Compiles, but not tested in detail.

File:

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

trunk/psphot/src/psphotSourceSize.c

@@ -2 +2 @@
 # include <gsl/gsl_sf_gamma.h>
 
-static float psphotModelContour(const psImage *image, const psImage *weight, const psImage *mask,
+static float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
                                 psImageMaskType maskVal, const pmModel *model, float Ro);
 
@@ -62 +62 @@
 
         psF32 **resid  = source->pixels->data.F32;
-        psF32 **weight = source->weight->data.F32;
+        psF32 **variance = source->variance->data.F32;
         psImageMaskType **mask    = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
 
         // check for extendedness: measure the delta flux significance at the 1 sigma contour
-        source->extNsigma = psphotModelContour(source->pixels, source->weight, source->maskObj, maskVal,
+        source->extNsigma = psphotModelContour(source->pixels, source->variance, source->maskObj, maskVal,
                                                source->modelPSF, 1.0);
 
@@ -103 +103 @@
         // Compare the central pixel with those on either side, for the four possible lines through it.
 
-        // Soften weights (add systematic error)
-        float softening = soft * PS_SQR(source->peak->flux); // Softening for weights
+        // Soften variances (add systematic error)
+        float softening = soft * PS_SQR(source->peak->flux); // Softening for variances
 
         // 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 dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1];
         float nX = cX / sqrtf(dcX + softening);
 
         // 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 dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0];
         float nY = cY / sqrtf(dcY + softening);
 
         // 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 dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1];
         float nL = cL / sqrtf(dcL + softening);
 
         // 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 dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1];
         float nR = cR / sqrtf(dcR + softening);
 
@@ -160 +160 @@
         // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
         if (source->crNsigma > CR_NSIGMA_LIMIT) {
-            // XXX still testing... : psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask);
-            psphotMaskCosmicRay_Old (source, maskVal, crMask);
+            // XXX still testing... : psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask);
+            psphotMaskCosmicRay_Old (source, maskVal, crMask);
         }
     }
@@ -190 +190 @@
 // deviation in sigmas.  This is measured on the residual image - should we ignore negative
 // deviations?
-static float psphotModelContour(const psImage *image, const psImage *weight, const psImage *mask,
+static float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
                                 psImageMaskType maskVal, const pmModel *model, float Ro)
 {
@@ -240 +240 @@
         if (yPixM >= 0 && yPixM < image->numRows &&
             !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) {
-            float dSigma = image->data.F32[yPixM][xPix] / sqrtf(weight->data.F32[yPixM][xPix]);
+            float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]);
             nSigma += dSigma;
             nPts++;
@@ -251 +251 @@
         if (yPixP >= 0 && yPixP < image->numRows &&
             !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) {
-            float dSigma = image->data.F32[yPixP][xPix] / sqrtf(weight->data.F32[yPixP][xPix]);
+            float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]);
             nSigma += dSigma;
             nPts++;
@@ -274 +274 @@
     pmFootprint *footprint = peak->footprint;
     if (!footprint) {
-        // if we have not footprint, use the old code to mask by isophot
-        psphotMaskCosmicRay_Old (source, maskVal, crMask);
-        return true;
+        // if we have not footprint, use the old code to mask by isophot
+        psphotMaskCosmicRay_Old (source, maskVal, crMask);
+        return true;
     }
 
     if (!footprint->spans) {
-        // if we have not footprint, use the old code to mask by isophot
-        psphotMaskCosmicRay_Old (source, maskVal, crMask);
-        return true;
+        // if we have not footprint, use the old code to mask by isophot
+        psphotMaskCosmicRay_Old (source, maskVal, crMask);
+        return true;
     }
 
     // mask all of the pixels covered by the spans of the footprint
     for (int j = 1; j < footprint->spans->n; j++) {
-        pmSpan *span1 = footprint->spans->data[j];
-
-        int iy = span1->y;
-        int xs = span1->x0;
-        int xe = span1->x1;
-        
-        for (int ix = xs; ix < xe; ix++) {
-            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
-        }
+        pmSpan *span1 = footprint->spans->data[j];
+
+        int iy = span1->y;
+        int xs = span1->x0;
+        int xe = span1->x1;
+
+        for (int ix = xs; ix < xe; ix++) {
+            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
+        }
     }
     return true;
@@ -308 +308 @@
     psImage *mask   = source->maskView;
     psImage *pixels = source->pixels;
-    psImage *weight = source->weight;
+    psImage *variance = source->variance;
 
     // XXX This should be a recipe variable
@@ -318 +318 @@
     // mark the pixels in this row to the left, then the right
     for (int ix = xo; ix >= 0; ix--) {
-        float SN = pixels->data.F32[yo][ix] / sqrt(weight->data.F32[yo][ix]);
-        if (SN > SN_LIMIT) {
-            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
-        }
+        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
+        if (SN > SN_LIMIT) {
+            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
+        }
     }
     for (int ix = xo + 1; ix < pixels->numCols; ix++) {
-        float SN = pixels->data.F32[yo][ix] / sqrt(weight->data.F32[yo][ix]);
-        if (SN > SN_LIMIT) {
-            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
-        }
+        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
+        if (SN > SN_LIMIT) {
+            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
+        }
     }
 
@@ -333 +333 @@
     // first go up:
     for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) {
-        // mark the pixels in this row to the left, then the right
-        for (int ix = 0; ix < pixels->numCols; ix++) {
-            float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]);
-            if (SN < SN_LIMIT) continue;
-
-            bool valid = false;
-            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask);
-            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0;
-            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0;
-
-            if (!valid) continue;
-            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
-        }
+        // mark the pixels in this row to the left, then the right
+        for (int ix = 0; ix < pixels->numCols; ix++) {
+            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
+            if (SN < SN_LIMIT) continue;
+
+            bool valid = false;
+            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask);
+            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0;
+            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0;
+
+            if (!valid) continue;
+            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
+        }
     }
     // next go down:
     for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) {
-        // mark the pixels in this row to the left, then the right
-        for (int ix = 0; ix < pixels->numCols; ix++) {
-            float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]);
-            if (SN < SN_LIMIT) continue;
-
-            bool valid = false;
-            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask);
-            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0;
-            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0;
-
-            if (!valid) continue;
-            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
-        }
+        // mark the pixels in this row to the left, then the right
+        for (int ix = 0; ix < pixels->numCols; ix++) {
+            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
+            if (SN < SN_LIMIT) continue;
+
+            bool valid = false;
+            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask);
+            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0;
+            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0;
+
+            if (!valid) continue;
+            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
+        }
     }
     return true;