Changeset 31451 for trunk/psModules/src/objects/pmGrowthCurveGenerate.c

Timestamp:

May 5, 2011, 11:02:53 AM (15 years ago)

Author:

eugene

Message:

merge updates from eam branch 20110404

File:

• trunk/psModules/src/objects/pmGrowthCurveGenerate.c (modified) (3 diffs)

trunk/psModules/src/objects/pmGrowthCurveGenerate.c

@@ -50 +50 @@
 
 #include "pmSourcePhotometry.h"
-
-pmGrowthCurve *pmGrowthCurveForPosition (psImage *image, pmPSF *psf, bool ignore, psImageMaskType maskVal, psImageMaskType markVal, float xc, float yc);
+#include "pmGrowthCurveGenerate.h"
 
 /*****************************************************************************/
@@ -197 +196 @@
         // mask the given aperture and measure the apMag
         psImageKeepCircle (mask, xc, yc, radius, "OR", markVal);
-        if (!pmSourcePhotometryAper (&apMag, NULL, NULL, model, pixels, NULL, mask, maskVal)) {
+        if (!pmSourcePhotometryAper (NULL, &apMag, NULL, NULL, model, pixels, NULL, mask, maskVal)) {
             psFree (growth);
             psFree (view);
@@ -226 +225 @@
     return growth;
 }
+
+# define DEBUG 0
+# if (DEBUG)
+static FILE *fgr = NULL;
+# endif
+
+// we generate the growth curve for the center of the image with the specified psf model
+bool pmGrowthCurveGenerateFromSources (pmReadout *readout, pmPSF *psf, psArray *sources, bool INTERPOLATE_AP, psImageMaskType maskVal, psImageMaskType markVal)
+{
+    PS_ASSERT_PTR_NON_NULL(readout, false);
+    PS_ASSERT_PTR_NON_NULL(readout->image, false);
+
+    // maskVal is used to test for rejected pixels, and must include markVal
+    maskVal |= markVal;
+
+    pmSourcePhotometryMode photMode = INTERPOLATE_AP ? PM_SOURCE_PHOT_INTERP : 0;
+    
+    // measure the growth curve for each PSF source and average them together
+    psArray *growths = psArrayAllocEmpty (100);
+
+# if (DEBUG)
+    fgr = fopen ("growth.mags.dat", "w");
+# endif
+
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+
+        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
+
+        pmGrowthCurve *growth = pmGrowthCurveForSource (source, psf, photMode, maskVal, markVal);
+        if (!growth) continue;
+        
+        psArrayAdd (growths, 100, growth);
+        psFree (growth);
+    }
+    psAssert (growths->n, "cannot build growth curve (no valid PSF stars?)");
+
+# if (DEBUG)
+    fclose (fgr);
+# endif
+
+    // just use a simple sample median to get the 'best' value from each growth curve...
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+    psVector *values = psVectorAlloc (growths->n, PS_DATA_F32);
+
+    // loop over a range of source fluxes
+    // no need to interpolate since we have forced the object center
+    // to 0.5, 0.5 above
+    for (int i = 0; i < psf->growth->radius->n; i++) {
+
+        // median the values for each radial bin
+        values->n = 0;
+        for (int j = 0; j < growths->n; j++) {
+            pmGrowthCurve *growth = growths->data[j];
+            if (!isfinite(growth->apMag->data.F32[i])) continue;
+            psVectorAppend (values, growth->apMag->data.F32[i] - growth->refMag);
+        }
+        if (values->n == 0) {
+            psf->growth->apMag->data.F32[i] = NAN;
+        } else {
+            if (!psVectorStats (stats, values, NULL, NULL, 0)) {
+                psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+                return false;
+            }
+            psf->growth->apMag->data.F32[i] = stats->sampleMedian;
+        }
+    }
+
+    psf->growth->fitMag = psf->growth->apMag->data.F32[psf->growth->radius->n-1];
+    psf->growth->apRef = psVectorInterpolate (psf->growth->radius, psf->growth->apMag, psf->growth->refRadius);
+    psf->growth->apLoss = psf->growth->fitMag - psf->growth->apRef;
+
+    psLogMsg ("psphot.growth", 4, "GrowthCurve : apLoss : %f (fitMag - apMag @ ref : %f - %f)\n", psf->growth->apLoss, psf->growth->fitMag, psf->growth->apRef);
+
+    psFree (growths);
+    psFree (stats);
+    psFree (values);
+
+    return true;
+}
+
+pmGrowthCurve *pmGrowthCurveForSource (pmSource *source, pmPSF *psf, pmSourcePhotometryMode photMode, psImageMaskType maskVal, psImageMaskType markVal) {
+
+    float radius;
+
+    assert (psf->growth);
+
+    float minRadius = psf->growth->radius->data.F32[0];
+    pmGrowthCurve *growth = pmGrowthCurveAlloc (minRadius, psf->growth->maxRadius, psf->growth->refRadius);
+
+    // measure the fitMag for this source (for normalization)
+    // pmSourcePhotometryModel (&fitMag, NULL, source->psfModel);
+    growth->fitMag = source->psfMag;
+
+    float xc = source->peak->xf;
+    float yc = source->peak->yf;
+
+    // Loop over the range of radii
+    for (int i = 0; i < growth->radius->n; i++) {
+
+        radius = growth->radius->data.F32[i];
+
+        // mask the given aperture and measure the apMag
+        psImageKeepCircle (source->maskObj, xc, yc, radius, "OR", markVal);
+
+        if (!pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, radius)) {
+            psFree (growth);
+            return NULL;
+        }
+
+        // if (!pmSourcePhotometryAper (NULL, &apMag, NULL, NULL, NULL, source->pixels, NULL, source->maskObj, maskVal)) {
+        //     psFree (growth);
+        //     return NULL;
+        // }
+
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
+
+        growth->apMag->data.F32[i] = source->apMag;
+    }
+    psAssert(growth->refBin >= 0, "invalid growth reference bin");
+    psAssert(growth->refBin < growth->apMag->n, "invalid growth reference bin");
+    growth->refMag = growth->apMag->data.F32[growth->refBin];
+
+    // Loop over the range of radii
+# if (DEBUG)
+    for (int i = 0; i < growth->radius->n; i++) {
+        fprintf (fgr, "%f %f  %f %f %f %f\n", xc, yc, growth->radius->data.F32[i], growth->apMag->data.F32[i], growth->fitMag, growth->refMag);
+    }
+# endif
+
+    return growth;
+}