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

Timestamp:

Mar 10, 2008, 2:49:14 PM (18 years ago)

Author:

Paul Price

Message:

The NANs in the pmPSFtry->metric vector were coming from
pmSourceMagnitudes with mode=PM_SOURCE_PHOT_INTERP, and some bad
pixels in the source postage stamp. The NANs were effectively grown
by the call to psImageShift, while the mask was not grown. This
results in NAN pixels being accumulated in pmSourcePhotometryAper,
resulting in a NAN aperMag, so that the metric is also NAN.

Another possible error I found while debugging is that
pmSourcePixelWeight and pmSourcePhotometryAper should be run with the
mask value set to "maskVal | mark", so that only the circular aperture
is used.

File:

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

trunk/psModules/src/objects/pmSourcePhotometry.c

@@ -3 +3 @@
  *  @author EAM, IfA; GLG, MHPCC
  *
- *  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-02-03 21:58:56 $
+ *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-03-11 00:49:14 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -85 +85 @@
     model = pmSourceGetModel (&isPSF, source);
     if (model == NULL) {
-        psTrace ("psModules.objects", 3, "fail mag : no valid model");
+        psTrace ("psModules.objects", 3, "fail mag : no valid model");
         return false;
     }
@@ -101 +101 @@
     // measure PSF model photometry
     if (psf->FluxScale) {
-        // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
-        double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y);
-        if (!isfinite(fluxScale)) {
-            // XXX objects on the edge can be slightly outside -- if we get an 
-            // error, use the full fit.
-            psErrorClear();
-            status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
-        } else {
-            if (isfinite(fluxScale) && (fluxScale > 0.0)) {
-                source->psfMag = -2.5*log10(fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]);
-            } else {
-                source->psfMag = NAN;
-            }
-        }
+        // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
+        double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y);
+        if (!isfinite(fluxScale)) {
+            // XXX objects on the edge can be slightly outside -- if we get an
+            // error, use the full fit.
+            psErrorClear();
+            status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
+        } else {
+            if (isfinite(fluxScale) && (fluxScale > 0.0)) {
+                source->psfMag = -2.5*log10(fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]);
+            } else {
+                source->psfMag = NAN;
+            }
+        }
     } else {
-        status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
+        status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
     }
 
@@ -124 +124 @@
     // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS
     if ((mode & PM_SOURCE_PHOT_APCORR) && isPSF && psf && psf->ApTrend) {
-        // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
-        double apTrend = pmTrend2DEval (psf->ApTrend, (float)source->peak->x, (float)source->peak->y);
+        // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
+        double apTrend = pmTrend2DEval (psf->ApTrend, (float)source->peak->x, (float)source->peak->y);
         source->psfMag += apTrend;
     }
 
     if (!isfinite(SN) || (SN < AP_MIN_SN)) {
-        psTrace ("psModules.objects", 3, "fail mag : bad SN: %f (limit: %f)", SN, AP_MIN_SN);
+        psTrace ("psModules.objects", 3, "fail mag : bad SN: %f (limit: %f)", SN, AP_MIN_SN);
         return false;
     }
@@ -145 +145 @@
     // we need to shift the flux in the selected pixels (but not the mask)
     // psImageShift ();
-    psImage *flux;
+    psImage *flux = NULL, *mask = NULL; // Star flux and mask images, to photometer
     if (mode & PM_SOURCE_PHOT_INTERP) {
         float dx = 0.5 - x + (int)x;
@@ -151 +151 @@
         x += dx;
         y += dy;
-        flux = psImageShift (NULL, source->pixels, dx, dy, 0.0, PS_INTERPOLATE_BICUBE);
-        // flux = psImageShift (NULL, source->pixels, dx, dy, 0.0, PS_INTERPOLATE_BILINEAR);
+
+        if (!psImageShiftMask(&flux, &mask, source->pixels, source->maskObj, maskVal, dx, dy,
+                              NAN, 0xff, PS_INTERPOLATE_BICUBE)) {
+            // Not much we can do about it
+            psErrorClear();
+            psTrace ("psModules.objects", 3, "fail shift");
+            return false;
+        }
 
         // XXX this is test code to verify the shift is doing the right thing (seems to be)
@@ -209 +215 @@
     } else {
         flux = source->pixels;
+        mask = source->maskObj;
     }
     // set aperture mask circle to model radius
@@ -219 +226 @@
     // XXX is this supposed to use the weight or the flux?
     if (mode & PM_SOURCE_PHOT_WEIGHT) {
-        pmSourcePixelWeight (&source->pixWeight, model, source->pixels, source->maskObj, maskVal);
+        pmSourcePixelWeight (&source->pixWeight, model, source->pixels, source->maskObj, maskVal | mark);
     }
 
     // measure object aperture photometry
-    status = pmSourcePhotometryAper  (&source->apMag, model, flux, source->maskObj, maskVal);
+    status = pmSourcePhotometryAper  (&source->apMag, model, flux, mask, maskVal | mark);
     if (!status) {
-        psTrace ("psModules.objects", 3, "fail mag : bad Ap Mag");
+        psTrace ("psModules.objects", 3, "fail mag : bad Ap Mag");
         psErrorCode last = psErrorCodeLast();
-        // XXX flag the object?
+        // XXX flag the object?
         if (last == PM_ERR_PHOTOM) {
             // the aper mag was undefined (flux < 0).  bad object, but
@@ -246 +253 @@
     }
     if (mode & PM_SOURCE_PHOT_INTERP) {
-        psFree (flux);
+        psFree(flux);
+        psFree(mask);
     }