Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 6204) +++ trunk/psLib/src/math/psStats.c (revision 6215) @@ -14,6 +14,14 @@ * stats->binsize * - * @version $Revision: 1.162 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-26 21:10:22 $ + * XXX: Must do + * nSubsample points + * use ->min and ->max (PS_STAT_USE_RANGE) + * use ->binsize (PS_STAT_USE_BINSIZE) + * + * + * + * + * @version $Revision: 1.163 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-01-26 23:49:11 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -97,16 +105,6 @@ return true; - case PS_STAT_ROBUST_MEAN: - *value = stats->robustMean; - psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); - return true; - case PS_STAT_ROBUST_MEDIAN: *value = stats->robustMedian; - psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); - return true; - - case PS_STAT_ROBUST_MODE: - *value = stats->robustMode; psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); return true; @@ -449,4 +447,6 @@ unmasked element within the specified min/max range). Otherwise, return "false". + +XXX: Can you use psVectorCountPixelMask here? *****************************************************************************/ bool p_psVectorCheckNonEmpty(const psVector* myVector, @@ -503,4 +503,6 @@ number of non-masked pixels in the vector that fall within the min/max range, if specified. + +XXX: Can you use psVectorCountPixelMask here? *****************************************************************************/ psS32 p_psVectorNValues(const psVector* myVector, @@ -766,5 +768,5 @@ PS_VECTOR_PRINT_F32(smooth); } - psTrace(__func__, 4, "---- %s(psVector) end ----\n", __func__); + psTrace(__func__, 4, "---- %s() end ----\n", __func__); return(smooth); } @@ -858,5 +860,5 @@ /****************************************************************************** -p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the +p_psVectorSampleStdevOLD(myVector, maskVector, maskVal, stats): calculates the stdev of the input vector. Inputs @@ -868,4 +870,5 @@ NULL +XXX: remove this *****************************************************************************/ void p_psVectorSampleStdevOLD(const psVector* myVector, @@ -1126,5 +1129,5 @@ p_psMemSetPersistent(statsTmp, true); } else { - // EAM : initialize structure if already allocated + // Initialize structure if already allocated statsTmp->sampleMean = NAN; statsTmp->sampleMedian = NAN; @@ -1132,20 +1135,21 @@ statsTmp->sampleUQ = NAN; statsTmp->sampleLQ = NAN; - statsTmp->robustMean = NAN; statsTmp->robustMedian = NAN; - statsTmp->robustMode = NAN; statsTmp->robustStdev = NAN; statsTmp->robustUQ = NAN; statsTmp->robustLQ = NAN; - statsTmp->robustN50 = -1; // XXX: This is never used - statsTmp->robustNfit = -1; + statsTmp->robustN50 = -1; + statsTmp->fittedMean = NAN; + statsTmp->fittedStdev = NAN; + statsTmp->fittedNfit = -1; statsTmp->clippedMean = NAN; statsTmp->clippedStdev = NAN; - statsTmp->clippedNvalues = -1; // XXX: This is never used + statsTmp->clippedNvalues = -1; statsTmp->clipSigma = 3.0; statsTmp->clipIter = 3; statsTmp->min = NAN; statsTmp->max = NAN; - statsTmp->binsize = NAN; // XXX: This is never used + statsTmp->binsize = NAN; + statsTmp->nSubsample = 100000; } @@ -1688,27 +1692,34 @@ } } - psBool iterate = true; + psS32 iterate = 1; psF32 sigma; psStats *tmpStatsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); - while (iterate) { + while (iterate > 0) { psTrace(__func__, 6, "Iterating on Bin size.\n"); - // - // Determine the bin size of the robust histogram. This is done - // by computing the total range of data values and dividing by 1000.0. - // - rc = p_psVectorMin(myVector, tmpMaskVec, 1, tmpStatsMinMax); - rc|= p_psVectorMax(myVector, tmpMaskVec, 1, tmpStatsMinMax); - if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); - psFree(tmpStatsMinMax); - psFree(tmpMaskVec); - psFree(tmpScalar); - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); - } - psTrace(__func__, 6, "Data min/man is (%.2f, %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max); - psF32 binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / 1000.0f; - psTrace(__func__, 6, "Robust bin size is %.2f\n", binSize); + psF32 binSize = 0.0; + if ((iterate == 1) && (stats->options & PS_STAT_USE_BINSIZE)) { + // Set initial bin size to the specified value. + binSize = stats->binsize; + psTrace(__func__, 6, "Setting initial robust bin size to %.2f\n", binSize); + } else { + // Determine the bin size of the robust histogram. This is done + // by computing the total range of data values and dividing by 1000.0. + + rc = p_psVectorMin(myVector, tmpMaskVec, 1, tmpStatsMinMax); + rc|= p_psVectorMax(myVector, tmpMaskVec, 1, tmpStatsMinMax); + if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); + psFree(tmpStatsMinMax); + psFree(tmpMaskVec); + psFree(tmpScalar); + psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); + return(1); + } + psTrace(__func__, 6, "Data min/man is (%.2f, %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max); + binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / 1000.0f; + psTrace(__func__, 6, "Robust bin size is %.2f\n", binSize); + } + psTrace(__func__, 6, "Initial robust bin size is %.2f\n", binSize); // @@ -1726,6 +1737,5 @@ } // XXX: Set these to the number of unmasked data points? - stats->robustNfit = 0.0; - stats->robustN50 = 0.0; + stats->robustN50 = 0; psFree(tmpStatsMinMax); psFree(tmpMaskVec); @@ -1900,4 +1910,5 @@ sigma = (binHiF32 - binLoF32) / 2.0; psTrace(__func__, 6, "The current sigma is %f.\n", sigma); + stats->robustStdev = sigma; // @@ -1921,8 +1932,8 @@ psFree(robustHistogram); psFree(cumulativeRobustHistogram); - + iterate++; } else { psTrace(__func__, 6, "*************: No more iteration. sigma is %f\n", sigma); - iterate = false; + iterate = 0; } } @@ -1981,11 +1992,7 @@ psTrace(__func__, 6, "The 25 and 75 percent data point exact positions are (%f, %f).\n", binLo25F32, binHi25F32); - // - // - // New algorithm for calculated mean and stdev. - // - // psS32 N50 = 0; for (psS32 i = 0 ; i < myVector->n ; i++) { + // XXX: use maskVal here? if ((0 == tmpMaskVec->data.U8[i]) && (binLo25F32 <= myVector->data.F32[i]) && @@ -1997,163 +2004,199 @@ psTrace(__func__, 6, "The robustN50 is %d.\n", N50); - psF32 dN = (psF32) (0.17 * N50); - if (dN < 1.0) { - dN = 1.0; - } else if (dN > 4.0) { - dN = 4.0; - } - psF32 newBinSize = sigma / dN; - - rc = p_psVectorMin(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax); - rc|= p_psVectorMax(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax); - if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); - psFree(tmpStatsMinMax); - psFree(robustHistogram); - psFree(cumulativeRobustHistogram); - psFree(tmpScalar); - psFree(tmpMaskVec); - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); - } - - numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / newBinSize); - psTrace(__func__, 6, "The new min/max values are (%f, %f).\n", tmpStatsMinMax->min, tmpStatsMinMax->max); - psTrace(__func__, 6, "The new bin size is %f.\n", newBinSize); - psTrace(__func__, 6, "The numBins is %d\n", numBins); - - psHistogram *newHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins); - newHistogram = psVectorHistogram(newHistogram, myVector, errors, tmpMaskVec, maskVal|1); - if (psTraceGetLevel(__func__) >= 8) { - PS_VECTOR_PRINT_F32(newHistogram->nums); - } - - // - // Smooth the resulting histogram with a Gaussian with sigma_s = 1 - // bin. - // - psF32 sigma_s = newBinSize; - - psVector *newHistogramSmoothed = p_psVectorSmoothHistGaussian(newHistogram, sigma_s); - if (psTraceGetLevel(__func__) >= 8) { - PS_VECTOR_PRINT_F32(newHistogramSmoothed); - } - - // - // Find the bin with the peak value in the range 2 sigma of the - // robust histogram median. - // - - psS32 binMin = 0; - psS32 binMax = 0; - tmpScalar->data.F32 = stats->robustMedian - (2.0 * sigma); - if (tmpScalar->data.F32 <= newHistogram->bounds->data.F32[0]) { - binMin = 0; - } else { - binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); - } - - tmpScalar->data.F32 = stats->robustMedian + (2.0 + sigma); - if (tmpScalar->data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) { - binMax = newHistogram->bounds->n-1; - } else { - binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); - } - if ((binMin < 0) || (binMax < 0)) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the +- 2.0 * sigma bins\n"); - psFree(tmpMaskVec); - psFree(robustHistogram); - psFree(cumulativeRobustHistogram); - psFree(tmpStatsMinMax); - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); - } - - psS32 binNum = binMin; - psF32 binMaxNums = newHistogramSmoothed->data.F32[binNum]; - for (psS32 i = binMin+1 ; i <= binMax ; i++) { - if (newHistogramSmoothed->data.F32[i] > binMaxNums) { - binNum = i; - binMaxNums = newHistogramSmoothed->data.F32[i]; - } - } - psTrace(__func__, 6, "The peak bin is %d, with %f data.n", binNum, binMaxNums); - - // - // Fit a Gaussian to the bins in the range 20 sigma of the robust - // histogram median. - // - tmpScalar->data.F32 = stats->robustMedian - (20.0 * sigma); - if (tmpScalar->data.F32 < tmpStatsMinMax->min) { - binMin = 0; - } else { - binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); - // XXX: check for errors here. - } - tmpScalar->data.F32 = stats->robustMedian + (20.0 * sigma); - if (tmpScalar->data.F32 > tmpStatsMinMax->max) { - binMax = newHistogramSmoothed->n - 1; - } else { - binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); - // XXX: check for errors here. - } - psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); - psVector *xTmp = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); - psArray *x = psArrayAlloc((1 + (binMax - binMin))); - stats->robustNfit = 0; - psS32 j = 0; - - for (psS32 i = binMin ; i <= binMax ; i++) { - y->data.F32[j] = newHistogramSmoothed->data.F32[i]; - x->data[j] = (psPtr *) psVectorAlloc(1, PS_TYPE_F32); - ((psVector *) x->data[j])->data.F32[0] = PS_BIN_MIDPOINT(newHistogram, i); - xTmp->data.F32[j] = PS_BIN_MIDPOINT(newHistogram, i); - - stats->robustNfit+= newHistogramSmoothed->data.F32[i]; - j++; - } - if (psTraceGetLevel(__func__) >= 8) { - // XXX: Print the x array somehow. - PS_VECTOR_PRINT_F32(y); - } - - // XXX: Use the min/max routines for this - psF32 minY = FLT_MAX; - psF32 maxY = -FLT_MAX; - for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) { - if (y->data.F32[i] > maxY) { - maxY = y->data.F32[i]; - } - if (y->data.F32[i] < minY) { - minY = y->data.F32[i]; - } - } - // - // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0]) - // XXX: Use the normalize routines for this. - // - for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) { - y->data.F32[i]= (y->data.F32[i] - minY) / (maxY - minY); - } - - // - psMinimization *min = psMinimizationAlloc(100, 0.01); - psVector *params = psVectorAlloc(2, PS_TYPE_F32); - // Initial guess for the mean ([0]) and standard dev. - params->data.F32[0] = stats->robustMedian; - params->data.F32[1] = sigma; - if (psTraceGetLevel(__func__) >= 8) { - PS_VECTOR_PRINT_F32(params); - PS_VECTOR_PRINT_F32(y); - } - psFree(xTmp); - rcBool = psMinimizeLMChi2(min, NULL, params, NULL, x, y, NULL, psMinimizeLMChi2Gauss1D); - - if (rcBool != true) { - psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); - psFree(tmpStatsMinMax); - psFree(robustHistogram); - psFree(cumulativeRobustHistogram); - psFree(tmpScalar); + // ************************************************************************ + if ((stats->options & PS_STAT_FITTED_MEAN) || + (stats->options & PS_STAT_FITTED_STDEV)) { + // + // New algorithm for calculated mean and stdev. + // + + // XXX: Where is this documented? + psF32 dN = (psF32) (0.17 * N50); + if (dN < 1.0) { + dN = 1.0; + } else if (dN > 4.0) { + dN = 4.0; + } + psF32 newBinSize = sigma / dN; + + rc = p_psVectorMin(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax); + rc|= p_psVectorMax(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax); + if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpScalar); + psFree(tmpMaskVec); + psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); + return(1); + } + + numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / newBinSize); + psTrace(__func__, 6, "The new min/max values are (%f, %f).\n", tmpStatsMinMax->min, tmpStatsMinMax->max); + psTrace(__func__, 6, "The new bin size is %f.\n", newBinSize); + psTrace(__func__, 6, "The numBins is %d\n", numBins); + + psHistogram *newHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins); + newHistogram = psVectorHistogram(newHistogram, myVector, errors, tmpMaskVec, maskVal|1); + if (psTraceGetLevel(__func__) >= 8) { + PS_VECTOR_PRINT_F32(newHistogram->nums); + } + + // + // FITTED STATISTICS HERE + // + // Smooth the resulting histogram with a Gaussian with sigma_s = 1 + // bin. + // + psF32 sigma_s = newBinSize; + + psVector *newHistogramSmoothed = p_psVectorSmoothHistGaussian(newHistogram, sigma_s); + if (psTraceGetLevel(__func__) >= 8) { + PS_VECTOR_PRINT_F32(newHistogramSmoothed); + } + + // + // Find the bin with the peak value in the range 2 sigma of the + // robust histogram median. + // + + psS32 binMin = 0; + psS32 binMax = 0; + tmpScalar->data.F32 = stats->robustMedian - (2.0 * sigma); + if (tmpScalar->data.F32 <= newHistogram->bounds->data.F32[0]) { + binMin = 0; + } else { + binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); + } + + tmpScalar->data.F32 = stats->robustMedian + (2.0 + sigma); + if (tmpScalar->data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) { + binMax = newHistogram->bounds->n-1; + } else { + binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); + } + if ((binMin < 0) || (binMax < 0)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the +- 2.0 * sigma bins\n"); + psFree(tmpMaskVec); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpStatsMinMax); + psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); + return(1); + } + + psS32 binNum = binMin; + psF32 binMaxNums = newHistogramSmoothed->data.F32[binNum]; + for (psS32 i = binMin+1 ; i <= binMax ; i++) { + if (newHistogramSmoothed->data.F32[i] > binMaxNums) { + binNum = i; + binMaxNums = newHistogramSmoothed->data.F32[i]; + } + } + psTrace(__func__, 6, "The peak bin is %d, with %f data.n", binNum, binMaxNums); + + // + // Fit a Gaussian to the bins in the range 20 sigma of the robust + // histogram median. + // + tmpScalar->data.F32 = stats->robustMedian - (20.0 * sigma); + if (tmpScalar->data.F32 < tmpStatsMinMax->min) { + binMin = 0; + } else { + binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); + // XXX: check for errors here. + } + tmpScalar->data.F32 = stats->robustMedian + (20.0 * sigma); + if (tmpScalar->data.F32 > tmpStatsMinMax->max) { + binMax = newHistogramSmoothed->n - 1; + } else { + binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); + // XXX: check for errors here. + } + psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); + psVector *xTmp = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); + psArray *x = psArrayAlloc((1 + (binMax - binMin))); + psS32 j = 0; + + for (psS32 i = binMin ; i <= binMax ; i++) { + y->data.F32[j] = newHistogramSmoothed->data.F32[i]; + x->data[j] = (psPtr *) psVectorAlloc(1, PS_TYPE_F32); + ((psVector *) x->data[j])->data.F32[0] = PS_BIN_MIDPOINT(newHistogram, i); + xTmp->data.F32[j] = PS_BIN_MIDPOINT(newHistogram, i); + j++; + } + if (psTraceGetLevel(__func__) >= 8) { + // XXX: Print the x array somehow. + PS_VECTOR_PRINT_F32(y); + } + + // XXX: Use the min/max routines for this + psF32 minY = FLT_MAX; + psF32 maxY = -FLT_MAX; + for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) { + if (y->data.F32[i] > maxY) { + maxY = y->data.F32[i]; + } + if (y->data.F32[i] < minY) { + minY = y->data.F32[i]; + } + } + // + // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0]) + // XXX: Use the normalize routines for this. + // + for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) { + y->data.F32[i]= (y->data.F32[i] - minY) / (maxY - minY); + } + + // + psMinimization *min = psMinimizationAlloc(100, 0.01); + psVector *params = psVectorAlloc(2, PS_TYPE_F32); + // Initial guess for the mean ([0]) and standard dev. + params->data.F32[0] = stats->robustMedian; + params->data.F32[1] = sigma; + if (psTraceGetLevel(__func__) >= 8) { + PS_VECTOR_PRINT_F32(params); + PS_VECTOR_PRINT_F32(y); + } + psFree(xTmp); + rcBool = psMinimizeLMChi2(min, NULL, params, NULL, x, y, NULL, psMinimizeLMChi2Gauss1D); + + if (rcBool != true) { + psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpScalar); + psFree(newHistogram); + psFree(x); + psFree(y); + psFree(min); + psFree(params); + psFree(tmpMaskVec); + psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); + return(1); + } + if (psTraceGetLevel(__func__) >= 8) { + PS_VECTOR_PRINT_F32(params); + } + + // + // The fitted mean mean_r is derived directly from the fitted + // Gaussian mean. + // + stats->fittedMean = params->data.F32[0]; + psTrace(__func__, 6, "The fitted mean is %f.\n", params->data.F32[0]); + + // + // The fitted standard deviation, SIGMA_r is determined by + // subtracting the smoothing scale in quadrature: + // SIGMA_r^2 = SIGMA^2 - sigma_s^2 + // + stats->fittedStdev = sqrt(PS_SQR(params->data.F32[1]) - PS_SQR(sigma_s)); + psTrace(__func__, 6, "The fitted stdev is %f.\n", stats->fittedStdev); + + psFree(newHistogramSmoothed); psFree(newHistogram); psFree(x); @@ -2161,36 +2204,9 @@ psFree(min); psFree(params); - psFree(tmpMaskVec); - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); - } - if (psTraceGetLevel(__func__) >= 8) { - PS_VECTOR_PRINT_F32(params); - } - - // - // The robust mean mean_r is derived directly from the fitted - // Gaussian mean. - // - stats->robustMean = params->data.F32[0]; - psTrace(__func__, 6, "The robust mean is %f.\n", params->data.F32[0]); - - // - // The robust standard deviation, SIGMA_r is determined by - // subtracting the smoothing scale in quadrature: - // SIGMA_r^2 = SIGMA^2 - sigma_s^2 - // - stats->robustStdev = sqrt(PS_SQR(params->data.F32[1]) - PS_SQR(sigma_s)); - psTrace(__func__, 6, "The robust stdev is %f.\n", stats->robustStdev); - - psFree(newHistogramSmoothed); + } + psFree(tmpStatsMinMax); psFree(cumulativeRobustHistogram); psFree(tmpScalar); - psFree(newHistogram); - psFree(x); - psFree(y); - psFree(min); - psFree(params); psFree(tmpMaskVec); psFree(robustHistogram); @@ -2236,16 +2252,15 @@ newStruct->sampleUQ = NAN; newStruct->sampleLQ = NAN; - newStruct->robustMean = NAN; newStruct->robustMedian = NAN; - newStruct->robustMode = NAN; newStruct->robustStdev = NAN; newStruct->robustUQ = NAN; newStruct->robustLQ = NAN; newStruct->robustN50 = -1; // XXX: This is never used - newStruct->robustNfit = -1; + newStruct->fittedMean = NAN; + newStruct->fittedStdev = NAN; + newStruct->fittedNfit = -1; newStruct->clippedMean = NAN; newStruct->clippedStdev = NAN; newStruct->clippedNvalues = -1; // XXX: This is never used - // XXX: Where do these values come from? newStruct->clipSigma = 3.0; newStruct->clipIter = 3; @@ -2253,7 +2268,8 @@ newStruct->max = NAN; newStruct->binsize = NAN; // XXX: This is never used + newStruct->nSubsample = 100000; newStruct->options = options; - psTrace(__func__, 3, "---- %s(psStats) end ----\n", __func__); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); return (newStruct); } @@ -2318,5 +2334,5 @@ newHist->uniform = true; - psTrace(__func__, 3, "---- %s(psHistogram) end ----\n", __func__); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); return (newHist); } @@ -2363,5 +2379,5 @@ newHist->uniform = false; - psTrace(__func__, 3, "---- %s(psHistogram) end ----\n", __func__); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); return (newHist); } @@ -2594,5 +2610,5 @@ } - psTrace(__func__, 3, "---- %s(psHistogram) end ----\n", __func__); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); return (out); } @@ -2674,5 +2690,5 @@ } - psTrace(__func__, 4,"---- %s(psVector) end ----\n", __func__); + psTrace(__func__, 4,"---- %s() end ----\n", __func__); return (tmp); } @@ -2708,4 +2724,5 @@ PS_ASSERT_VECTOR_TYPE(errors, in->type.type, stats); } + // XXX: Assert that "in" is F64, F32, U16, or S8 psVector* inF32 = NULL; @@ -2729,47 +2746,51 @@ } + if ((stats->options & PS_STAT_USE_BINSIZE) && (stats->min >= stats->max)) { + PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(stats->binsize, 0.0, stats); + } + // ************************************************************************ if (stats->options & PS_STAT_SAMPLE_MEAN) { if (0 != p_psVectorSampleMean(inF32, errorsF32, mask, maskVal, stats)) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: psVectorStats(): p_psVectorSampleMean() returned an error.\n"); - } - } + psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMean() returned an error.\n"); + stats->sampleMean = NAN; + } + } + // ************************************************************************ if (stats->options & PS_STAT_SAMPLE_MEDIAN) { if (false == p_psVectorSampleMedian(inF32, mask, maskVal, stats)) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: psVectorStats(): p_psVectorSampleMedian() returned an error.\n"); - } - } + psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMedian() returned an error.\n"); + stats->sampleMedian = NAN; + } + } + // ************************************************************************ if (stats->options & PS_STAT_SAMPLE_STDEV) { if (0 != p_psVectorSampleMean(inF32, errorsF32, mask, maskVal, stats)) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: psVectorStats(): p_psVectorSampleMean() returned an error.\n"); - } - p_psVectorSampleStdev(inF32, errorsF32, mask, maskVal, stats); - } + psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMean() returned an error.\n"); + stats->sampleMean = NAN; + } else { + p_psVectorSampleStdev(inF32, errorsF32, mask, maskVal, stats); + } + } + // ************************************************************************ if (stats->options & PS_STAT_SAMPLE_QUARTILE) { if (false == p_psVectorSampleQuartiles(inF32, mask, maskVal, stats)) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: psVectorStats(): p_psVectorSampleQuartiles() returned an error.\n"); - } - } - // Since the various robust stats quantities share much computation, they - // are grouped together in a single private function: - // p_psVectorRobustStats() - if ((stats->options & PS_STAT_ROBUST_MEAN) || - (stats->options & PS_STAT_ROBUST_MEDIAN) || - (stats->options & PS_STAT_ROBUST_MODE) || + psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleQuartiles() returned an error.\n"); + stats->sampleLQ = NAN; + stats->sampleUQ = NAN; + } + } + + // ************************************************************************ + if ((stats->options & PS_STAT_ROBUST_MEDIAN) || (stats->options & PS_STAT_ROBUST_STDEV) || - (stats->options & PS_STAT_ROBUST_QUARTILE)) { + (stats->options & PS_STAT_ROBUST_QUARTILE) || + (stats->options & PS_STAT_FITTED_MEAN) || + (stats->options & PS_STAT_FITTED_STDEV)) { if (0 != p_psVectorRobustStats(inF32, errorsF32, mask, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, - PS_ERRORTEXT_psStats_STATS_FAILED); - // XXX: Set to NAN - // XXX: Is this the right thing to do? - // XXX: If so, do it for other funcs? + psError(PS_ERR_UNKNOWN, false, PS_ERRORTEXT_psStats_STATS_FAILED); psFree(stats); psTrace(__func__, 3,"---- %s(NULL) end ----\n", __func__); @@ -2778,31 +2799,26 @@ } - // XXX: Different conditions for return -1 and -2? + // ************************************************************************ if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) { psS32 rc = p_psVectorClippedStats(inF32, errorsF32, mask, maskVal, stats); - if (-1 == rc) { - psError(PS_ERR_UNKNOWN, false, - "Failed to calculate clipped statistics for input psVector.\n"); + if (rc < 0) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate clipped statistics for input psVector.\n"); stats->clippedMean = NAN; stats->clippedStdev = NAN; - } else if (-2 == rc) { - psLogMsg(__func__, PS_LOG_WARN, "Failed to calculate clipped statistics for input psVector."); - stats->clippedMean = NAN; - stats->clippedStdev = NAN; - } - } + } + } + // ************************************************************************ if (stats->options & PS_STAT_MAX) { if (0 != p_psVectorMax(inF32, mask, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, - "Failed to calculate vector maximum"); + psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector maximum"); stats->max = NAN; } } + // ************************************************************************ if (stats->options & PS_STAT_MIN) { if (0 != p_psVectorMin(inF32, mask, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, - "Failed to calculate vector minimum"); + psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector minimum"); stats->min = NAN; } @@ -2815,5 +2831,5 @@ psFree(errorsF32); } - psTrace(__func__, 3,"---- %s(psStats) end ----\n", __func__); + psTrace(__func__, 3,"---- %s() end ----\n", __func__); return (stats); }