IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25476


Ignore:
Timestamp:
Sep 21, 2009, 8:37:01 PM (17 years ago)
Author:
eugene
Message:

re-organize pmPSFtry (split into clearer files); modify 2D psf variation analysis to use psf-aper mags as the metric

Location:
branches/eam_branches/20090715/psModules/src
Files:
2 added
9 edited
1 moved

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20090715/psModules/src/extras/pmVisual.c

    r23989 r25476  
    2222#include "pmSubtractionStamps.h"
    2323#include "pmTrend2D.h"
     24#include "pmPSF.h"
     25#include "pmPSFtry.h"
    2426#include "pmFPAExtent.h"
    2527

  • branches/eam_branches/20090715/psModules/src/objects/Makefile.am

    r25407 r25476  
    5050        pmPSF_IO.c \
    5151        pmPSFtry.c \
     52        pmPSFtryModel.c \
     53        pmPSFtryFitEXT.c \
     54        pmPSFtryMakePSF.c \
     55        pmPSFtryFitPSF.c \
     56        pmPSFtryMetric.c \
    5257        pmTrend2D.c \
    5358        pmGrowthCurveGenerate.c \

  • branches/eam_branches/20090715/psModules/src/objects/pmPSF.h

    r24206 r25476  
    3838 *
    3939 */
    40 typedef struct
    41 {
     40typedef struct {
    4241    pmModelType type;                   ///< PSF Model in use
    4342    psArray *params;                    ///< Model parameters (psPolynomial2D)
    … …  
    6564    pmGrowthCurve *growth;              ///< apMag vs Radius
    6665    pmResiduals *residuals;             ///< normalized residual image (no spatial variation)
    67 }
    68 pmPSF;
     66} pmPSF;
    6967
    7068typedef struct {

  • branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.c

    r25455 r25476  
    6363    psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree);
    6464
    65     test->psf       = pmPSFAlloc (options);
     65    test->psf       = NULL;
    6666    test->metric    = psVectorAlloc (sources->n, PS_TYPE_F32);
    6767    test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F32);
    … …  
    8787    PS_ASSERT_PTR(ptr, false);
    8888    return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmPSFtryFree);
    89 }
    90 
    91 
    92 // build a pmPSFtry for the given model:
    93 // - fit each source with the free-floating model
    94 // - construct the pmPSF from the collection of models
    95 // - fit each source with the PSF-parameter models
    96 // - measure the pmPSF quality metric (dApResid)
    97 
    98 // sources used in for pmPSFtry may be masked by the analysis
    99 // mask values indicate the reason the source was rejected:
    100 
    101 // generate a pmPSFtry with a copy of the test PSF sources
    102 pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal)
    103 {
    104     bool status;
    105     int Next = 0;
    106     int Npsf = 0;
    107 
    108     // validate the requested model name
    109     options->type = pmModelClassGetType (modelName);
    110     if (options->type == -1) {
    111         psError (PS_ERR_UNKNOWN, true, "invalid model name %s", modelName);
    112         return NULL;
    113     }
    114 
    115     pmPSFtry *psfTry = pmPSFtryAlloc (sources, options);
    116     if (psfTry == NULL) {
    117         psError (PS_ERR_UNKNOWN, false, "failed to allocate psf model");
    118         return NULL;
    119     }
    120 
    121     // maskVal is used to test for rejected pixels, and must include markVal
    122     maskVal |= markVal;
    123 
    124     // stage 1:  fit an EXT model to all candidates PSF sources -- this is independent of the modeled 2D variations in the PSF
    125     if (!pmPSFtryFitEXT(psfTry, options, maskVal, markVal)) {
    126         psError(PS_ERR_UNKNOWN, false, "failed to fit EXT models to sources for psf model");
    127         psFree(psfTry);
    128         return NULL;
    129     }     
    130 
    131     for (int i = 0; i < Norder; i++) {
    132         // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation
    133         if (!pmPSFFromPSFtry (psfTry, Nx, Ny)) {
    134             psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources");
    135             psFree(psfTry);
    136             return NULL;
    137         }
    138 
    139         // stage 3: refit with fixed shape parameters, measure pmPSFtryMetric
    140         if (!pmPSFtryFitPSF (psfTry, Nx, Ny)) {
    141             psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources");
    142             psFree(psfTry);
    143             return NULL;
    144         }
    145     }
    146     // XXXXX this is probably not used any more.  Are the chisq of the fits so bad? can we
    147     // fix them by softening the errors on the brightest pixels?
    148 
    149     // measure the chi-square trend as a function of flux (PAR[PM_PAR_I0])
    150     // this should be linear for Poisson errors and quadratic for constant sky errors
    151     psVector *flux  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    152     psVector *chisq = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    153     psVector *mask  = psVectorAlloc (psfTry->sources->n, PS_TYPE_VECTOR_MASK);
    154 
    155     // generate the x and y vectors, and mask missing models
    156     for (int i = 0; i < psfTry->sources->n; i++) {
    157         pmSource *source = psfTry->sources->data[i];
    158         if (source->modelPSF == NULL) {
    159             flux->data.F32[i] = 0.0;
    160             chisq->data.F32[i] = 0.0;
    161             mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
    162         } else {
    163             flux->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_I0];
    164             chisq->data.F32[i] = source->modelPSF->chisq / source->modelPSF->nDOF;
    165             mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
    166         }
    167     }
    168 
    169     // use 3hi/3lo sigma clipping on the chisq fit
    170     psStats *stats = options->stats;
    171 
    172     // linear clipped fit of chisq trend vs flux
    173     if (options->chiFluxTrend) {
    174         bool result = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, options->stats,
    175                                                   mask, 0xff, chisq, NULL, flux);
    176         psStatsOptions meanStat = psStatsMeanOption(options->stats->options); // Statistic for mean
    177         psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev
    178 
    179         psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n",
    180                   psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat));
    181 
    182         psFree(flux);
    183         psFree(mask);
    184         psFree(chisq);
    185 
    186         if (!result) {
    187             psError(PS_ERR_UNKNOWN, false, "Failed to fit psf->ChiTrend");
    188             psFree(psfTry);
    189             return NULL;
    190         }
    191     }
    192 
    193     for (int i = 0; i < psfTry->psf->ChiTrend->nX + 1; i++) {
    194         psLogMsg ("pmPSFtry", 4, "chisq vs flux fit term %d: %f +/- %f\n", i,
    195                   psfTry->psf->ChiTrend->coeff[i]*pow(10000, i),
    196                   psfTry->psf->ChiTrend->coeffErr[i]*pow(10000,i));
    197     }
    198 
    199     // XXX this function wants aperture radius for pmSourcePhotometry
    200     if (!pmPSFtryMetric (psfTry, options)) {
    201         psError(PS_ERR_UNKNOWN, false, "Attempt to fit PSF with model %s failed.", modelName);
    202         psFree (psfTry);
    203         return NULL;
    204     }
    205 
    206     psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n",
    207               modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias);
    208 
    209     return (psfTry);
    210 }
    211 
    212 bool pmPSFtryMetric (pmPSFtry *psfTry, pmPSFOptions *options)
    213 {
    214     PS_ASSERT_PTR_NON_NULL(psfTry, false);
    215     PS_ASSERT_PTR_NON_NULL(options, false);
    216     PS_ASSERT_PTR_NON_NULL(psfTry->sources, false);
    217 
    218     float RADIUS = options->radius;
    219 
    220     // the measured (aperture - fit) magnitudes (dA == psfTry->metric)
    221     //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
    222     //     dA = dAo + dsky/flux
    223     //   where flux is the flux of the star
    224     // we fit this trend to find the infinite flux aperture correction (dAo),
    225     //   the nominal sky bias (dsky), and the error on dAo
    226     // the values of dA are contaminated by stars with close neighbors in the aperture
    227     //   we use an outlier rejection to avoid this bias
    228 
    229     // r2rflux = radius^2 * ten(0.4*fitMag);
    230     psVector *r2rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    231 
    232     for (int i = 0; i < psfTry->sources->n; i++) {
    233         if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL)
    234             continue;
    235         r2rflux->data.F32[i] = PS_SQR(RADIUS) * pow(10.0, 0.4*psfTry->fitMag->data.F32[i]);
    236     }
    237 
    238     // XXX test dump of aperture residual data
    239     if (psTraceGetLevel("psModules.objects") >= 5) {
    240         FILE *f = fopen ("apresid.dat", "w");
    241         for (int i = 0; i < psfTry->sources->n; i++) {
    242             int keep = (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL);
    243 
    244             pmSource *source = psfTry->sources->data[i];
    245             float x = source->peak->x;
    246             float y = source->peak->y;
    247 
    248             fprintf (f, "%d  %d, %f %f %f  %f %f %f \n",
    249                      i, keep, x, y,
    250                      psfTry->fitMag->data.F32[i],
    251                      r2rflux->data.F32[i],
    252                      psfTry->metric->data.F32[i],
    253                      psfTry->metricErr->data.F32[i]);
    254         }
    255         fclose (f);
    256     }
    257 
    258     // This analysis of the apResid statistics is only approximate.  The fitted magnitudes
    259     // measured at this point (in the PSF fit) use Poisson errors, and are thus biased as a
    260     // function of magnitude.  We re-measure the apResid statistics later in psphot using the
    261     // linear, constant-error fitting.  Do not reject outliers with excessive vigor here.
    262 
    263     // fit ApTrend only to r2rflux, ignore x,y,flux variations for now
    264     // linear clipped fit of ApResid to r2rflux
    265     psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
    266     poly->coeffMask[1] = PS_POLY_MASK_SET; // fit only a constant offset (no SKYBIAS)
    267 
    268     // XXX replace this with a psVectorStats call?  since we are not fitting the trend
    269     bool result = psVectorClipFitPolynomial1D(poly, options->stats, psfTry->mask, PSFTRY_MASK_ALL,
    270                                               psfTry->metric, psfTry->metricErr, r2rflux);
    271     if (!result) {
    272         psError(PS_ERR_UNKNOWN, false, "Failed to fit clipped poly");
    273 
    274         psFree(poly);
    275         psFree(r2rflux);
    276 
    277         return false;
    278     }
    279     psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev
    280     psLogMsg ("pmPSFtryMetric", 4, "apresid: %f +/- %f; from statistics of %ld psf stars\n", poly->coeff[0],
    281               psStatsGetValue(options->stats, stdevStat), psfTry->sources->n);
    282 
    283     float dSys = psVectorSystematicError (psfTry->metric, psfTry->metricErr, 0.1);
    284     fprintf (stderr, "systematic error: %f\n", dSys);
    285 
    286     int n = 0;
    287     psVector *bright = psVectorAllocEmpty (psfTry->metric->n, PS_TYPE_F32);
    288     psVector *brightErr = psVectorAllocEmpty (psfTry->metric->n, PS_TYPE_F32);
    289     for (int i = 0; i < psfTry->metric->n; i++) {
    290         if (!isfinite(psfTry->metric->data.F32[i])) continue;
    291         if (!isfinite(psfTry->metricErr->data.F32[i])) continue;
    292         if (psfTry->metricErr->data.F32[i] <= 0.0) continue;
    293         if (psfTry->metricErr->data.F32[i] > 0.005) continue;
    294         bright->data.F32[n] = psfTry->metric->data.F32[i];
    295         brightErr->data.F32[n] = psfTry->metricErr->data.F32[i];
    296         n++;
    297     }
    298     bright->n = brightErr->n = n;
    299 
    300     float dSysBright = psVectorSystematicError (bright, brightErr, 0.1);
    301     fprintf (stderr, "bright systematic error: %f\n", dSysBright);
    302     psFree(bright);
    303     psFree(brightErr);
    304 
    305     // XXX test dump of fitted model (dump when tracing?)
    306     if (psTraceGetLevel("psModules.objects") >= 4) {
    307         FILE *f = fopen ("resid.dat", "w");
    308         psVector *apfit = psPolynomial1DEvalVector (poly, r2rflux);
    309         for (int i = 0; i < psfTry->sources->n; i++) {
    310             int keep = (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL);
    311 
    312             pmSource *source = psfTry->sources->data[i];
    313             float x = source->peak->x;
    314             float y = source->peak->y;
    315 
    316             fprintf (f, "%d  %d, %f %f %f  %f %f %f  %f\n",
    317                      i, keep, x, y,
    318                      psfTry->fitMag->data.F32[i],
    319                      r2rflux->data.F32[i],
    320                      psfTry->metric->data.F32[i],
    321                      psfTry->metricErr->data.F32[i],
    322                      apfit->data.F32[i]);
    323         }
    324         fclose (f);
    325         psFree (apfit);
    326     }
    327 
    328     // XXX drop the skyBias value, or include above??
    329     psfTry->psf->skyBias  = poly->coeff[1];
    330     psfTry->psf->ApResid  = poly->coeff[0];
    331     psfTry->psf->dApResid = psStatsGetValue(options->stats, stdevStat);
    332 
    333     psFree (r2rflux);
    334     psFree (poly);
    335 
    336     return true;
    337 }
    338 
    339 /*
    340   (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
    341   (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
    342   (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
    343 */
    344 
    345 bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) {
    346 
    347     // we are doing a robust fit.  after each pass, we drop points which are more deviant than
    348     // three sigma.  mask is currently updated for each pass.
    349 
    350     // The shape parameters (SXX, SXY, SYY) are strongly coupled.  We have to handle them very
    351     // carefully.  First, we convert them to the Ellipse Polarization terms (E0, E1, E2) for
    352     // each source and fit this set of parameters.  These values are less tightly coupled, but
    353     // are still inter-related.  The fitted values do a good job of constraining the major axis
    354     // and the position angle, but the minor axis is weakly measured.  When we apply the PSF
    355     // model to construct a source model, we convert the fitted values of E0,E1,E2 to the shape
    356     // parameters, with the constraint that the minor axis must be greater than a minimum
    357     // threshold.
    358 
    359     // XXX re-read the sextractor manual on handling 'infinitely thin' sources...
    360 
    361     // convert the measured source shape paramters to polarization terms
    362     psVector *e0   = psVectorAlloc (sources->n, PS_TYPE_F32);
    363     psVector *e1   = psVectorAlloc (sources->n, PS_TYPE_F32);
    364     psVector *e2   = psVectorAlloc (sources->n, PS_TYPE_F32);
    365     psVector *mag  = psVectorAlloc (sources->n, PS_TYPE_F32);
    366 
    367     for (int i = 0; i < sources->n; i++) {
    368         // skip any masked sources (failed to fit one of the model steps or get a magnitude)
    369         if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
    370 
    371         pmSource *source = sources->data[i];
    372         assert (source->modelEXT); // all unmasked sources should have modelEXT
    373 
    374         psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32);
    375 
    376         e0->data.F32[i] = pol.e0;
    377         e1->data.F32[i] = pol.e1;
    378         e2->data.F32[i] = pol.e2;
    379 
    380         float flux = source->modelEXT->params->data.F32[PM_PAR_I0];
    381         mag->data.F32[i] = (flux > 0.0) ? -2.5*log(flux) : -100.0;
    382     }
    383 
    384     if (psf->psfTrendMode == PM_TREND_MAP) {
    385         float scatterTotal = 0.0;
    386         float scatterTotalMin = FLT_MAX;
    387         int entryMin = -1;
    388 
    389         psVector *dz = NULL;
    390         psVector *mask = psVectorAlloc (sources->n, PS_TYPE_VECTOR_MASK);
    391 
    392         // check the fit residuals and increase Nx,Ny until the error is minimized
    393         // pmPSFParamTrend increases the number along the longer of x or y
    394         for (int i = 1; i <= PS_MAX (psf->trendNx, psf->trendNy); i++) {
    395 
    396             // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky)
    397             for (int i = 0; i < mask->n; i++) {
    398                 mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    399             }
    400             if (!pmPSFFitShapeParamsMap (psf, i, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) {
    401                 break;
    402             }
    403 
    404             // store the resulting scatterTotal values and the scales, redo the best
    405             if (scatterTotal < scatterTotalMin) {
    406                 scatterTotalMin = scatterTotal;
    407                 entryMin = i;
    408             }
    409         }
    410         if (entryMin == -1) {
    411             psError (PS_ERR_UNKNOWN, false, "failed to find image map for shape params");
    412             return false;
    413         }
    414 
    415         // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky)
    416         for (int i = 0; i < mask->n; i++) {
    417             mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    418         }
    419         if (!pmPSFFitShapeParamsMap (psf, entryMin, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) {
    420             psAbort ("failed pmPSFFitShapeParamsMap on second pass?");
    421         }
    422 
    423         pmTrend2D *trend = psf->params->data[PM_PAR_E0];
    424         psf->trendNx = trend->map->map->numCols;
    425         psf->trendNy = trend->map->map->numRows;
    426 
    427         // copy mask back to srcMask
    428         for (int i = 0; i < mask->n; i++) {
    429             srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    430         }
    431 
    432         psFree (mask);
    433         psFree (dz);
    434     } else {
    435 
    436         // XXX iterate Nx, Ny based on scatter?
    437         // XXX we force the x & y order to be the same
    438         // modify the order to correspond to the actual number of matched stars:
    439         int order = PS_MAX (psf->trendNx, psf->trendNy);
    440         if ((sources->n < 15) && (order >= 3)) order = 2;
    441         if ((sources->n < 11) && (order >= 2)) order = 1;
    442         if ((sources->n <  8) && (order >= 1)) order = 0;
    443         if ((sources->n <  3)) {
    444             psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params");
    445             return false;
    446         }
    447         psf->trendNx = order;
    448         psf->trendNy = order;
    449 
    450         // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
    451         // This way, the parameters masked by one of the fits will be applied to the others
    452         for (int i = 0; i < psf->psfTrendStats->clipIter; i++) {
    453 
    454             psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
    455             psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
    456 
    457             pmTrend2D *trend = NULL;
    458             float mean, stdev;
    459 
    460             // XXX we are using the same stats structure on each pass: do we need to re-init it?
    461             bool status = true;
    462 
    463             trend = psf->params->data[PM_PAR_E0];
    464             status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e0, NULL);
    465             mean = psStatsGetValue (trend->stats, meanOption);
    466             stdev = psStatsGetValue (trend->stats, stdevOption);
    467             psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n);
    468             pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask);
    469 
    470             trend = psf->params->data[PM_PAR_E1];
    471             status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e1, NULL);
    472             mean = psStatsGetValue (trend->stats, meanOption);
    473             stdev = psStatsGetValue (trend->stats, stdevOption);
    474             psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n);
    475             pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask);
    476 
    477             trend = psf->params->data[PM_PAR_E2];
    478             status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e2, NULL);
    479             mean = psStatsGetValue (trend->stats, meanOption);
    480             stdev = psStatsGetValue (trend->stats, stdevOption);
    481             psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n);
    482             pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask);
    483 
    484             if (!status) {
    485                 psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params");
    486                 return false;
    487             }
    488         }
    489     }
    490 
    491     // test dump of the psf parameters
    492     if (psTraceGetLevel("psModules.objects") >= 4) {
    493         FILE *f = fopen ("pol.dat", "w");
    494         fprintf (f, "# x y  :  e0obs e1obs e2obs  : e0fit e1fit e2fit : mask\n");
    495         for (int i = 0; i < e0->n; i++) {
    496             fprintf (f, "%f %f  :  %f %f %f  : %f %f %f  : %d\n",
    497                      x->data.F32[i], y->data.F32[i],
    498                      e0->data.F32[i], e1->data.F32[i], e2->data.F32[i],
    499                      pmTrend2DEval (psf->params->data[PM_PAR_E0], x->data.F32[i], y->data.F32[i]),
    500                      pmTrend2DEval (psf->params->data[PM_PAR_E1], x->data.F32[i], y->data.F32[i]),
    501                      pmTrend2DEval (psf->params->data[PM_PAR_E2], x->data.F32[i], y->data.F32[i]),
    502                      srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]);
    503         }
    504         fclose (f);
    505     }
    506 
    507     psFree (e0);
    508     psFree (e1);
    509     psFree (e2);
    510     psFree (mag);
    511     return true;
    512 }
    513 
    514 // fit the shape variations as a psImageMap for the given scale factor
    515 bool pmPSFFitShapeParamsMap (pmPSF *psf, int scale, float *scatterTotal, psVector *mask, psVector *x, psVector *y, psVector *mag, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz) {
    516 
    517     int Nx, Ny;
    518 
    519     // set the map scale to match the aspect ratio : for a scale of 1, we guarantee
    520     // that we have a single cell
    521     if (psf->fieldNx > psf->fieldNy) {
    522         Nx = scale;
    523         float AR = psf->fieldNy / (float) psf->fieldNx;
    524         Ny = (int) (Nx * AR + 0.5);
    525         Ny = PS_MAX (1, Ny);
    526     } else {
    527         Ny = scale;
    528         float AR = psf->fieldNx / (float) psf->fieldNy;
    529         Nx = (int) (Ny * AR + 0.5);
    530         Nx = PS_MAX (1, Nx);
    531     }
    532 
    533     // do we have enough sources for this fine of a grid?
    534     if (x->n < 10*Nx*Ny) {
    535         return false;
    536     }
    537 
    538     // XXX check this against the exising type
    539     pmTrend2DMode psfTrendMode = PM_TREND_MAP;
    540 
    541     psImageBinning *binning = psImageBinningAlloc();
    542     binning->nXruff = Nx;
    543     binning->nYruff = Ny;
    544     binning->nXfine = psf->fieldNx;
    545     binning->nYfine = psf->fieldNy;
    546     psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
    547     psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo);
    548 
    549     psFree (psf->params->data[PM_PAR_E0]);
    550     psFree (psf->params->data[PM_PAR_E1]);
    551     psFree (psf->params->data[PM_PAR_E2]);
    552 
    553     int nIter = psf->psfTrendStats->clipIter;
    554     psf->psfTrendStats->clipIter = 1;
    555     psf->params->data[PM_PAR_E0] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats);
    556     psf->params->data[PM_PAR_E1] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats);
    557     psf->params->data[PM_PAR_E2] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats);
    558     psFree (binning);
    559 
    560     // if the map is 1x1 (a single value), we measure the resulting ensemble scatter
    561 
    562     // if the map is more finely sampled, divide the values into two sets: measure the fit from
    563     // one set and the scatter from the other set.
    564     psVector *x_fit = NULL;
    565     psVector *y_fit = NULL;
    566     psVector *x_tst = NULL;
    567     psVector *y_tst = NULL;
    568 
    569     psVector *e0obs_fit = NULL;
    570     psVector *e1obs_fit = NULL;
    571     psVector *e2obs_fit = NULL;
    572     psVector *e0obs_tst = NULL;
    573     psVector *e1obs_tst = NULL;
    574     psVector *e2obs_tst = NULL;
    575 
    576     if (scale == 1) {
    577         x_fit  = psMemIncrRefCounter (x);
    578         y_fit  = psMemIncrRefCounter (y);
    579         x_tst  = psMemIncrRefCounter (x);
    580         y_tst  = psMemIncrRefCounter (y);
    581         e0obs_fit = psMemIncrRefCounter (e0obs);
    582         e1obs_fit = psMemIncrRefCounter (e1obs);
    583         e2obs_fit = psMemIncrRefCounter (e2obs);
    584         e0obs_tst = psMemIncrRefCounter (e0obs);
    585         e1obs_tst = psMemIncrRefCounter (e1obs);
    586         e2obs_tst = psMemIncrRefCounter (e2obs);
    587     } else {
    588         x_fit  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    589         y_fit  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    590         x_tst  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    591         y_tst  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    592         e0obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    593         e1obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    594         e2obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    595         e0obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    596         e1obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    597         e2obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    598         for (int i = 0; i < e0obs_fit->n; i++) {
    599             // e0obs->n ==  8 or 9:
    600             // i = 0, 1, 2, 3 : 2i+0 = 0, 2, 4, 6
    601             // i = 0, 1, 2, 3 : 2i+1 = 1, 3, 5, 7
    602             x_fit->data.F32[i] = x->data.F32[2*i+0];
    603             x_tst->data.F32[i] = x->data.F32[2*i+1];
    604             y_fit->data.F32[i] = y->data.F32[2*i+0];
    605             y_tst->data.F32[i] = y->data.F32[2*i+1];
    606 
    607             e0obs_fit->data.F32[i] = e0obs->data.F32[2*i+0];
    608             e0obs_tst->data.F32[i] = e0obs->data.F32[2*i+1];
    609             e1obs_fit->data.F32[i] = e1obs->data.F32[2*i+0];
    610             e1obs_tst->data.F32[i] = e1obs->data.F32[2*i+1];
    611             e2obs_fit->data.F32[i] = e2obs->data.F32[2*i+0];
    612             e2obs_tst->data.F32[i] = e2obs->data.F32[2*i+1];
    613         }
    614     }
    615 
    616     // the mask marks the values not used to calculate the ApTrend
    617     psVector *fitMask = psVectorAlloc (x_fit->n, PS_TYPE_VECTOR_MASK);
    618     // copy mask values to fitMask as a starting point
    619     for (int i = 0; i < fitMask->n; i++) {
    620         fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    621     }
    622 
    623     // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
    624     // This way, the parameters masked by one of the fits will be applied to the others
    625     for (int i = 0; i < nIter; i++) {
    626         // XXX we are using the same stats structure on each pass: do we need to re-init it?
    627         psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
    628         psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
    629 
    630         pmTrend2D *trend = NULL;
    631         float mean, stdev;
    632 
    633         // XXX we are using the same stats structure on each pass: do we need to re-init it?
    634         bool status = true;
    635 
    636         trend = psf->params->data[PM_PAR_E0];
    637         status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e0obs_fit, NULL);
    638         mean = psStatsGetValue (trend->stats, meanOption);
    639         stdev = psStatsGetValue (trend->stats, stdevOption);
    640         psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs_fit->n);
    641         // pmTrend2DPrintMap (trend);
    642         psImageMapCleanup (trend->map);
    643         // pmTrend2DPrintMap (trend);
    644         pmSourceVisualPSFModelResid (trend, x, y, e0obs, mask);
    645 
    646         trend = psf->params->data[PM_PAR_E1];
    647         status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e1obs_fit, NULL);
    648         mean = psStatsGetValue (trend->stats, meanOption);
    649         stdev = psStatsGetValue (trend->stats, stdevOption);
    650         psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs_fit->n);
    651         // pmTrend2DPrintMap (trend);
    652         psImageMapCleanup (trend->map);
    653         // pmTrend2DPrintMap (trend);
    654         pmSourceVisualPSFModelResid (trend, x, y, e1obs, mask);
    655 
    656         trend = psf->params->data[PM_PAR_E2];
    657         status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e2obs_fit, NULL);
    658         mean = psStatsGetValue (trend->stats, meanOption);
    659         stdev = psStatsGetValue (trend->stats, stdevOption);
    660         psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n);
    661         // pmTrend2DPrintMap (trend);
    662         psImageMapCleanup (trend->map);
    663         // pmTrend2DPrintMap (trend);
    664         pmSourceVisualPSFModelResid (trend, x, y, e2obs, mask);
    665     }
    666     psf->psfTrendStats->clipIter = nIter; // restore default setting
    667 
    668     // construct the fitted values and the residuals
    669     psVector *e0fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E0], fitMask, 0xff, x_tst, y_tst);
    670     psVector *e1fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E1], fitMask, 0xff, x_tst, y_tst);
    671     psVector *e2fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E2], fitMask, 0xff, x_tst, y_tst);
    672 
    673     psVector *e0res = (psVector *) psBinaryOp (NULL, (void *) e0obs_tst, "-", (void *) e0fit);
    674     psVector *e1res = (psVector *) psBinaryOp (NULL, (void *) e1obs_tst, "-", (void *) e1fit);
    675     psVector *e2res = (psVector *) psBinaryOp (NULL, (void *) e2obs_tst, "-", (void *) e2fit);
    676 
    677     // measure scatter for the unfitted points
    678     // psTraceSetLevel ("psLib.math.vectorSampleStdev", 10);
    679     // psTraceSetLevel ("psLib.math.vectorClippedStats", 10);
    680     pmPSFShapeParamsScatter (scatterTotal, e0res, e1res, e2res, fitMask, 0xff, psStatsStdevOption(psf->psfTrendStats->options));
    681     // psTraceSetLevel ("psLib.math.vectorSampleStdev", 0);
    682     // psTraceSetLevel ("psLib.math.vectorClippedStats", 0);
    683 
    684     psLogMsg ("psphot.psftry", PS_LOG_INFO, "result of %d x %d grid\n", Nx, Ny);
    685     psLogMsg ("psphot.psftry", PS_LOG_INFO, "systematic scatter: %f\n", *scatterTotal);
    686 
    687     psFree (x_fit);
    688     psFree (y_fit);
    689     psFree (x_tst);
    690     psFree (y_tst);
    691 
    692     psFree (e0obs_fit);
    693     psFree (e1obs_fit);
    694     psFree (e2obs_fit);
    695     psFree (e0obs_tst);
    696     psFree (e1obs_tst);
    697     psFree (e2obs_tst);
    698 
    699     psFree (e0fit);
    700     psFree (e1fit);
    701     psFree (e2fit);
    702 
    703     psFree (e0res);
    704     psFree (e1res);
    705     psFree (e2res);
    706 
    707     // XXX copy fitMask values back to mask
    708     for (int i = 0; i < fitMask->n; i++) {
    709         mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    710     }
    711     psFree (fitMask);
    712 
    713     return true;
    714 }
    715 
    716 // calculate the scatter of the parameters
    717 bool pmPSFShapeParamsScatter(float *scatterTotal, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, psVectorMaskType maskValue, psStatsOptions stdevOpt)
    718 {
    719 
    720     // psStats *stats = psStatsAlloc(stdevOpt);
    721     psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_STDEV);
    722 
    723     // calculate the root-mean-square of E0, E1, E2
    724     float dEsquare = 0.0;
    725     psStatsInit (stats);
    726     if (!psVectorStats (stats, e0res, NULL, mask, maskValue)) {
    727         psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    728         return false;
    729     }
    730     dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
    731 
    732     psStatsInit (stats);
    733     if (!psVectorStats (stats, e1res, NULL, mask, maskValue)) {
    734         psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    735         return false;
    736     }
    737     dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
    738 
    739     psStatsInit (stats);
    740     if (!psVectorStats (stats, e2res, NULL, mask, maskValue)) {
    741         psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    742         return false;
    743     }
    744     dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
    745 
    746     *scatterTotal = sqrtf(dEsquare);
    747 
    748     psFree(stats);
    749     return true;
    750 }
    751 
    752 // calculate the minimum scatter of the parameters
    753 bool pmPSFShapeParamsErrors(float *errorFloor, psVector *mag, psVector *e0res, psVector *e1res,
    754                             psVector *e2res, psVector *mask, int nGroup, psStatsOptions stdevOpt)
    755 {
    756 
    757     psStats *statsS = psStatsAlloc(stdevOpt);
    758 
    759     // measure the trend in bins with 10 values each; if < 10 total, use them all
    760     int nBin = PS_MAX (mag->n / nGroup, 1);
    761 
    762     // use mag to group parameters in sequence
    763     psVector *index = psVectorSortIndex (NULL, mag);
    764 
    765     // subset vectors for mag and dap values within the given range
    766     psVector *dE0subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
    767     psVector *dE1subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
    768     psVector *dE2subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
    769     psVector *mkSubset  = psVectorAllocEmpty (nGroup, PS_TYPE_VECTOR_MASK);
    770 
    771     int n = 0;
    772     float min = INFINITY;               // Minimum error
    773     for (int i = 0; i < nBin; i++) {
    774         int j;
    775         int nValid = 0;
    776         for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) {
    777             int N = index->data.U32[n];
    778             dE0subset->data.F32[j] = e0res->data.F32[N];
    779             dE1subset->data.F32[j] = e1res->data.F32[N];
    780             dE2subset->data.F32[j] = e2res->data.F32[N];
    781 
    782             mkSubset->data.PS_TYPE_VECTOR_MASK_DATA[j]   = mask->data.PS_TYPE_VECTOR_MASK_DATA[N];
    783             if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[N]) nValid ++;
    784         }
    785         if (nValid < 3) continue;
    786 
    787         dE0subset->n = j;
    788         dE1subset->n = j;
    789         dE2subset->n = j;
    790         mkSubset->n = j;
    791 
    792         // calculate the root-mean-square of E0, E1, E2
    793         float dEsquare = 0.0;
    794         psStatsInit (statsS);
    795         if (!psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff)) {
    796         }
    797         dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
    798 
    799         psStatsInit (statsS);
    800         if (!psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff)) {
    801             psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    802             return false;
    803         }
    804         dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
    805 
    806         psStatsInit (statsS);
    807         if (!psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff)) {
    808             psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    809             return false;
    810         }
    811         dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
    812 
    813         if (isfinite(dEsquare)) {
    814             float err = sqrtf(dEsquare);
    815             if (err < min) {
    816                 min = err;
    817             }
    818         }
    819     }
    820     *errorFloor = min;
    821 
    822     psFree (dE0subset);
    823     psFree (dE1subset);
    824     psFree (dE2subset);
    825     psFree (mkSubset);
    826 
    827     psFree(index);
    828 
    829     psFree(statsS);
    830 
    831     return true;
    83289}
    83390
    … …  
    886143        if (mask->data.PS_TYPE_VECTOR_MASK_DATA[n]) continue;
    887144        res2mean += PS_SQR(residuals->data.F32[n]);
    888         ChiSq += PS_SQR(residuals->data.F32[n] / errors->data.F32[n]);
     145        ChiSq += PS_SQR(residuals->data.F32[n]) / PS_SQR(errors->data.F32[n]);
    889146        nPts += 1.0;
    890147    }
    … …  
    914171        float dS = (ChiSq - 1.0) / dRdS;
    915172        S2guess += dS;
     173        S2guess = PS_MAX(0.0, S2guess);
    916174
    917175        psLogMsg ("psModules", 3, "ChiSquare: %f, dS: %f, S2 guess: %f\n", ChiSq, dS, S2guess);

  • branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.h

    r25353 r25476  
    8989 *
    9090 */
    91 pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType mark);
     91pmPSFtry *pmPSFtryModel (
     92    const psArray *sources,             ///< PSF sources to use in the pmPSF model analysis
     93    const char *modelName,              ///< human-readable name of desired model
     94    pmPSFOptions *options,
     95    psImageMaskType maskVal,
     96    psImageMaskType mark
     97    );
     98
     99/** fit EXT models to all possible psf sources */
     100bool pmPSFtryFitEXT (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal);
     101
     102bool pmPSFtryMakePSF (pmPSFtry *psfTry);
     103
     104bool pmPSFtryFitPSF (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal);
    92105
    93106/** pmPSFtryMetric()
    … …  
    97110 *
    98111 */
    99 bool pmPSFtryMetric(
    100     pmPSFtry *psfTry,                  ///< Add comment.
    101     pmPSFOptions *options              ///< PSF fitting options
    102 );
     112bool pmPSFtryMetric(pmPSFtry *psfTry);
    103113
    104114/** pmPSFtryMetric_Alt()
    … …  
    112122    float RADIUS                        ///< Add comment.
    113123);
     124
     125bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask);
     126
     127float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction);
     128
     129/// @}
     130# endif
    114131
    115132/**
    … …  
    125142 *
    126143 */
    127 bool pmPSFFromPSFtry (pmPSFtry *psfTry);
    128144
    129 bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask);
    130 bool pmPSFFitShapeParamsMap (pmPSF *psf, int scale, float *scatterTotal, psVector *mask, psVector *x, psVector *y, psVector *mag, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz);
    131 bool pmPSFShapeParamsScatter(float *scatterTotal, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, psVectorMaskType maskValue, psStatsOptions stdevOpt);
    132 bool pmPSFShapeParamsErrors (float *errorFloor, psVector *mag, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, int nGroup, psStatsOptions stdevOpt);
    133 
    134 float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction);
    135 
    136 /// @}
    137 # endif

  • branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitEXT.c

    r25455 r25476  
    6161        }
    6262
    63         source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type);
     63        source->modelEXT = pmSourceModelGuess (source, options->type);
    6464        if (source->modelEXT == NULL) {
    6565            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    … …  
    8686        Next ++;
    8787    }
    88     psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit ext:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Next, sources->n);
    89     psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, sources->n);
     88    psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit ext:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Next, psfTry->sources->n);
     89    psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, psfTry->sources->n);
    9090
    9191    if (Next == 0) {

  • branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitPSF.c

    r25455 r25476  
    3939// stage 3: Refit with fixed shape parameters.  This function uses the LMM fitting, but could
    4040// be re-written to use the simultaneous linear fitting (see psphotFitSourcesLinear.c)
    41 bool pmPSFtryFitPSF (pmPSFtry *psfTry) {
     41bool pmPSFtryFitPSF (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) {
     42
     43    bool status;
    4244
    4345    psTimerStart ("psf.fit");
     46
     47    // maskVal is used to test for rejected pixels, and must include markVal
     48    maskVal |= markVal;
     49
     50    int Npsf = 0;
    4451    for (int i = 0; i < psfTry->sources->n; i++) {
    4552
    … …  
    98105    psfTry->psf->nPSFstars = Npsf;
    99106
    100     psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Npsf, sources->n);
    101     psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, sources->n);
     107    psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Npsf, psfTry->sources->n);
     108    psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, psfTry->sources->n);
    102109
    103110    if (Npsf == 0) {
    104111        psError(PS_ERR_UNKNOWN, false, "No sources with good PSF fits after model is built.");
    105         psFree(psfTry);
    106         return NULL;
     112        return false;
    107113    }
    108114
     115    return true;
     116}

  • branches/eam_branches/20090715/psModules/src/objects/pmPSFtryMakePSF.c

    • Property svn:mergeinfo set to (toggle deleted branches)
      /branches/czw_branch/cleanup/psModules/src/objects/pmPSFFromPSFtry.c24713-25285
      /branches/eam_branches/20090522/psModules/src/objects/pmPSFFromPSFtry.c24238-24573
      /branches/pap/psModules/src/objects/pmPSFFromPSFtry.c25238-25382
      /branches/pap_mops/psModules/src/objects/pmPSFFromPSFtry.c25137-25255
      /trunk/psModules/src/objects/pmPSFFromPSFtry.c24799-25395
    r25455 r25476  
    11/** @file  pmPSFtry.c
    2  *
    3  *  XXX: need description of file purpose
     2 *  @brief generate a pmPSF from a collection of EXT measurments of likely PSF stars.
    43 *
    54 *  @author EAM, IfA
    … …  
    4443Note: some of the array entries may be NULL (failed fits); ignore them.
    4544 *****************************************************************************/
    46 bool pmPSFFromPSFtry (pmPSFtry *psfTry, int Nx, int Ny)
     45bool pmPSFtryMakePSF (pmPSFtry *psfTry)
    4746{
    4847    PS_ASSERT_PTR_NON_NULL(psfTry, false);
    … …  
    5049
    5150    pmPSF *psf = psfTry->psf;
     51    psVector *srcMask = psfTry->mask;
    5252
    5353    // construct the fit vectors from the collection of objects
    5454    psVector *x  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    5555    psVector *y  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    56     psVector *z  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    5756
    5857    // construct the x,y terms
    5958    for (int i = 0; i < psfTry->sources->n; i++) {
     59        if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
     60
    6061        pmSource *source = psfTry->sources->data[i];
    61         if (source->modelEXT == NULL)
    62             continue;
     62        assert (source->modelEXT); // all unmasked sources should have modelEXT
    6363
    6464        x->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS];
    … …  
    6666    }
    6767
    68     if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, psfTry->mask)) {
     68    // fit the shape parameters (SXX, SYY, SXY) as a function of position
     69    if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, srcMask)) {
    6970        psFree(x);
    7071        psFree(y);
    71         psFree(z);
    7272        return false;
    7373    }
    7474
    75     // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY)
    76     // apply the values of Nx, Ny determined above for E0,E1,E2 to the remaining parameters
     75    // vector to store the other parameter values
     76    psVector *z  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
     77
     78    // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY);
     79    // fit the remaining parameters.
    7780    for (int i = 0; i < psf->params->n; i++) {
    7881        switch (i) {
    … …  
    9194        // select the per-object fitted data for this parameter
    9295        for (int j = 0; j < psfTry->sources->n; j++) {
     96            // skip any masked sources (failed to fit one of the model steps or get a magnitude)
     97            if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
     98
    9399            pmSource *source = psfTry->sources->data[j];
    94             if (source->modelEXT == NULL) continue;
     100            assert (source->modelEXT); // all unmasked sources should have modelEXT
     101
    95102            z->data.F32[j] = source->modelEXT->params->data.F32[i];
    96103        }
    97 
    98         psImageBinning *binning = psImageBinningAlloc();
    99         binning->nXruff = psf->trendNx;
    100         binning->nYruff = psf->trendNy;
    101         binning->nXfine = psf->fieldNx;
    102         binning->nYfine = psf->fieldNy;
    103 
    104         if (psf->psfTrendMode == PM_TREND_MAP) {
    105             psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
    106             psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo);
    107         }
    108 
    109         // free existing trend, re-alloc
    110         psFree (psf->params->data[i]);
    111         psf->params->data[i] = pmTrend2DNoImageAlloc (psf->psfTrendMode, binning, psf->psfTrendStats);
    112         psFree (binning);
    113104
    114105        // fit the collection of measured parameters to the PSF 2D model
    115106        // the mask is carried from previous steps and updated with this operation
    116107        // the weight is either the flux error or NULL, depending on 'psf->poissonErrorParams'
    117         if (!pmTrend2DFit (psf->params->data[i], psfTry->mask, 0xff, x, y, z, NULL)) {
     108        if (!pmTrend2DFit (psf->params->data[i], srcMask, 0xff, x, y, z, NULL)) {
    118109            psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i);
    119110            psFree(x);
    … …  
    154145}
    155146
     147// fit the shape parameters using the supplied order (pmPSF->trendNx,trendNy)
     148bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) {
     149
     150    // we are doing a robust fit.  after each pass, we drop points which are more deviant than
     151    // three sigma.  the source mask (srcMask) is updated for each pass. 
     152
     153    // The shape parameters (SXX, SXY, SYY) are strongly coupled.  We have to handle them very
     154    // carefully.  First, we convert them to the Ellipse Polarization terms (E0, E1, E2) for
     155    // each source and fit this set of parameters.  These values are less tightly coupled, but
     156    // are still inter-related.  The fitted values do a good job of constraining the major axis
     157    // and the position angle, but the minor axis is weakly measured.  When we apply the PSF
     158    // model to construct a source model, we convert the fitted values of E0,E1,E2 to the shape
     159    // parameters, with the constraint that the minor axis must be greater than a minimum
     160    // threshold.
     161
     162    // XXX re-read the sextractor manual on handling 'infinitely thin' sources...
     163
     164    // storage vectors for the polarization terms & mags
     165    psVector *e0   = psVectorAlloc (sources->n, PS_TYPE_F32);
     166    psVector *e1   = psVectorAlloc (sources->n, PS_TYPE_F32);
     167    psVector *e2   = psVectorAlloc (sources->n, PS_TYPE_F32);
     168
     169    // convert the measured source shape paramters to polarization terms
     170    for (int i = 0; i < sources->n; i++) {
     171        // skip any masked sources (failed to fit one of the model steps or get a magnitude)
     172        if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
     173
     174        pmSource *source = sources->data[i];
     175        assert (source->modelEXT); // all unmasked sources should have modelEXT
     176
     177        psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32);
     178
     179        e0->data.F32[i] = pol.e0;
     180        e1->data.F32[i] = pol.e1;
     181        e2->data.F32[i] = pol.e2;
     182    }
     183
     184    // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
     185    // This way, the parameters masked by one of the fits will be applied to the others
     186    for (int i = 0; i < psf->psfTrendStats->clipIter; i++) {
     187        // XXX we are using the same stats structure on each pass: do we need to re-init it?
     188        // XXX we hardwire this to SAMPLE stats above (psphotChoosePSF.c), hardwire here instead?
     189        psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
     190        psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
     191
     192        pmTrend2D *trend = NULL;
     193        float mean, stdev;
     194
     195        // XXX we are using the same stats structure on each pass: do we need to re-init it?
     196        bool status = true;
     197
     198        trend = psf->params->data[PM_PAR_E0];
     199        trend->stats->clipIter = 1; // in allocation, this value is set to the value of nIter, but we should use 1 here
     200        status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e0, NULL);
     201        mean = psStatsGetValue (trend->stats, meanOption);
     202        stdev = psStatsGetValue (trend->stats, stdevOption);
     203        psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n);
     204        if (psf->psfTrendMode == PM_TREND_MAP) psImageMapCleanup (trend->map);
     205        // pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask);
     206
     207        trend = psf->params->data[PM_PAR_E1];
     208        trend->stats->clipIter = 1; // in allocation, this value is set to the value of nIter, but we should use 1 here
     209        status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e1, NULL);
     210        mean = psStatsGetValue (trend->stats, meanOption);
     211        stdev = psStatsGetValue (trend->stats, stdevOption);
     212        psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n);
     213        if (psf->psfTrendMode == PM_TREND_MAP) psImageMapCleanup (trend->map);
     214        // pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask);
     215
     216        trend = psf->params->data[PM_PAR_E2];
     217        trend->stats->clipIter = 1; // in allocation, this value is set to the value of nIter, but we should use 1 here
     218        status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e2, NULL);
     219        mean = psStatsGetValue (trend->stats, meanOption);
     220        stdev = psStatsGetValue (trend->stats, stdevOption);
     221        psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n);
     222        if (psf->psfTrendMode == PM_TREND_MAP) psImageMapCleanup (trend->map);
     223        // pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask);
     224
     225        if (!status) {
     226            psError (PS_ERR_UNKNOWN, true, "failed to fit PSF shape params");
     227            return false;
     228        }
     229    }
     230
     231    // test dump of the psf parameters
     232    if (psTraceGetLevel("psModules.objects") >= 4) {
     233        FILE *f = fopen ("pol.dat", "w");
     234        fprintf (f, "# x y  :  e0obs e1obs e2obs  : e0fit e1fit e2fit : mask\n");
     235        for (int i = 0; i < e0->n; i++) {
     236            fprintf (f, "%f %f  :  %f %f %f  : %f %f %f  : %d\n",
     237                     x->data.F32[i], y->data.F32[i],
     238                     e0->data.F32[i], e1->data.F32[i], e2->data.F32[i],
     239                     pmTrend2DEval (psf->params->data[PM_PAR_E0], x->data.F32[i], y->data.F32[i]),
     240                     pmTrend2DEval (psf->params->data[PM_PAR_E1], x->data.F32[i], y->data.F32[i]),
     241                     pmTrend2DEval (psf->params->data[PM_PAR_E2], x->data.F32[i], y->data.F32[i]),
     242                     srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]);
     243        }
     244        fclose (f);
     245    }
     246
     247    psFree (e0);
     248    psFree (e1);
     249    psFree (e2);
     250    return true;
     251}
     252

  • branches/eam_branches/20090715/psModules/src/objects/pmSourceVisual.c

    r25268 r25476  
    55#include <pslib.h>
    66#include "pmTrend2D.h"
     7#include "pmPSF.h"
     8#include "pmPSFtry.h"
    79#include "pmSourceVisual.h"
    810
    … …  
    2729bool pmSourcePlotPoints3D (int myKapa, Graphdata *graphdata, psVector *xn, psVector *yn, psVector *zn, float theta, float phi);
    2830
     31bool pmSourceVisualPlotPSFMetric (pmPSFtry *psfTry) {
     32
     33    Graphdata graphdata;
     34
     35    if (!pmVisualIsVisual()) return true;
     36
     37    if (kapa1 == -1) {
     38        kapa1 = KapaOpenNamedSocket ("kapa", "pmSource:plots");
     39        if (kapa1 == -1) {
     40            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     41            pmVisualSetVisual(false);
     42            return false;
     43        }
     44    }
     45
     46    KapaClearPlots (kapa1);
     47    KapaInitGraph (&graphdata);
     48
     49    psVector *x = psVectorAllocEmpty (psfTry->sources->n, PS_TYPE_F32);
     50    psVector *y = psVectorAllocEmpty (psfTry->sources->n, PS_TYPE_F32);
     51    psVector *dy = psVectorAllocEmpty(psfTry->sources->n, PS_TYPE_F32);
     52
     53    graphdata.xmin = +32.0;
     54    graphdata.xmax = -32.0;
     55    graphdata.ymin = +32.0;
     56    graphdata.ymax = -32.0;
     57
     58    // construct the plot vectors
     59    int n = 0;
     60    for (int i = 0; i < psfTry->sources->n; i++) {
     61        if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue;
     62        x->data.F32[n] = psfTry->fitMag->data.F32[i];
     63        y->data.F32[n] = psfTry->metric->data.F32[i];
     64        dy->data.F32[n] = psfTry->metricErr->data.F32[i];
     65        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
     66        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
     67        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]);
     68        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]);
     69        n++;
     70    }
     71    x->n = y->n = dy->n = n;
     72
     73    float range;
     74    range = graphdata.xmax - graphdata.xmin;
     75    graphdata.xmax += 0.05*range;
     76    graphdata.xmin -= 0.05*range;
     77    range = graphdata.ymax - graphdata.ymin;
     78    graphdata.ymax += 0.05*range;
     79    graphdata.ymin -= 0.05*range;
     80
     81    KapaSetLimits (kapa1, &graphdata);
     82
     83    KapaSetFont (kapa1, "helvetica", 14);
     84    KapaBox (kapa1, &graphdata);
     85    KapaSendLabel (kapa1, "PSF Mag", KAPA_LABEL_XM);
     86    KapaSendLabel (kapa1, "Ap Mag - PSF Mag", KAPA_LABEL_YM);
     87
     88    graphdata.color = KapaColorByName ("black");
     89    graphdata.ptype = 2;
     90    graphdata.size = 0.5;
     91    graphdata.style = 2;
     92    graphdata.etype |= 0x01;
     93
     94    KapaPrepPlot (kapa1, n, &graphdata);
     95    KapaPlotVector (kapa1, n, x->data.F32, "x");
     96    KapaPlotVector (kapa1, n, y->data.F32, "y");
     97    KapaPlotVector (kapa1, n, dy->data.F32, "dym");
     98    KapaPlotVector (kapa1, n, dy->data.F32, "dyp");
     99
     100    psFree (x);
     101    psFree (y);
     102    psFree (dy);
     103
     104    // pause and wait for user input:
     105    // continue, save (provide name), ??
     106    char key[10];
     107    fprintf (stdout, "[c]ontinue? ");
     108    if (!fgets(key, 8, stdin)) {
     109        psWarning("Unable to read option");
     110    }
     111    return true;
     112}
    29113
    30114bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask) {

  • branches/eam_branches/20090715/psModules/src/objects/pmSourceVisual.h

    r23242 r25476  
    1818
    1919bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask);
     20bool pmSourceVisualPlotPSFMetric (pmPSFtry *try);
    2021
    2122/// @}
Note: See TracChangeset for help on using the changeset viewer.