Changeset 31153 for trunk/psModules/src/objects/pmSourcePhotometry.c

Timestamp:

Apr 4, 2011, 1:04:41 PM (15 years ago)

Author:

eugene

Message:

updates to pmPeak to better distinguish peak flux versions; updates to visualization; add bits for substantial suspect masking; consolidate assignment of source position and flux based on peak, moments, etc; improve footprint culling process; fix PSF_QF and PSF_QF_PERFECT calculations; fix source model chisq values

Location:

trunk/psModules

Files:

• . (modified) (1 prop)
• src/objects/pmSourcePhotometry.c (modified) (11 diffs)

trunk/psModules

@@ -28 +28 @@
 ChangeLog
 psmodules-*.tar.*
+a.out.dSYM

@@ -28 +28 @@
 ChangeLog
 psmodules-*.tar.*
+a.out.dSYM

trunk/psModules/src/objects/pmSourcePhotometry.c

@@ -23 +23 @@
 #include "pmFPAMaskWeight.h"
 
+#include "pmConfigMask.h"
 #include "pmTrend2D.h"
 #include "pmResiduals.h"
@@ -49 +50 @@
 static float AP_MIN_SN = 0.0;
 
-bool pmSourceMagnitudesInit (psMetadata *config)
-{
-    PS_ASSERT_PTR_NON_NULL(config, false);
+// make this a bit more clever and dynamic
+static psImageMaskType maskSuspect  = 0;
+static psImageMaskType maskSpike    = 0;
+static psImageMaskType maskStarCore = 0;
+static psImageMaskType maskBurntool = 0;
+static psImageMaskType maskConvPoor = 0;
+
+bool pmSourceMagnitudesInit (pmConfig *config, psMetadata *recipe)
+{
+    PS_ASSERT_PTR_NON_NULL(recipe, false);
     bool status;
 
-    float limit = psMetadataLookupF32 (&status, config, "AP_MIN_SN");
+    // we are going to test specially against these poor values
+    if (config) {
+        maskSpike    = pmConfigMaskGet("SPIKE", config);
+        maskStarCore = pmConfigMaskGet("STARCORE", config);
+        maskBurntool = pmConfigMaskGet("BURNTOOL", config);
+        maskConvPoor = pmConfigMaskGet("CONV.POOR", config);
+        maskSuspect  = maskSpike | maskStarCore | maskBurntool | maskConvPoor;
+    }
+
+    float limit = psMetadataLookupF32 (&status, recipe, "AP_MIN_SN");
     if (status) {
         AP_MIN_SN = limit;
@@ -77 +94 @@
 // XXX masked region should be (optionally) elliptical
 // if mode is PM_SOURCE_PHOT_PSFONLY, we skip all other magnitudes
-bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal, psImageMaskType markVal)
+bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal, psImageMaskType markVal, float radius)
 {
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -166 +183 @@
     // measure the contribution of included pixels
     if (mode & PM_SOURCE_PHOT_WEIGHT) {
-        pmSourcePixelWeight (&source->pixWeightNotBad, &source->pixWeightNotPoor, model, source->maskObj, maskVal, markVal);
+        pmSourcePixelWeight (source, model, source->maskObj, maskVal, radius);
     }
 
@@ -342 +359 @@
 
 // return source aperture magnitude
-bool pmSourcePixelWeight (float *pixWeightNotBad, float *pixWeightNotPoor, pmModel *model, psImage *mask, psImageMaskType maskVal, psImageMaskType markVal)
-{
-    PS_ASSERT_PTR_NON_NULL(pixWeightNotBad, false);
-    PS_ASSERT_PTR_NON_NULL(pixWeightNotPoor, false);
+bool pmSourcePixelWeight (pmSource *source, pmModel *model, psImage *mask, psImageMaskType maskVal, float radius)
+{
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    source->pixWeightNotBad = NAN;
+    source->pixWeightNotPoor = NAN;
+
     PS_ASSERT_PTR_NON_NULL(mask, false);
     PS_ASSERT_PTR_NON_NULL(model, false);
@@ -355 +374 @@
     float value;
 
+    float spikeSum = 0;
+    float starcoreSum = 0;
+    float burntoolSum = 0;
+    float convpoorSum = 0;
+
     int Xo, Yo, dP;
     int dX, DX, NX;
     int dY, DY, NY;
 
-    *pixWeightNotBad = 0.0;
-    *pixWeightNotPoor = 0.0;
+    float radius2 = PS_SQR(radius);
 
     // we only care about the value of the object model, not the local sky
@@ -387 +410 @@
     NY = mask->numRows;
 
+    psImageMaskType maskBad = maskVal;
+    maskBad &= ~maskSuspect;
+
     // measure modelSum and validSum.  this function is applied to a sources' subimage.  the
     // value of DX is chosen (see above) to cover the full possible size of the subimage if it
     // were not by an edge; ie, if the source is cut in half by an image edge, we correctly
     // count the virtual pixels off the edge in normalizing the value of the pixWeight
+
+    // we skip any pixels [real or virtual] outside of the specified radius (nominally the aperture radius)
     for (int ix = -DX; ix < DX + 1; ix++) {
+        if (ix > radius) continue;
         int mx = ix + dX;
         for (int iy = -DY; iy < DY + 1; iy++) {
+            if (iy > radius) continue;
+            if (ix*ix + iy*iy > radius2) continue;
             int my = iy + dY;
 
@@ -409 +440 @@
             if (my >= NY) continue;
 
-            if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskVal)) {
+            // count pixels which are masked only with bad pixels
+            if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskBad)) {
                 notBadSum += value;
             }
-            if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & ~markVal)) {
+
+            // count pixels which are masked with an mask bit (bad or poor)
+            if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskVal)) {
                 notPoorSum += value;
             }
