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

Timestamp:

Apr 25, 2007, 3:20:29 PM (19 years ago)

Author:

magnier

Message:

incorporating updates from eam_02_branch (cached models, pmPSF FITS IO)

File:

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

trunk/psModules/src/objects/pmSourcePhotometry.c

@@ -3 +3 @@
  *  @author EAM, IfA; GLG, MHPCC
  *
- *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -106 +106 @@
     // replace source flux
     // XXX need to be certain which model and size of mask for prior subtraction
+    // XXX full model or just analytical?
+    // XXX use pmSourceAdd instead?
     if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
-        pmModelAdd (source->pixels, source->mask, model, false, false);
+        pmModelAdd (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL);
     }
 
@@ -180 +182 @@
     // set aperture mask circle to model radius
     // XXX use a different radius for apertures and fits...
-    psImageKeepCircle (source->mask, x, y, model->radiusFit, "OR", PM_MASK_MARK);
+    // XXX can I remove this?  the source should have the mask defined when it is constructed or
+    // when the fit / aperture radius is changed...
+    psImageKeepCircle (source->maskObj, x, y, model->radiusFit, "OR", PM_MASK_MARK);
 
     // measure the weight of included pixels
     // XXX is this supposed to use the weight or the flux?
     if (mode & PM_SOURCE_PHOT_WEIGHT) {
-        pmSourcePixelWeight (&source->pixWeight, model, source->pixels, source->mask);
+        pmSourcePixelWeight (&source->pixWeight, model, source->pixels, source->maskObj);
     }
 
     // measure object aperture photometry
-    status = pmSourcePhotometryAper  (&source->apMag, model, flux, source->mask);
+    status = pmSourcePhotometryAper  (&source->apMag, model, flux, source->maskObj);
 
     // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS
@@ -206 +210 @@
 
     // unmask aperture
-    psImageKeepCircle (source->mask, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));
+    // XXX can I remove this?  this will probably break things downstream...
+    psImageKeepCircle (source->maskObj, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));
 
     // if source was originally subtracted, re-subtract object, leave local sky
+    // XXX replace with pmSourceSub...
     if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
-        pmModelSub (source->pixels, source->mask, model, false, false);
+        pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL);
     }
 
@@ -363 +369 @@
 }
 
+# if (0)
 double pmSourceCrossProduct (const pmSource *Mi,
                              const pmSource *Mj,
@@ -379 +386 @@
     const psImage *Wi = Mi->weight;
 
-    const psImage *Ti = Mi->mask;
-    const psImage *Tj = Mj->mask;
+    const psImage *Ti = Mi->maskObj;
+    const psImage *Tj = Mj->maskObj;
 
     Xs = PS_MAX (Pi->col0, Pj->col0);
@@ -435 +442 @@
     const psImage *Wi = Mi->weight;
 
-    const psImage *Ti = Mi->mask;
-    const psImage *Tj = Mj->mask;
+    const psImage *Ti = Mi->maskObj;
+    const psImage *Tj = Mj->maskObj;
 
     Xs = PS_MAX (Pi->col0, Pj->col0);
@@ -483 +490 @@
     const psImage *Pi = Mi->pixels;
     const psImage *Wi = Mi->weight;
-    const psImage *Ti = Mi->mask;
+    const psImage *Ti = Mi->maskObj;
 
     // note that this is addressing the same image pixels,
@@ -521 +528 @@
     return flux;
 }
+# endif
 
 bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight)
@@ -543 +551 @@
 }
 
-// given a source, which model is currently appropriate?
-// choose PSF or EXT based on source->type, but fall back on PSF
-// if the EXT model is NULL
-pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source)
-{
-
-    pmModel *model;
-
-    switch (source->type) {
-    case PM_SOURCE_TYPE_STAR:
-        model = source->modelPSF;
-        if (model == NULL)
-            return NULL;
-        if (isPSF)
-            *isPSF = true;
-        return model;
-
-    case PM_SOURCE_TYPE_EXTENDED:
-        model = source->modelEXT;
-        if (isPSF)
-            *isPSF = false;
-        if (model == NULL) {
-            model = source->modelPSF;
-            if (isPSF)
-                *isPSF = true;
-        }
-        if (model == NULL) {
-            if (isPSF)
-                *isPSF = FALSE;
-            return NULL;
-        }
-        return (model);
-        break;
-
-    default:
-        if (isPSF)
-            *isPSF = false;
-        return NULL;
-    }
-    return NULL;
-}
+
+double pmSourceModelWeight(const pmSource *Mi,
+                      int term,
+                      const bool unweighted_sum) // should the cross product be weighted?
+{
+    double flux = 0, wt = 0, factor = 0;
+
+    const psImage *Pi = Mi->modelFlux;
+    const psImage *Wi = Mi->weight;
+    const psImage *Ti = Mi->maskObj;
+
+    for (int yi = 0; yi < Pi->numRows; yi++) {
+        for (int xi = 0; xi < Pi->numCols; xi++) {
+            if (Ti->data.U8[yi][xi])
+                continue;
+            if (!unweighted_sum) {
+                wt = 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");
+            }
+
+            if (unweighted_sum) {
+                flux += (factor * Pi->data.F32[yi][xi]);
+            } else {
+                flux += (factor * Pi->data.F32[yi][xi]) / wt;
+            }
+        }
+    }
+    return flux;
+}
+
+double pmSourceModelDotModel (const pmSource *Mi,
+                              const pmSource *Mj,
+                              const bool unweighted_sum) // should the cross product be weighted?
+{
+    int Xs, Xe, Ys, Ye;
+    int xi, xj, yi, yj;
+    int xIs, xJs, yIs, yJs;
+    int xIe, yIe;
+    double flux, wt;
+
+    const psImage *Pi = Mi->modelFlux;
+    const psImage *Pj = Mj->modelFlux;
+
+    const psImage *Wi = Mi->weight;
+
+    const psImage *Ti = Mi->maskObj;
+    const psImage *Tj = Mj->maskObj;
+
+    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.U8[yi][xi])
+                continue;
+            if (Tj->data.U8[yj][xj])
+                continue;
+
+            // XXX skip the nonsense weight pixels?
+            if (unweighted_sum) {
+                flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
+            } else {
+                wt = Wi->data.F32[yi][xi];
+                if (wt > 0) {
+                    flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
+                }
+            }
+        }
+    }
+    return flux;
+}
+
+double pmSourceDataDotModel (const pmSource *Mi,
+                             const pmSource *Mj,
+                             const bool unweighted_sum) // should the cross product be weighted?
+{
+
+    int Xs, Xe, Ys, Ye;
+    int xi, xj, yi, yj;
+    int xIs, xJs, yIs, yJs;
+    int xIe, yIe;
+    double flux, wt;
+
+    const psImage *Pi = Mi->pixels;
+    const psImage *Pj = Mj->modelFlux;
+
+    const psImage *Wi = Mi->weight;
+
+    const psImage *Ti = Mi->maskObj;
+    const psImage *Tj = Mj->maskObj;
+
+    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.U8[yi][xi])
+                continue;
+            if (Tj->data.U8[yj][xj])
+                continue;
+
+            // XXX skip the nonsense weight pixels?
+            if (unweighted_sum) {
+                flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
+            } else {
+                wt = Wi->data.F32[yi][xi];
+                if (wt > 0) {
+                    flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
+                }
+            }
+        }
+    }
+    return flux;
+}
+