Changeset 31361

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

Location:

branches/eam_branches/ipp-20110404/psModules/src

Files:

• camera/pmFPAMaskWeight.c (modified) (1 diff)
• imcombine/pmPSFEnvelope.c (modified) (1 diff)
• objects/Makefile.am (modified) (1 diff)
• objects/pmGrowthCurve.c (modified) (1 diff)
• objects/pmGrowthCurveGenerate.c (modified) (2 diffs)
• objects/pmGrowthCurveGenerate.h (added)
• objects/pmPSF.h (modified) (1 diff)
• objects/pmPSF_IO.c (modified) (4 diffs)
• objects/pmSource.h (modified) (3 diffs)
• objects/pmSourceIO_CMF_PS1_V3.c (modified) (4 diffs)
• objects/pmSourceMoments.c (modified) (21 diffs)
• objects/pmSourcePhotometry.c (modified) (9 diffs)
• objects/pmSourceUtils.h (modified) (1 diff)
• psmodules.h (modified) (1 diff)

branches/eam_branches/ipp-20110404/psModules/src/camera/pmFPAMaskWeight.c

@@ -587 +587 @@
                 double imageValue, varianceValue; // Image and variance value from interpolation
                 psImageMaskType maskValue = 0; // Mask value from interpolation
-                psImageInterpolateStatus status = psImageInterpolate(&imageValue, &varianceValue, &maskValue, x, y, interp);
+
+                // interpolate to pixel center (index + 0.5)
+                psImageInterpolateStatus status = psImageInterpolate(&imageValue, &varianceValue, &maskValue, x + 0.5, y + 0.5, interp);
                 if (status == PS_INTERPOLATE_STATUS_ERROR || status == PS_INTERPOLATE_STATUS_OFF) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to interpolate readout at %d,%d", x, y);

branches/eam_branches/ipp-20110404/psModules/src/imcombine/pmPSFEnvelope.c

@@ -380 +380 @@
 
         // measure the source moments: tophat windowing, no pixel S/N cutoff
