IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 24206


Ignore:
Timestamp:
May 15, 2009, 5:05:25 PM (17 years ago)
Author:
Paul Price
Message:

Fitting a trend of chi2 as a function of flux was often failing when fitting a PSF for the fake sources in pmPSFEnvelope (because all sources have similar flux). Made this part optional via an entry in pmPSFOptions (default true).

Location:
trunk/psModules/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/imcombine/pmPSFEnvelope.c

    r23975 r24206  
    3333
    3434
    35 //#define TESTING                         // Enable test output
     35#define TESTING                         // Enable test output
    3636#define PEAK_FLUX 1.0e4                 // Peak flux for each source
    3737#define SKY_VALUE 0.0e0                 // Sky value for fake image
    … …  
    327327    options->psfFieldXo = 0;
    328328    options->psfFieldYo = 0;
     329    options->chiFluxTrend = false;      // All sources have similar flux, so fitting a trend often fails
    329330
    330331    pmSourceFitModelInit(SOURCE_FIT_ITERATIONS, 0.01, VARIANCE_VAL, true);

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

    r23487 r24206  
    7575    options->poissonErrorsParams  = true;
    7676
     77    options->chiFluxTrend = true;
     78
    7779    return options;
    7880}

  • trunk/psModules/src/objects/pmPSF.h

    r23487 r24206  
    8282    bool          poissonErrorsParams; ///< use poission errors for model parameter fitting
    8383    float         radius;
     84    bool          chiFluxTrend;         // Fit a trend in Chi2 as a function of flux?
    8485} pmPSFOptions;
    8586

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

    r23989 r24206  
    4343
    4444    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");
     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");
    5353    }
    5454    return true;
    … …  
    6363    float Wt = 0.0;
    6464    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         }
     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        }
    7070    }
    7171
    … …  
    7575    // XXX for now, we are just replacing bad pixels with the Mean
    7676    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         }
     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        }
    8282    }
    8383    return true;
    … …  
    174174
    175175        pmSource *source = psfTry->sources->data[i];
    176         if (!source->moments) {
     176        if (!source->moments) {
    177177            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    178             continue;
    179         }
    180         if (!source->moments->nPixels) {
     178            continue;
     179        }
     180        if (!source->moments->nPixels) {
    181181            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    182             continue;
    183         }
     182            continue;
     183        }
    184184
    185185        source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type);
    … …  
    311311
    312312    // linear clipped fit of chisq trend vs flux
    313     bool result = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, options->stats, mask, 0xff, chisq, NULL, flux);
    314     psStatsOptions meanStat = psStatsMeanOption(options->stats->options); // Statistic for mean
    315     psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev
    316 
    317     psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n",
    318               psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat));
    319 
    320     psFree(flux);
    321     psFree(mask);
    322     psFree(chisq);
    323 
    324     if (!result) {
    325         psError(PS_ERR_UNKNOWN, false, "Failed to fit psf->ChiTrend");
    326         psFree(psfTry);
    327         return NULL;
     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        }
    328331    }
    329332
    … …  
    600603
    601604    for (int i = 0; i < sources->n; i++) {
    602         // skip any masked sources (failed to fit one of the model steps or get a magnitude)
    603         if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
    604        
     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
    605608        pmSource *source = sources->data[i];
    606609        assert (source->modelEXT); // all unmasked sources should have modelEXT
    … …  
    628631        for (int i = 1; i <= PS_MAX (psf->trendNx, psf->trendNy); i++) {
    629632
    630             // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky)
    631             for (int i = 0; i < mask->n; i++) {
    632                 mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    633             }
     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            }
    634637            if (!pmPSFFitShapeParamsMap (psf, i, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) {
    635638                break;
    … …  
    647650        }
    648651
    649         // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky)
    650         for (int i = 0; i < mask->n; i++) {
    651             mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    652         }
     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        }
    653656        if (!pmPSFFitShapeParamsMap (psf, entryMin, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) {
    654657            psAbort ("failed pmPSFFitShapeParamsMap on second pass?");
    … …  
    659662        psf->trendNy = trend->map->map->numRows;
    660663
    661         // copy mask back to srcMask
    662         for (int i = 0; i < mask->n; i++) {
    663             srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    664         }
     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        }
    665668
    666669        psFree (mask);
    … …  
    686689        for (int i = 0; i < psf->psfTrendStats->clipIter; i++) {
    687690
    688             psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
    689             psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
    690 
    691             pmTrend2D *trend = NULL;
    692             float mean, stdev;
     691            psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
     692            psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
     693
     694            pmTrend2D *trend = NULL;
     695            float mean, stdev;
    693696
    694697            // XXX we are using the same stats structure on each pass: do we need to re-init it?
    695698            bool status = true;
    696699
    697             trend = psf->params->data[PM_PAR_E0];
     700            trend = psf->params->data[PM_PAR_E0];
    698701            status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e0, NULL);
    699             mean = psStatsGetValue (trend->stats, meanOption);
    700             stdev = psStatsGetValue (trend->stats, stdevOption);
     702            mean = psStatsGetValue (trend->stats, meanOption);
     703            stdev = psStatsGetValue (trend->stats, stdevOption);
    701704            psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n);
    702             pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask);
    703 
    704             trend = psf->params->data[PM_PAR_E1];
     705            pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask);
     706
     707            trend = psf->params->data[PM_PAR_E1];
    705708            status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e1, NULL);
    706             mean = psStatsGetValue (trend->stats, meanOption);
    707             stdev = psStatsGetValue (trend->stats, stdevOption);
     709            mean = psStatsGetValue (trend->stats, meanOption);
     710            stdev = psStatsGetValue (trend->stats, stdevOption);
    708711            psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n);
    709             pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask);
    710 
    711             trend = psf->params->data[PM_PAR_E2];
     712            pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask);
     713
     714            trend = psf->params->data[PM_PAR_E2];
    712715            status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e2, NULL);
    713             mean = psStatsGetValue (trend->stats, meanOption);
    714             stdev = psStatsGetValue (trend->stats, stdevOption);
     716            mean = psStatsGetValue (trend->stats, meanOption);
     717            stdev = psStatsGetValue (trend->stats, stdevOption);
    715718            psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n);
    716             pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask);
     719            pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask);
    717720
    718721            if (!status) {
    … …  
    755758    if (psf->fieldNx > psf->fieldNy) {
    756759        Nx = scale;
    757         float AR = psf->fieldNy / (float) psf->fieldNx;
    758         Ny = (int) (Nx * AR + 0.5);
     760        float AR = psf->fieldNy / (float) psf->fieldNx;
     761        Ny = (int) (Nx * AR + 0.5);
    759762        Ny = PS_MAX (1, Ny);
    760763    } else {
    761764        Ny = scale;
    762         float AR = psf->fieldNx / (float) psf->fieldNy;
    763         Nx = (int) (Ny * AR + 0.5);
     765        float AR = psf->fieldNx / (float) psf->fieldNy;
     766        Nx = (int) (Ny * AR + 0.5);
    764767        Nx = PS_MAX (1, Nx);
    765768    }
    … …  
    809812
    810813    if (scale == 1) {
    811         x_fit  = psMemIncrRefCounter (x);
    812         y_fit  = psMemIncrRefCounter (y);
    813         x_tst  = psMemIncrRefCounter (x);
    814         y_tst  = psMemIncrRefCounter (y);
    815         e0obs_fit = psMemIncrRefCounter (e0obs);
    816         e1obs_fit = psMemIncrRefCounter (e1obs);
    817         e2obs_fit = psMemIncrRefCounter (e2obs);
    818         e0obs_tst = psMemIncrRefCounter (e0obs);
    819         e1obs_tst = psMemIncrRefCounter (e1obs);
    820         e2obs_tst = psMemIncrRefCounter (e2obs);
     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);
    821824    } else {
    822         x_fit  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    823         y_fit  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    824         x_tst  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    825         y_tst  = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    826         e0obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    827         e1obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    828         e2obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    829         e0obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    830         e1obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    831         e2obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32);
    832         for (int i = 0; i < e0obs_fit->n; i++) {
    833             // e0obs->n ==  8 or 9:
    834             // i = 0, 1, 2, 3 : 2i+0 = 0, 2, 4, 6
    835             // i = 0, 1, 2, 3 : 2i+1 = 1, 3, 5, 7
    836             x_fit->data.F32[i] = x->data.F32[2*i+0]; 
    837             x_tst->data.F32[i] = x->data.F32[2*i+1]; 
    838             y_fit->data.F32[i] = y->data.F32[2*i+0]; 
    839             y_tst->data.F32[i] = y->data.F32[2*i+1]; 
    840 
    841             e0obs_fit->data.F32[i] = e0obs->data.F32[2*i+0]; 
    842             e0obs_tst->data.F32[i] = e0obs->data.F32[2*i+1]; 
    843             e1obs_fit->data.F32[i] = e1obs->data.F32[2*i+0];
    844             e1obs_tst->data.F32[i] = e1obs->data.F32[2*i+1];
    845             e2obs_fit->data.F32[i] = e2obs->data.F32[2*i+0];
    846             e2obs_tst->data.F32[i] = e2obs->data.F32[2*i+1];
    847         }
     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        }
    848851    }
    849852
    … …  
    852855    // copy mask values to fitMask as a starting point
    853856    for (int i = 0; i < fitMask->n; i++) {
    854         fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i];
     857        fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    855858    }
    856859
    … …  
    859862    for (int i = 0; i < nIter; i++) {
    860863        // XXX we are using the same stats structure on each pass: do we need to re-init it?
    861         psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
    862         psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
    863 
    864         pmTrend2D *trend = NULL;
    865         float mean, stdev;
    866 
    867         // XXX we are using the same stats structure on each pass: do we need to re-init it?
    868         bool status = true;
    869 
    870         trend = psf->params->data[PM_PAR_E0];
    871         status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e0obs_fit, NULL);
    872         mean = psStatsGetValue (trend->stats, meanOption);
    873         stdev = psStatsGetValue (trend->stats, stdevOption);
    874         psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs_fit->n);
    875         // printTrendMap (trend);
    876         psImageMapCleanup (trend->map);
    877         // printTrendMap (trend);
    878         pmSourceVisualPSFModelResid (trend, x, y, e0obs, mask);
    879 
    880         trend = psf->params->data[PM_PAR_E1];
    881         status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e1obs_fit, NULL);
    882         mean = psStatsGetValue (trend->stats, meanOption);
    883         stdev = psStatsGetValue (trend->stats, stdevOption);
    884         psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs_fit->n);
    885         // printTrendMap (trend);
    886         psImageMapCleanup (trend->map);
    887         // printTrendMap (trend);
    888         pmSourceVisualPSFModelResid (trend, x, y, e1obs, mask);
    889 
    890         trend = psf->params->data[PM_PAR_E2];
    891         status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e2obs_fit, NULL);
    892         mean = psStatsGetValue (trend->stats, meanOption);
    893         stdev = psStatsGetValue (trend->stats, stdevOption);
    894         psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n);
    895         // printTrendMap (trend);
    896         psImageMapCleanup (trend->map);
    897         // printTrendMap (trend);
    898         pmSourceVisualPSFModelResid (trend, x, y, e2obs, mask);
     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);
    899902    }
    900903    psf->psfTrendStats->clipIter = nIter; // restore default setting
    … …  
    941944    // XXX copy fitMask values back to mask
    942945    for (int i = 0; i < fitMask->n; i++) {
    943         mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
     946        mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i];
    944947    }
    945948    psFree (fitMask);
    … …  
    959962    psStatsInit (stats);
    960963    if (!psVectorStats (stats, e0res, NULL, mask, maskValue)) {
    961         psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    962         return false;
     964        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     965        return false;
    963966    }
    964967    dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
    … …  
    966969    psStatsInit (stats);
    967970    if (!psVectorStats (stats, e1res, NULL, mask, maskValue)) {
    968         psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    969         return false;
     971        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     972        return false;
    970973    }
    971974    dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
    … …  
    973976    psStatsInit (stats);
    974977    if (!psVectorStats (stats, e2res, NULL, mask, maskValue)) {
    975         psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    976         return false;
     978        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     979        return false;
    977980    }
    978981    dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt));
    … …  
    10071010    for (int i = 0; i < nBin; i++) {
    10081011        int j;
    1009         int nValid = 0;
     1012        int nValid = 0;
    10101013        for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) {
    10111014            int N = index->data.U32[n];
    … …  
    10151018
    10161019            mkSubset->data.PS_TYPE_VECTOR_MASK_DATA[j]   = mask->data.PS_TYPE_VECTOR_MASK_DATA[N];
    1017             if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[N]) nValid ++;
    1018         }
    1019         if (nValid < 3) continue;
     1020            if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[N]) nValid ++;
     1021        }
     1022        if (nValid < 3) continue;
    10201023
    10211024        dE0subset->n = j;
    … …  
    10281031        psStatsInit (statsS);
    10291032        if (!psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff)) {
    1030         }
     1033        }
    10311034        dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
    10321035
    10331036        psStatsInit (statsS);
    10341037        if (!psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff)) {
    1035             psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    1036             return false;
    1037         }       
    1038         dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
     1038            psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     1039            return false;
     1040        }
     1041        dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
    10391042
    10401043        psStatsInit (statsS);
    10411044        if (!psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff)) {
    1042             psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    1043             return false;
    1044         }
     1045            psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     1046            return false;
     1047        }
    10451048        dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));
    10461049
Note: See TracChangeset for help on using the changeset viewer.