Changeset 27840 for branches/simtest_nebulous_branches/psModules/src/objects/pmPSFtry.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/objects/pmPSFtry.c (modified) (3 diffs)

branches/simtest_nebulous_branches/psModules/src/objects/pmPSFtry.c

@@ -37 +37 @@
 #include "pmSourceVisual.h"
 
-bool printTrendMap (pmTrend2D *trend) {
-
-    if (!trend->map) return false;
-    if (!trend->map->map) return false;
-
-    for (int j = 0; j < trend->map->map->numRows; j++) {
-        for (int i = 0; i < trend->map->map->numCols; i++) {
-            fprintf (stderr, "%5.2f  ", trend->map->map->data.F32[j][i]);
-        }
-        fprintf (stderr, "\t\t\t");
-        for (int i = 0; i < trend->map->map->numCols; i++) {
-            fprintf (stderr, "%5.2f  ", trend->map->error->data.F32[j][i]);
-        }
-        fprintf (stderr, "\n");
-    }
-    return true;
-}
-
-bool psImageMapCleanup (psImageMap *map) {
-
-    if ((map->map->numRows == 1) && (map->map->numCols == 1)) return true;
-
-    // find the weighted average of all pixels
-    float Sum = 0.0;
-    float Wt = 0.0;
-    for (int j = 0; j < map->map->numRows; j++) {
-        for (int i = 0; i < map->map->numCols; i++) {
-            if (!isfinite(map->map->data.F32[j][i])) continue;
-            Sum += map->map->data.F32[j][i] * map->error->data.F32[j][i];
-            Wt += map->error->data.F32[j][i];
-        }
-    }
-
-    float Mean = Sum / Wt;
-
-    // do any of the pixels in the map need to be repaired?
-    // XXX for now, we are just replacing bad pixels with the Mean
-    for (int j = 0; j < map->map->numRows; j++) {
-        for (int i = 0; i < map->map->numCols; i++) {
-            if (isfinite(map->map->data.F32[j][i]) &&
-                (map->error->data.F32[j][i] > 0.0)) continue;
-            map->map->data.F32[j][i] = Mean;
-        }
-    }
-    return true;
-}
-
 // ********  pmPSFtry functions  **************************************************
 // * pmPSFtry holds a single pmPSF model test, with the input sources, the freely
@@ -110 +63 @@
     psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree);
 
-    test->psf       = pmPSFAlloc (options);
+    test->psf       = NULL; 
     test->metric    = psVectorAlloc (sources->n, PS_TYPE_F32);
     test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F32);
@@ -136 +89 @@
 }
 