-        if (!pmSourceMoments(source, maxRadius, 0.0, 0.0, maskVal)) {
+        if (!pmSourceMoments(source, maxRadius, 0.0, 0.0, 0.0, maskVal)) {
             // Can't do anything about it; limp along as best we can
             psErrorClear();

branches/eam_branches/ipp-20110404/psModules/src/objects/Makefile.am

@@ -114 +114 @@
         pmTrend2D.h \
         pmGrowthCurve.h \
+        pmGrowthCurveGenerate.h \
         pmSourceMatch.h \
         pmDetEff.h \

branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurve.c

@@ -81 +81 @@
     while (radius < Rlin) {
         // fprintf (stderr, "r: %f\n", radius);
-        psVectorAppend (growth->radius, radius);
+        psVectorAppend (growth->radius, radius - 0.001);
         radius += 1.0;
     }    
     while (radius < maxRadius) {
         // fprintf (stderr, "r: %f\n", radius);
-        psVectorAppend (growth->radius, radius);
+        psVectorAppend (growth->radius, radius - 0.001);
         radius *= fR;
         radius = (int) (radius + 0.5);

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

branches/eam_branches/ipp-20110404/psModules/src/objects/pmPSF.h

@@ -94 +94 @@
 double pmPSF_SXYtoModel (psF32 *fittedPar);
 
-bool pmGrowthCurveGenerate (pmReadout *readout, pmPSF *psf, bool ignore, psImageMaskType maskVal, psImageMaskType mark);
 pmPSF *pmPSFBuildSimple (char *typeName, float sxx, float syy, float sxy, ...);
 

branches/eam_branches/ipp-20110404/psModules/src/objects/pmPSF_IO.c

@@ -65 +65 @@
 bool pmPSFmodelRead_ApTrend (pmPSF *psf, pmFPAfile *file);
 bool pmPSFmodelWrite_ApTrend (pmFPAfile *file, pmPSF *psf);
+
+bool pmPSFmodelRead_GrowthCurve (pmPSF *psf, pmFPAfile *file);
+bool pmPSFmodelWrite_GrowthCurve (pmFPAfile *file, pmPSF *psf);
 
 bool pmPSFmodelCheckDataStatusForView (const pmFPAview *view, const pmFPAfile *file)
@@ -600 +603 @@
     if (!pmPSFmodelWrite_ApTrend(file, psf)) {
         psError(psErrorCodeLast(), false, "Unable to write PSF ApTrend");
+        return false;
+    }
+
+    if (!pmPSFmodelWrite_GrowthCurve(file, psf)) {
+        psError(psErrorCodeLast(), false, "Unable to write PSF Growth Curve");
         return false;
     }
@@ -1104 +1112 @@
     }
 
+    if (!pmPSFmodelRead_GrowthCurve(psf, file)) {
+        psError(psErrorCodeLast(), false, "Unable to read PSF Growth Curve");
+        return false;
+    }
+
     psMetadataAdd (chipAnalysis, PS_LIST_TAIL, "PSPHOT.PSF",     PS_DATA_UNKNOWN,  "psphot psf", psf);
     psFree (psf);
@@ -1230 +1243 @@
 }
 
+// write aperture trend to a FITS table
+bool pmPSFmodelWrite_GrowthCurve (pmFPAfile *file, pmPSF *psf) {
+
+    pmGrowthCurve *growth = psf->growth;
+    if (growth == NULL) { 
+        psWarning ("no PSF Growth Curve to write out, skipping");
+        return true; 
+    }
+
+    // we need to write a header for the table,
+    psMetadata *header = psMetadataAlloc();
+
+    psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_MIN_RAD", 0, "", growth->radius->data.F32[0]);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_MAX_RAD", 0, "", growth->maxRadius);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_REF_RAD", 0, "", growth->refRadius);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_AP_LOSS", 0, "", growth->apLoss);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_AP_REF",  0, "", growth->apRef);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_FIT_MAG", 0, "", growth->fitMag);
+
+    // build a FITS table of the ApTrend (only 1)
+    psArray *table = psArrayAllocEmpty (100);
+    for (int i = 0; i < growth->apMag->n; i++) {
+        psMetadata *row = psMetadataAlloc ();
+        psMetadataAddF32 (row, PS_LIST_TAIL, "RADIUS", 0, "", growth->radius->data.F32[i]);
+        psMetadataAddF32 (row, PS_LIST_TAIL, "AP_MAG", 0, "", growth->apMag->data.F32[i]);
+        psArrayAdd (table, 100, row);
+        psFree (row);
+    }
+
+    // write an empty FITS segment if we have no PSF information
+    if (table->n == 0) {
+        psError(PM_ERR_PROG, true, "No PSF data to write.");
+        psFree(table);
+        psFree(header);
+        return false;
+    } 
+
+    psTrace ("pmFPAfile", 5, "writing psf Growth Curve data %s\n", "GROWTH_CURVE");
+    if (!psFitsWriteTable(file->fits, header, table, "GROWTH_CURVE")) {
+        psError(psErrorCodeLast(), false, "Error writing psf table data %s\n", "GROWTH_CURVE");
+        psFree(table);
+        psFree(header);
+        return false;
+    }
+
+    psFree (table);
+    psFree (header);
+    return true;
+}
+
+// read aperture trend to a FITS table
+bool pmPSFmodelRead_GrowthCurve (pmPSF *psf, pmFPAfile *file) {
+
+    bool status;
+
+    // move fits pointer to AP_TREND section
+    // advance to the table data extension
+    if (!psFitsMoveExtNameClean (file->fits, "GROWTH_CURVE")) {
+        psWarning ("no Growth Curve data in PSF file, skipping");
+        return true;
+    }
+
+    psMetadata *header = psFitsReadHeader (NULL, file->fits);
+    if (!header) {
+        psError(psErrorCodeLast(), false, "Unable to read GROWTH_CURVE header.");
+        return false;
+    }
+        
+    // read the raw table data
+    psArray *table = psFitsReadTable (file->fits);
+    if (!table) {
+        psError(psErrorCodeLast(), false, "Unable to read GROWTH_CURVE table.");
+        psFree(header);
+        return false;
+    }
+
+    float minRadius = psMetadataLookupF32 (&status, header, "GROWTH_MIN_RAD"); if (!status) return false;
+    float maxRadius = psMetadataLookupF32 (&status, header, "GROWTH_MAX_RAD"); if (!status) return false;
+    float refRadius = psMetadataLookupF32 (&status, header, "GROWTH_REF_RAD"); if (!status) return false;
+
+    psf->growth = pmGrowthCurveAlloc(minRadius, maxRadius, refRadius);
+
+    psf->growth->apLoss = psMetadataLookupF32 (&status, header, "GROWTH_AP_LOSS"); if (!status) return false;
+    psf->growth->apRef  = psMetadataLookupF32 (&status, header, "GROWTH_AP_REF"); if (!status) return false;
+    psf->growth->fitMag = psMetadataLookupF32 (&status, header, "GROWTH_FIT_MAG"); if (!status) return false;
+
+    // fill in the matching psf->params entries
+    for (int i = 0; i < table->n; i++) {
+        psMetadata *row = table->data[i];
+        psf->growth->apMag->data.F32[i] = psMetadataLookupF32 (&status, row, "AP_MAG"); if (!status) return false;
+    }
+
+    psFree (header);
+    psFree (table);
+    return true;
+}
+
 bool pmPSFmodelReadPSFClump (psMetadata *analysis, psMetadata *header) {
 

branches/eam_branches/ipp-20110404/psModules/src/objects/pmSource.h

@@ -82 +82 @@
     psArray *blends;                    ///< collection of sources thought to be confused with object
     float psfMag;                       ///< calculated from flux in modelPSF
+    float psfMagErr;                    ///< error in psfMag
     float psfFlux;                      ///< calculated from flux in modelPSF
-    float psfFluxErr;                   ///< calculated from flux in modelPSF
-    float extMag;                       ///< calculated from flux in modelEXT
-    float errMag;                       ///< error in psfMag OR extMag (depending on type)
+    float psfFluxErr;                   ///< error in psfFlux
+    float extMag;                       ///< calculated from flux in modelEXT -- NOTE this is not actually used
     float apMag;                        ///< apMag corresponding to psfMag or extMag (depending on type)
     float apMagRaw;                     ///< raw mag in given aperture
@@ -254 +254 @@
     float sigma,      ///< size of Gaussian window function (<= 0.0 -> skip window)
     float minSN,              ///< minimum pixel significance
+    float minKronRadius,      ///< minimum pixel significance
     psImageMaskType maskVal
 );
