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

Timestamp:

Oct 2, 2009, 3:11:32 PM (17 years ago)

Author:

eugene

Message:

check in changes from genes development branch : extensive changes to moments calculation, psf model generation, aperture residuals

File:

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

trunk/psModules/src/objects/pmSourcePhotometry.c

@@ -84 +84 @@
 
     // we must have a valid model
+    // XXX allow aperture magnitudes for sources without a model
     model = pmSourceGetModel (&isPSF, source);
     if (model == NULL) {
@@ -90 +91 @@
     }
 
+    // XXX handle negative flux, low-significance
     if (model->dparams->data.F32[PM_PAR_I0] > 0) {
         SN = model->params->data.F32[PM_PAR_I0] / model->dparams->data.F32[PM_PAR_I0];
@@ -101 +103 @@
 
     // measure PSF model photometry
-    if (psf->FluxScale) {
+    // XXX TEST: do not use flux scale
+    if (0 && 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;
-            }
-        }
+        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]);
     } else {
         status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
     }
 
-    // if we have a collection of model fits, one of them is a pointer to modelEXT?
-    // XXX not guaranteed
+    // if we have a collection of model fits, check if one of them is a pointer to modelEXT
     if (source->modelFits) {
-      for (int i = 0; i < source->modelFits->n; i++) {
-        pmModel *model = source->modelFits->data[i];
-        status = pmSourcePhotometryModel (&model->mag, model);
-      }
-      if (source->modelEXT) {
-        source->extMag = source->modelEXT->mag;
-      }
+        bool foundEXT = false;
+        for (int i = 0; i < source->modelFits->n; i++) {
+            pmModel *model = source->modelFits->data[i];
+            status = pmSourcePhotometryModel (&model->mag, model);
+            if (model == source->modelEXT) foundEXT = true;
+        }
+        if (foundEXT) {
+            source->extMag = source->modelEXT->mag;
+        } else {
+            status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
+        }
     } else {
-      if (source->modelEXT) {
-        status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
-      }
+        if (source->modelEXT) {
+            status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
+        }
     }
 
@@ -160 +157 @@
     }
 
-    // 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->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
-        pmModelAdd (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
-    }
+    // if we measure aperture magnitudes, the source must not currently be subtracted!
+    psAssert (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED), "cannot measure ap mags if source is subtracted!");
 
     // if we are measuring aperture photometry and applying the growth correction,
@@ -201 +193 @@
     if (isfinite (source->apMag) && isPSF && psf) {
         if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) {
-            source->apMag += pmGrowthCurveCorrect (psf->growth, model->radiusFit);
+            source->apMag += pmGrowthCurveCorrect (psf->growth, source->apRadius);
         }
         if (mode & PM_SOURCE_PHOT_APCORR) {
+            // XXX this should be removed -- we no longer fit for the 'sky bias'
             rflux   = pow (10.0, 0.4*source->psfMag);
-            source->apMag -= PS_SQR(model->radiusFit)*rflux * psf->skyBias + psf->skySat / rflux;
+            source->apMag -= PS_SQR(source->apRadius)*rflux * psf->skyBias + psf->skySat / rflux;
         }
     }
@@ -211 +204 @@
         psFree(flux);
         psFree(mask);
-    }
-
-    // if source was originally subtracted, re-subtract object, leave local sky
-    // XXX replace with pmSourceSub...
-    if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
-        pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
     }
 
@@ -763 +750 @@
     return flux;
 }
-
-// XXX this is test code to verify the shift is doing the right thing (seems to be)
-# if (0)
-// measure centroid of unshifted gaussian (should be 16.0,16.0)
-        {
-          psImage *image = source->pixels;
-          float xo = 0.0;
-          float yo = 0.0;
-          float xo2 = 0.0;
-          float yo2 = 0.0;
-          float no = 0.0;
-          for (int j = 0; j < image->numRows; j++)
-          {
-            for (int i = 0; i < image->numCols; i++) {
-              xo += i*image->data.F32[j][i];
-              yo += j*image->data.F32[j][i];
-              xo2 += i*i*image->data.F32[j][i];
-              yo2 += j*j*image->data.F32[j][i];
-              no += image->data.F32[j][i];
-            }
-          }
-          xo /= no;
-          yo /= no;
-          xo2 = sqrt (xo2/no - xo*xo);
-          yo2 = sqrt (yo2/no - yo*yo);
-          fprintf (stderr, "pre-shift centroid: %f,%f, sigma: %f,%f: flux: %f\n", xo, yo, xo2, yo2, no);
-        }
-
-// measure centroid of unshifted gaussian (should be 16.0,16.0)
-        {
-          psImage *image = flux;
-          float xo = 0.0;
-          float yo = 0.0;
-          float xo2 = 0.0;
-          float yo2 = 0.0;
-          float no = 0.0;
-          for (int j = 0; j < image->numRows; j++)
-          {
-            for (int i = 0; i < image->numCols; i++) {
-              xo += i*image->data.F32[j][i];
-              yo += j*image->data.F32[j][i];
-              xo2 += i*i*image->data.F32[j][i];
-              yo2 += j*j*image->data.F32[j][i];
-              no += image->data.F32[j][i];
-            }
-          }
-          xo /= no;
-          yo /= no;
-          xo2 = sqrt (xo2/no - xo*xo);
-          yo2 = sqrt (yo2/no - yo*yo);
-          fprintf (stderr, "pre-shift centroid: %f,%f, sigma: %f,%f: flux: %f\n", xo, yo, xo2, yo2, no);
-        }
-# endif