+float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction) {
 
-// build a pmPSFtry for the given model:
-// - fit each source with the free-floating model
-// - construct the pmPSF from the collection of models
-// - fit each source with the PSF-parameter models
-// - measure the pmPSF quality metric (dApResid)
+    psAssert(residuals, "residuals cannot be NULL");
+    psAssert(errors, "errors cannot be NULL");
+    psAssert(residuals->n == errors->n, "residuals and errors must be the same length");
 
-// sources used in for pmPSFtry may be masked by the analysis
-// mask values indicate the reason the source was rejected:
+    // given a vector of residuals and their formal errors, calculated the necessary systematic
+    // error needed to yield a reduced chisq of 1.0, after first tossing out the clipFraction
+    // highest chi-square contributors (allowed outliers)
 
-// generate a pmPSFtry with a copy of the test PSF sources
-pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal)
-{
-    bool status;
-    int Next = 0;
-    int Npsf = 0;
+    psVector *mask  = psVectorAlloc(residuals->n, PS_TYPE_VECTOR_MASK);
+    psVector *chisq = psVectorAlloc(residuals->n, PS_TYPE_F32);
 
-    // validate the requested model name
-    options->type = pmModelClassGetType (modelName);
-    if (options->type == -1) {
-        psError (PS_ERR_UNKNOWN, true, "invalid model name %s", modelName);
-        return NULL;
+    // calculate the chisq vector:
+    int Ngood = 0;
+    for (int i = 0; i < residuals->n; i++) {
+        chisq->data.F32[i] = PS_MAX_F32;
+        if (!isfinite(residuals->data.F32[i])) continue;
+        if (!isfinite(errors->data.F32[i])) continue;
+        if (errors->data.F32[i] <= 0.0) continue;
+        chisq->data.F32[i] = PS_SQR(residuals->data.F32[i] / errors->data.F32[i]);
+        Ngood ++;
     }
 
-    pmPSFtry *psfTry = pmPSFtryAlloc (sources, options);
-    if (psfTry == NULL) {
-        psError (PS_ERR_UNKNOWN, false, "failed to allocate psf model");
-        return NULL;
+    psVector *index = psVectorSortIndex(NULL, chisq);
+
+    // toss out the clipFraction highest chisq values
+    for (int i = 0; i < residuals->n; i++) {
+        int n = index->data.S32[i];
+        if (i < (1.0 - clipFraction)*Ngood) {
+            mask->data.PS_TYPE_VECTOR_MASK_DATA[n] = 0;
+        } else {
+            mask->data.PS_TYPE_VECTOR_MASK_DATA[n] = 1;
+        }
     }
 
-    // maskVal is used to test for rejected pixels, and must include markVal
-    maskVal |= markVal;
+    // Ndof ~= Ngood
+    // Chisq_Ndof = sum(residuals_i^2 / error_i^2) / Ndof
+    // choose S2 such than Chisq^sys_Ndof = sum(residuals_i^2 / (error_i^2 + S2)) / Ndof = 1.0
+    
+    // use Newton-Raphson to solve for S2:
 
-    // stage 1:  fit an EXT model to all candidates PSF sources -- this is independent of the modeled 2D variations in the PSF
-    psTimerStart ("psf.fit");
-    for (int i = 0; i < psfTry->sources->n; i++) {
+    // use median sigma to calculate the initial guess for S2:
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
+    psVectorStats (stats, errors, NULL, mask, 1);
+    float errorMedian = stats->sampleMedian;
+    
+    float nPts = 0.0;
+    float res2mean = 0.0;
+    float ChiSq = 0.0;
+    for (int i = 0; i < residuals->n; i++) {
+        int n = index->data.S32[i];
+        if (mask->data.PS_TYPE_VECTOR_MASK_DATA[n]) continue;
+        res2mean += PS_SQR(residuals->data.F32[n]);
+        ChiSq += PS_SQR(residuals->data.F32[n]) / PS_SQR(errors->data.F32[n]);
+        nPts += 1.0;
+    }
+    res2mean /= nPts;
+    ChiSq /= nPts;
+    
+    float S2guess = res2mean - PS_SQR(errorMedian);
 
-        pmSource *source = psfTry->sources->data[i];
-        if (!source->moments) {
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
-            continue;
-        }
-        if (!source->moments->nPixels) {
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
-            continue;
-        }
+    psLogMsg ("psModules", 10, "ChiSquare: %f, Ntotal: %ld, Ngood: %d, Nkeep: %.0f, S2 guess: %f\n", 
+              ChiSq, residuals->n, Ngood, nPts, S2guess);
 
-        source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type);
-        if (source->modelEXT == NULL) {
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
-            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y);
-            continue;
-        }
+    for (int iter = 0; iter < 10; iter++) {
 
-        // set object mask to define valid pixels
-        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal);
+        ChiSq = 0.0;
+        float dRdS = 0.0;
+        for (int i = 0; i < residuals->n; i++) {
+            int n = index->data.S32[i];
+            if (mask->data.PS_TYPE_VECTOR_MASK_DATA[n]) continue;
+            float error2 = PS_SQR(errors->data.F32[n]) + S2guess;
+            ChiSq += PS_SQR(residuals->data.F32[n]) / error2;
+            dRdS += PS_SQR(residuals->data.F32[n]) / PS_SQR(error2);
+        }
+        ChiSq /= nPts;
+        dRdS /= nPts;
 
-        // fit model as EXT, not PSF
-        status = pmSourceFitModel (source, source->modelEXT, PM_SOURCE_FIT_EXT, maskVal);
+        // Note the sign on dS: dRdS above is -1 * dR/dS formally
+        float dS = (ChiSq - 1.0) / dRdS;
+        S2guess += dS;
+        S2guess = PS_MAX(0.0, S2guess);
 
-        // clear object mask to define valid pixels
-        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
-
-        // exclude the poor fits
-        if (!status) {
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
-            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y);
-            continue;
-        }
-        Next ++;
-    }
-    psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit ext:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Next, sources->n);
-    psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, sources->n);
-
-    if (Next == 0) {
-        psError(PS_ERR_UNKNOWN, false, "No sources with good extended fits from which to determine PSF.");
-        psFree(psfTry);
-        return NULL;
+        psLogMsg ("psModules", 10, "ChiSquare: %f, dS: %f, S2 guess: %f\n", ChiSq, dS, S2guess);
     }
 
-    // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation
-    if (!pmPSFFromPSFtry (psfTry)) {
-        psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources");
-        psFree(psfTry);
-        return NULL;
-    }
+    // free local allocations
+    psFree (mask);
+    psFree (chisq);
+    psFree (stats);
+    psFree (index);
 