@@ -273 +274 @@
   float xGuess, float yGuess);
 
+float pmSourceMinKronRadius(psArray *sources, float PSF_SN_LIM);
+
 pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source);
 

branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c

@@ -116 +116 @@
         psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->psfMagErr);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental flux (counts)",         source->psfFlux);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux",             source->psfFluxErr);
@@ -161 +161 @@
         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kinner);
         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kouter);
+
+        // Do NOT write these : not consistent with the definition of PS1_V3 in Ohana/src/libautocode/dev/cmf-ps1-v3.d
         // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_FLUX",   PS_DATA_F32, "Kron Flux (in 1.0 R1)",                      moments.KronCore);
-        // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_ERROR",  PS_DATA_F32, "Kron Error (in 1.0 R1)",                     moments.KronCoreErr);
+        // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_ERROR",  PS_DATA_F32, "Kron Error (in 1.0 R1)",                     moments.KronCoreErr);
         psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
         psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                      source->mode2);
@@ -270 +272 @@
         // XXX use these to determine PAR[PM_PAR_I0]?
         source->psfMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG");
-        source->errMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
+        source->psfMagErr    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
         source->apMag     = psMetadataLookupF32 (&status, row, "AP_MAG");
 
         // XXX this scaling is incorrect: does not include the 2 \pi AREA factor
         PAR[PM_PAR_I0]    = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN;
-        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN;
+        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
 
         pmPSF_AxesToModel (PAR, axes);
@@ -292 +294 @@
 
         // we no longer sort by S/N, only flux
-        // if (isfinite (source->errMag) && (source->errMag > 0.0)) {
-        //   source->peak->SN = 1.0 / source->errMag;
+        // if (isfinite (source->psfMagErr) && (source->psfMagErr > 0.0)) {
+        //   source->peak->SN = 1.0 / source->psfMagErr;
         // } else {
         //   source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N

branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourceMoments.c

@@ -66 +66 @@
 // if mode & EXTERNAL or mode2 & MATCHED, do not re-calculate the centroid (use peak as centroid)
 
-bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal)
+bool pmSourceMoments(pmSource *source, float radius, float sigma, float minSN, float minKronRadius, psImageMaskType maskVal)
 {
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -74 +74 @@
 
     // this function assumes the sky has been well-subtracted for the image
-    psF32 sky = 0.0;
+    float sky = 0.0;
 
     if (source->moments == NULL) {
@@ -80 +80 @@
     }
 
-    psF32 Sum = 0.0;
-    psF32 Var = 0.0;
-    psF32 SumCore = 0.0;
-    psF32 VarCore = 0.0;
-    psF32 R2 = PS_SQR(radius);
-    psF32 minSN2 = PS_SQR(minSN);
-    psF32 rsigma2 = 0.5 / PS_SQR(sigma);
+    float Sum = 0.0;
+    float Var = 0.0;
+    float SumCore = 0.0;
+    float VarCore = 0.0;
+    float R2 = PS_SQR(radius);
+    float minSN2 = PS_SQR(minSN);
+    float rsigma2 = 0.5 / PS_SQR(sigma);
 
     // a note about coordinates: coordinates of objects throughout psphot refer to the primary
@@ -109 +109 @@
     // Xn  = SUM (x - xc)^n * (z - sky)
 
-    psF32 RFW = 0.0;
-    psF32 RHW = 0.0;
-
-    psF32 RF = 0.0;
-    psF32 RH = 0.0;
-    psF32 RS = 0.0;
-    psF32 XX = 0.0;
-    psF32 XY = 0.0;
-    psF32 YY = 0.0;
-    psF32 XXX = 0.0;
-    psF32 XXY = 0.0;
-    psF32 XYY = 0.0;
-    psF32 YYY = 0.0;
-    psF32 XXXX = 0.0;
-    psF32 XXXY = 0.0;
-    psF32 XXYY = 0.0;
-    psF32 XYYY = 0.0;
-    psF32 YYYY = 0.0;
+    float RFW = 0.0;
+    float RHW = 0.0;
+
+    float RF = 0.0;
+    float RH = 0.0;
+    float RS = 0.0;
+    float XX = 0.0;
+    float XY = 0.0;
+    float YY = 0.0;
+    float XXX = 0.0;
+    float XXY = 0.0;
+    float XYY = 0.0;
+    float YYY = 0.0;
+    float XXXX = 0.0;
+    float XXXY = 0.0;
+    float XXYY = 0.0;
+    float XYYY = 0.0;
+    float YYYY = 0.0;
 
     Sum = 0.0;  // the second pass may include slightly different pixels, re-determine Sum
@@ -140 +140 @@
     // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
     // xCM, yCM from pixel coords to pixel index here.
-    psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
-    psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
+    float xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
+    float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
 
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
 
-        psF32 yDiff = row - yCM;
+        float yDiff = row - yCM;
         if (fabs(yDiff) > radius) continue;
 
-        psF32 *vPix = source->pixels->data.F32[row];
-        psF32 *vWgt = source->variance->data.F32[row];
+        float *vPix = source->pixels->data.F32[row];
+        float *vWgt = source->variance->data.F32[row];
 
         psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
@@ -164 +164 @@
             if (isnan(*vPix)) continue;
 
-            psF32 xDiff = col - xCM;
+            float xDiff = col - xCM;
             if (fabs(xDiff) > radius) continue;
 
             // radius is just a function of (xDiff, yDiff)
-            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
             if (r2 > R2) continue;
 
-            psF32 fDiff = *vPix - sky;
-            psF32 pDiff = fDiff;
-            psF32 wDiff = *vWgt;
+            float fDiff = *vPix - sky;
+            float pDiff = fDiff;
+            float wDiff = *vWgt;
 
             // skip pixels below specified significance level.  this is allowed, but should be
@@ -184 +184 @@
             // weighting over weights the sky for faint sources
             if (sigma > 0.0) {
-                psF32 z = r2 * rsigma2;
+                float z = r2 * rsigma2;
                 assert (z >= 0.0);
-                psF32 weight  = exp(-z);
+                float weight  = exp(-z);
 
                 wDiff *= weight;
@@ -195 +195 @@
 
             // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
-            psF32 r = sqrt(r2);
-            psF32 rf = r * fDiff;
-            psF32 rh = sqrt(r) * fDiff;
-            psF32 rs = fDiff;
-
-            psF32 rfw = r * pDiff;
-            psF32 rhw = sqrt(r) * pDiff;
-
-            psF32 x = xDiff * pDiff;
-            psF32 y = yDiff * pDiff;
-
-            psF32 xx = xDiff * x;
-            psF32 xy = xDiff * y;
-            psF32 yy = yDiff * y;
-
-            psF32 xxx = xDiff * xx / r;
-            psF32 xxy = xDiff * xy / r;
-            psF32 xyy = xDiff * yy / r;
-            psF32 yyy = yDiff * yy / r;
-
-            psF32 xxxx = xDiff * xxx / r2;
-            psF32 xxxy = xDiff * xxy / r2;
-            psF32 xxyy = xDiff * xyy / r2;
-            psF32 xyyy = xDiff * yyy / r2;
-            psF32 yyyy = yDiff * yyy / r2;
+            float r = sqrt(r2);
+            float rf = r * fDiff;
+            float rh = sqrt(r) * fDiff;
+            float rs = fDiff;
+
+            float rfw = r * pDiff;
+            float rhw = sqrt(r) * pDiff;
+
+            float x = xDiff * pDiff;
+            float y = yDiff * pDiff;
+
+            float xx = xDiff * x;
+            float xy = xDiff * y;
+            float yy = yDiff * y;
+
+            float xxx = xDiff * xx / r;
+            float xxy = xDiff * xy / r;
+            float xyy = xDiff * yy / r;
+            float yyy = yDiff * yy / r;
+
+            float xxxx = xDiff * xxx / r2;
+            float xxxy = xDiff * xxy / r2;
+            float xxyy = xDiff * xyy / r2;
+            float xyyy = xDiff * yyy / r2;
+            float yyyy = yDiff * yyy / r2;
 
             RF  += rf;
@@ -245 +245 @@
     }
 
-    // if Mrf (first radial moment) is << sigma, we are getting into low-significance
-    // territory.  saturate at 0.75*sigma.  conversely, if Mrf is > radius, we are clearly
-    // making an error.  saturate at radius.
-    source->moments->Mrf = MIN(radius, MAX(0.75*sigma, RF/RS));
-    source->moments->Mrh = MIN(radius, MAX(0.75*sigma, RH/RS));
+    source->moments->Mrf = RF/RS;
+    source->moments->Mrh = RH/RS;
 
     source->moments->Mxx = XX/Sum;
@@ -266 +263 @@
     source->moments->Myyyy = YYYY/Sum;
 
-    // XXX TEST:
-    // source->moments->KronFinner = RFW/Sum;
-    // source->moments->KronFouter = sigma;
-
-    // Calculate the Kron magnitude (make this block optional?)
-    float radKinner = 1.0*source->moments->Mrf;
-    float radKron   = 2.5*source->moments->Mrf;
-    float radKouter = 4.0*source->moments->Mrf;
+    // if Mrf (first radial moment) is very small, we are getting into low-significance
+    // territory.  saturate at minKronRadius.  conversely, if Mrf is > radius, we are clearly
+    // making an error.  saturate at radius.
+    float kronRefRadius = MIN(radius, MAX(minKronRadius, source->moments->Mrf));
+
+    float radKinner = 1.0*kronRefRadius;
+    float radKron   = 2.5*kronRefRadius;
+    float radKouter = 4.0*kronRefRadius;
 
     int nKronPix = 0;
@@ -285 +282 @@
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
 
-        psF32 yDiff = row - yCM;
+        float yDiff = row - yCM;
         if (fabs(yDiff) > radKouter) continue;
 
-        psF32 *vPix = source->pixels->data.F32[row];
-        psF32 *vWgt = source->variance->data.F32[row];
+        float *vPix = source->pixels->data.F32[row];
+        float *vWgt = source->variance->data.F32[row];
 
         psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
@@ -304 +301 @@
             if (isnan(*vPix)) continue;
 
-            psF32 xDiff = col - xCM;
+            float xDiff = col - xCM;
             if (fabs(xDiff) > radKouter) continue;
 
             // radKron is just a function of (xDiff, yDiff)
-            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
-
-            psF32 pDiff = *vPix - sky;
-            psF32 wDiff = *vWgt;
+            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+
+            float pDiff = *vPix - sky;
+            float wDiff = *vWgt;
 
             // skip pixels below specified significance level.  this is allowed, but should be
@@ -318 +315 @@
             if (PS_SQR(pDiff) < minSN2*wDiff) continue;
 
-            psF32 r  = sqrt(r2);
+            float r  = sqrt(r2);
             if (r < radKron) {
                 Sum += pDiff;
@@ -363 +360 @@
 }
 
-bool pmSourceMomentsGetCentroid(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal, float xGuess, float yGuess) { 
+bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) { 
 
     // First Pass: calculate the first moments (these are subtracted from the coordinates below)
@@ -370 +367 @@
     // .. etc
 
-    psF32 sky = 0.0;
-
-    psF32 peakPixel = -PS_MAX_F32;
+    float sky = 0.0;
+
+    float peakPixel = -PS_MAX_F32;
     psS32 numPixels = 0;
-    psF32 Sum = 0.0;
-    psF32 Var = 0.0;
-    psF32 X1 = 0.0;
-    psF32 Y1 = 0.0;
-    psF32 R2 = PS_SQR(radius);
-    psF32 minSN2 = PS_SQR(minSN);
-    psF32 rsigma2 = 0.5 / PS_SQR(sigma);
+    float Sum = 0.0;
+    float Var = 0.0;
+    float X1 = 0.0;
+    float Y1 = 0.0;
+    float R2 = PS_SQR(radius);
+    float minSN2 = PS_SQR(minSN);
+    float rsigma2 = 0.5 / PS_SQR(sigma);
 
     float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
@@ -386 +383 @@
 
     // we are guaranteed to have a valid pixel and variance at this location (right? right?)
-    // psF32 weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
+    // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
     // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
@@ -398 +395 @@
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
 
-        psF32 yDiff = row + 0.5 - yPeak;
+        float yDiff = row + 0.5 - yPeak;
         if (fabs(yDiff) > radius) continue;
 
-        psF32 *vPix = source->pixels->data.F32[row];
-        psF32 *vWgt = source->variance->data.F32[row];
+        float *vPix = source->pixels->data.F32[row];
+        float *vWgt = source->variance->data.F32[row];
 
         psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
@@ -417 +414 @@
             if (isnan(*vPix)) continue;
 
-            psF32 xDiff = col + 0.5 - xPeak;
+            float xDiff = col + 0.5 - xPeak;
             if (fabs(xDiff) > radius) continue;
 
             // radius is just a function of (xDiff, yDiff)
-            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
             if (r2 > R2) continue;
 
-            psF32 pDiff = *vPix - sky;
-            psF32 wDiff = *vWgt;
+            float pDiff = *vPix - sky;
+            float wDiff = *vWgt;
 
             // skip pixels below specified significance level.  for a PSFs, this
@@ -438 +435 @@
             // weighting over weights the sky for faint sources
             if (sigma > 0.0) {
-                psF32 z  = r2*rsigma2;
+                float z  = r2*rsigma2;
                 assert (z >= 0.0);
-                psF32 weight  = exp(-z);
+                float weight  = exp(-z);
 
                 wDiff *= weight;
@@ -449 +446 @@
             Sum += pDiff;
 
-            psF32 xWght = xDiff * pDiff;
-            psF32 yWght = yDiff * pDiff;
+            float xWght = xDiff * pDiff;
+            float yWght = yDiff * pDiff;
 
             X1  += xWght;
@@ -507 +504 @@
     return true;
 }
+
+float pmSourceMinKronRadius(psArray *sources, float PSF_SN_LIM) {
+
+    psVector *radii = psVectorAllocEmpty(100, PS_TYPE_F32);
+
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *src = sources->data[i]; // Source of interest
+        if (!src || !src->moments) {
+            continue;
+        }
+
+        if (src->mode & PM_SOURCE_MODE_BLEND) {
+            continue;
+        }
+
+        if (!src->moments->nPixels) continue;
+
+        if (src->moments->SN < PSF_SN_LIM) continue;
+
+        // XXX put in Mxx,Myy cut based on clump location
+
+        psVectorAppend(radii, src->moments->Mrf);
+    }
+
+    // find the peak in this image
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+    if (!psVectorStats (stats, radii, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
+        psFree(stats);
+        return NAN;
+    }
+
+    float minRadius = stats->sampleMedian;
+
+    psFree(radii);
+    psFree(stats);
+    return minRadius;
+}
+

branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourcePhotometry.c

@@ -82 +82 @@
    - apMag  : only if S/N > AP_MIN_SN
    : is optionally corrected for curve-of-growth if:
-   - the source is a STAR (PSF)
    - the option is selected (mode & PM_SOURCE_PHOT_GROWTH)
    - psfMag : all sources with non-NULL modelPSF
    : is optionally corrected for aperture residual if:
-   - the source is a STAR (PSF)
    - the option is selected (mode & PM_SOURCE_PHOT_APCORR)
-
    - extMag : all sources with non-NULL modelEXT
 **/
@@ -100 +97 @@
 
     int status = false;
-    bool isPSF;
     float x, y;
     float SN;
-    pmModel *model;
 
     source->psfMag    = NAN;
     source->extMag    = NAN;
-    source->errMag    = NAN;
+    source->psfMagErr = NAN;
     source->apMag     = NAN;
     source->apMagRaw  = NAN;
@@ -113 +108 @@
     source->apFluxErr = NAN;
 
-    // we must have a valid model
-    // XXX allow aperture magnitudes for sources without a model
-    model = pmSourceGetModel (&isPSF, source);
-    if (model == NULL) {
-        psTrace ("psModules.objects", 3, "fail mag : no valid model");
+    // XXXXXX review:
+    // Select the 'best' model -- this is used for PSF_QF,_PERFECT & ???. isPSF is true if this
+    // object is a PSF (not extended).  We must have a valid model.  XXX NOTE: allow aperture
+    // magnitudes for sources without a model
+
+    // select the psf model
+    pmModel *modelPSF = source->modelPSF;
+    if (modelPSF == NULL) {
+        psTrace ("psModules.objects", 3, "fail mag : no valid PSF model");
         return false;
     }
 
-    // XXX handle negative flux, low-significance
-    if (model->dparams->data.F32[PM_PAR_I0] > 0) {
-        SN = fabs(model->params->data.F32[PM_PAR_I0] / model->dparams->data.F32[PM_PAR_I0]);
-        source->errMag = 1.0 / SN;
+    // get the error on the PSF model magnitude
+    if (modelPSF->dparams->data.F32[PM_PAR_I0] > 0) {
+        SN = fabs(modelPSF->params->data.F32[PM_PAR_I0] / modelPSF->dparams->data.F32[PM_PAR_I0]);
+        source->psfMagErr = 1.0 / SN;
     } else {
         SN = NAN;
-        source->errMag = NAN;
-    }
-    x = model->params->data.F32[PM_PAR_XPOS];
-    y = model->params->data.F32[PM_PAR_YPOS];
+        source->psfMagErr = NAN;
+    }
+    // the source position is used to recenter the aperture for ap photometry
+    x = modelPSF->params->data.F32[PM_PAR_XPOS];
+    y = modelPSF->params->data.F32[PM_PAR_YPOS];
 
     // measure PSF model photometry
-    // XXX TEST: do not use flux scale
-    // XXX NOTE: turn this back on?
-    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);
-        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->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0];
-        source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0];
-        source->psfMag = -2.5*log10(source->psfFlux);
-    } else {
-        status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, source->modelPSF);
-        source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]);
-    }
+    status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, modelPSF);
+    source->psfFluxErr = source->psfFlux * source->psfMagErr;
+
+# if (0)
+    // XXX NOTE: old code to use the flux scale.  test & turn this back on?  if so, need to save with psf model
+    // 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);
+    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->psfFlux = fluxScale * modelPSF->params->data.F32[PM_PAR_I0];
+    source->psfFluxErr = fluxScale * modelPSF->dparams->data.F32[PM_PAR_I0];
+    source->psfMag = -2.5*log10(source->psfFlux);
+# endif
 
     if (mode == PM_SOURCE_PHOT_PSFONLY) {
@@ -152 +151 @@
     }
 
