Index: /trunk/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- /trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 24205) +++ /trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 24206) @@ -33,5 +33,5 @@ -//#define TESTING // Enable test output +#define TESTING // Enable test output #define PEAK_FLUX 1.0e4 // Peak flux for each source #define SKY_VALUE 0.0e0 // Sky value for fake image @@ -327,4 +327,5 @@ options->psfFieldXo = 0; options->psfFieldYo = 0; + options->chiFluxTrend = false; // All sources have similar flux, so fitting a trend often fails pmSourceFitModelInit(SOURCE_FIT_ITERATIONS, 0.01, VARIANCE_VAL, true); Index: /trunk/psModules/src/objects/pmPSF.c =================================================================== --- /trunk/psModules/src/objects/pmPSF.c (revision 24205) +++ /trunk/psModules/src/objects/pmPSF.c (revision 24206) @@ -75,4 +75,6 @@ options->poissonErrorsParams = true; + options->chiFluxTrend = true; + return options; } Index: /trunk/psModules/src/objects/pmPSF.h =================================================================== --- /trunk/psModules/src/objects/pmPSF.h (revision 24205) +++ /trunk/psModules/src/objects/pmPSF.h (revision 24206) @@ -82,4 +82,5 @@ bool poissonErrorsParams; ///< use poission errors for model parameter fitting float radius; + bool chiFluxTrend; // Fit a trend in Chi2 as a function of flux? } pmPSFOptions; Index: /trunk/psModules/src/objects/pmPSFtry.c =================================================================== --- /trunk/psModules/src/objects/pmPSFtry.c (revision 24205) +++ /trunk/psModules/src/objects/pmPSFtry.c (revision 24206) @@ -43,12 +43,12 @@ 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"); + 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; @@ -63,9 +63,9 @@ 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]; - } + 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]; + } } @@ -75,9 +75,9 @@ // 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; - } + 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; @@ -174,12 +174,12 @@ pmSource *source = psfTry->sources->data[i]; - if (!source->moments) { + if (!source->moments) { psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - continue; - } - if (!source->moments->nPixels) { + continue; + } + if (!source->moments->nPixels) { psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - continue; - } + continue; + } source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type); @@ -311,19 +311,22 @@ // linear clipped fit of chisq trend vs flux - 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; + 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; + } } @@ -600,7 +603,7 @@ 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; - + // 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 @@ -628,8 +631,8 @@ 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]; - } + // 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; @@ -647,8 +650,8 @@ } - // 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]; - } + // 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?"); @@ -659,8 +662,8 @@ 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]; - } + // 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); @@ -686,33 +689,33 @@ 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; + 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]; + 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); + 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]; + 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); + 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]; + 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); + 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); + pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask); if (!status) { @@ -755,11 +758,11 @@ if (psf->fieldNx > psf->fieldNy) { Nx = scale; - float AR = psf->fieldNy / (float) psf->fieldNx; - Ny = (int) (Nx * AR + 0.5); + 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); + float AR = psf->fieldNx / (float) psf->fieldNy; + Nx = (int) (Ny * AR + 0.5); Nx = PS_MAX (1, Nx); } @@ -809,41 +812,41 @@ 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); + 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]; - } + 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]; + } } @@ -852,5 +855,5 @@ // 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]; + fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i]; } @@ -859,42 +862,42 @@ 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); + 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 @@ -941,5 +944,5 @@ // 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]; + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i]; } psFree (fitMask); @@ -959,6 +962,6 @@ psStatsInit (stats); if (!psVectorStats (stats, e0res, NULL, mask, maskValue)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; } dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt)); @@ -966,6 +969,6 @@ psStatsInit (stats); if (!psVectorStats (stats, e1res, NULL, mask, maskValue)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; } dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt)); @@ -973,6 +976,6 @@ psStatsInit (stats); if (!psVectorStats (stats, e2res, NULL, mask, maskValue)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; } dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt)); @@ -1007,5 +1010,5 @@ for (int i = 0; i < nBin; i++) { int j; - int nValid = 0; + int nValid = 0; for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) { int N = index->data.U32[n]; @@ -1015,7 +1018,7 @@ 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; + if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[N]) nValid ++; + } + if (nValid < 3) continue; dE0subset->n = j; @@ -1028,19 +1031,19 @@ 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)); + 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; - } + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt));