Changeset 12495 for trunk/psLib/src/math/psStats.c
- Timestamp:
- Mar 19, 2007, 12:48:11 PM (19 years ago)
- File:
-
- 1 edited
-
trunk/psLib/src/math/psStats.c (modified) (25 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/math/psStats.c
r12458 r12495 13 13 * use ->min and ->max (PS_STAT_USE_RANGE) 14 14 * 15 * @version $Revision: 1.20 3$ $Name: not supported by cvs2svn $16 * @date $Date: 2007-03-1 6 00:33:55$15 * @version $Revision: 1.204 $ $Name: not supported by cvs2svn $ 16 * @date $Date: 2007-03-19 22:48:11 $ 17 17 * 18 18 * Copyright 2006 IfA, University of Hawaii … … 67 67 (0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1])) 68 68 69 // set the bin closest to the corresponding value. if USE_END is +/- 1, 69 // set the bin closest to the corresponding value. if USE_END is +/- 1, 70 70 // out-of-range saturates on lower/upper bin REGARDLESS of actual value 71 71 #define PS_BIN_FOR_VALUE(RESULT, VECTOR, VALUE, USE_END) { \ 72 psVectorBinaryDisectResult result; \73 psScalar tmpScalar; \74 tmpScalar.type.type = PS_TYPE_F32; \75 tmpScalar.data.F32 = (VALUE); \76 RESULT = psVectorBinaryDisect (&result, VECTOR, &tmpScalar); \77 switch (result) { \78 case PS_BINARY_DISECT_PASS: \79 break; \80 case PS_BINARY_DISECT_OUTSIDE_RANGE: \81 psTrace ("psLib.math", 6, "selected bin outside range"); \82 if (USE_END == -1) { RESULT = 0; } \83 if (USE_END == +1) { RESULT = VECTOR->n - 1; } \84 break; \85 case PS_BINARY_DISECT_INVALID_INPUT: \86 case PS_BINARY_DISECT_INVALID_TYPE: \87 psAbort ("programming error"); \88 break; \72 psVectorBinaryDisectResult result; \ 73 psScalar tmpScalar; \ 74 tmpScalar.type.type = PS_TYPE_F32; \ 75 tmpScalar.data.F32 = (VALUE); \ 76 RESULT = psVectorBinaryDisect (&result, VECTOR, &tmpScalar); \ 77 switch (result) { \ 78 case PS_BINARY_DISECT_PASS: \ 79 break; \ 80 case PS_BINARY_DISECT_OUTSIDE_RANGE: \ 81 psTrace ("psLib.math", 6, "selected bin outside range"); \ 82 if (USE_END == -1) { RESULT = 0; } \ 83 if (USE_END == +1) { RESULT = VECTOR->n - 1; } \ 84 break; \ 85 case PS_BINARY_DISECT_INVALID_INPUT: \ 86 case PS_BINARY_DISECT_INVALID_TYPE: \ 87 psAbort ("programming error"); \ 88 break; \ 89 89 } } 90 90 … … 92 92 float dX, dY, Xo, Yo, Xt; \ 93 93 if (BIN == BOUNDS->n - 1) { \ 94 dX = 0.5*(BOUNDS->data.F32[BIN+1] - BOUNDS->data.F32[BIN-1]); \95 dY = VECTOR->data.F32[BIN] - VECTOR->data.F32[BIN-1]; \96 Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \97 Yo = VECTOR->data.F32[BIN]; \94 dX = 0.5*(BOUNDS->data.F32[BIN+1] - BOUNDS->data.F32[BIN-1]); \ 95 dY = VECTOR->data.F32[BIN] - VECTOR->data.F32[BIN-1]; \ 96 Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \ 97 Yo = VECTOR->data.F32[BIN]; \ 98 98 } else { \ 99 dX = 0.5*(BOUNDS->data.F32[BIN+2] - BOUNDS->data.F32[BIN]); \100 dY = VECTOR->data.F32[BIN+1] - VECTOR->data.F32[BIN]; \101 Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \102 Yo = VECTOR->data.F32[BIN]; \103 } \104 if (dY != 0) { \105 Xt = (VALUE - Yo)*dX/dY + Xo; \106 } else { \107 Xt = Xo; \108 } \109 Xt = PS_MIN (BOUNDS->data.F32[BIN+1], PS_MAX(BOUNDS->data.F32[BIN], Xt)); \99 dX = 0.5*(BOUNDS->data.F32[BIN+2] - BOUNDS->data.F32[BIN]); \ 100 dY = VECTOR->data.F32[BIN+1] - VECTOR->data.F32[BIN]; \ 101 Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \ 102 Yo = VECTOR->data.F32[BIN]; \ 103 } \ 104 if (dY != 0) { \ 105 Xt = (VALUE - Yo)*dX/dY + Xo; \ 106 } else { \ 107 Xt = Xo; \ 108 } \ 109 Xt = PS_MIN (BOUNDS->data.F32[BIN+1], PS_MAX(BOUNDS->data.F32[BIN], Xt)); \ 110 110 psTrace("psLib.math", 6, "(Xo, Yo, dX, dY, Xt, Yt) is (%.2f %.2f %.2f %.2f %.2f %.2f)\n", \ 111 111 Xo, Yo, dX, dY, Xt, VALUE); \ … … 166 166 *****************************************************************************/ 167 167 static bool vectorSampleMean(const psVector* myVector, 168 const psVector* errors,169 const psVector* maskVector,170 psMaskType maskVal,171 psStats* stats)168 const psVector* errors, 169 const psVector* maskVector, 170 psMaskType maskVal, 171 psStats* stats) 172 172 { 173 173 psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); … … 234 234 *****************************************************************************/ 235 235 static long vectorMinMax(const psVector* myVector, 236 const psVector* maskVector,237 psMaskType maskVal,238 psStats* stats239 )236 const psVector* maskVector, 237 psMaskType maskVal, 238 psStats* stats 239 ) 240 240 { 241 241 psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); … … 377 377 *****************************************************************************/ 378 378 static bool vectorSampleStdev(const psVector* myVector, 379 const psVector* errors,380 const psVector* maskVector,381 psMaskType maskVal,382 psStats* stats)379 const psVector* errors, 380 const psVector* maskVector, 381 psMaskType maskVal, 382 psStats* stats) 383 383 { 384 384 psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); … … 536 536 for (long j = 0; j < myVector->n; j++) { 537 537 if (!tmpMask->data.U8[j] && 538 fabsf(myVector->data.F32[j] - clippedMean) > stats->clipSigma * errors->data.F32[j]) {538 fabsf(myVector->data.F32[j] - clippedMean) > stats->clipSigma * errors->data.F32[j]) { 539 539 tmpMask->data.U8[j] = 0xff; 540 540 psTrace("psLib.math", 10, "Clipped %ld: %f +/- %f\n", j, … … 547 547 for (long j = 0; j < myVector->n; j++) { 548 548 if (!tmpMask->data.U8[j] && 549 fabsf(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) {549 fabsf(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) { 550 550 tmpMask->data.U8[j] = 0xff; 551 551 psTrace("psLib.math", 10, "Clipped %ld: %f\n", j, myVector->data.F32[j]); … … 667 667 // Data range calculation failed 668 668 psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); 669 goto escape;669 goto escape; 670 670 } 671 671 psTrace("psLib.math", 6, "Data min/max is (%.2f, %.2f)\n", min, max); … … 731 731 psTrace("psLib.math", 6, "Total data points is %ld\n", totalDataPoints); 732 732 733 PS_BIN_FOR_VALUE(binMedian, cumulative->nums, totalDataPoints/2.0, 0);733 PS_BIN_FOR_VALUE(binMedian, cumulative->nums, totalDataPoints/2.0, 0); 734 734 psTrace("psLib.math", 6, "The median bin is %ld (%.2f to %.2f)\n", binMedian, 735 735 cumulative->bounds->data.F32[binMedian], cumulative->bounds->data.F32[binMedian+1]); … … 737 737 // ADD step 3: Interpolate to the exact 50% position: this is the robust histogram median. 738 738 stats->robustMedian = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, 739 binMedian, totalDataPoints/2.0);739 binMedian, totalDataPoints/2.0); 740 740 if (isnan(stats->robustMedian)) { 741 741 psError(PS_ERR_UNKNOWN, false, "Failed to fit a quadratic and calculate the 50-percent position.\n"); 742 goto escape;742 goto escape; 743 743 } 744 744 psTrace("psLib.math", 6, "Current robust median is %f\n", stats->robustMedian); 745 745 746 746 // ADD step 4: Find the bins which contains the 15.8655% (-1 sigma) and 84.1345% (+1 sigma) data points 747 // also find the 30.8538% (-0.5 sigma) and 69.1462% (+0.5 sigma) points748 PS_BIN_FOR_VALUE(binLo, cumulative->nums, totalDataPoints * 0.158655f, -1);749 PS_BIN_FOR_VALUE(binHi, cumulative->nums, totalDataPoints * 0.841345f, +1);750 751 PS_BIN_FOR_VALUE(binL2, cumulative->nums, totalDataPoints * 0.308538f, -1);752 PS_BIN_FOR_VALUE(binH2, cumulative->nums, totalDataPoints * 0.691462f, +1);747 // also find the 30.8538% (-0.5 sigma) and 69.1462% (+0.5 sigma) points 748 PS_BIN_FOR_VALUE(binLo, cumulative->nums, totalDataPoints * 0.158655f, -1); 749 PS_BIN_FOR_VALUE(binHi, cumulative->nums, totalDataPoints * 0.841345f, +1); 750 751 PS_BIN_FOR_VALUE(binL2, cumulative->nums, totalDataPoints * 0.308538f, -1); 752 PS_BIN_FOR_VALUE(binH2, cumulative->nums, totalDataPoints * 0.691462f, +1); 753 753 754 754 psTrace("psLib.math", 6, "The 15.8655%% and 84.1345%% data point bins are (%ld, %ld).\n", … … 759 759 if ((binLo < 0) || (binHi < 0)) { 760 760 psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 15.8655%% and 84.1345%% data points.\n"); 761 goto escape;761 goto escape; 762 762 } 763 763 … … 769 769 770 770 // find the +0.5 and -0.5 sigma points with linear interpolation. binLo and binHi are the bins 771 // containing the value of interest. constrain the answer to be within the current bin.772 // (extrapolation should not be needed and will result in errors)771 // containing the value of interest. constrain the answer to be within the current bin. 772 // (extrapolation should not be needed and will result in errors) 773 773 float binLoF32, binHiF32, binL2F32, binH2F32; 774 PS_BIN_INTERPOLATE (binLoF32, cumulative->nums, cumulative->bounds, binL2, totalDataPoints * 0.158655f);775 PS_BIN_INTERPOLATE (binHiF32, cumulative->nums, cumulative->bounds, binH2, totalDataPoints * 0.841345f);776 PS_BIN_INTERPOLATE (binL2F32, cumulative->nums, cumulative->bounds, binL2, totalDataPoints * 0.308538f);777 PS_BIN_INTERPOLATE (binH2F32, cumulative->nums, cumulative->bounds, binH2, totalDataPoints * 0.691462f);774 PS_BIN_INTERPOLATE (binLoF32, cumulative->nums, cumulative->bounds, binL2, totalDataPoints * 0.158655f); 775 PS_BIN_INTERPOLATE (binHiF32, cumulative->nums, cumulative->bounds, binH2, totalDataPoints * 0.841345f); 776 PS_BIN_INTERPOLATE (binL2F32, cumulative->nums, cumulative->bounds, binL2, totalDataPoints * 0.308538f); 777 PS_BIN_INTERPOLATE (binH2F32, cumulative->nums, cumulative->bounds, binH2, totalDataPoints * 0.691462f); 778 778 779 779 // report +/- 1 sigma points … … 810 810 // Free the histograms; they will be recreated on the next iteration, with new bounds 811 811 psFree(histogram); 812 histogram = NULL;812 histogram = NULL; 813 813 814 814 psFree(cumulative); 815 cumulative = NULL;815 cumulative = NULL; 816 816 } else { 817 817 // We've got the bin size correct now … … 823 823 // XXX test lines while studying algorithm errors 824 824 // fprintf (stderr, "robust stats test %7.1f +/- %7.1f : %4ld %4ld %4ld %4ld %4ld : %f %f %f\n", 825 // stats->robustMedian, stats->robustStdev, 825 // stats->robustMedian, stats->robustStdev, 826 826 // binLo, binL2, binMedian, binH2, binHi, binSize, max, min); 827 827 … … 835 835 // positions. 836 836 psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, 837 binLo25, totalDataPoints * 0.25f);837 binLo25, totalDataPoints * 0.25f); 838 838 psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, 839 binHi25, totalDataPoints * 0.75f);839 binHi25, totalDataPoints * 0.75f); 840 840 if (isnan(binLo25F32) || isnan(binHi25F32)) { 841 841 psError(PS_ERR_UNKNOWN, false, 842 842 "could not determine the robustUQ: fitQuadraticSearchForYThenReturnX() returned a NAN.\n"); 843 goto escape;843 goto escape; 844 844 } 845 845 … … 851 851 for (long i = 0 ; i < myVector->n ; i++) { 852 852 if (!mask->data.U8[i] && 853 (binLo25F32 <= myVector->data.F32[i]) && (binHi25F32 >= myVector->data.F32[i])) {853 (binLo25F32 <= myVector->data.F32[i]) && (binHi25F32 >= myVector->data.F32[i])) { 854 854 N50++; 855 855 } … … 986 986 } 987 987 988 // select the min and max bins, saturating on the lower and upper end-points989 long binMin, binMax;990 PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, -1);991 PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, +1);988 // select the min and max bins, saturating on the lower and upper end-points 989 long binMin, binMax; 990 PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, -1); 991 PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, +1); 992 992 993 993 // Generate the variables that will be used in the Gaussian fitting … … 1164 1164 } 1165 1165 1166 // select the min and max bins, saturating on the lower and upper end-points1167 long binMin, binMax;1168 PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, -1);1169 PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, +1);1166 // select the min and max bins, saturating on the lower and upper end-points 1167 long binMin, binMax; 1168 PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, -1); 1169 PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, +1); 1170 1170 1171 1171 // Generate the variables that will be used in the Gaussian fitting … … 1357 1357 } 1358 1358 1359 // select the min and max bins, saturating on the lower and upper end-points1360 long binMin, binMax;1361 PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, -1);1362 PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, +1);1359 // select the min and max bins, saturating on the lower and upper end-points 1360 long binMin, binMax; 1361 PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, -1); 1362 PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, +1); 1363 1363 1364 1364 // search for mode (peak of histogram within range mean-2sigma - mean+2sigma … … 1494 1494 1495 1495 // calculate upper mean & stdev from parabolic fit -- use this as the result 1496 #ifndef PS_NO_TRACE 1496 1497 float upperStdev = sqrt(-0.5/poly->coeff[2]); 1497 1498 float upperMean = poly->coeff[1]*PS_SQR(upperStdev); … … 1499 1500 psTrace("psLib.math", 6, "The upper mean is %f.\n", upperMean); 1500 1501 psTrace("psLib.math", 6, "The upper stdev is %f.\n", upperStdev); 1502 #endif 1501 1503 1502 1504 psFree (poly); … … 2043 2045 *****************************************************************************/ 2044 2046 static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, 2045 psVector *yVec,2046 psS32 binNum,2047 psF32 yVal2047 psVector *yVec, 2048 psS32 binNum, 2049 psF32 yVal 2048 2050 ) 2049 2051 { … … 2084 2086 // 2085 2087 if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) || 2086 ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) {2088 ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) { 2087 2089 psError(PS_ERR_BAD_PARAMETER_VALUE, true, 2088 2090 _("Specified yVal, %g, is not within y-range, %g to %g."), … … 2094 2096 // 2095 2097 if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) || 2096 ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) {2098 ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) { 2097 2099 psError(PS_ERR_UNKNOWN, true, 2098 2100 "This routine must be called with monotonically increasing or decreasing data points.\n");
Note:
See TracChangeset
for help on using the changeset viewer.
