Changeset 6859 for branches/rel10_ifa/psModules/src/objects/pmSourcePhotometry.c

Timestamp:

Apr 14, 2006, 11:43:59 AM (20 years ago)

Author:

magnier

Message:

added JPEG and PSF I/O methods to pmFPAfile, added pmSourcePixelWeights

File:

• branches/rel10_ifa/psModules/src/objects/pmSourcePhotometry.c (modified) (5 diffs)

branches/rel10_ifa/psModules/src/objects/pmSourcePhotometry.c

@@ -3 +3 @@
  *  @author EAM, IfA; GLG, MHPCC
  *
- *  @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-01 02:47:20 $
+ *  @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-14 21:43:59 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -25 +25 @@
 #include "pmSourcePhotometry.h"
 
+# define DO_SKY 0
+
 static float AP_MIN_SN = 0.0;
 
@@ -105 +107 @@
     status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
 
+    // measure the weight of included pixels
+    status = pmSourcePixelWeight (&source->pixWeight, model, source->pixels, source->mask);
+
     // measure object aperture photometry
     if (SN > AP_MIN_SN) {
@@ -169 +174 @@
 }
 
-# define DO_SKY 0
 // return source aperture magnitude
 bool pmSourcePhotometryAper (float *apMag, pmModel *model, psImage *image, psImage *mask)
@@ -202 +206 @@
 }
 
-# if (0)
-    // XXX split into ap, psf, ext mags?
-    bool pmSourcePhotometry (float *fitMag, float *apMag, pmModel *model, psImage *image, psImage *mask)
-{
-
-    float obsSum = 0;
-    float fitSum = 0;
-    float sky = model->params->data.F32[0];
-
-    *fitMag = 99.0;
-    *obsMag = 99.0;
-
-    // measure fitMag
-    pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);
-    fitSum = modelFluxFunc (model->params);
-    if (fitSum <= 0)
-        return false;
-    if (!isfinite(fitSum))
-        return false;
-    *fitMag = -2.5*log10(fitSum);
-
-    // measure apMag
-    for (int ix = 0; ix < image->numCols; ix++) {
-        for (int iy = 0; iy < image->numRows; iy++) {
-            if (mask->data.U8[iy][ix])
-                continue;
-            obsSum += image->data.F32[iy][ix] - sky;
-        }
-    }
-    if (obsSum <= 0)
-        return false;
-    *obsMag = -2.5*log10(obsSum);
-
+// return source aperture magnitude
+bool pmSourcePixelWeight (float *pixWeight, pmModel *model, psImage *image, psImage *mask)
+{
+    float modelSum = 0;
+    float validSum = 0;
+    float sky = 0;
+    float value;
+
+    int Xo, Yo, dP;
+    int dX, DX, NX;
+    int dY, DY, NY;
+
+    *pixWeight = 0.0;
+
+    if (model == NULL)
+        return false;
+
+    // we only care about the value of the object model, not the local sky
+    if (DO_SKY) {
+        sky = model->params->data.F32[0];
+    } else {
+        sky = 0;
+    }
+
+    // the model function returns the source flux at a position
+    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
+    psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
+    psVector *params = model->params;
+
+    Xo = params->data.F32[2];
+    Yo = params->data.F32[3];
+
+    dX = Xo - image->col0;
+    dP = image->numCols - dX;
+    DX = PS_MAX(dX, dP);
+    NX = image->numCols;
+
+    dY = Yo - image->row0;
+    dP = image->numRows - dY;
+    DY = PS_MAX(dY, dP);
+    NY = image->numRows;
+
+    // measure modelSum and validSum
+    // XXX this does not work for sources near the edge: we need to measure for
+    // a full square region
+    for (int ix = -DX; ix < DX + 1; ix++) {
+        int mx = ix + dX;
+        for (int iy = -DY; iy < DY + 1; iy++) {
+            int my = iy + dY;
+
+            coord->data.F32[0] = (psF32) (ix + Xo);
+            coord->data.F32[1] = (psF32) (iy + Yo);
+
+            // for the full model, add all points
+            value = modelFunc (NULL, params, coord) - sky;
+            modelSum += value;
+
+            // include count only the unmasked pixels within the image area
+            if (mx < 0)
+                continue;
+            if (my < 0)
+                continue;
+            if (mx >= NX)
+                continue;
+            if (my >= NY)
+                continue;
+            if (mask->data.U8[my][mx])
+                continue;
+
+            validSum += value;
+        }
+    }
+    if (validSum <= 0)
+        return false;
+
+    *pixWeight = validSum / modelSum;
     return (true);
 }
-# endif
 
 float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj)