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

Timestamp:

Mar 30, 2010, 1:29:50 PM (16 years ago)

Author:

eugene

Message:

add diff stats; add flux to CMF_PS1_DV1; output PSF residual mask; more accurate calculation of radius of source model

File:

• trunk/psModules/src/objects/pmSourcePhotometry.c (modified) (6 diffs)

trunk/psModules/src/objects/pmSourcePhotometry.c

@@ -109 +109 @@
         psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y);
         psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y);
-        source->psfMag = -2.5*log10(fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]);
+        source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0];
+        source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0];
+        source->psfMag = -2.5*log10(source->psfFlux);
     } else {
-        status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
+        status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, source->modelPSF);
+        source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]);
     }
 
@@ -119 +122 @@
         for (int i = 0; i < source->modelFits->n; i++) {
             pmModel *model = source->modelFits->data[i];
-            status = pmSourcePhotometryModel (&model->mag, model);
+            status = pmSourcePhotometryModel (&model->mag, NULL, model);
             if (model == source->modelEXT) foundEXT = true;
         }
@@ -125 +128 @@
             source->extMag = source->modelEXT->mag;
         } else {
-            status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
+            status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
         }
     } else {
         if (source->modelEXT) {
-            status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
+            status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
         }
     }
@@ -143 +146 @@
     if (mode & PM_SOURCE_PHOT_WEIGHT) {
         pmSourcePixelWeight (&source->pixWeight, model, source->maskObj, maskVal);
+    }
+
+    // measure the contribution of included pixels
+    if (mode & PM_SOURCE_PHOT_DIFFSTATS) {
+        pmSourceMeasureDiffStats (source, maskVal);
     }
 
@@ -217 +225 @@
 
 // return source model magnitude
-bool pmSourcePhotometryModel (float *fitMag, pmModel *model)
-{
-    PS_ASSERT_PTR_NON_NULL(fitMag, false);
-    if (model == NULL) {
-        return false;
-    }
-
-    float fitSum = 0;
-    *fitMag = NAN;
+bool pmSourcePhotometryModel (float *fitMag, float *fitFlux, pmModel *model)
+{
+    psAssert (fitMag || fitFlux, "at least one of magnitude or flux must be requested (not NULL)");
+    if (model == NULL) return false;
+
+    float mag  = NAN;
+    float flux = NAN;
 
     // measure fitMag
-    fitSum = model->modelFlux (model->params);
-    if (fitSum <= 0)
-        return false;
-    if (!isfinite(fitSum))
-        return false;
-    *fitMag = -2.5*log10(fitSum);
+    flux = model->modelFlux (model->params);
+    if (flux > 0) {
+        mag = -2.5*log10(flux);
+    }
+    if (fitMag) {
+        *fitMag = mag;
+    }
+    if (fitFlux) {
+        *fitFlux = flux;
+    }
+
+    if (flux <= 0) return false;
+    if (!isfinite(flux)) return false;
 
     return (true);
@@ -356 +369 @@
 
     *pixWeight = validSum / modelSum;
+    return (true);
+}
+
+# define FLUX_LIMIT 3.0
+
+// return source aperture magnitude
+bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal)
+{
+    PS_ASSERT_PTR_NON_NULL(source, false);
+
+    if (source->diffStats == NULL) {
+        source->diffStats = pmSourceDiffStatsAlloc();
+    }
+
+    float fGood = 0.0;
+    float fBad  = 0.0;
+    int   nGood = 0;
+    int   nMask = 0;
+    int   nBad  = 0;
+    
+    psImage *flux     = source->pixels;
+    psImage *variance = source->variance;
+    psImage *mask     = source->maskObj;
+
+    for (int iy = 0; iy < flux->numRows; iy++) {
+        for (int ix = 0; ix < flux->numCols; ix++) {
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) {
+                nMask ++;
+                continue;
+            }
+
+            float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
+
+            if (SN > +FLUX_LIMIT) { 
+                nGood ++;
+                fGood += fabs(flux->data.F32[iy][ix]);
+            }
+
+            if (SN < -FLUX_LIMIT) { 
+                nBad ++;
+                fBad += fabs(flux->data.F32[iy][ix]);
+            }
+        }
+    }
+
+    source->diffStats->nGood      = nGood;
+    source->diffStats->fRatio     = (fGood + fBad         == 0.0) ? NAN : fGood / (fGood + fBad);          
+    source->diffStats->nRatioBad  = (nGood + nBad         == 0)   ? NAN : nGood / (float) (nGood + nBad);          
+    source->diffStats->nRatioMask = (nGood + nMask        == 0)   ? NAN : nGood / (float) (nGood + nMask);         
+    source->diffStats->nRatioAll  = (nGood + nMask + nBad == 0)   ? NAN : nGood / (float) (nGood + nMask + nBad);
+
     return (true);
 }