IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 2, 2009, 3:11:32 PM (17 years ago)
Author:
eugene
Message:

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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects/pmPSFtry.c

    r24206 r25754  
    3737#include "pmSourceVisual.h"
    3838
    39 bool printTrendMap (pmTrend2D *trend) {
    40 
    41     if (!trend->map) return false;
    42     if (!trend->map->map) return false;
    43 
    44     for (int j = 0; j < trend->map->map->numRows; j++) {
    45         for (int i = 0; i < trend->map->map->numCols; i++) {
    46             fprintf (stderr, "%5.2f  ", trend->map->map->data.F32[j][i]);
    47         }
    48         fprintf (stderr, "\t\t\t");
    49         for (int i = 0; i < trend->map->map->numCols; i++) {
    50             fprintf (stderr, "%5.2f  ", trend->map->error->data.F32[j][i]);
    51         }
    52         fprintf (stderr, "\n");
    53     }
    54     return true;
    55 }
    56 
    57 bool psImageMapCleanup (psImageMap *map) {
    58 
    59     if ((map->map->numRows == 1) && (map->map->numCols == 1)) return true;
    60 
    61     // find the weighted average of all pixels
    62     float Sum = 0.0;
    63     float Wt = 0.0;
    64     for (int j = 0; j < map->map->numRows; j++) {
    65         for (int i = 0; i < map->map->numCols; i++) {
    66             if (!isfinite(map->map->data.F32[j][i])) continue;
    67             Sum += map->map->data.F32[j][i] * map->error->data.F32[j][i];
    68             Wt += map->error->data.F32[j][i];
    69         }
    70     }
    71 
    72     float Mean = Sum / Wt;
    73 
    74     // do any of the pixels in the map need to be repaired?
    75     // XXX for now, we are just replacing bad pixels with the Mean
    76     for (int j = 0; j < map->map->numRows; j++) {
    77         for (int i = 0; i < map->map->numCols; i++) {
    78             if (isfinite(map->map->data.F32[j][i]) &&
    79                 (map->error->data.F32[j][i] > 0.0)) continue;
    80             map->map->data.F32[j][i] = Mean;
    81         }
    82     }
    83     return true;
    84 }
    85 
    8639// ********  pmPSFtry functions  **************************************************
    8740// * pmPSFtry holds a single pmPSF model test, with the input sources, the freely
    … …  
    11063    psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree);
    11164
    112     test->psf       = pmPSFAlloc (options);
     65    test->psf       = NULL;
    11366    test->metric    = psVectorAlloc (sources->n, PS_TYPE_F32);
    11467    test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F32);
    … …  
    13689}
    13790
     91float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction) {
    13892
    139 // build a pmPSFtry for the given model:
    140 // - fit each source with the free-floating model
    141 // - construct the pmPSF from the collection of models
    142 // - fit each source with the PSF-parameter models
    143 // - measure the pmPSF quality metric (dApResid)
     93    psAssert(residuals, "residuals cannot be NULL");
     94    psAssert(errors, "errors cannot be NULL");
     95    psAssert(residuals->n == errors->n, "residuals and errors must be the same length");
    14496
    145 // sources used in for pmPSFtry may be masked by the analysis
    146 // mask values indicate the reason the source was rejected:
     97    // given a vector of residuals and their formal errors, calculated the necessary systematic
     98    // error needed to yield a reduced chisq of 1.0, after first tossing out the clipFraction
     99    // highest chi-square contributors (allowed outliers)
    147100
    148 // generate a pmPSFtry with a copy of the test PSF sources
    149 pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal)
    150 {
    151     bool status;
    152     int Next = 0;
    153     int Npsf = 0;
     101    psVector *mask  = psVectorAlloc(residuals->n, PS_TYPE_VECTOR_MASK);
     102    psVector *chisq = psVectorAlloc(residuals->n, PS_TYPE_F32);
    154103
    155     // validate the requested model name
    156     options->type = pmModelClassGetType (modelName);
    157     if (options->type == -1) {
    158         psError (PS_ERR_UNKNOWN, true, "invalid model name %s", modelName);
    159         return NULL;
     104    // calculate the chisq vector:
     105    int Ngood = 0;
     106    for (int i = 0; i < residuals->n; i++) {
     107        chisq->data.F32[i] = PS_MAX_F32;
     108        if (!isfinite(residuals->data.F32[i])) continue;
     109        if (!isfinite(errors->data.F32[i])) continue;
     110        if (errors->data.F32[i] <= 0.0) continue;
     111        chisq->data.F32[i] = PS_SQR(residuals->data.F32[i] / errors->data.F32[i]);
     112        Ngood ++;
    160113    }
    161114
    162     pmPSFtry *psfTry = pmPSFtryAlloc (sources, options);
    163     if (psfTry == NULL) {
    164         psError (PS_ERR_UNKNOWN, false, "failed to allocate psf model");
    165         return NULL;
     115    psVector *index = psVectorSortIndex(NULL, chisq);
     116
     117    // toss out the clipFraction highest chisq values
     118    for (int i = 0; i < residuals->n; i++) {
     119        int n = index->data.S32[i];
     120        if (i < (1.0 - clipFraction)*Ngood) {
     121            mask->data.PS_TYPE_VECTOR_MASK_DATA[n] = 0;
     122        } else {
     123            mask->data.PS_TYPE_VECTOR_MASK_DATA[n] = 1;
     124        }
    166125    }
    167126
    168     // maskVal is used to test for rejected pixels, and must include markVal
    169     maskVal |= markVal;
     127    // Ndof ~= Ngood
     128    // Chisq_Ndof = sum(residuals_i^2 / error_i^2) / Ndof
     129    // choose S2 such than Chisq^sys_Ndof = sum(residuals_i^2 / (error_i^2 + S2)) / Ndof = 1.0
     130   
     131    // use Newton-Raphson to solve for S2:
    170132
    171     // stage 1:  fit an EXT model to all candidates PSF sources -- this is independent of the modeled 2D variations in the PSF
    172     psTimerStart ("psf.fit");
    173     for (int i = 0; i < psfTry->sources->n; i++) {
     133    // use median sigma to calculate the initial guess for S2:
     134    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
     135    psVectorStats (stats, errors, NULL, mask, 1);
     136    float errorMedian = stats->sampleMedian;
     137   
     138    float nPts = 0.0;
     139    float res2mean = 0.0;
     140    float ChiSq = 0.0;
     141    for (int i = 0; i < residuals->n; i++) {
     142        int n = index->data.S32[i];
     143        if (mask->data.PS_TYPE_VECTOR_MASK_DATA[n]) continue;
     144        res2mean += PS_SQR(residuals->data.F32[n]);
     145        ChiSq += PS_SQR(residuals->data.F32[n]) / PS_SQR(errors->data.F32[n]);
     146        nPts += 1.0;
     147    }
     148    res2mean /= nPts;
     149    ChiSq /= nPts;
     150   
     151    float S2guess = res2mean - PS_SQR(errorMedian);
    174152
    175         pmSource *source = psfTry->sources->data[i];
    176         if (!source->moments) {
    177             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    178             continue;
    179         }
    180         if (!source->moments->nPixels) {
    181             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    182             continue;
    183         }
     153    psLogMsg ("psModules", 10, "ChiSquare: %f, Ntotal: %ld, Ngood: %d, Nkeep: %.0f, S2 guess: %f\n",
     154              ChiSq, residuals->n, Ngood, nPts, S2guess);
    184155
    185         source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type);
    186         if (source->modelEXT == NULL) {
    187             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    188             psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y);
    189             continue;
    190         }
     156    for (int iter = 0; iter < 10; iter++) {
    191157
    192         // set object mask to define valid pixels
    193         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal);
     158        ChiSq = 0.0;
     159        float dRdS = 0.0;
     160        for (int i = 0; i < residuals->n; i++) {
     161            int n = index->data.S32[i];
     162            if (mask->data.PS_TYPE_VECTOR_MASK_DATA[n]) continue;
     163            float error2 = PS_SQR(errors->data.F32[n]) + S2guess;
     164            ChiSq += PS_SQR(residuals->data.F32[n]) / error2;
     165            dRdS += PS_SQR(residuals->data.F32[n]) / PS_SQR(error2);
     166        }
     167        ChiSq /= nPts;
     168        dRdS /= nPts;
    194169
    195         // fit model as EXT, not PSF
    196         status = pmSourceFitModel (source, source->modelEXT, PM_SOURCE_FIT_EXT, maskVal);
     170        // Note the sign on dS: dRdS above is -1 * dR/dS formally
     171        float dS = (ChiSq - 1.0) / dRdS;
     172        S2guess += dS;
     173        S2guess = PS_MAX(0.0, S2guess);
    197174
    198         // clear object mask to define valid pixels
    199         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
    200 
    201         // exclude the poor fits
    202         if (!status) {
    203             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    204             psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y);
    205             continue;
    206         }
    207         Next ++;
    208     }
    209     psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit ext:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Next, sources->n);
    210     psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, sources->n);
    211 
    212     if (Next == 0) {
    213         psError(PS_ERR_UNKNOWN, false, "No sources with good extended fits from which to determine PSF.");
    214         psFree(psfTry);
    215         return NULL;
     175        psLogMsg ("psModules", 10, "ChiSquare: %f, dS: %f, S2 guess: %f\n", ChiSq, dS, S2guess);
    216176    }
    217177
    218     // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation
    219     if (!pmPSFFromPSFtry (psfTry)) {
    220         psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources");
    221         psFree(psfTry);
    222         return NULL;
    223     }
     178    // free local allocations
     179    psFree (mask);
     180    psFree (chisq);
     181    psFree (stats);
     182    psFree (index);
    224183
    225     // stage 3: refit with fixed shape parameters
    226     psTimerStart ("psf.fit");
    227     for (int i = 0; i < psfTry->sources->n; i++) {
    228 
    229         pmSource *source = psfTry->sources->data[i];
    230 
    231         // masked for: bad model fit, outlier in parameters
    232         if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) {
    233             psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : source is masked\n", i, source->peak->x, source->peak->y);
    234             continue;
    235         }
    236 
    237         // set shape for this model based on PSF
    238         source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf);
    239         if (source->modelPSF == NULL) {
    240             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL;
    241             abort();
    242             continue;
    243         }
    244         source->modelPSF->radiusFit = options->radius;
    245 
    246         // set object mask to define valid pixels
    247         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal);
    248 
    249         // fit the PSF model to the source
    250         status = pmSourceFitModel (source, source->modelPSF, PM_SOURCE_FIT_PSF, maskVal);
    251 
    252         // skip poor fits
    253         if (!status) {
    254             psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
    255             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL;
    256             psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y);
    257             continue;
    258         }
    259 
    260         status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal);
    261         if (!status || isnan(source->apMag)) {
    262             psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
    263             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT;
    264             psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y);
    265             continue;
    266         }
    267 
    268         // clear object mask to define valid pixels
    269         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
    270 
    271         psfTry->fitMag->data.F32[i] = source->psfMag;
    272         psfTry->metric->data.F32[i] = source->apMag - source->psfMag;
    273         psfTry->metricErr->data.F32[i] = source->errMag;
    274 
    275         psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n);
    276         Npsf ++;
    277     }
    278     psfTry->psf->nPSFstars = Npsf;
    279 
    280     psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Npsf, sources->n);
    281     psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, sources->n);
    282 
    283     if (Npsf == 0) {
    284         psError(PS_ERR_UNKNOWN, false, "No sources with good PSF fits after model is built.");
    285         psFree(psfTry);
    286         return NULL;
    287     }
    288 
    289     // measure the chi-square trend as a function of flux (PAR[PM_PAR_I0])
    290     // this should be linear for Poisson errors and quadratic for constant sky errors
    291     psVector *flux  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    292     psVector *chisq = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    293     psVector *mask  = psVectorAlloc (psfTry->sources->n, PS_TYPE_VECTOR_MASK);
    294 
    295     // generate the x and y vectors, and mask missing models
    296     for (int i = 0; i < psfTry->sources->n; i++) {
    297         pmSource *source = psfTry->sources->data[i];
    298         if (source->modelPSF == NULL) {
    299             flux->data.F32[i] = 0.0;
    300             chisq->data.F32[i] = 0.0;
    301             mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
    302         } else {
    303             flux->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_I0];
    304             chisq->data.F32[i] = source->modelPSF->chisq / source->modelPSF->nDOF;
    305             mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
    306         }
    307     }
    308 
    309     // use 3hi/3lo sigma clipping on the chisq fit
    310     psStats *stats = options->stats;
    311 
    312     // linear clipped fit of chisq trend vs flux
    313     if (options->chiFluxTrend) {
    314         bool result = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, options->stats,
    315                                                   mask, 0xff, chisq, NULL, flux);
    316         psStatsOptions meanStat = psStatsMeanOption(options->stats->options); // Statistic for mean
    317         psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev
    318 
    319         psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n",
    320                   psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat));
    321 
    322         psFree(flux);
    323         psFree(mask);
    324         psFree(chisq);
    325 
    326         if (!result) {
    327             psError(PS_ERR_UNKNOWN, false, "Failed to fit psf->ChiTrend");
    328             psFree(psfTry);
    329             return NULL;
    330         }
    331     }
    332 
    333     for (int i = 0; i < psfTry->psf->ChiTrend->nX + 1; i++) {
    334         psLogMsg ("pmPSFtry", 4, "chisq vs flux fit term %d: %f +/- %f\n", i,
    335                   psfTry->psf->ChiTrend->coeff[i]*pow(10000, i),
    336                   psfTry->psf->ChiTrend->coeffErr[i]*pow(10000,i));
    337     }
    338 
    339     // XXX this function wants aperture radius for pmSourcePhotometry
    340     if (!pmPSFtryMetric (psfTry, options)) {
    341         psError(PS_ERR_UNKNOWN, false, "Attempt to fit PSF with model %s failed.", modelName);
    342         psFree (psfTry);
    343         return NULL;
    344     }
    345 
    346     psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n",
    347               modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias);
    348 
    349     return (psfTry);
     184    return (sqrt(S2guess));
    350185}
    351186
    352 bool pmPSFtryMetric (pmPSFtry *psfTry, pmPSFOptions *options)
    353 {
    354     PS_ASSERT_PTR_NON_NULL(psfTry, false);
    355     PS_ASSERT_PTR_NON_NULL(options, false);
    356     PS_ASSERT_PTR_NON_NULL(psfTry->sources, false);
    357 
    358     float RADIUS = options->radius;
    359 
    360     // the measured (aperture - fit) magnitudes (dA == psfTry->metric)
    361     //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
    362     //     dA = dAo + dsky/flux
    363     //   where flux is the flux of the star
    364     // we fit this trend to find the infinite flux aperture correction (dAo),
    365     //   the nominal sky bias (dsky), and the error on dAo
    366     // the values of dA are contaminated by stars with close neighbors in the aperture
    367     //   we use an outlier rejection to avoid this bias
    368 
    369     // r2rflux = radius^2 * ten(0.4*fitMag);
    370     psVector *r2rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    371 
    372     for (int i = 0; i < psfTry->sources->n; i++) {
    373         if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL)
    374             continue;
    375         r2rflux->data.F32[i] = PS_SQR(RADIUS) * pow(10.0, 0.4*psfTry->fitMag->data.F32[i]);
    376     }
    377 
    378     // XXX test dump of aperture residual data
    379     if (psTraceGetLevel("psModules.objects") >= 5) {
    380         FILE *f = fopen ("apresid.dat", "w");
    381         for (int i = 0; i < psfTry->sources->n; i++) {
    382             int keep = (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL);
    383 
    384             pmSource *source = psfTry->sources->data[i];
    385             float x = source->peak->x;
    386             float y = source->peak->y;
    387 
    388             fprintf (f, "%d  %d, %f %f %f  %f %f %f \n",
    389                      i, keep, x, y,
    390                      psfTry->fitMag->data.F32[i],
    391                      r2rflux->data.F32[i],
    392                      psfTry->metric->data.F32[i],
    393                      psfTry->metricErr->data.F32[i]);
    394         }
    395         fclose (f);
    396     }
    397 
    398     // This analysis of the apResid statistics is only approximate.  The fitted magnitudes
    399     // measured at this point (in the PSF fit) use Poisson errors, and are thus biased as a
    400     // function of magnitude.  We re-measure the apResid statistics later in psphot using the
    401     // linear, constant-error fitting.  Do not reject outliers with excessive vigor here.
    402 
    403     // fit ApTrend only to r2rflux, ignore x,y,flux variations for now
    404     // linear clipped fit of ApResid to r2rflux
    405     psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
    406     poly->coeffMask[1] = PS_POLY_MASK_SET; // fit only a constant offset (no SKYBIAS)
    407 
    408     // XXX replace this with a psVectorStats call?  since we are not fitting the trend
    409     bool result = psVectorClipFitPolynomial1D(poly, options->stats, psfTry->mask, PSFTRY_MASK_ALL,
    410                                               psfTry->metric, psfTry->metricErr, r2rflux);
    411     if (!result) {
    412         psError(PS_ERR_UNKNOWN, false, "Failed to fit clipped poly");
    413 
    414         psFree(poly);
    415         psFree(r2rflux);
    416 
    417         return false;
    418     }
    419     psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev
    420     psLogMsg ("pmPSFtryMetric", 4, "apresid: %f +/- %f; from statistics of %ld psf stars\n", poly->coeff[0],
    421               psStatsGetValue(options->stats, stdevStat), psfTry->sources->n);
    422 
    423 
    424 
    425     // XXX test dump of fitted model (dump when tracing?)
    426     if (psTraceGetLevel("psModules.objects") >= 4) {
    427         FILE *f = fopen ("resid.dat", "w");
    428         psVector *apfit = psPolynomial1DEvalVector (poly, r2rflux);
    429         for (int i = 0; i < psfTry->sources->n; i++) {
    430             int keep = (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL);
    431 
    432             pmSource *source = psfTry->sources->data[i];
    433             float x = source->peak->x;
    434             float y = source->peak->y;
    435 
    436             fprintf (f, "%d  %d, %f %f %f  %f %f %f  %f\n",
    437                      i, keep, x, y,
    438                      psfTry->fitMag->data.F32[i],
    439                      r2rflux->data.F32[i],
    440                      psfTry->metric->data.F32[i],
    441                      psfTry->metricErr->data.F32[i],
    442                      apfit->data.F32[i]);
    443         }
    444         fclose (f);
    445         psFree (apfit);
    446     }
    447 
    448     // XXX drop the skyBias value, or include above??
    449     psfTry->psf->skyBias  = poly->coeff[1];
    450     psfTry->psf->ApResid  = poly->coeff[0];
    451     psfTry->psf->dApResid = psStatsGetValue(options->stats, stdevStat);
    452 
    453     psFree (r2rflux);
    454     psFree (poly);
    455 
    456     return true;
    457 }
    458 
    459 /*
    460   (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
    461   (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
    462   (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
    463 */
    464 
    465 /*****************************************************************************
    466 pmPSFFromPSFtry (psfTry): build a PSF model from a collection of
    467 source->modelEXT entries.  The PSF ignores the first 4 (independent) model
    468 parameters and constructs a polynomial fit to the remaining as a function of
    469 image coordinate.
    470     input: psfTry with fitted source->modelEXT collection, pre-allocated psf
    471 Note: some of the array entries may be NULL (failed fits); ignore them.
    472  *****************************************************************************/
    473 bool pmPSFFromPSFtry (pmPSFtry *psfTry)
    474 {
    475     PS_ASSERT_PTR_NON_NULL(psfTry, false);
    476     PS_ASSERT_PTR_NON_NULL(psfTry->sources, false);
    477 
    478     pmPSF *psf = psfTry->psf;
    479 
    480     // construct the fit vectors from the collection of objects
    481     psVector *x  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    482     psVector *y  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    483     psVector *z  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    484 
    485     // construct the x,y terms
    486     for (int i = 0; i < psfTry->sources->n; i++) {
    487         pmSource *source = psfTry->sources->data[i];
    488         if (source->modelEXT == NULL)
    489             continue;
    490 
    491         x->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS];
    492         y->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS];
    493     }
    494 
    495     if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, psfTry->mask)) {
    496         psFree(x);
    497         psFree(y);
    498         psFree(z);
    499         return false;
    500     }
    501 
    502     // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY)
    503     // apply the values of Nx, Ny determined above for E0,E1,E2 to the remaining parameters
    504     for (int i = 0; i < psf->params->n; i++) {
    505         switch (i) {
    506           case PM_PAR_SKY:
    507           case PM_PAR_I0:
    508           case PM_PAR_XPOS:
    509           case PM_PAR_YPOS:
    510           case PM_PAR_SXX:
    511           case PM_PAR_SYY:
    512           case PM_PAR_SXY:
    513             continue;
    514           default:
    515             break;
    516         }
    517 
    518         // select the per-object fitted data for this parameter
    519         for (int j = 0; j < psfTry->sources->n; j++) {
    520             pmSource *source = psfTry->sources->data[j];
    521             if (source->modelEXT == NULL) continue;
    522             z->data.F32[j] = source->modelEXT->params->data.F32[i];
    523         }
    524 
    525         psImageBinning *binning = psImageBinningAlloc();
    526         binning->nXruff = psf->trendNx;
    527         binning->nYruff = psf->trendNy;
    528         binning->nXfine = psf->fieldNx;
    529         binning->nYfine = psf->fieldNy;
    530 
    531         if (psf->psfTrendMode == PM_TREND_MAP) {
    532             psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
    533             psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo);
    534         }
    535 
    536         // free existing trend, re-alloc
    537         psFree (psf->params->data[i]);
    538         psf->params->data[i] = pmTrend2DNoImageAlloc (psf->psfTrendMode, binning, psf->psfTrendStats);
    539         psFree (binning);
    540 
    541         // fit the collection of measured parameters to the PSF 2D model
    542         // the mask is carried from previous steps and updated with this operation
    543         // the weight is either the flux error or NULL, depending on 'psf->poissonErrorParams'
    544         if (!pmTrend2DFit (psf->params->data[i], psfTry->mask, 0xff, x, y, z, NULL)) {
    545             psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i);
    546             psFree(x);
    547             psFree(y);
    548             psFree(z);
    549             return false;
    550         }
    551     }
    552 
    553     // test dump of star parameters vs position (compare with fitted values)
    554     if (psTraceGetLevel("psModules.objects") >= 4) {
    555         FILE *f = fopen ("params.dat", "w");
    556 
    557         for (int j = 0; j < psfTry->sources->n; j++) {
    558             pmSource *source = psfTry->sources->data[j];
    559             if (source == NULL) continue;
    560             if (source->modelEXT == NULL) continue;
    561 
    562             pmModel *modelPSF = pmModelFromPSF (source->modelEXT, psf);
    563 
    564             fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[PM_PAR_XPOS], source->modelEXT->params->data.F32[PM_PAR_YPOS]);
    565 
    566             for (int i = 0; i < psf->params->n; i++) {
    567                 if (psf->params->data[i] == NULL) continue;
    568                 fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[i], modelPSF->params->data.F32[i]);
    569             }
    570             fprintf (f, "%f %d\n", source->modelEXT->chisq, source->modelEXT->nIter);
    571 
    572             psFree(modelPSF);
    573         }
    574         fclose (f);
    575     }
    576 
    577     psFree (x);
    578     psFree (y);
    579     psFree (z);
    580     return true;
    581 }
    582 
    583 
    584 bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) {
    585 
    586     // we are doing a robust fit.  after each pass, we drop points which are more deviant than
    587     // three sigma.  mask is currently updated for each pass.
    588 
    589     // The shape parameters (SXX, SXY, SYY) are strongly coupled.  We have to handle them very
    590     // carefully.  First, we convert them to the Ellipse Polarization terms (E0, E1, E2) for
    591     // each source and fit this set of parameters.  These values are less tightly coupled, but
    592     // are still inter-related.  The fitted values do a good job of constraining the major axis
    593     // and the position angle, but the minor axis is weakly measured.  When we apply the PSF
    594     // model to construct a source model, we convert the fitted values of E0,E1,E2 to the shape
    595     // parameters, with the constraint that the minor axis must be greater than a minimum
    596     // threshold.
    597 
    598     // convert the measured source shape paramters to polarization terms
    599     psVector *e0   = psVectorAlloc (sources->n, PS_TYPE_F32);
    600     psVector *e1   = psVectorAlloc (sources->n, PS_TYPE_F32);
    601     psVector *e2   = psVectorAlloc (sources->n, PS_TYPE_F32);
    602     psVector *mag  = psVectorAlloc (sources->n, PS_TYPE_F32);
    603 
    604     for (int i = 0; i < sources->n; i++) {
    605         // skip any masked sources (failed to fit one of the model steps or get a magnitude)
    606         if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
    607 
    608         pmSource *source = sources->data[i];
    609         assert (source->modelEXT); // all unmasked sources should have modelEXT
    610 
    611         psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32);
    612 
    613         e0->data.F32[i] = pol.e0;
    614         e1->data.F32[i] = pol.e1;
    615         e2->data.F32[i] = pol.e2;
    616 
    617         float flux = source->modelEXT->params->data.F32[PM_PAR_I0];
    618         mag->data.F32[i] = (flux > 0.0) ? -2.5*log(flux) : -100.0;
    619     }
    620 
    621     if (psf->psfTrendMode == PM_TREND_MAP) {
    622         float scatterTotal = 0.0;
    623         float scatterTotalMin = FLT_MAX;
    624         int entryMin = -1;
    625 
    626         psVector *dz = NULL;
    627         psVector *mask = psVectorAlloc (sources->n, PS_TYPE_VECTOR_MASK);
    628 
    629         // check the fit residuals and increase Nx,Ny until the error is minimized
    630         // pmPSFParamTrend increases the number along the longer of x or y
    631         for (int i = 1; i <= PS_MAX (psf->trendNx, psf->trendNy); i++) {
    632 
    633             // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky)
    634             for (int i = 0; i < mask->n; i++) {
    635                 mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    636             }
    637             if (!pmPSFFitShapeParamsMap (psf, i, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) {
    638                 break;
    639             }
    640 
    641             // store the resulting scatterTotal values and the scales, redo the best
    642             if (scatterTotal < scatterTotalMin) {
    643                 scatterTotalMin = scatterTotal;
    644                 entryMin = i;
    645             }
    646         }
    647         if (entryMin == -1) {
    648             psError (PS_ERR_UNKNOWN, false, "failed to find image map for shape params");
    649             return false;
    650         }
    651 
    652         // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky)
    653         for (int i = 0; i < mask->n; i++) {
    654             mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    655         }
    656         if (!pmPSFFitShapeParamsMap (psf, entryMin, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) {
    657             psAbort ("failed pmPSFFitShapeParamsMap on second pass?");
    658         }
    659 
    660         pmTrend2D *trend = psf->params->data[PM_PAR_E0];
    661         psf->trendNx = trend->map->map->numCols;
    662         psf->trendNy = trend->map->map->numRows;
    663 
    664         // copy mask back to srcMask
    665         for (int i = 0; i < mask->n; i++) {
    666             srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    667         }
    668 
    669         psFree (mask);
    670         psFree (dz);
    671     } else {
    672 
    673         // XXX iterate Nx, Ny based on scatter?
    674         // XXX we force the x & y order to be the same
    675         // modify the order to correspond to the actual number of matched stars:
    676         int order = PS_MAX (psf->trendNx, psf->trendNy);
    677         if ((sources->n < 15) && (order >= 3)) order = 2;
    678         if ((sources->n < 11) && (order >= 2)) order = 1;
    679         if ((sources->n <  8) && (order >= 1)) order = 0;
    680         if ((sources->n <  3)) {
    681             psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params");
    682             return false;
    683         }
    684         psf->trendNx = order;
    685         psf->trendNy = order;
    686 
    687         // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
    688         // This way, the parameters masked by one of the fits will be applied to the others
    689         for (int i = 0; i < psf->psfTrendStats->clipIter; i++) {
    690 
    691             psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
    692             psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
    693 
    694             pmTrend2D *trend = NULL;
    695             float mean, stdev;
    696 
    697             // XXX we are using the same stats structure on each pass: do we need to re-init it?
    698             bool status = true;
    699 
    700             trend = psf->params->data[PM_PAR_E0];
    701             status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e0, NULL);
    702             mean = psStatsGetValue (trend->stats, meanOption);
    703             stdev = psStatsGetValue (trend->stats, stdevOption);
    704             psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n);
    705             pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask);
    706 
    707             trend = psf->params->data[PM_PAR_E1];
    708             status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e1, NULL);
    709             mean = psStatsGetValue (trend->stats, meanOption);
    710             stdev = psStatsGetValue (trend->stats, stdevOption);
    711             psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n);
    712             pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask);
    713 
    714             trend = psf->params->data[PM_PAR_E2];
    715             status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e2, NULL);
    716             mean = psStatsGetValue (trend->stats, meanOption);
    717             stdev = psStatsGetValue (trend->stats, stdevOption);
    718             psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n);
    719             pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask);
    720 
    721             if (!status) {
    722                 psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params");
    723                 return false;
    724             }
    725         }
    726     }
    727 
    728     // test dump of the psf parameters
    729     if (psTraceGetLevel("psModules.objects") >= 4) {
    730         FILE *f = fopen ("pol.dat", "w");
    731         fprintf (f, "# x y  :  e0obs e1obs e2obs  : e0fit e1fit e2fit : mask\n");
    732         for (int i = 0; i < e0->n; i++) {
    733             fprintf (f, "%f %f  :  %f %f %f  : %f %f %f  : %d\n",
    734                      x->data.F32[i], y->data.F32[i],
    735                      e0->data.F32[i], e1->data.F32[i], e2->data.F32[i],
    736                      pmTrend2DEval (psf->params->data[PM_PAR_E0], x->data.F32[i], y->data.F32[i]),
    737                      pmTrend2DEval (psf->params->data[PM_PAR_E1], x->data.F32[i], y->data.F32[i]),
    738                      pmTrend2DEval (psf->params->data[PM_PAR_E2], x->data.F32[i], y->data.F32[i]),
    739                      srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]);
    740         }
    741         fclose (f);
    742     }
    743 
    744     psFree (e0);
    745     psFree (e1);
    746     psFree (e2);
    747     psFree (mag);
    748     return true;
    749 }
    750 
    751 // fit the shape variations as a psImageMap for the given scale factor
    752 bool pmPSFFitShapeParamsMap (pmPSF *psf, int scale, float *scatterTotal, psVector *mask, psVector *x, psVector *y, psVector *mag, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz) {
    753 
    754     int Nx, Ny;
    755 
    756     // set the map scale to match the aspect ratio : for a scale of 1, we guarantee
    757     // that we have a single cell
    758     if (psf->fieldNx > psf->fieldNy) {
    759         Nx = scale;
    760         float AR = psf->fieldNy / (float) psf->fieldNx;
    761         Ny = (int) (Nx * AR + 0.5);
    762         Ny = PS_MAX (1, Ny);
    763     } else {
    764         Ny = scale;
    765         float AR = psf->fieldNx / (float) psf->fieldNy;
    766         Nx = (int) (Ny * AR + 0.5);
    767         Nx = PS_MAX (1, Nx);
    768     }
    769 
    770     // do we have enough sources for this fine of a grid?
    771     if (x->n < 10*Nx*Ny) {
    772         return false;
    773     }
    774 
    775     // XXX check this against the exising type
    776     pmTrend2DMode psfTrendMode = PM_TREND_MAP;
    777 
    778     psImageBinning *binning = psImageBinningAlloc();
    779     binning->nXruff = Nx;
    780     binning->nYruff = Ny;
    781     binning->nXfine = psf->fieldNx;
    782     binning->nYfine = psf->fieldNy;
    783     psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
    784     psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo);
    785 
    786     psFree (psf->params->data[PM_PAR_E0]);
    787     psFree (psf->params->data[PM_PAR_E1]);
    788     psFree (psf->params->data[PM_PAR_E2]);
    789 
    790     int nIter = psf->psfTrendStats->clipIter;
    791     psf->psfTrendStats->clipIter = 1;
    792     psf->params->data[PM_PAR_E0] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats);
    793     psf->params->data[PM_PAR_E1] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats);
    794     psf->params->data[PM_PAR_E2] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats);
    795     psFree (binning);
    796 
    797     // if the map is 1x1 (a single value), we measure the resulting ensemble scatter
    798 
    799     // if the map is more finely sampled, divide the values into two sets: measure the fit from
    800     // one set and the scatter from the other set.
    801     psVector *x_fit = NULL;
    802     psVector *y_fit = NULL;
    803     psVector *x_tst = NULL;
    804     psVector *y_tst = NULL;
    805 
    806     psVector *e0obs_fit = NULL;
    807     psVector *e1obs_fit = NULL;
    808     psVector *e2obs_fit = NULL;
    809     psVector *e0obs_tst = NULL;
    810     psVector *e1obs_tst = NULL;
    811     psVector *e2obs_tst = NULL;
    812 
    813     if (scale == 1) {
    814         x_fit  = psMemIncrRefCounter (x);
    815         y_fit  = psMemIncrRefCounter (y);
    816         x_tst  = psMemIncrRefCounter (x);
    817         y_tst  = psMemIncrRefCounter (y);
    818         e0obs_fit = psMemIncrRefCounter (e0obs);
    819         e1obs_fit = psMemIncrRefCounter (e1obs);
    820         e2obs_fit = psMemIncrRefCounter (e2obs);
    821         e0obs_tst = psMemIncrRefCounter (e0obs);
    822         e1obs_tst = psMemIncrRefCounter (e1obs);
    823         e2obs_tst = psMemIncrRefCounter (e2obs);
    824     } else {
    825         x_fit  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    826         y_fit  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    827         x_tst  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    828         y_tst  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    829         e0obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    830         e1obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    831         e2obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    832         e0obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    833         e1obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    834         e2obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    835         for (int i = 0; i < e0obs_fit->n; i++) {
    836             // e0obs->n ==  8 or 9:
    837             // i = 0, 1, 2, 3 : 2i+0 = 0, 2, 4, 6
    838             // i = 0, 1, 2, 3 : 2i+1 = 1, 3, 5, 7
    839             x_fit->data.F32[i] = x->data.F32[2*i+0];
    840             x_tst->data.F32[i] = x->data.F32[2*i+1];
    841             y_fit->data.F32[i] = y->data.F32[2*i+0];
    842             y_tst->data.F32[i] = y->data.F32[2*i+1];
    843 
    844             e0obs_fit->data.F32[i] = e0obs->data.F32[2*i+0];
    845             e0obs_tst->data.F32[i] = e0obs->data.F32[2*i+1];
    846             e1obs_fit->data.F32[i] = e1obs->data.F32[2*i+0];
    847             e1obs_tst->data.F32[i] = e1obs->data.F32[2*i+1];
    848             e2obs_fit->data.F32[i] = e2obs->data.F32[2*i+0];
    849             e2obs_tst->data.F32[i] = e2obs->data.F32[2*i+1];
    850         }
    851     }
    852 
    853     // the mask marks the values not used to calculate the ApTrend
    854     psVector *fitMask = psVectorAlloc (x_fit->n, PS_TYPE_VECTOR_MASK);
    855     // copy mask values to fitMask as a starting point
    856     for (int i = 0; i < fitMask->n; i++) {
    857         fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    858     }
    859 
    860     // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
    861     // This way, the parameters masked by one of the fits will be applied to the others
    862     for (int i = 0; i < nIter; i++) {
    863         // XXX we are using the same stats structure on each pass: do we need to re-init it?
    864         psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
    865         psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
    866 
    867         pmTrend2D *trend = NULL;
    868         float mean, stdev;
    869 
    870         // XXX we are using the same stats structure on each pass: do we need to re-init it?
    871         bool status = true;
    872 
    873         trend = psf->params->data[PM_PAR_E0];
    874         status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e0obs_fit, NULL);
    875         mean = psStatsGetValue (trend->stats, meanOption);
    876         stdev = psStatsGetValue (trend->stats, stdevOption);
    877         psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs_fit->n);
    878         // printTrendMap (trend);
    879         psImageMapCleanup (trend->map);
    880         // printTrendMap (trend);
    881         pmSourceVisualPSFModelResid (trend, x, y, e0obs, mask);
    882 
    883         trend = psf->params->data[PM_PAR_E1];
    884         status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e1obs_fit, NULL);
    885         mean = psStatsGetValue (trend->stats, meanOption);
    886         stdev = psStatsGetValue (trend->stats, stdevOption);
    887         psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs_fit->n);
    888         // printTrendMap (trend);
    889         psImageMapCleanup (trend->map);
    890         // printTrendMap (trend);
    891         pmSourceVisualPSFModelResid (trend, x, y, e1obs, mask);
    892 
    893         trend = psf->params->data[PM_PAR_E2];
    894         status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e2obs_fit, NULL);
    895         mean = psStatsGetValue (trend->stats, meanOption);
    896         stdev = psStatsGetValue (trend->stats, stdevOption);
    897         psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n);
    898         // printTrendMap (trend);
    899         psImageMapCleanup (trend->map);
    900         // printTrendMap (trend);
    901         pmSourceVisualPSFModelResid (trend, x, y, e2obs, mask);
    902     }
    903     psf->psfTrendStats->clipIter = nIter; // restore default setting
    904 
    905     // construct the fitted values and the residuals
    906     psVector *e0fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E0], fitMask, 0xff, x_tst, y_tst);
    907     psVector *e1fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E1], fitMask, 0xff, x_tst, y_tst);
    908     psVector *e2fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E2], fitMask, 0xff, x_tst, y_tst);
    909 
    910     psVector *e0res = (psVector *) psBinaryOp (NULL, (void *) e0obs_tst, "-", (void *) e0fit);
    911     psVector *e1res = (psVector *) psBinaryOp (NULL, (void *) e1obs_tst, "-", (void *) e1fit);
    912     psVector *e2res = (psVector *) psBinaryOp (NULL, (void *) e2obs_tst, "-", (void *) e2fit);
    913 
    914     // measure scatter for the unfitted points
    915     // psTraceSetLevel ("psLib.math.vectorSampleStdev", 10);
    916     // psTraceSetLevel ("psLib.math.vectorClippedStats", 10);
    917     pmPSFShapeParamsScatter (scatterTotal, e0res, e1res, e2res, fitMask, 0xff, psStatsStdevOption(psf->psfTrendStats->options));
    918     // psTraceSetLevel ("psLib.math.vectorSampleStdev", 0);
    919     // psTraceSetLevel ("psLib.math.vectorClippedStats", 0);
    920 
    921     psLogMsg ("psphot.psftry", PS_LOG_INFO, "result of %d x %d grid\n", Nx, Ny);
    922     psLogMsg ("psphot.psftry", PS_LOG_INFO, "systematic scatter: %f\n", *scatterTotal);
    923 
    924     psFree (x_fit);
    925     psFree (y_fit);
    926     psFree (x_tst);
    927     psFree (y_tst);
    928 
    929     psFree (e0obs_fit);
    930     psFree (e1obs_fit);
    931     psFree (e2obs_fit);
    932     psFree (e0obs_tst);
    933     psFree (e1obs_tst);
    934     psFree (e2obs_tst);
    935 
    936     psFree (e0fit);
    937     psFree (e1fit);
    938     psFree (e2fit);
    939 
    940     psFree (e0res);
    941     psFree (e1res);
    942     psFree (e2res);
    943 
    944     // XXX copy fitMask values back to mask
    945     for (int i = 0; i < fitMask->n; i++) {
    946         mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    947     }
    948     psFree (fitMask);
    949 
    950     return true;
    951 }
    952 
    953 // calculate the scatter of the parameters
    954 bool pmPSFShapeParamsScatter(float *scatterTotal, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, psVectorMaskType maskValue, psStatsOptions stdevOpt)
    955 {
    956 
    957     // psStats *stats = psStatsAlloc(stdevOpt);
    958     psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_STDEV);
    959 
    960     // calculate the root-mean-square of E0, E1, E2
    961     float dEsquare = 0.0;
    962     psStatsInit (stats);
    963     if (!psVectorStats (stats, e0res, NULL, mask, maskValue)) {
    964         psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    965         return false;
    966     }
    967     dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
    968 
    969     psStatsInit (stats);
    970     if (!psVectorStats (stats, e1res, NULL, mask, maskValue)) {
    971         psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    972         return false;
    973     }
    974     dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
    975 
    976     psStatsInit (stats);
    977     if (!psVectorStats (stats, e2res, NULL, mask, maskValue)) {
    978         psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    979         return false;
    980     }
    981     dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
    982 
    983     *scatterTotal = sqrtf(dEsquare);
    984 
    985     psFree(stats);
    986     return true;
    987 }
    988 
    989 // calculate the minimum scatter of the parameters
    990 bool pmPSFShapeParamsErrors(float *errorFloor, psVector *mag, psVector *e0res, psVector *e1res,
    991                             psVector *e2res, psVector *mask, int nGroup, psStatsOptions stdevOpt)
    992 {
    993 
    994     psStats *statsS = psStatsAlloc(stdevOpt);
    995 
    996     // measure the trend in bins with 10 values each; if < 10 total, use them all
    997     int nBin = PS_MAX (mag->n / nGroup, 1);
    998 
    999     // use mag to group parameters in sequence
    1000     psVector *index = psVectorSortIndex (NULL, mag);
    1001 
    1002     // subset vectors for mag and dap values within the given range
    1003     psVector *dE0subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
    1004     psVector *dE1subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
    1005     psVector *dE2subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
    1006     psVector *mkSubset  = psVectorAllocEmpty (nGroup, PS_TYPE_VECTOR_MASK);
    1007 
    1008     int n = 0;
    1009     float min = INFINITY;               // Minimum error
    1010     for (int i = 0; i < nBin; i++) {
    1011         int j;
    1012         int nValid = 0;
    1013         for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) {
    1014             int N = index->data.U32[n];
    1015             dE0subset->data.F32[j] = e0res->data.F32[N];
    1016             dE1subset->data.F32[j] = e1res->data.F32[N];
    1017             dE2subset->data.F32[j] = e2res->data.F32[N];
    1018 
    1019             mkSubset->data.PS_TYPE_VECTOR_MASK_DATA[j]   = mask->data.PS_TYPE_VECTOR_MASK_DATA[N];
    1020             if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[N]) nValid ++;
    1021         }
    1022         if (nValid < 3) continue;
    1023 
    1024         dE0subset->n = j;
    1025         dE1subset->n = j;
    1026         dE2subset->n = j;
    1027         mkSubset->n = j;
    1028 
    1029         // calculate the root-mean-square of E0, E1, E2
    1030         float dEsquare = 0.0;
    1031         psStatsInit (statsS);
    1032         if (!psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff)) {
    1033         }
    1034         dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
    1035 
    1036         psStatsInit (statsS);
    1037         if (!psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff)) {
    1038             psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    1039             return false;
    1040         }
    1041         dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
    1042 
    1043         psStatsInit (statsS);
    1044         if (!psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff)) {
    1045             psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    1046             return false;
    1047         }
    1048         dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
    1049 
    1050         if (isfinite(dEsquare)) {
    1051             float err = sqrtf(dEsquare);
    1052             if (err < min) {
    1053                 min = err;
    1054             }
    1055         }
    1056     }
    1057     *errorFloor = min;
    1058 
    1059     psFree (dE0subset);
    1060     psFree (dE1subset);
    1061     psFree (dE2subset);
    1062     psFree (mkSubset);
    1063 
    1064     psFree(index);
    1065 
    1066     psFree(statsS);
    1067 
    1068     return true;
    1069 }
Note: See TracChangeset for help on using the changeset viewer.