-    // stage 3: refit with fixed shape parameters
-    psTimerStart ("psf.fit");
-    for (int i = 0; i < psfTry->sources->n; i++) {
-
-        pmSource *source = psfTry->sources->data[i];
-
-        // masked for: bad model fit, outlier in parameters
-        if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) {
-            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : source is masked\n", i, source->peak->x, source->peak->y);
-            continue;
-        }
-
-        // set shape for this model based on PSF
-        source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf);
-        if (source->modelPSF == NULL) {
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL;
-            abort();
-            continue;
-        }
-        source->modelPSF->radiusFit = options->radius;
-
-        // set object mask to define valid pixels
-        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal);
-
-        // fit the PSF model to the source
-        status = pmSourceFitModel (source, source->modelPSF, PM_SOURCE_FIT_PSF, maskVal);
-
-        // skip poor fits
-        if (!status) {
-            psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL;
-            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y);
-            continue;
-        }
-
-        status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal);
-        if (!status || isnan(source->apMag)) {
-            psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT;
-            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y);
-            continue;
-        }
-
-        // clear object mask to define valid pixels
-        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
-
-        psfTry->fitMag->data.F32[i] = source->psfMag;
-        psfTry->metric->data.F32[i] = source->apMag - source->psfMag;
-        psfTry->metricErr->data.F32[i] = source->errMag;
-
-        psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n);
-        Npsf ++;
-    }
-    psfTry->psf->nPSFstars = Npsf;
-
-    psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Npsf, sources->n);
-    psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, sources->n);
-
-    if (Npsf == 0) {
-        psError(PS_ERR_UNKNOWN, false, "No sources with good PSF fits after model is built.");
-        psFree(psfTry);
-        return NULL;
-    }
-
-    // measure the chi-square trend as a function of flux (PAR[PM_PAR_I0])
-    // this should be linear for Poisson errors and quadratic for constant sky errors
-    psVector *flux  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
-    psVector *chisq = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
-    psVector *mask  = psVectorAlloc (psfTry->sources->n, PS_TYPE_VECTOR_MASK);
-
-    // generate the x and y vectors, and mask missing models
-    for (int i = 0; i < psfTry->sources->n; i++) {
-        pmSource *source = psfTry->sources->data[i];
-        if (source->modelPSF == NULL) {
-            flux->data.F32[i] = 0.0;
-            chisq->data.F32[i] = 0.0;
-            mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
-        } else {
-            flux->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_I0];
-            chisq->data.F32[i] = source->modelPSF->chisq / source->modelPSF->nDOF;
-            mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
-        }
-    }
-
-    // use 3hi/3lo sigma clipping on the chisq fit
-    psStats *stats = options->stats;
-
-    // linear clipped fit of chisq trend vs flux
-    if (options->chiFluxTrend) {
-        bool result = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, options->stats,
-                                                  mask, 0xff, chisq, NULL, flux);
-        psStatsOptions meanStat = psStatsMeanOption(options->stats->options); // Statistic for mean
-        psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev
-
-        psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n",
-                  psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat));
-
-        psFree(flux);
-        psFree(mask);
-        psFree(chisq);
-
-        if (!result) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to fit psf->ChiTrend");
-            psFree(psfTry);
-            return NULL;
-        }
-    }
-
-    for (int i = 0; i < psfTry->psf->ChiTrend->nX + 1; i++) {
-        psLogMsg ("pmPSFtry", 4, "chisq vs flux fit term %d: %f +/- %f\n", i,
-                  psfTry->psf->ChiTrend->coeff[i]*pow(10000, i),
-                  psfTry->psf->ChiTrend->coeffErr[i]*pow(10000,i));
-    }
-
-    // XXX this function wants aperture radius for pmSourcePhotometry
-    if (!pmPSFtryMetric (psfTry, options)) {
-        psError(PS_ERR_UNKNOWN, false, "Attempt to fit PSF with model %s failed.", modelName);
-        psFree (psfTry);
-        return NULL;
-    }
-
-    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n",
-              modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias);
-
-    return (psfTry);
+    return (sqrt(S2guess));
 }
 
-bool pmPSFtryMetric (pmPSFtry *psfTry, pmPSFOptions *options)
-{
-    PS_ASSERT_PTR_NON_NULL(psfTry, false);
-    PS_ASSERT_PTR_NON_NULL(options, false);
-    PS_ASSERT_PTR_NON_NULL(psfTry->sources, false);
-
-    float RADIUS = options->radius;
-
-    // the measured (aperture - fit) magnitudes (dA == psfTry->metric)
-    //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
-    //     dA = dAo + dsky/flux
-    //   where flux is the flux of the star
-    // we fit this trend to find the infinite flux aperture correction (dAo),
-    //   the nominal sky bias (dsky), and the error on dAo
-    // the values of dA are contaminated by stars with close neighbors in the aperture
-    //   we use an outlier rejection to avoid this bias
-
-    // r2rflux = radius^2 * ten(0.4*fitMag);
-    psVector *r2rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
-
-    for (int i = 0; i < psfTry->sources->n; i++) {
-        if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL)
-            continue;
-        r2rflux->data.F32[i] = PS_SQR(RADIUS) * pow(10.0, 0.4*psfTry->fitMag->data.F32[i]);
-    }
-
-    // XXX test dump of aperture residual data
-    if (psTraceGetLevel("psModules.objects") >= 5) {
-        FILE *f = fopen ("apresid.dat", "w");
-        for (int i = 0; i < psfTry->sources->n; i++) {
-            int keep = (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL);
-
-            pmSource *source = psfTry->sources->data[i];
-            float x = source->peak->x;
-            float y = source->peak->y;
-
-            fprintf (f, "%d  %d, %f %f %f  %f %f %f \n",
-                     i, keep, x, y,
-                     psfTry->fitMag->data.F32[i],
-                     r2rflux->data.F32[i],
-                     psfTry->metric->data.F32[i],
-                     psfTry->metricErr->data.F32[i]);
-        }
-        fclose (f);
-    }
-
-    // This analysis of the apResid statistics is only approximate.  The fitted magnitudes
-    // measured at this point (in the PSF fit) use Poisson errors, and are thus biased as a
-    // function of magnitude.  We re-measure the apResid statistics later in psphot using the
-    // linear, constant-error fitting.  Do not reject outliers with excessive vigor here.
-
-    // fit ApTrend only to r2rflux, ignore x,y,flux variations for now
-    // linear clipped fit of ApResid to r2rflux
-    psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
-    poly->coeffMask[1] = PS_POLY_MASK_SET; // fit only a constant offset (no SKYBIAS)
-
-    // XXX replace this with a psVectorStats call?  since we are not fitting the trend
-    bool result = psVectorClipFitPolynomial1D(poly, options->stats, psfTry->mask, PSFTRY_MASK_ALL,
-                                              psfTry->metric, psfTry->metricErr, r2rflux);
-    if (!result) {
-        psError(PS_ERR_UNKNOWN, false, "Failed to fit clipped poly");
-
-        psFree(poly);
-        psFree(r2rflux);
-
-        return false;
-    }
-    psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev
-    psLogMsg ("pmPSFtryMetric", 4, "apresid: %f +/- %f; from statistics of %ld psf stars\n", poly->coeff[0],
-              psStatsGetValue(options->stats, stdevStat), psfTry->sources->n);
-
-
-
-    // XXX test dump of fitted model (dump when tracing?)
-    if (psTraceGetLevel("psModules.objects") >= 4) {
-        FILE *f = fopen ("resid.dat", "w");
-        psVector *apfit = psPolynomial1DEvalVector (poly, r2rflux);
-        for (int i = 0; i < psfTry->sources->n; i++) {
-            int keep = (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL);
-
-            pmSource *source = psfTry->sources->data[i];
-            float x = source->peak->x;
-            float y = source->peak->y;
-
-            fprintf (f, "%d  %d, %f %f %f  %f %f %f  %f\n",
-                     i, keep, x, y,
-                     psfTry->fitMag->data.F32[i],
-                     r2rflux->data.F32[i],
-                     psfTry->metric->data.F32[i],
-                     psfTry->metricErr->data.F32[i],
-                     apfit->data.F32[i]);
-        }
-        fclose (f);
-        psFree (apfit);
-    }
-
-    // XXX drop the skyBias value, or include above??
-    psfTry->psf->skyBias  = poly->coeff[1];
-    psfTry->psf->ApResid  = poly->coeff[0];
-    psfTry->psf->dApResid = psStatsGetValue(options->stats, stdevStat);
-
-    psFree (r2rflux);
-    psFree (poly);
-
-    return true;
-}
-
-/*
-  (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
-  (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
-  (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
-*/
-
-/*****************************************************************************
-pmPSFFromPSFtry (psfTry): build a PSF model from a collection of
-source->modelEXT entries.  The PSF ignores the first 4 (independent) model
-parameters and constructs a polynomial fit to the remaining as a function of
-image coordinate.
-    input: psfTry with fitted source->modelEXT collection, pre-allocated psf
-Note: some of the array entries may be NULL (failed fits); ignore them.
- *****************************************************************************/
-bool pmPSFFromPSFtry (pmPSFtry *psfTry)
-{
-    PS_ASSERT_PTR_NON_NULL(psfTry, false);
-    PS_ASSERT_PTR_NON_NULL(psfTry->sources, false);
-
-    pmPSF *psf = psfTry->psf;
-
-    // construct the fit vectors from the collection of objects
-    psVector *x  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
-    psVector *y  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
-    psVector *z  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
-
-    // construct the x,y terms
-    for (int i = 0; i < psfTry->sources->n; i++) {
-        pmSource *source = psfTry->sources->data[i];
-        if (source->modelEXT == NULL)
-            continue;
-
-        x->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS];
-        y->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS];
-    }
-
-    if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, psfTry->mask)) {
-        psFree(x);
-        psFree(y);
-        psFree(z);
-        return false;
-    }
-
-    // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY)
-    // apply the values of Nx, Ny determined above for E0,E1,E2 to the remaining parameters
-    for (int i = 0; i < psf->params->n; i++) {
-        switch (i) {
-          case PM_PAR_SKY:
-          case PM_PAR_I0:
-          case PM_PAR_XPOS:
-          case PM_PAR_YPOS:
-          case PM_PAR_SXX:
-          case PM_PAR_SYY:
-          case PM_PAR_SXY:
-            continue;
-          default:
-            break;
-        }
-
-        // select the per-object fitted data for this parameter
-        for (int j = 0; j < psfTry->sources->n; j++) {
-            pmSource *source = psfTry->sources->data[j];
-            if (source->modelEXT == NULL) continue;
-            z->data.F32[j] = source->modelEXT->params->data.F32[i];
-        }
-
-        psImageBinning *binning = psImageBinningAlloc();
-        binning->nXruff = psf->trendNx;
-        binning->nYruff = psf->trendNy;
-        binning->nXfine = psf->fieldNx;
-        binning->nYfine = psf->fieldNy;
-
-        if (psf->psfTrendMode == PM_TREND_MAP) {
-            psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
-            psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo);
-        }
-
-        // free existing trend, re-alloc
-        psFree (psf->params->data[i]);
-        psf->params->data[i] = pmTrend2DNoImageAlloc (psf->psfTrendMode, binning, psf->psfTrendStats);
-        psFree (binning);
-
-        // fit the collection of measured parameters to the PSF 2D model
-        // the mask is carried from previous steps and updated with this operation
-        // the weight is either the flux error or NULL, depending on 'psf->poissonErrorParams'
-        if (!pmTrend2DFit (psf->params->data[i], psfTry->mask, 0xff, x, y, z, NULL)) {
-            psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i);
-            psFree(x);
-            psFree(y);
-            psFree(z);
-            return false;
-        }
-    }
-
-    // test dump of star parameters vs position (compare with fitted values)
-    if (psTraceGetLevel("psModules.objects") >= 4) {
-        FILE *f = fopen ("params.dat", "w");
-
-        for (int j = 0; j < psfTry->sources->n; j++) {
-            pmSource *source = psfTry->sources->data[j];
-            if (source == NULL) continue;
-            if (source->modelEXT == NULL) continue;
-
-            pmModel *modelPSF = pmModelFromPSF (source->modelEXT, psf);
-
-            fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[PM_PAR_XPOS], source->modelEXT->params->data.F32[PM_PAR_YPOS]);
-
-            for (int i = 0; i < psf->params->n; i++) {
-                if (psf->params->data[i] == NULL) continue;
-                fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[i], modelPSF->params->data.F32[i]);
-            }
-            fprintf (f, "%f %d\n", source->modelEXT->chisq, source->modelEXT->nIter);
-
-            psFree(modelPSF);
-        }
-        fclose (f);
-    }
-
-    psFree (x);
-    psFree (y);
-    psFree (z);
-    return true;
-}
-
-
-bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) {
-
-    // we are doing a robust fit.  after each pass, we drop points which are more deviant than
-    // three sigma.  mask is currently updated for each pass.
-
-    // The shape parameters (SXX, SXY, SYY) are strongly coupled.  We have to handle them very
-    // carefully.  First, we convert them to the Ellipse Polarization terms (E0, E1, E2) for
-    // each source and fit this set of parameters.  These values are less tightly coupled, but
-    // are still inter-related.  The fitted values do a good job of constraining the major axis
-    // and the position angle, but the minor axis is weakly measured.  When we apply the PSF
-    // model to construct a source model, we convert the fitted values of E0,E1,E2 to the shape
-    // parameters, with the constraint that the minor axis must be greater than a minimum
-    // threshold.
-
-    // convert the measured source shape paramters to polarization terms
-    psVector *e0   = psVectorAlloc (sources->n, PS_TYPE_F32);
-    psVector *e1   = psVectorAlloc (sources->n, PS_TYPE_F32);
-    psVector *e2   = psVectorAlloc (sources->n, PS_TYPE_F32);
-    psVector *mag  = psVectorAlloc (sources->n, PS_TYPE_F32);
-
-    for (int i = 0; i < sources->n; i++) {
-        // skip any masked sources (failed to fit one of the model steps or get a magnitude)
-        if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
-
-        pmSource *source = sources->data[i];
-        assert (source->modelEXT); // all unmasked sources should have modelEXT
-
-        psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32);
-
-        e0->data.F32[i] = pol.e0;
-        e1->data.F32[i] = pol.e1;
-        e2->data.F32[i] = pol.e2;
-
-        float flux = source->modelEXT->params->data.F32[PM_PAR_I0];
-        mag->data.F32[i] = (flux > 0.0) ? -2.5*log(flux) : -100.0;
-    }
-
-    if (psf->psfTrendMode == PM_TREND_MAP) {
-        float scatterTotal = 0.0;
-        float scatterTotalMin = FLT_MAX;
-        int entryMin = -1;
-
-        psVector *dz = NULL;
-        psVector *mask = psVectorAlloc (sources->n, PS_TYPE_VECTOR_MASK);
-
-        // check the fit residuals and increase Nx,Ny until the error is minimized
-        // pmPSFParamTrend increases the number along the longer of x or y
-        for (int i = 1; i <= PS_MAX (psf->trendNx, psf->trendNy); i++) {
-
-            // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky)
-            for (int i = 0; i < mask->n; i++) {
-                mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
-            }
-            if (!pmPSFFitShapeParamsMap (psf, i, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) {
-                break;
-            }
-
-            // store the resulting scatterTotal values and the scales, redo the best
-            if (scatterTotal < scatterTotalMin) {
-                scatterTotalMin = scatterTotal;
-                entryMin = i;
-            }
-        }
-        if (entryMin == -1) {
-            psError (PS_ERR_UNKNOWN, false, "failed to find image map for shape params");
-            return false;
-        }
-
-        // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky)
-        for (int i = 0; i < mask->n; i++) {
-            mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
-        }
-        if (!pmPSFFitShapeParamsMap (psf, entryMin, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) {
-            psAbort ("failed pmPSFFitShapeParamsMap on second pass?");
-        }
-
-        pmTrend2D *trend = psf->params->data[PM_PAR_E0];
-        psf->trendNx = trend->map->map->numCols;
-        psf->trendNy = trend->map->map->numRows;
-
-        // copy mask back to srcMask
-        for (int i = 0; i < mask->n; i++) {
-            srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i];
-        }
-
-        psFree (mask);
-        psFree (dz);
-    } else {
-
-        // XXX iterate Nx, Ny based on scatter?
-        // XXX we force the x & y order to be the same
-        // modify the order to correspond to the actual number of matched stars:
-        int order = PS_MAX (psf->trendNx, psf->trendNy);
-        if ((sources->n < 15) && (order >= 3)) order = 2;
-        if ((sources->n < 11) && (order >= 2)) order = 1;
-        if ((sources->n <  8) && (order >= 1)) order = 0;
-        if ((sources->n <  3)) {
-            psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params");
-            return false;
-        }
-        psf->trendNx = order;
-        psf->trendNy = order;
-
-        // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
-        // This way, the parameters masked by one of the fits will be applied to the others
-        for (int i = 0; i < psf->psfTrendStats->clipIter; i++) {
-
-            psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
-            psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
-
-            pmTrend2D *trend = NULL;
-            float mean, stdev;
-
-            // XXX we are using the same stats structure on each pass: do we need to re-init it?
-            bool status = true;
-
-            trend = psf->params->data[PM_PAR_E0];
-            status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e0, NULL);
-            mean = psStatsGetValue (trend->stats, meanOption);
-            stdev = psStatsGetValue (trend->stats, stdevOption);
-            psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n);
-            pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask);
-
-            trend = psf->params->data[PM_PAR_E1];
-            status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e1, NULL);
-            mean = psStatsGetValue (trend->stats, meanOption);
-            stdev = psStatsGetValue (trend->stats, stdevOption);
-            psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n);
-            pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask);
-
-            trend = psf->params->data[PM_PAR_E2];
-            status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e2, NULL);
-            mean = psStatsGetValue (trend->stats, meanOption);
-            stdev = psStatsGetValue (trend->stats, stdevOption);
-            psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n);
-            pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask);
-
-            if (!status) {
-                psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params");
-                return false;
-            }
-        }
-    }
-
-    // test dump of the psf parameters
-    if (psTraceGetLevel("psModules.objects") >= 4) {
-        FILE *f = fopen ("pol.dat", "w");
-        fprintf (f, "# x y  :  e0obs e1obs e2obs  : e0fit e1fit e2fit : mask\n");
-        for (int i = 0; i < e0->n; i++) {
-            fprintf (f, "%f %f  :  %f %f %f  : %f %f %f  : %d\n",
-                     x->data.F32[i], y->data.F32[i],
-                     e0->data.F32[i], e1->data.F32[i], e2->data.F32[i],
-                     pmTrend2DEval (psf->params->data[PM_PAR_E0], x->data.F32[i], y->data.F32[i]),
-                     pmTrend2DEval (psf->params->data[PM_PAR_E1], x->data.F32[i], y->data.F32[i]),
-                     pmTrend2DEval (psf->params->data[PM_PAR_E2], x->data.F32[i], y->data.F32[i]),
-                     srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]);
-        }
-        fclose (f);
-    }
-
-    psFree (e0);
-    psFree (e1);
-    psFree (e2);
-    psFree (mag);
-    return true;
-}
-
-// fit the shape variations as a psImageMap for the given scale factor
-bool pmPSFFitShapeParamsMap (pmPSF *psf, int scale, float *scatterTotal, psVector *mask, psVector *x, psVector *y, psVector *mag, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz) {
-
-    int Nx, Ny;
-
-    // set the map scale to match the aspect ratio : for a scale of 1, we guarantee
-    // that we have a single cell
-    if (psf->fieldNx > psf->fieldNy) {
-        Nx = scale;
-        float AR = psf->fieldNy / (float) psf->fieldNx;
-        Ny = (int) (Nx * AR + 0.5);
-        Ny = PS_MAX (1, Ny);
-    } else {
-        Ny = scale;
-        float AR = psf->fieldNx / (float) psf->fieldNy;
-        Nx = (int) (Ny * AR + 0.5);
-        Nx = PS_MAX (1, Nx);
-    }
-
-    // do we have enough sources for this fine of a grid?
-    if (x->n < 10*Nx*Ny) {
-        return false;
-    }
-
-    // XXX check this against the exising type
-    pmTrend2DMode psfTrendMode = PM_TREND_MAP;
-
-    psImageBinning *binning = psImageBinningAlloc();
-    binning->nXruff = Nx;
-    binning->nYruff = Ny;
-    binning->nXfine = psf->fieldNx;
-    binning->nYfine = psf->fieldNy;
-    psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
-    psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo);
-
-    psFree (psf->params->data[PM_PAR_E0]);
-    psFree (psf->params->data[PM_PAR_E1]);
-    psFree (psf->params->data[PM_PAR_E2]);
-
-    int nIter = psf->psfTrendStats->clipIter;
-    psf->psfTrendStats->clipIter = 1;
-    psf->params->data[PM_PAR_E0] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats);
-    psf->params->data[PM_PAR_E1] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats);
-    psf->params->data[PM_PAR_E2] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats);
-    psFree (binning);
-
-    // if the map is 1x1 (a single value), we measure the resulting ensemble scatter
-
-    // if the map is more finely sampled, divide the values into two sets: measure the fit from
-    // one set and the scatter from the other set.
-    psVector *x_fit = NULL;
-    psVector *y_fit = NULL;
-    psVector *x_tst = NULL;
-    psVector *y_tst = NULL;
-
-    psVector *e0obs_fit = NULL;
-    psVector *e1obs_fit = NULL;
-    psVector *e2obs_fit = NULL;
-    psVector *e0obs_tst = NULL;
-    psVector *e1obs_tst = NULL;
-    psVector *e2obs_tst = NULL;
-
-    if (scale == 1) {
-        x_fit  = psMemIncrRefCounter (x);
-        y_fit  = psMemIncrRefCounter (y);
-        x_tst  = psMemIncrRefCounter (x);
-        y_tst  = psMemIncrRefCounter (y);
-        e0obs_fit = psMemIncrRefCounter (e0obs);
-        e1obs_fit = psMemIncrRefCounter (e1obs);
-        e2obs_fit = psMemIncrRefCounter (e2obs);
-        e0obs_tst = psMemIncrRefCounter (e0obs);
-        e1obs_tst = psMemIncrRefCounter (e1obs);
-        e2obs_tst = psMemIncrRefCounter (e2obs);
-    } else {
-        x_fit  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
-        y_fit  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
-        x_tst  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
-        y_tst  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
-        e0obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
-        e1obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
-        e2obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
-        e0obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
-        e1obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
-        e2obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
-        for (int i = 0; i < e0obs_fit->n; i++) {
-            // e0obs->n ==  8 or 9:
-            // i = 0, 1, 2, 3 : 2i+0 = 0, 2, 4, 6
-            // i = 0, 1, 2, 3 : 2i+1 = 1, 3, 5, 7
-            x_fit->data.F32[i] = x->data.F32[2*i+0];
-            x_tst->data.F32[i] = x->data.F32[2*i+1];
-            y_fit->data.F32[i] = y->data.F32[2*i+0];
-            y_tst->data.F32[i] = y->data.F32[2*i+1];
-
-            e0obs_fit->data.F32[i] = e0obs->data.F32[2*i+0];
-            e0obs_tst->data.F32[i] = e0obs->data.F32[2*i+1];
-            e1obs_fit->data.F32[i] = e1obs->data.F32[2*i+0];
-            e1obs_tst->data.F32[i] = e1obs->data.F32[2*i+1];
-            e2obs_fit->data.F32[i] = e2obs->data.F32[2*i+0];
-            e2obs_tst->data.F32[i] = e2obs->data.F32[2*i+1];
-        }
-    }
-
-    // the mask marks the values not used to calculate the ApTrend
-    psVector *fitMask = psVectorAlloc (x_fit->n, PS_TYPE_VECTOR_MASK);
-    // copy mask values to fitMask as a starting point
-    for (int i = 0; i < fitMask->n; i++) {
-        fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i];
-    }
-
-    // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
-    // This way, the parameters masked by one of the fits will be applied to the others
-    for (int i = 0; i < nIter; i++) {
-        // XXX we are using the same stats structure on each pass: do we need to re-init it?
-        psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
-        psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
-
-        pmTrend2D *trend = NULL;
-        float mean, stdev;
-
-        // XXX we are using the same stats structure on each pass: do we need to re-init it?
-        bool status = true;
-
-        trend = psf->params->data[PM_PAR_E0];
-        status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e0obs_fit, NULL);
-        mean = psStatsGetValue (trend->stats, meanOption);
-        stdev = psStatsGetValue (trend->stats, stdevOption);
-        psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs_fit->n);
-        // printTrendMap (trend);
-        psImageMapCleanup (trend->map);
-        // printTrendMap (trend);
-        pmSourceVisualPSFModelResid (trend, x, y, e0obs, mask);
-
-        trend = psf->params->data[PM_PAR_E1];
-        status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e1obs_fit, NULL);
-        mean = psStatsGetValue (trend->stats, meanOption);
-        stdev = psStatsGetValue (trend->stats, stdevOption);
-        psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs_fit->n);
-        // printTrendMap (trend);
-        psImageMapCleanup (trend->map);
-        // printTrendMap (trend);
-        pmSourceVisualPSFModelResid (trend, x, y, e1obs, mask);
-
-        trend = psf->params->data[PM_PAR_E2];
-        status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e2obs_fit, NULL);
-        mean = psStatsGetValue (trend->stats, meanOption);
-        stdev = psStatsGetValue (trend->stats, stdevOption);
-        psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n);
-        // printTrendMap (trend);
-        psImageMapCleanup (trend->map);
-        // printTrendMap (trend);
-        pmSourceVisualPSFModelResid (trend, x, y, e2obs, mask);
-    }
-    psf->psfTrendStats->clipIter = nIter; // restore default setting
-
-    // construct the fitted values and the residuals
-    psVector *e0fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E0], fitMask, 0xff, x_tst, y_tst);
-    psVector *e1fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E1], fitMask, 0xff, x_tst, y_tst);
-    psVector *e2fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E2], fitMask, 0xff, x_tst, y_tst);
-
-    psVector *e0res = (psVector *) psBinaryOp (NULL, (void *) e0obs_tst, "-", (void *) e0fit);
-    psVector *e1res = (psVector *) psBinaryOp (NULL, (void *) e1obs_tst, "-", (void *) e1fit);
-    psVector *e2res = (psVector *) psBinaryOp (NULL, (void *) e2obs_tst, "-", (void *) e2fit);
-
-    // measure scatter for the unfitted points
-    // psTraceSetLevel ("psLib.math.vectorSampleStdev", 10);
-    // psTraceSetLevel ("psLib.math.vectorClippedStats", 10);
-    pmPSFShapeParamsScatter (scatterTotal, e0res, e1res, e2res, fitMask, 0xff, psStatsStdevOption(psf->psfTrendStats->options));
-    // psTraceSetLevel ("psLib.math.vectorSampleStdev", 0);
-    // psTraceSetLevel ("psLib.math.vectorClippedStats", 0);
-
-    psLogMsg ("psphot.psftry", PS_LOG_INFO, "result of %d x %d grid\n", Nx, Ny);
-    psLogMsg ("psphot.psftry", PS_LOG_INFO, "systematic scatter: %f\n", *scatterTotal);
-
-    psFree (x_fit);
-    psFree (y_fit);
-    psFree (x_tst);
-    psFree (y_tst);
-
-    psFree (e0obs_fit);
-    psFree (e1obs_fit);
-    psFree (e2obs_fit);
-    psFree (e0obs_tst);
-    psFree (e1obs_tst);
-    psFree (e2obs_tst);
-
-    psFree (e0fit);
-    psFree (e1fit);
-    psFree (e2fit);
-
-    psFree (e0res);
-    psFree (e1res);
-    psFree (e2res);
-
-    // XXX copy fitMask values back to mask
-    for (int i = 0; i < fitMask->n; i++) {
-        mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
-    }
-    psFree (fitMask);
-
-    return true;
-}
-
-// calculate the scatter of the parameters
-bool pmPSFShapeParamsScatter(float *scatterTotal, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, psVectorMaskType maskValue, psStatsOptions stdevOpt)
-{
-
-    // psStats *stats = psStatsAlloc(stdevOpt);
-    psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_STDEV);
-
-    // calculate the root-mean-square of E0, E1, E2
-    float dEsquare = 0.0;
-    psStatsInit (stats);
-    if (!psVectorStats (stats, e0res, NULL, mask, maskValue)) {
-        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
-        return false;
-    }
-    dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
-
-    psStatsInit (stats);
-    if (!psVectorStats (stats, e1res, NULL, mask, maskValue)) {
-        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
-        return false;
-    }
-    dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
-
-    psStatsInit (stats);
-    if (!psVectorStats (stats, e2res, NULL, mask, maskValue)) {
-        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
-        return false;
-    }
-    dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
-
-    *scatterTotal = sqrtf(dEsquare);
-
-    psFree(stats);
-    return true;
-}
-
-// calculate the minimum scatter of the parameters
-bool pmPSFShapeParamsErrors(float *errorFloor, psVector *mag, psVector *e0res, psVector *e1res,
-                            psVector *e2res, psVector *mask, int nGroup, psStatsOptions stdevOpt)
-{
-
-    psStats *statsS = psStatsAlloc(stdevOpt);
-
-    // measure the trend in bins with 10 values each; if < 10 total, use them all
-    int nBin = PS_MAX (mag->n / nGroup, 1);
-
-    // use mag to group parameters in sequence
-    psVector *index = psVectorSortIndex (NULL, mag);
-
-    // subset vectors for mag and dap values within the given range
-    psVector *dE0subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
-    psVector *dE1subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
-    psVector *dE2subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
-    psVector *mkSubset  = psVectorAllocEmpty (nGroup, PS_TYPE_VECTOR_MASK);
-
-    int n = 0;
-    float min = INFINITY;               // Minimum error
-    for (int i = 0; i < nBin; i++) {
-        int j;
-        int nValid = 0;
-        for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) {
-            int N = index->data.U32[n];
-            dE0subset->data.F32[j] = e0res->data.F32[N];
-            dE1subset->data.F32[j] = e1res->data.F32[N];
-            dE2subset->data.F32[j] = e2res->data.F32[N];
-
-            mkSubset->data.PS_TYPE_VECTOR_MASK_DATA[j]   = mask->data.PS_TYPE_VECTOR_MASK_DATA[N];
-            if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[N]) nValid ++;
-        }
-        if (nValid < 3) continue;
-
-        dE0subset->n = j;
-        dE1subset->n = j;
-        dE2subset->n = j;
-        mkSubset->n = j;
-
-        // calculate the root-mean-square of E0, E1, E2
-        float dEsquare = 0.0;
-        psStatsInit (statsS);
-        if (!psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff)) {
-        }
-        dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
-
-        psStatsInit (statsS);
-        if (!psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff)) {
-            psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
-            return false;
-        }
-        dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
-
-        psStatsInit (statsS);
-        if (!psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff)) {
-            psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
-            return false;
-        }
-        dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
-
-        if (isfinite(dEsquare)) {
-            float err = sqrtf(dEsquare);
-            if (err < min) {
-                min = err;
-            }
-        }
-    }
-    *errorFloor = min;
-
-    psFree (dE0subset);
-    psFree (dE1subset);
-    psFree (dE2subset);
-    psFree (mkSubset);
-
-    psFree(index);
-
-    psFree(statsS);
-
-    return true;
-}