Changeset 31361 for branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurveGenerate.c

Timestamp:

Apr 24, 2011, 10:09:38 PM (15 years ago)

Author:

eugene

Message:

determine the min kron radius = radius for bright PSF stars; mask interpolation used index instead of pixel coordinates in psImageInterpolate; add minKronRadius to pmSource; change errMag to psfMagErr; in growth curve, avoid perfect int radii (they can be inconsistent on + and - due to float precision; add pmGrowthCurveForSources; psfMagErr is always calculated from psfModel; apply ApTrend to all source psf mags and fluxes (not just psf sources); always use the psfModel for PSF_QF,_PERFECT; use BILINEAR to interpolate for aperture mags (BIQUAD implementation is wrong; apply growth curve to apMag & apFlux for all sources; save GrowthCurve with PSF model and read from PSF model

File:

• branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurveGenerate.c (modified) (2 diffs)

branches/eam_branches/ipp-20110404/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"
 
 /*****************************************************************************/
@@ -226 +225 @@
     return growth;
 }
+
+// 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);
+
+    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 (0)
+    FILE *f = fopen ("growth.dat", "w");
+    assert (f);
+
+    for (int j = 0; j < psf->growth->radius->n; j++) {
+
+        fprintf (f, "%f : ", psf->growth->radius->data.F32[j]);
+
+        // XXX dump the growth curves:
+        for (int i = 0; i < growths->n; i++) {
+            
+            pmGrowthCurve *growth = growths->data[i];
+            if (!growth) continue;
+
+            fprintf (f, "%8.4f ", growth->apMag->data.F32[j]);
+        }
+        fprintf (f, "\n");
+    }
+    fclose (f);
+# 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->fitMag);
+        }
+        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;
+    }
+    return growth;
+}
+
+# if (0)
+    // solve for the best growth and the offsets per star to minimize the scatter
+
+    // generate the per-star and per-radius sums
+    int Nradius = psf->growth->radius->n;
+    int Nstar = growths->n;
+
+    psVector *Mstar   = psVectorAlloc (Nstar, PS_DATA_F32);
+    psVector *Mradius = psVectorAlloc (Nradius, PS_DATA_F32);
+
+    psVectorInit (Mstar, 0.0);
+    psVectorInit (Mradius, 0.0);
+
+    for (int i = 0; i < Nstar; i++) {
+        pmGrowthCurve *growth = growths->data[i];
+        for (int j = 0; j < Nradius; j++) {
+            Mstar->data.F32[i] += growth->apMag->data.F32[j];
+            Mradius->data.F32[j] += growth->apMag->data.F32[j];
+        }
+    }
+
+    psImage *A = psImageAlloc(Nstar + Nradius, Nstar + Nradius, PS_DATA_F32);
+    psVector *B = psVectorAlloc(Nstar + Nradius, PS_DATA_F32);
+
+    psImageInit (A, 0.0);
+    psVectorInit (B, 0.0);
+
+    for (int i = 0; i < Nstar; i++) {
+        A->data.F32[i][i] = Nradius;
+        B->data.F32[i] = Mstar->data.F32[i];
+    }
+    for (int i = 0; i < Nradius; i++) {
+        int j = i + Nstar;
+        A->data.F32[j][j] = Nstar;
+        B->data.F32[j] = Mradius->data.F32[i];
+    }
+    
+    if (!psMatrixGJSolve(A, B)) {
+        psAbort("GJ failure");
+    }
+
+    for (int i = 0; i < Nradius; i++) {
+        psf->growth->apMag->data.F32[i] = B->data.F32[Nstar+i];
+    }
+
+    // XXX this analysis produces a biased growth curve: points with small radius are more
+    // likely to err on the side of being too faint (bacause of mis-aligned aperture), while on
+    // the bright side the tend to be too bright (because of contaminating neighbors).  Bright
+    // stars are less likely to be biased, 
+
+    // XXX at this point, we could iterate and exclude some of the stars with low data quality
+    // XXX we could also weight the above by flux or something
+
+# endif