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

Timestamp:

Jun 26, 2012, 11:31:38 AM (14 years ago)

Author:

eugene

Message:

re-enable MODEL_VAR option for linear photometry fit

Location:

trunk/psModules

Files:

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

trunk/psModules/src/objects/pmSourcePhotometry.c

@@ -899 +899 @@
 }
 
-# if (HAVE_MODEL_VAR) 
 double pmSourceModelWeight(const pmSource *Mi, int term, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
-# else
-double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
-# endif
 {
     PS_ASSERT_PTR_NON_NULL(Mi, NAN);
-# if (HAVE_MODEL_VAR) 
-    double flux = 0, wt = 1.0, factor = 0;
-# else
-    double flux = 0, wt = 0, factor = 0;
-# endif
+    double flux = 0;
+    double wt = 1.0;
+    double factor = 0;
 
     const psImage *Pi = Mi->modelFlux;
     assert (Pi != NULL);
 
-# if (HAVE_MODEL_VAR)
     const psImage *Wi = NULL;
     switch (fitVarMode) {
@@ -921 +914 @@
         break;
       case PM_SOURCE_PHOTFIT_IMAGE_VAR:
+      case PM_SOURCE_PHOTFIT_MODEL_SKY:
+        Wi = Mi->variance;
+        psAssert (Wi, "programming error");
+        break;
+      case PM_SOURCE_PHOTFIT_MODEL_VAR:
+        Wi = Mi->modelVar;
+        psAssert (Wi, "programming error");
+        break;
+      case PM_SOURCE_PHOTFIT_NONE:
+        psAbort("programming error");
+    }   
+    const psImage *Ti = Mi->maskObj;
+    assert (Ti != NULL);
+
+    for (int yi = 0; yi < Pi->numRows; yi++) {
+        for (int xi = 0; xi < Pi->numCols; xi++) {
+            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
+                continue;
+            if (fitVarMode != PM_SOURCE_PHOTFIT_CONST) {
+                wt = covarFactor * Wi->data.F32[yi][xi];
+                if (wt == 0) continue;
+            }
+            switch (term) {
+              case 0:
+                factor = 1;
+                break;
+              case 1:
+                factor = xi + Pi->col0;
+                break;
+              case 2:
+                factor = yi + Pi->row0;
+                break;
+              default:
+                psAbort("invalid term for pmSourceWeight");
+            }
+
+            // wt is 1.0 for CONST
+            flux += (factor * Pi->data.F32[yi][xi]) / wt;
+        }
+    }
+    return flux;
+}
+
+double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
+{
+    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
+    PS_ASSERT_PTR_NON_NULL(Mj, NAN);
+    int Xs, Xe, Ys, Ye;
+    int xi, xj, yi, yj;
+    int xIs, xJs, yIs, yJs;
+    int xIe, yIe;
+    double flux;
+    double wt = 1.0;
+
+    const psImage *Pi = Mi->modelFlux;
+    assert (Pi != NULL);
+    const psImage *Pj = Mj->modelFlux;
+    assert (Pj != NULL);
+
+    const psImage *Wi = NULL;
+    switch (fitVarMode) {
+      case PM_SOURCE_PHOTFIT_CONST:
+        break;
+      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
+      case PM_SOURCE_PHOTFIT_MODEL_SKY:
+        Wi = Mi->variance;
+        psAssert (Wi, "programming error");
+        break;
+      case PM_SOURCE_PHOTFIT_MODEL_VAR:
+        Wi = Mi->modelVar;
+        psAssert (Wi, "programming error");
+        break;
+      case PM_SOURCE_PHOTFIT_NONE:
+        psAbort("programming error");
+    }   
+
+    const psImage *Ti = Mi->maskObj;
+    assert (Ti != NULL);
+    const psImage *Tj = Mj->maskObj;
+    assert (Tj != NULL);
+
+    Xs = PS_MAX (Pi->col0, Pj->col0);
+    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
+
+    Ys = PS_MAX (Pi->row0, Pj->row0);
+    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
+
+    xIs = Xs - Pi->col0;
+    xJs = Xs - Pj->col0;
+    yIs = Ys - Pi->row0;
+    yJs = Ys - Pj->row0;
+
+    xIe = Xe - Pi->col0;
+    yIe = Ye - Pi->row0;
+
+    // note that weight is addressing the same image pixels
+    flux = 0;
+    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
+        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
+            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
+                continue;
+            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj] & maskVal)
+                continue;
+
+            float value = (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
+            switch (fitVarMode) {
+              case PM_SOURCE_PHOTFIT_CONST:
+                wt = 1.0;
+                break;
+              case PM_SOURCE_PHOTFIT_IMAGE_VAR:
+              case PM_SOURCE_PHOTFIT_MODEL_SKY:
+              case PM_SOURCE_PHOTFIT_MODEL_VAR:
+                wt = covarFactor * Wi->data.F32[yi][xi];
+                break;
+              case PM_SOURCE_PHOTFIT_NONE:
+                psAbort("programming error");
+            }
+            // skip pixels with nonsense weight values
+            if (wt <= 0) continue;
+
+            flux += value / wt;
+        }
+    }
+    return flux;
+}
+
+double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
+{
+    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
+    PS_ASSERT_PTR_NON_NULL(Mj, NAN);
+    int Xs, Xe, Ys, Ye;
+    int xi, xj, yi, yj;
+    int xIs, xJs, yIs, yJs;
+    int xIe, yIe;
+    double flux;
+    double wt = 1.0;
+
+    const psImage *Pi = Mi->pixels;
+    assert (Pi != NULL);
+    const psImage *Pj = Mj->modelFlux;
+    assert (Pj != NULL);
+
+    const psImage *Wi = NULL;
+    switch (fitVarMode) {
+      case PM_SOURCE_PHOTFIT_CONST:
+        break;
+      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
+      case PM_SOURCE_PHOTFIT_MODEL_SKY:
         Wi = Mi->variance;
         psAssert (Wi, "programming error");
@@ -931 +1072 @@
         psAbort("programming error");
     }   
-# else
-    const psImage *Wi = Mi->variance;
-    if (!unweighted_sum) {
-        assert (Wi != NULL);
-    }
-# endif
-    const psImage *Ti = Mi->maskObj;
-    assert (Ti != NULL);
-
-    for (int yi = 0; yi < Pi->numRows; yi++) {
-        for (int xi = 0; xi < Pi->numCols; xi++) {
-            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
-                continue;
-# if (HAVE_MODEL_VAR)
-            if (fitVarMode != PM_SOURCE_PHOTFIT_CONST) {
-                wt = covarFactor * Wi->data.F32[yi][xi];
-                if (wt == 0) continue;
-            }
-# else
-            if (!unweighted_sum) {
-                wt = covarFactor * Wi->data.F32[yi][xi];
-                if (wt == 0)
-                    continue;
-            }
-# endif
-
-            switch (term) {
-              case 0:
-                factor = 1;
-                break;
-              case 1:
-                factor = xi + Pi->col0;
-                break;
-              case 2:
-                factor = yi + Pi->row0;
-                break;
-              default:
-                psAbort("invalid term for pmSourceWeight");
-            }
-
-# if (HAVE_MODEL_VAR)
-            // wt is 1.0 for CONST
-            flux += (factor * Pi->data.F32[yi][xi]) / wt;
-# else
-            if (unweighted_sum) {
-                flux += (factor * Pi->data.F32[yi][xi]);
-            } else {
-                flux += (factor * Pi->data.F32[yi][xi]) / wt;
-            }
-# endif
-        }
-    }
-    return flux;
-}
-
-# if (HAVE_MODEL_VAR)
-double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
-# else
-double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
-# endif
-{
-    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
-    PS_ASSERT_PTR_NON_NULL(Mj, NAN);
-    int Xs, Xe, Ys, Ye;
-    int xi, xj, yi, yj;
-    int xIs, xJs, yIs, yJs;
-    int xIe, yIe;
-# if (HAVE_MODEL_VAR)
-    double flux;
-    double wt = 1.0;
-# else
-    double flux, wt;
-# endif
-
-    const psImage *Pi = Mi->modelFlux;
-    assert (Pi != NULL);
-    const psImage *Pj = Mj->modelFlux;
-    assert (Pj != NULL);
-
-# if (HAVE_MODEL_VAR)
-    const psImage *Wi = NULL;
-    switch (fitVarMode) {
-      case PM_SOURCE_PHOTFIT_CONST:
-        break;
-      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
-        Wi = Mi->variance;
-        psAssert (Wi, "programming error");
-        break;
-      case PM_SOURCE_PHOTFIT_MODEL_VAR:
-        Wi = Mi->modelVar;
-        psAssert (Wi, "programming error");
-        break;
-      case PM_SOURCE_PHOTFIT_NONE:
-        psAbort("programming error");
-    }   
-# else
-    const psImage *Wi = Mi->variance;
-    if (!unweighted_sum) {
-        assert (Wi != NULL);
-    }
-# endif
 
     const psImage *Ti = Mi->maskObj;
@@ -1052 +1092 @@
     yIe = Ye - Pi->row0;
 
-    // note that weight is addressing the same image pixels
+    // note that weight is addressing the same image pixels,
     flux = 0;
     for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
@@ -1061 +1101 @@
                 continue;
 
-# if (HAVE_MODEL_VAR)
             float value = (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
             switch (fitVarMode) {
@@ -1068 +1107 @@
                 break;
               case PM_SOURCE_PHOTFIT_IMAGE_VAR:
+              case PM_SOURCE_PHOTFIT_MODEL_SKY:
               case PM_SOURCE_PHOTFIT_MODEL_VAR:
                 wt = covarFactor * Wi->data.F32[yi][xi];
@@ -1078 +1118 @@
 
             flux += value / wt;
-# else
-            // XXX skip the nonsense weight pixels?
-            if (unweighted_sum) {
-                flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
-            } else {
-                wt = covarFactor * Wi->data.F32[yi][xi];
-                if (wt > 0) {
-                    flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
-                }
-            }
-# endif
         }
     }
     return flux;
 }