+    // get the EXT model photometry (all EXT models)
     // if we have a collection of model fits, check if one of them is a pointer to modelEXT
     if (source->modelFits) {
@@ -172 +172 @@
     }
 
-    // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS
-    // XXX add a flag for "ap_mag is corrected?"
-    if ((mode & PM_SOURCE_PHOT_APCORR) && isPSF && psf && psf->ApTrend) {
+    // Correct psfMag to match aperture magnitude system (NOTE : Growth curve is already applied to ApTrend)
+    if ((mode & PM_SOURCE_PHOT_APCORR) && 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);
@@ -182 +181 @@
     }
 
-    // measure the contribution of included pixels
+    // measure the contribution of included pixels to the PSF model fit
     if (mode & PM_SOURCE_PHOT_WEIGHT) {
-        pmSourcePixelWeight (source, model, source->maskObj, maskVal, radius);
+        pmSourcePixelWeight (source, modelPSF, source->maskObj, maskVal, radius);
     }
 
@@ -218 +217 @@
         y += dy;
 
-        if (!psImageShiftMask(&flux, &mask, source->pixels, source->maskObj, maskVal, dx, dy,
-                              NAN, 0xff, PS_INTERPOLATE_BIQUADRATIC)) {
+        // if (!psImageShiftMask(&flux, &mask, source->pixels, source->maskObj, maskVal, dx, dy, NAN, 0xff, PS_INTERPOLATE_LANCZOS2)) {
+        // if (!psImageShiftMask(&flux, &mask, source->pixels, source->maskObj, maskVal, dx, dy, NAN, 0xff, PS_INTERPOLATE_BIQUADRATIC)) {
+        if (!psImageShiftMask(&flux, &mask, source->pixels, source->maskObj, maskVal, dx, dy, NAN, 0xff, PS_INTERPOLATE_BILINEAR)) {
             // Not much we can do about it
             psErrorClear();
@@ -233 +233 @@
 
     // measure object aperture photometry
-    status = pmSourcePhotometryAperSource (source, model, flux, variance, mask, maskVal);
+    status = pmSourcePhotometryAperSource (source, modelPSF, flux, variance, mask, maskVal);
     if (!status) {
         psTrace ("psModules.objects", 3, "fail mag : bad Ap Mag");
@@ -242 +242 @@
     // detection limits (esp near bright neighbors)
     source->apMag = source->apMagRaw;
-    if (isfinite (source->apMag) && isPSF && psf) {
+    if (isfinite (source->apMag) && psf) {
         if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) {
-            source->apMag = source->apMagRaw + pmGrowthCurveCorrect (psf->growth, source->apRadius);
-            // XXX correct the apFlux?
+            float apOffset = pmGrowthCurveCorrect (psf->growth, source->apRadius);
+            source->apMag = source->apMagRaw + apOffset;
+            source->apFlux *= pow(10.0, -0.4*apOffset);
+            source->apFluxErr *= pow(10.0, -0.4*apOffset);
         }
     }

branches/eam_branches/ipp-20110404/psModules/src/objects/pmSourceUtils.h

@@ -38 +38 @@
 
 /// @}
-# endif /* PM_MODEL_CLASS_H */
+# endif /* PM_SOURCE_UTILS_H */

branches/eam_branches/ipp-20110404/psModules/src/psmodules.h

@@ -148 +148 @@
 #include <pmModelUtils.h>
 #include <pmSourcePhotometry.h>
+#include <pmGrowthCurveGenerate.h>
 #include <pmSourceVisual.h>
 #include <pmSourceMatch.h>