Changeset 27838 for branches/tap_branches/psModules/src/objects/pmSourcePhotometry.c

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

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

branches/tap_branches/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);
 }
@@ -545 +609 @@
 
 // determine chisq, etc for linear normalization-only fit
-bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance,
-                    psImageMaskType maskVal)
+bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal, const float covarFactor)
 {
     PS_ASSERT_PTR_NON_NULL(model, false);
@@ -561 +624 @@
             if (variance->data.F32[j][i] <= 0)
                 continue;
-            dC += PS_SQR (image->data.F32[j][i]) / variance->data.F32[j][i];
+            dC += PS_SQR (image->data.F32[j][i]) / (covarFactor * variance->data.F32[j][i]);
             Npix ++;
         }
@@ -573 +636 @@
 
 
-double pmSourceModelWeight(const pmSource *Mi,
-                      int term,
-                      const bool unweighted_sum) // should the cross product be weighted?
+double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor)
 {
     PS_ASSERT_PTR_NON_NULL(Mi, NAN);
@@ -594 +655 @@
                 continue;
             if (!unweighted_sum) {
-                wt = Wi->data.F32[yi][xi];
+                wt = covarFactor * Wi->data.F32[yi][xi];
                 if (wt == 0)
                     continue;
@@ -623 +684 @@
 }
 
-double pmSourceModelDotModel (const pmSource *Mi,
-                              const pmSource *Mj,
-                              const bool unweighted_sum) // should the cross product be weighted?
+double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor)
 {
     PS_ASSERT_PTR_NON_NULL(Mi, NAN);
@@ -677 +736 @@
                 flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
             } else {
-                wt = Wi->data.F32[yi][xi];
+                wt = covarFactor * Wi->data.F32[yi][xi];
                 if (wt > 0) {
                     flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
@@ -687 +746 @@
 }
 
-double pmSourceDataDotModel (const pmSource *Mi,
-                             const pmSource *Mj,
-                             const bool unweighted_sum) // should the cross product be weighted?
+double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor)
 {
     PS_ASSERT_PTR_NON_NULL(Mi, NAN);
@@ -741 +798 @@
                 flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
             } else {
-                wt = Wi->data.F32[yi][xi];
+                wt = covarFactor * Wi->data.F32[yi][xi];
                 if (wt > 0) {
                     flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;