-
-# if (HAVE_MODEL_VAR)
-double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
-# else
-double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
-# endif
-{
-    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
-    PS_ASSERT_PTR_NON_NULL(Mj, NAN);
-    int Xs, Xe, Ys, Ye;
-    int xi, xj, yi, yj;
-    int xIs, xJs, yIs, yJs;
-    int xIe, yIe;
-# if (HAVE_MODEL_VAR)
-    double flux;
-    double wt = 1.0;
-# else
-    double flux, wt;
-# endif
-
-    const psImage *Pi = Mi->pixels;
-    assert (Pi != NULL);
-    const psImage *Pj = Mj->modelFlux;
-    assert (Pj != NULL);
-
-# if (HAVE_MODEL_VAR)
-    const psImage *Wi = NULL;
-    switch (fitVarMode) {
-      case PM_SOURCE_PHOTFIT_CONST:
-        break;
-      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
-        Wi = Mi->variance;
-        psAssert (Wi, "programming error");
-        break;
-      case PM_SOURCE_PHOTFIT_MODEL_VAR:
-        Wi = Mi->modelVar;
-        psAssert (Wi, "programming error");
-        break;
-      case PM_SOURCE_PHOTFIT_NONE:
-        psAbort("programming error");
-    }   
-# else
-    const psImage *Wi = Mi->variance;
-    if (!unweighted_sum) {
-        assert (Wi != NULL);
-    }
-# endif
-
-    const psImage *Ti = Mi->maskObj;
-    assert (Ti != NULL);
-    const psImage *Tj = Mj->maskObj;
-    assert (Tj != NULL);
-
-    Xs = PS_MAX (Pi->col0, Pj->col0);
-    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
-
-    Ys = PS_MAX (Pi->row0, Pj->row0);
-    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
-
-    xIs = Xs - Pi->col0;
-    xJs = Xs - Pj->col0;
-    yIs = Ys - Pi->row0;
-    yJs = Ys - Pj->row0;
-
-    xIe = Xe - Pi->col0;
-    yIe = Ye - Pi->row0;
-
-    // note that weight is addressing the same image pixels,
-    flux = 0;
-    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
-        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
-            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
-                continue;
-            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj] & maskVal)
-                continue;
-
-# if (HAVE_MODEL_VAR)
-            float value = (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
-            switch (fitVarMode) {
-              case PM_SOURCE_PHOTFIT_CONST:
-                wt = 1.0;
-                break;
-              case PM_SOURCE_PHOTFIT_IMAGE_VAR:
-              case PM_SOURCE_PHOTFIT_MODEL_VAR:
-                wt = covarFactor * Wi->data.F32[yi][xi];
-                break;
-              case PM_SOURCE_PHOTFIT_NONE:
-                psAbort("programming error");
-            }
-            // skip pixels with nonsense weight values
-            if (wt <= 0) continue;
-
-            flux += value / wt;
-
-# else
-            // XXX skip the nonsense weight pixels?
-            if (unweighted_sum) {
-                flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
-            } else {
-                wt = covarFactor * Wi->data.F32[yi][xi];
-                if (wt > 0) {
-                    flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
-                }
-            }
-# endif
-        }
-    }
-    return flux;
-}