+
+            // count pixels which are masked with an mask bit (bad or poor)
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskSpike) {
+                spikeSum += value;
+            }
+            // count pixels which are masked with an mask bit (bad or poor)
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskStarCore) {
+                starcoreSum += value;
+            }
+            // count pixels which are masked with an mask bit (bad or poor)
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskBurntool) {
+                burntoolSum += value;
+            }
+            // count pixels which are masked with an mask bit (bad or poor)
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskConvPoor) {
+                convpoorSum += value;
+            }
         }
     }
     psFree (coord);
 
-    *pixWeightNotBad  = notBadSum  / modelSum;
-    *pixWeightNotPoor = notPoorSum / modelSum;
+    source->pixWeightNotBad  = notBadSum  / modelSum;
+    source->pixWeightNotPoor = notPoorSum / modelSum;
+
+    if ((spikeSum/modelSum) > 0.25) {
+        source->mode2 |= PM_SOURCE_MODE2_ON_SPIKE;
+    }
+    if ((starcoreSum/modelSum) > 0.25) {
+        source->mode2 |= PM_SOURCE_MODE2_ON_STARCORE;
+    }
+    if ((burntoolSum/modelSum) > 0.25) {
+        source->mode2 |= PM_SOURCE_MODE2_ON_BURNTOOL;
+    }
+    if ((convpoorSum/modelSum) > 0.25) {
+        source->mode2 |= PM_SOURCE_MODE2_ON_CONVPOOR;
+    }
+
+    if (isfinite(source->pixWeightNotBad) && isfinite(source->pixWeightNotPoor)) {
+        psAssert (source->pixWeightNotBad <= source->pixWeightNotPoor, "error: all bad pixels should also be poor");
+    }
 
     return (true);
@@ -427 +495 @@
 # define FLUX_LIMIT 3.0
 
-// measure stats that may be using in difference images for distinguishing real sources from bad residuals
+// measure stats that may be used in difference images for distinguishing real sources from bad residuals
 bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
 {
@@ -673 +741 @@
 # endif
 
-// determine chisq, etc for linear normalization-only fit
-bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal, const float covarFactor)
+// determine chisq, nPix, nDOF, chisqNorm : model->nPar must be set
+bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal)
 {
     PS_ASSERT_PTR_NON_NULL(model, false);
@@ -689 +757 @@
             if (variance->data.F32[j][i] <= 0)
                 continue;
-            dC += PS_SQR (image->data.F32[j][i]) / (covarFactor * variance->data.F32[j][i]);
+            dC += PS_SQR (image->data.F32[j][i]) / variance->data.F32[j][i];
             Npix ++;
         }
     }
     model->nPix = Npix;
-    model->nDOF = Npix - 1;
+    model->nDOF = Npix - model->nPar;
     model->chisq = dC;
+    model->chisqNorm = dC / model->nDOF;
 
     return (true);
 }
 
+
+// return source aperture magnitude
+bool pmSourceChisqUnsubtracted (pmSource *source, pmModel *model, psImageMaskType maskVal)
+{
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(model, false);
+
+    float dC = 0.0;
+    int Npix = 0;
+
+    // the model function returns the source flux at a position
+    psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
+
+    psVector *params = model->params;
+    psImage  *image = source->pixels;
+    psImage  *mask = source->maskObj;
+    psImage  *variance = source->variance;
+
+    int dX = image->col0;
+    int dY = image->row0;
+
+    for (int iy = 0; iy < image->numRows; iy++) {
+        for (int ix = 0; ix < image->numCols; ix++) {
+
+            // skip pixels which are masked
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue;
+
+            if (variance->data.F32[iy][ix] <= 0) continue;
+
+            coord->data.F32[0] = (psF32) ix + dX + 0.5;
+            coord->data.F32[1] = (psF32) iy + dY + 0.5;
+
+            // for the full model, add all points
+            float value = model->modelFunc (NULL, params, coord);
+
+            // fprintf (stderr, "%d, %d : %f, %f : %f - %f : %f\n", 
+            // ix, iy, coord->data.F32[0], coord->data.F32[1], image->data.F32[iy][ix], value, dC);
+
+            dC += PS_SQR (image->data.F32[iy][ix] - value) / variance->data.F32[iy][ix];
+            Npix ++;
+        }
+    }
+    model->nPix = Npix;
+    model->nDOF = Npix - model->nPar;
+    model->chisq = dC;
+    model->chisqNorm = dC / model->nDOF;
+
+    psFree (coord);
+    return (true);
+}
 
 double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)

@@ -28 +28 @@
 ChangeLog
 psmodules-*.tar.*
+a.out.dSYM