Index: /trunk/psLib/src/math/psStats.c =================================================================== --- /trunk/psLib/src/math/psStats.c (revision 12436) +++ /trunk/psLib/src/math/psStats.c (revision 12437) @@ -13,6 +13,6 @@ * use ->min and ->max (PS_STAT_USE_RANGE) * - * @version $Revision: 1.201 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-02-13 03:08:17 $ + * @version $Revision: 1.202 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-03-14 02:43:47 $ * * Copyright 2006 IfA, University of Hawaii @@ -67,4 +67,46 @@ (0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1])) +// set the bin closest to the corresponding value. if USE_END is +/- 1, +// out-of-range saturates on lower/upper bin REGARDLESS of actual value +#define PS_BIN_FOR_VALUE(RESULT, VECTOR, VALUE, USE_END) \ + tmpScalar.data.F32 = (VALUE); \ + RESULT = psVectorBinaryDisect (&result, VECTOR, &tmpScalar); \ + switch (result) { \ + case PS_BINARY_DISECT_PASS: \ + break; \ + case PS_BINARY_DISECT_OUTSIDE_RANGE: \ + psTrace ("psLib.math", 6, "selected bin outside range"); \ + if (USE_END == -1) { RESULT = 0; } \ + if (USE_END == +1) { RESULT = VECTOR->n - 1; } \ + break; \ + case PS_BINARY_DISECT_INVALID_INPUT: \ + case PS_BINARY_DISECT_INVALID_TYPE: \ + psAbort ("programming error"); \ + break; \ + } + +# define PS_BIN_INTERPOLATE(RESULT, VECTOR, BOUNDS, BIN, VALUE) { \ + float dX, dY, Xo, Yo, Xt; \ + if (BIN == BOUNDS->n - 1) { \ + dX = 0.5*(BOUNDS->data.F32[BIN+1] - BOUNDS->data.F32[BIN-1]); \ + dY = VECTOR->data.F32[BIN] - VECTOR->data.F32[BIN-1]; \ + Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \ + Yo = VECTOR->data.F32[BIN]; \ + } else { \ + dX = 0.5*(BOUNDS->data.F32[BIN+2] - BOUNDS->data.F32[BIN]); \ + dY = VECTOR->data.F32[BIN+1] - VECTOR->data.F32[BIN]; \ + Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \ + Yo = VECTOR->data.F32[BIN]; \ + } \ + if (dY != 0) { \ + Xt = (VALUE - Yo)*dX/dY + Xo; \ + } else { \ + Xt = Xo; \ + } \ + Xt = PS_MIN (BOUNDS->data.F32[BIN+1], PS_MAX(BOUNDS->data.F32[BIN], Xt)); \ + psTrace("psLib.math", 6, "(Xo, Yo, dX, dY, Xt, Yt) is (%.2f %.2f %.2f %.2f %.2f %.2f)\n", \ + Xo, Yo, dX, dY, Xt, VALUE); \ + RESULT = Xt; } + /*****************************************************************************/ /* TYPE DEFINITIONS */ @@ -120,9 +162,9 @@ To optmize this, use a macro and ifdef in or out the three states (errors, mask, range) *****************************************************************************/ -static bool vectorSampleMean(const psVector* myVector, - const psVector* errors, - const psVector* maskVector, - psMaskType maskVal, - psStats* stats) + static bool vectorSampleMean(const psVector* myVector, + const psVector* errors, + const psVector* maskVector, + psMaskType maskVal, + psStats* stats) { psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); @@ -188,9 +230,9 @@ (mask: 1, range: 0): 0.200 sec 0.244 sec *****************************************************************************/ -static long vectorMinMax(const psVector* myVector, - const psVector* maskVector, - psMaskType maskVal, - psStats* stats - ) + static long vectorMinMax(const psVector* myVector, + const psVector* maskVector, + psMaskType maskVal, + psStats* stats + ) { psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); @@ -331,9 +373,9 @@ *****************************************************************************/ -static bool vectorSampleStdev(const psVector* myVector, - const psVector* errors, - const psVector* maskVector, - psMaskType maskVal, - psStats* stats) + static bool vectorSampleStdev(const psVector* myVector, + const psVector* errors, + const psVector* maskVector, + psMaskType maskVal, + psStats* stats) { psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); @@ -422,5 +464,5 @@ psMaskType maskValInput, psStats* stats - ) + ) { psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); @@ -491,5 +533,5 @@ for (long j = 0; j < myVector->n; j++) { if (!tmpMask->data.U8[j] && - fabsf(myVector->data.F32[j] - clippedMean) > stats->clipSigma * errors->data.F32[j]) { + fabsf(myVector->data.F32[j] - clippedMean) > stats->clipSigma * errors->data.F32[j]) { tmpMask->data.U8[j] = 0xff; psTrace("psLib.math", 10, "Clipped %ld: %f +/- %f\n", j, @@ -502,5 +544,5 @@ for (long j = 0; j < myVector->n; j++) { if (!tmpMask->data.U8[j] && - fabsf(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) { + fabsf(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) { tmpMask->data.U8[j] = 0xff; psTrace("psLib.math", 10, "Clipped %ld: %f\n", j, myVector->data.F32[j]); @@ -584,5 +626,6 @@ } - psScalar tmpScalar; // Temporary scalar variable, for p_psVectorBinDisect + psVectorBinaryDisectResult result; + psScalar tmpScalar; // Temporary scalar variable, for psVectorBinaryDisect tmpScalar.type.type = PS_TYPE_F32; @@ -609,4 +652,9 @@ long totalDataPoints = 0; // Total number of (unmasked) data points + float binSize = 0.0; // Size of bins for histogram + long binLo, binHi; + long binL2, binH2; + long binMedian; + // Iterate to get the best bin size for (int iterate = 1; iterate > 0; iterate++) { @@ -620,8 +668,5 @@ // Data range calculation failed psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); - psFree(statsMinMax); - psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); - psFree(mask); - return false; + goto escape; } psTrace("psLib.math", 6, "Data min/max is (%.2f, %.2f)\n", min, max); @@ -629,5 +674,4 @@ // If all data points have the same value, then we set the appropriate members of stats and return. if (fabs(max - min) <= FLT_EPSILON) { - psTrace("psLib.math", 5, "All data points have the same value: %f.\n", min); stats->robustMedian = min; stats->robustStdev = 0.0; @@ -635,16 +679,14 @@ stats->robustLQ = min; stats->robustN50 = numValid; - psFree(statsMinMax); - stats->results |= PS_STAT_ROBUST_MEDIAN; stats->results |= PS_STAT_ROBUST_STDEV; stats->results |= PS_STAT_ROBUST_QUARTILE; - + psTrace("psLib.math", 5, "All data points have the same value: %f.\n", min); psTrace("psLib.math", 4, "---- %s(0) end ----\n", __func__); psFree(mask); + psFree(statsMinMax); return true; } - float binSize = 0.0; // Size of bins for histogram if ((iterate == 1) && (stats->options & PS_STAT_USE_BINSIZE)) { // Set initial bin size to the specified value. @@ -668,7 +710,5 @@ if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for robust statistics.\n"); - psFree(histogram); - psFree(statsMinMax); - return false; + goto escape; } if (psTraceGetLevel("psLib.math") >= 8) { @@ -691,23 +731,6 @@ totalDataPoints = cumulative->nums->data.F32[numBins - 1]; psTrace("psLib.math", 6, "Total data points is %ld\n", totalDataPoints); - long binMedian; - if (totalDataPoints/2.0 < cumulative->nums->data.F32[0]) { - // Special case: the median is in the first bin. Perhaps we should recode this. - // XXX: Must try a special case where the median in in the last bin. - binMedian = 0; - } else { - tmpScalar.data.F32 = totalDataPoints/2.0; - binMedian = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar); - if (binMedian < 0) { - psError(PS_ERR_UNKNOWN, false, - "Failed to calculate the 50 precent data point (%ld).\n", binMedian); - psFree(statsMinMax); - psFree(histogram); - psFree(cumulative); - psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); - psFree(mask); - return false; - } - } + + PS_BIN_FOR_VALUE(binMedian, cumulative->nums, totalDataPoints/2.0, 0); psTrace("psLib.math", 6, "The median bin is %ld (%.2f to %.2f)\n", binMedian, cumulative->bounds->data.F32[binMedian], cumulative->bounds->data.F32[binMedian+1]); @@ -715,37 +738,19 @@ // ADD step 3: Interpolate to the exact 50% position: this is the robust histogram median. stats->robustMedian = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, - binMedian, totalDataPoints/2.0); + binMedian, totalDataPoints/2.0); if (isnan(stats->robustMedian)) { - psError(PS_ERR_UNKNOWN, false, - "Failed to fit a quadratic and calculate the 50-percent position.\n"); - psFree(statsMinMax); - psFree(histogram); - psFree(cumulative); - psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); - psFree(mask); - return false; + psError(PS_ERR_UNKNOWN, false, "Failed to fit a quadratic and calculate the 50-percent position.\n"); + goto escape; } psTrace("psLib.math", 6, "Current robust median is %f\n", stats->robustMedian); - // ADD step 4: Find the bins which contains the 15.8655% and 84.1345% data points. - long binLo = 0; // Bin containing the 15.8655% mark - tmpScalar.data.F32 = totalDataPoints * 0.158655f; - if (tmpScalar.data.F32 > cumulative->nums->data.F32[0]) { - // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined. - binLo = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar); - if (binLo > cumulative->nums->n-1) { - binLo = cumulative->nums->n-1; - } - } - - long binHi = 0; // Bin containing the 84.1345% mark - tmpScalar.data.F32 = totalDataPoints * 0.841345f; - if (tmpScalar.data.F32 > cumulative->nums->data.F32[0]) { - // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined. - binHi = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar); - if (binHi > cumulative->nums->n-1) { - binHi = cumulative->nums->n-1; - } - } + // ADD step 4: Find the bins which contains the 15.8655% (-1 sigma) and 84.1345% (+1 sigma) data points + // also find the 30.8538% (-0.5 sigma) and 69.1462% (+0.5 sigma) points + PS_BIN_FOR_VALUE(binLo, cumulative->nums, totalDataPoints * 0.158655f, -1); + PS_BIN_FOR_VALUE(binHi, cumulative->nums, totalDataPoints * 0.841345f, +1); + + PS_BIN_FOR_VALUE(binL2, cumulative->nums, totalDataPoints * 0.308538f, -1); + PS_BIN_FOR_VALUE(binH2, cumulative->nums, totalDataPoints * 0.691462f, +1); + psTrace("psLib.math", 6, "The 15.8655%% and 84.1345%% data point bins are (%ld, %ld).\n", binLo, binHi); @@ -755,10 +760,5 @@ if ((binLo < 0) || (binHi < 0)) { psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 15.8655%% and 84.1345%% data points.\n"); - psFree(statsMinMax); - psFree(histogram); - psFree(cumulative); - psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); - psFree(mask); - return false; + goto escape; } @@ -769,38 +769,12 @@ binHi, cumulative->nums->data.F32[binHi], cumulative->nums->data.F32[binHi+1]); - float deltaBounds, deltaNums, prevPixels; - float percentNums, base, binLoF32; - - // find the -1 sigma points with linear interpolation - deltaBounds = cumulative->bounds->data.F32[binLo+1] - cumulative->bounds->data.F32[binLo]; - if (binLo == 0) { - deltaNums = cumulative->nums->data.F32[0]; - prevPixels = 0; - } else { - deltaNums = cumulative->nums->data.F32[binLo] - cumulative->nums->data.F32[binLo-1]; - prevPixels = cumulative->nums->data.F32[binLo-1]; - } - percentNums = (totalDataPoints * 0.158655f) - prevPixels; - base = cumulative->bounds->data.F32[binLo]; - binLoF32 = base + (deltaBounds / deltaNums) * percentNums; // Value for the 15.8655% mark - psTrace("psLib.math", 6, - "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n", - base, deltaBounds, deltaNums, prevPixels, percentNums); - - // find the +1 sigma points with linear interpolation - deltaBounds = cumulative->bounds->data.F32[binHi+1] - cumulative->bounds->data.F32[binHi]; - if (binHi == 0) { - deltaNums = cumulative->nums->data.F32[0]; - prevPixels = 0; - } else { - deltaNums = cumulative->nums->data.F32[binHi] - cumulative->nums->data.F32[binHi-1]; - prevPixels = cumulative->nums->data.F32[binHi-1]; - } - percentNums = (totalDataPoints * 0.841345f) - prevPixels; - base = cumulative->bounds->data.F32[binHi]; - float binHiF32 = base + (deltaBounds / deltaNums) * percentNums; // Value for the 84.1345% mark - psTrace("psLib.math", 6, - "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n", - base, deltaBounds, deltaNums, prevPixels, percentNums); + // find the +0.5 and -0.5 sigma points with linear interpolation. binLo and binHi are the bins + // containing the value of interest. constrain the answer to be within the current bin. + // (extrapolation should not be needed and will result in errors) + float binLoF32, binHiF32, binL2F32, binH2F32; + PS_BIN_INTERPOLATE (binLoF32, cumulative->nums, cumulative->bounds, binL2, totalDataPoints * 0.158655f); + PS_BIN_INTERPOLATE (binHiF32, cumulative->nums, cumulative->bounds, binH2, totalDataPoints * 0.841345f); + PS_BIN_INTERPOLATE (binL2F32, cumulative->nums, cumulative->bounds, binL2, totalDataPoints * 0.308538f); + PS_BIN_INTERPOLATE (binH2F32, cumulative->nums, cumulative->bounds, binH2, totalDataPoints * 0.691462f); // report +/- 1 sigma points @@ -809,6 +783,8 @@ binLoF32, binHiF32); - // ADD step 5: Determine SIGMA as 1/2 of the distance between these positions. - sigma = (binHiF32 - binLoF32) / 2.0; + // ADD step 5: Determine SIGMA as (1/2 of) the distance between these positions. + // sigma = (binHiF32 - binLoF32) / 2.0; + sigma = (binH2F32 - binL2F32); + psTrace("psLib.math", 6, "The current sigma is %f.\n", sigma); stats->robustStdev = sigma; @@ -835,5 +811,8 @@ // Free the histograms; they will be recreated on the next iteration, with new bounds psFree(histogram); + histogram = NULL; + psFree(cumulative); + cumulative = NULL; } else { // We've got the bin size correct now @@ -842,32 +821,14 @@ } } - psFree(histogram); + + // XXX test lines while studying algorithm errors + // fprintf (stderr, "robust stats test %7.1f +/- %7.1f : %4ld %4ld %4ld %4ld %4ld : %f %f %f\n", + // stats->robustMedian, stats->robustStdev, + // binLo, binL2, binMedian, binH2, binHi, binSize, max, min); // ADD step 7: Find the bins which contains the 25% and 75% data points. - long binLo25 = 0; // Bin containing the 25% value - long binHi25 = 0; // Bin containing the 75% value - - tmpScalar.data.F32 = totalDataPoints * 0.25f; - if (tmpScalar.data.F32 <= cumulative->nums->data.F32[0]) { - // Special case where its in first bin. The last bin should be handled automatically. - binLo25 = 0; - } else { - binLo25 = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar); - } - tmpScalar.data.F32 = totalDataPoints * 0.75f; - if (tmpScalar.data.F32 <= cumulative->nums->data.F32[0]) { - // Special case where its in first bin. The last bin should be handled automatically. - binHi25 = 0; - } else { - binHi25 = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar); - } - if ((binLo25 < 0) || (binHi25 < 0)) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 25 and 75 percent data points\n"); - psFree(statsMinMax); - psFree(cumulative); - psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); - psFree(mask); - return false; - } + long binLo25, binHi25; + PS_BIN_FOR_VALUE (binLo25, cumulative->nums, totalDataPoints * 0.25f, 0); + PS_BIN_FOR_VALUE (binHi25, cumulative->nums, totalDataPoints * 0.75f, 0); psTrace("psLib.math", 6, "The 25-percent and 75-precent data point bins are (%ld, %ld).\n", binLo25, binHi25); @@ -875,19 +836,13 @@ // positions. psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, - binLo25, totalDataPoints * 0.25f); + binLo25, totalDataPoints * 0.25f); psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, - binHi25, totalDataPoints * 0.75f); - psFree(cumulative); - + binHi25, totalDataPoints * 0.75f); if (isnan(binLo25F32) || isnan(binHi25F32)) { psError(PS_ERR_UNKNOWN, false, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnX() returned a NAN.\n"); - psFree(statsMinMax); - psTrace("psLib.math", 4, "---- %s(1) end ----\n", __func__); - psFree(mask); - return false; - } - - stats->results |= PS_STAT_ROBUST_QUARTILE; + goto escape; + } + stats->robustLQ = binLo25F32; stats->robustUQ = binHi25F32; @@ -897,5 +852,5 @@ for (long i = 0 ; i < myVector->n ; i++) { if (!mask->data.U8[i] && - (binLo25F32 <= myVector->data.F32[i]) && (binHi25F32 >= myVector->data.F32[i])) { + (binLo25F32 <= myVector->data.F32[i]) && (binHi25F32 >= myVector->data.F32[i])) { N50++; } @@ -905,4 +860,6 @@ // Clean up + psFree(histogram); + psFree(cumulative); psFree(statsMinMax); psFree(mask); @@ -910,7 +867,27 @@ stats->results |= PS_STAT_ROBUST_MEDIAN; stats->results |= PS_STAT_ROBUST_STDEV; + stats->results |= PS_STAT_ROBUST_QUARTILE; psTrace("psLib.math", 4, "---- %s(0) end ----\n", __func__); return true; + +escape: + stats->robustMedian = NAN; + stats->robustStdev = NAN; + stats->robustUQ = NAN; + stats->robustLQ = NAN; + stats->robustN50 = 0; + stats->results |= PS_STAT_ROBUST_MEDIAN; + stats->results |= PS_STAT_ROBUST_STDEV; + stats->results |= PS_STAT_ROBUST_QUARTILE; + + psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); + + psFree(histogram); + psFree(cumulative); + psFree(statsMinMax); + psFree(mask); + + return false; } @@ -943,4 +920,6 @@ float guessMean = stats->robustMedian; // pass the guess mean + psVectorBinaryDisectResult result; + psTrace("psLib.math", 6, "The guess mean is %f.\n", guessMean); psTrace("psLib.math", 6, "The guess stdev is %f.\n", guessStdev); @@ -950,5 +929,5 @@ psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max - psScalar tmpScalar; // Temporary scalar variable, for p_psVectorBinDisect + psScalar tmpScalar; // Temporary scalar variable, for psVectorBinaryDisect tmpScalar.type.type = PS_TYPE_F32; @@ -995,7 +974,4 @@ } - long binMin = 0; // Low bin to check - long binMax = 0; // High bin to check - // Fit a Gaussian to the bins in the requested range about the guess mean // min,max overrides clipSigma @@ -1016,23 +992,8 @@ } - tmpScalar.data.F32 = guessMean - minFitSigma*guessStdev; - if (tmpScalar.data.F32 < min) { - binMin = 0; - } else { - binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar); - } - tmpScalar.data.F32 = guessMean + maxFitSigma*guessStdev; - if (tmpScalar.data.F32 > max) { - binMax = histogram->bounds->n - 1; - } else { - binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar); - } - if (binMin < 0 || binMax < 0) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the -%f sigma : +%f sigma bins\n", minFitSigma, maxFitSigma); - psFree(statsMinMax); - psFree(histogram); - psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); - return false; - } + // select the min and max bins, saturating on the lower and upper end-points + long binMin, binMax; + PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, -1); + PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, +1); // Generate the variables that will be used in the Gaussian fitting @@ -1142,4 +1103,5 @@ float guessMean = stats->robustMedian; // pass the guess mean + psVectorBinaryDisectResult result; psTrace("psLib.math", 6, "The ** starting ** guess mean is %f.\n", guessMean); psTrace("psLib.math", 6, "The ** starting ** guess stdev is %f.\n", guessStdev); @@ -1149,5 +1111,5 @@ psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max - psScalar tmpScalar; // Temporary scalar variable, for p_psVectorBinDisect + psScalar tmpScalar; // Temporary scalar variable, for psVectorBinaryDisect tmpScalar.type.type = PS_TYPE_F32; @@ -1194,7 +1156,4 @@ } - long binMin = 0; // Low bin to check - long binMax = 0; // High bin to check - // Fit a Gaussian to the bins in the requested range about the guess mean // min,max overrides clipSigma @@ -1215,23 +1174,8 @@ } - tmpScalar.data.F32 = guessMean - minFitSigma*guessStdev; - if (tmpScalar.data.F32 < min) { - binMin = 0; - } else { - binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar); - } - tmpScalar.data.F32 = guessMean + maxFitSigma*guessStdev; - if (tmpScalar.data.F32 > max) { - binMax = histogram->bounds->n - 1; - } else { - binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar); - } - if (binMin < 0 || binMax < 0) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the -%f sigma : +%f sigma bins\n", minFitSigma, maxFitSigma); - psFree(statsMinMax); - psFree(histogram); - psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); - return false; - } + // select the min and max bins, saturating on the lower and upper end-points + long binMin, binMax; + PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, -1); + PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, +1); // Generate the variables that will be used in the Gaussian fitting @@ -1355,4 +1299,5 @@ float guessMean = stats->robustMedian; // pass the guess mean + psVectorBinaryDisectResult result; psTrace("psLib.math", 6, "The ** starting ** guess mean is %f.\n", guessMean); psTrace("psLib.math", 6, "The ** starting ** guess stdev is %f.\n", guessStdev); @@ -1362,5 +1307,5 @@ psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max - psScalar tmpScalar; // Temporary scalar variable, for p_psVectorBinDisect + psScalar tmpScalar; // Temporary scalar variable, for psVectorBinaryDisect tmpScalar.type.type = PS_TYPE_F32; @@ -1409,7 +1354,4 @@ // now fit a Gaussian to the upper and lower halves about the peak independently - long binMin, binMax; - float upperMean, upperStdev; - // set the full-range upper and lower limits psF32 maxFitSigma = 2.0; @@ -1429,23 +1371,8 @@ } - tmpScalar.data.F32 = guessMean - minFitSigma*guessStdev; - if (tmpScalar.data.F32 < min) { - binMin = 0; - } else { - binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar); - } - tmpScalar.data.F32 = guessMean + maxFitSigma*guessStdev; - if (tmpScalar.data.F32 > max) { - binMax = histogram->bounds->n; - } else { - binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar) + 1; - } - if (binMin < 0 || binMax < 0) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the -%f sigma : +%f sigma bins\n", minFitSigma, maxFitSigma); - psFree(statsMinMax); - psFree(histogram); - psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); - return false; - } + // select the min and max bins, saturating on the lower and upper end-points + long binMin, binMax; + PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, -1); + PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, +1); // search for mode (peak of histogram within range mean-2sigma - mean+2sigma @@ -1463,5 +1390,5 @@ psTrace("psLib.math", 6, "The clipped numBins is %ld\n", binMax - binMin); psTrace("psLib.math", 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin); - psTrace("psLib.math", 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax), binMax); + psTrace("psLib.math", 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax - 1), binMax - 1); psTrace("psLib.math", 6, "The clipped peak is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binPeak), binPeak); @@ -1581,6 +1508,6 @@ // calculate upper mean & stdev from parabolic fit -- use this as the result - upperStdev = sqrt(-0.5/poly->coeff[2]); - upperMean = poly->coeff[1]*PS_SQR(upperStdev); + float upperStdev = sqrt(-0.5/poly->coeff[2]); + float upperMean = poly->coeff[1]*PS_SQR(upperStdev); psTrace("psLib.math", 6, "Parabolic Upper fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); psTrace("psLib.math", 6, "The upper mean is %f.\n", upperMean); @@ -1680,8 +1607,8 @@ // took way too long. // - #if 0 +#if 0 gaussWidth = 10; - #endif +#endif // // We determine the bin numbers (jMin:jMax) corresponding to a @@ -1926,5 +1853,5 @@ psStatsOptions psStatsOptionFromString(const char *string) { - #define READ_STAT(NAME, SYMBOL) \ +#define READ_STAT(NAME, SYMBOL) \ if (strcasecmp(string, NAME) == 0) { \ return SYMBOL; \ @@ -1964,5 +1891,5 @@ psString string = NULL; // String to return - #define WRITE_STAT(NAME, SYMBOL) \ +#define WRITE_STAT(NAME, SYMBOL) \ if (option & SYMBOL) { \ psStringAppend(&string, "%s ", NAME); \ @@ -2023,23 +1950,23 @@ { switch (option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE)) { - case PS_STAT_SAMPLE_MEAN: - case PS_STAT_SAMPLE_MEDIAN: - case PS_STAT_SAMPLE_STDEV: - case PS_STAT_SAMPLE_QUARTILE: - case PS_STAT_ROBUST_MEDIAN: - case PS_STAT_ROBUST_STDEV: - case PS_STAT_ROBUST_QUARTILE: - case PS_STAT_FITTED_MEAN: - case PS_STAT_FITTED_STDEV: - case PS_STAT_FITTED_MEAN_V2: - case PS_STAT_FITTED_STDEV_V2: - case PS_STAT_FITTED_MEAN_V3: - case PS_STAT_FITTED_STDEV_V3: - case PS_STAT_CLIPPED_MEAN: - case PS_STAT_CLIPPED_STDEV: - case PS_STAT_MAX: - case PS_STAT_MIN: + case PS_STAT_SAMPLE_MEAN: + case PS_STAT_SAMPLE_MEDIAN: + case PS_STAT_SAMPLE_STDEV: + case PS_STAT_SAMPLE_QUARTILE: + case PS_STAT_ROBUST_MEDIAN: + case PS_STAT_ROBUST_STDEV: + case PS_STAT_ROBUST_QUARTILE: + case PS_STAT_FITTED_MEAN: + case PS_STAT_FITTED_STDEV: + case PS_STAT_FITTED_MEAN_V2: + case PS_STAT_FITTED_STDEV_V2: + case PS_STAT_FITTED_MEAN_V3: + case PS_STAT_FITTED_STDEV_V3: + case PS_STAT_CLIPPED_MEAN: + case PS_STAT_CLIPPED_STDEV: + case PS_STAT_MAX: + case PS_STAT_MIN: return option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE); - default: + default: return 0; } @@ -2052,40 +1979,40 @@ // We could call psStatsSingle to check, but it would be a waste since we effectively do it anyway switch (option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE)) { - case PS_STAT_SAMPLE_MEAN: + case PS_STAT_SAMPLE_MEAN: return stats->sampleMean; - case PS_STAT_SAMPLE_MEDIAN: + case PS_STAT_SAMPLE_MEDIAN: return stats->sampleMedian; - case PS_STAT_SAMPLE_STDEV: + case PS_STAT_SAMPLE_STDEV: return stats->sampleStdev; - case PS_STAT_ROBUST_MEDIAN: + case PS_STAT_ROBUST_MEDIAN: return stats->robustMedian; - case PS_STAT_ROBUST_STDEV: + case PS_STAT_ROBUST_STDEV: return stats->robustStdev; - case PS_STAT_FITTED_MEAN: + case PS_STAT_FITTED_MEAN: return stats->fittedMean; - case PS_STAT_FITTED_STDEV: + case PS_STAT_FITTED_STDEV: return stats->fittedStdev; - case PS_STAT_FITTED_MEAN_V2: + case PS_STAT_FITTED_MEAN_V2: return stats->fittedMean; - case PS_STAT_FITTED_STDEV_V2: + case PS_STAT_FITTED_STDEV_V2: return stats->fittedStdev; - case PS_STAT_FITTED_MEAN_V3: + case PS_STAT_FITTED_MEAN_V3: return stats->fittedMean; - case PS_STAT_FITTED_STDEV_V3: + case PS_STAT_FITTED_STDEV_V3: return stats->fittedStdev; - case PS_STAT_CLIPPED_MEAN: + case PS_STAT_CLIPPED_MEAN: return stats->clippedMean; - case PS_STAT_CLIPPED_STDEV: + case PS_STAT_CLIPPED_STDEV: return stats->clippedStdev; - case PS_STAT_MAX: + case PS_STAT_MAX: return stats->max; - case PS_STAT_MIN: + case PS_STAT_MIN: return stats->min; - case PS_STAT_SAMPLE_QUARTILE: - case PS_STAT_ROBUST_QUARTILE: + case PS_STAT_SAMPLE_QUARTILE: + case PS_STAT_ROBUST_QUARTILE: psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cannot return a single quartile value; " "get them yourself.\n"); return NAN; - default: + default: return NAN; } @@ -2102,5 +2029,5 @@ psF32 xLo, psF32 xHi - ) + ) { psF64 tmp = sqrt((y - c)/a + (b*b)/(4.0 * a * a)); @@ -2130,8 +2057,8 @@ *****************************************************************************/ static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, - psVector *yVec, - psS32 binNum, - psF32 yVal - ) + psVector *yVec, + psS32 binNum, + psF32 yVal + ) { psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); @@ -2171,5 +2098,5 @@ // if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) || - ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) { + ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, _("Specified yVal, %g, is not within y-range, %g to %g."), @@ -2181,5 +2108,5 @@ // if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) || - ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) { + ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) { psError(PS_ERR_UNKNOWN, true, "This routine must be called with monotonically increasing or decreasing data points.\n"); @@ -2252,5 +2179,5 @@ const psVector *params, const psVector *coords - ) + ) { psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);