Changeset 4958 for trunk/psLib/src/math/psStats.c
- Timestamp:
- Sep 7, 2005, 11:35:50 AM (21 years ago)
- File:
-
- 1 edited
-
trunk/psLib/src/math/psStats.c (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/math/psStats.c
r4898 r4958 14 14 * stats->binsize 15 15 * 16 * @version $Revision: 1.1 39$ $Name: not supported by cvs2svn $17 * @date $Date: 2005-0 8-30 01:14:13$16 * @version $Revision: 1.140 $ $Name: not supported by cvs2svn $ 17 * @date $Date: 2005-09-07 21:35:50 $ 18 18 * 19 19 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 46 46 /* DEFINE STATEMENTS */ 47 47 /*****************************************************************************/ 48 #define PS_GAUSS_WIDTH 5 // The width of the Gaussian or boxcar smoothing. 48 #define PS_GAUSS_WIDTH 5 // The width of the Gaussian smoothing. 49 // This corresponds to N in the ADD. 49 50 #define PS_CLIPPED_NUM_ITER_LB 1 50 51 #define PS_CLIPPED_NUM_ITER_UB 10 … … 292 293 max of the input vector. If there was a problem with the max calculation, 293 294 this routine sets stats->max to NAN. 295 296 XXX: Do we need to factor errors into it? 294 297 *****************************************************************************/ 295 298 psS32 p_psVectorMax(const psVector* myVector, … … 614 617 /****************************************************************************** 615 618 p_psVectorSmoothHistGaussian(): This routine smoothes the data in the input 616 robustHistogram with a Gaussian of width sigma. 619 robustHistogram with a Gaussian of width sigma. It returns a psVector of the 620 smoothed data. 617 621 618 622 XXX: Only PS_TYPE_F32 is supported. … … 621 625 call that. Is that possible? 622 626 *****************************************************************************/ 623 psVector * p_psVectorSmoothHistGaussian(psHistogram* robustHistogram,627 psVector *p_psVectorSmoothHistGaussian(psHistogram *histogram, 624 628 psF32 sigma) 625 629 { 626 PS_ASSERT_PTR_NON_NULL(robustHistogram, NULL); 627 PS_ASSERT_PTR_NON_NULL(robustHistogram->bounds, NULL); 628 629 psS32 i = 0; // Loop index variable 630 psS32 j = 0; // Loop index variable 631 psF32 iMid; 632 psF32 jMid; 633 psS32 numBins = robustHistogram->nums->n; 634 psS32 numBounds = robustHistogram->bounds->n; 635 psVector* smooth = psVectorAlloc(numBins, PS_TYPE_F32); 630 PS_ASSERT_PTR_NON_NULL(histogram, NULL); 631 PS_ASSERT_PTR_NON_NULL(histogram->bounds, NULL); 632 PS_ASSERT_PTR_NON_NULL(histogram->nums, NULL); 633 634 psS32 numBins = histogram->nums->n; 635 psS32 numBounds = histogram->bounds->n; 636 psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32); 637 psF32 firstBound = histogram->bounds->data.F32[0]; 638 psF32 lastBound = histogram->bounds->data.F32[numBounds-1]; 639 psScalar x; 640 x.type.type = PS_TYPE_F32; 636 641 psS32 jMin = 0; 637 642 psS32 jMax = 0; 638 psF32 firstBound = robustHistogram->bounds->data.F32[0]; 639 psF32 lastBound = robustHistogram->bounds->data.F32[numBounds-1]; 643 644 if (histogram->uniform == false) { 645 // 646 // We get here if the histogram is non-uniform. 647 // 648 649 for (psS32 i = 0; i < numBins; i++) { 650 // Determine the midpoint of bin i. 651 psS32 iMid = PS_BIN_MIDPOINT(histogram, i); 652 653 // 654 // We determine the bin numbers (jMin:jMax) corresponding to a 655 // range of data values surrounding iMid. The range is of size: 656 // 2*PS_GAUSS_WIDTH*sigma 657 // 658 x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma); 659 if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { 660 jMin = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x); 661 if (jMin < 0) { 662 psError(PS_ERR_UNEXPECTED_NULL, 663 false, 664 PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); 665 return(NULL); 666 } 667 } else if (x.data.F32 <= firstBound) { 668 jMin = 0; 669 } else if (x.data.F32 >= lastBound) { 670 jMin = histogram->bounds->n - 1; 671 } 672 673 x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma); 674 if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { 675 jMax = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x); 676 if (jMax < 0) { 677 psError(PS_ERR_UNEXPECTED_NULL, 678 false, 679 PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); 680 return(NULL); 681 } 682 } else if (x.data.F32 <= firstBound) { 683 jMax = 0; 684 } else if (x.data.F32 >= lastBound) { 685 jMax = histogram->bounds->n - 1; 686 } 687 688 // 689 // Loop from jMin to jMax, computing the gaussian of data i. 690 // 691 smooth->data.F32[i] = 0.0; 692 for (psS32 j = jMin ; j <= jMax ; j++) { 693 psS32 jMid = PS_BIN_MIDPOINT(histogram, j); 694 smooth->data.F32[i] += 695 histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); 696 } 697 } 698 } else { 699 // 700 // We get here if the histogram is uniform. 701 // 702 703 for (psS32 i = 0; i < numBins; i++) { 704 psF32 binSize = histogram->bounds->data.F32[1] - histogram->bounds->data.F32[0]; 705 psS32 gaussWidth = (psS32) ((PS_GAUSS_WIDTH * sigma) / binSize); 706 707 // 708 // We determine the bin numbers (jMin:jMax) corresponding to a 709 // range of data values surrounding iMid. The range is of size: 710 // 2*PS_GAUSS_WIDTH*sigma 711 // 712 psS32 jMin = i - gaussWidth; 713 if (jMin < 0 ) { 714 jMin = 0; 715 } 716 psS32 jMax = i + gaussWidth; 717 if (jMax > (histogram->bounds->n - 1)) { 718 jMax = (histogram->bounds->n - 1); 719 } 720 721 // 722 // Loop from jMin to jMax, computing the gaussian of data i. 723 // 724 smooth->data.F32[i] = 0.0; 725 psS32 iMid = PS_BIN_MIDPOINT(histogram, i); 726 for (psS32 j = jMin ; j <= jMax ; j++) { 727 psS32 jMid = PS_BIN_MIDPOINT(histogram, j); 728 smooth->data.F32[i] += 729 histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); 730 } 731 } 732 } 733 734 return(smooth); 735 } 736 /****************************************************************************** 737 p_psVectorSmoothHistGaussianNEW(): This routine smoothes the data in the input 738 robustHistogram with a Gaussian of width sigma. It returns a psVector of the 739 smoothed data. 740 741 XXX: Only PS_TYPE_F32 is supported. 742 743 XXX: Write a general routine which smoothes a psVector. This routine should 744 call that. Is that possible? 745 *****************************************************************************/ 746 psVector *p_psVectorSmoothHistGaussianNEW(psHistogram *histogram, 747 psF32 sigma) 748 { 749 PS_ASSERT_PTR_NON_NULL(histogram, NULL); 750 PS_ASSERT_PTR_NON_NULL(histogram->bounds, NULL); 751 PS_ASSERT_PTR_NON_NULL(histogram->nums, NULL); 752 753 psS32 numBins = histogram->nums->n; 754 psS32 numBounds = histogram->bounds->n; 755 psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32); 756 psF32 firstBound = histogram->bounds->data.F32[0]; 757 psF32 lastBound = histogram->bounds->data.F32[numBounds-1]; 640 758 psScalar x; 641 642 759 x.type.type = PS_TYPE_F32; 643 for (i = 0; i < numBins; i++) { 644 // Determine the midpoint of bin i. 645 iMid = (robustHistogram->bounds->data.F32[i] + 646 robustHistogram->bounds->data.F32[i+1]) / 2.0; 647 648 649 // We determine the bin numbers corresponding to a range of data 650 // values surrounding iMid. The ranges is of size 651 // s*PS_GAUSS_WIDTH*sigma 652 653 // YYY: The p_psVectorBinDisect() routine does much of the work of 654 // the following conditionals, however, it also reports a warning 655 // message. I don't want the warning message so I reproduce the 656 // conditionals here. Maybe p_psVectorBinDisect() should not produce 657 // warnings? 658 659 x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma); 660 if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { 661 jMin = p_psVectorBinDisect( *(psVector* *)&robustHistogram->bounds, &x); 662 if (jMin < 0) { 663 psError(PS_ERR_UNEXPECTED_NULL, 664 false, 665 PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); 666 return(NULL); 667 } 668 } else if (x.data.F32 <= firstBound) { 669 jMin = 0; 670 } else if (x.data.F32 >= lastBound) { 671 jMin = robustHistogram->bounds->n - 1; 672 } 673 674 x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma); 675 if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { 676 jMax = p_psVectorBinDisect( *(psVector* *)&robustHistogram->bounds, &x); 677 if (jMax < 0) { 678 psError(PS_ERR_UNEXPECTED_NULL, 679 false, 680 PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); 681 return(NULL); 682 } 683 } else if (x.data.F32 <= firstBound) { 684 jMax = 0; 685 } else if (x.data.F32 >= lastBound) { 686 jMax = robustHistogram->bounds->n - 1; 687 } 688 689 smooth->data.F32[i] = 0.0; 690 for (j = jMin ; j <= jMax ; j++) { 691 jMid = (robustHistogram->bounds->data.F32[j] + 692 robustHistogram->bounds->data.F32[j+1]) / 2.0; 693 smooth->data.F32[i] += 694 robustHistogram->nums->data.F32[j] * 695 psGaussian(jMid, iMid, sigma, true); 760 psS32 jMin = 0; 761 psS32 jMax = 0; 762 763 if (histogram->uniform == false) { 764 // 765 // We get here if the histogram is non-uniform. 766 // 767 768 for (psS32 i = 0; i < numBins; i++) { 769 // Determine the midpoint of bin i. 770 psS32 iMid = PS_BIN_MIDPOINT(histogram, i); 771 772 // 773 // We determine the bin numbers (jMin:jMax) corresponding to a 774 // range of data values surrounding iMid. The range is of size: 775 // 2*PS_GAUSS_WIDTH*sigma 776 // 777 x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma); 778 if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { 779 jMin = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x); 780 if (jMin < 0) { 781 psError(PS_ERR_UNEXPECTED_NULL, 782 false, 783 PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); 784 return(NULL); 785 } 786 } else if (x.data.F32 <= firstBound) { 787 jMin = 0; 788 } else if (x.data.F32 >= lastBound) { 789 jMin = histogram->bounds->n - 1; 790 } 791 792 x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma); 793 if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { 794 jMax = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x); 795 if (jMax < 0) { 796 psError(PS_ERR_UNEXPECTED_NULL, 797 false, 798 PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); 799 return(NULL); 800 } 801 } else if (x.data.F32 <= firstBound) { 802 jMax = 0; 803 } else if (x.data.F32 >= lastBound) { 804 jMax = histogram->bounds->n - 1; 805 } 806 807 // 808 // Loop from jMin to jMax, computing the gaussian of data i. 809 // 810 smooth->data.F32[i] = 0.0; 811 for (psS32 j = jMin ; j <= jMax ; j++) { 812 psS32 jMid = PS_BIN_MIDPOINT(histogram, j); 813 smooth->data.F32[i] += 814 histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); 815 } 816 } 817 } else { 818 // 819 // We get here if the histogram is uniform. 820 // 821 822 for (psS32 i = 0; i < numBins; i++) { 823 psF32 binSize = histogram->bounds->data.F32[1] - histogram->bounds->data.F32[0]; 824 psS32 gaussWidth = (psS32) ((PS_GAUSS_WIDTH * sigma) / binSize); 825 826 // 827 // We determine the bin numbers (jMin:jMax) corresponding to a 828 // range of data values surrounding iMid. The range is of size: 829 // 2*PS_GAUSS_WIDTH*sigma 830 // 831 psS32 jMin = i - gaussWidth; 832 if (jMin < 0 ) { 833 jMin = 0; 834 } 835 psS32 jMax = i + gaussWidth; 836 if (jMax > (histogram->bounds->n - 1)) { 837 jMax = (histogram->bounds->n - 1); 838 } 839 840 // 841 // Loop from jMin to jMax, computing the gaussian of data i. 842 // 843 smooth->data.F32[i] = 0.0; 844 psS32 iMid = PS_BIN_MIDPOINT(histogram, i); 845 for (psS32 j = jMin ; j <= jMax ; j++) { 846 psS32 jMid = PS_BIN_MIDPOINT(histogram, j); 847 smooth->data.F32[i] += 848 histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); 849 } 696 850 } 697 851 } … … 1142 1296 These macros and functions define the following functions: 1143 1297 1144 <p_psNormalizeVectorRange(myData, low, high)1298 p_psNormalizeVectorRange(myData, low, high) 1145 1299 1146 1300 That assumes that the low/high arguments are PS_TYPE_F64; the vector myData … … 1251 1405 1252 1406 /****************************************************************************** 1253 p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes 1254 as input a 1-D polynomial of arbitrary order (though we are using 2nd-order 1255 polynomials here) and a range of x-values for which it is defined: 1256 [rangeLow, rangeHigh]. It determines the x-value of that polynomial such 1257 that f(x) == midpoint. This functions uses a binary-search algorithm on the 1258 range and assumes that the polynomial is monotonically increasing or 1259 decreasing within that range. 1407 p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, getThisValue): This routine 1408 takes as input a 1-D polynomial of arbitrary order and a range of x-values for 1409 which it is defined: [rangeLow, rangeHigh]. It determines the x-value of 1410 that polynomial such that f(x) == getThisValue. This function uses a 1411 binary-search algorithm on the range and assumes that the polynomial is 1412 monotonically increasing or decreasing within that range. 1260 1413 1261 1414 XXX: Terminate when f(x)-getThisValue is within some error tolerance. … … 1269 1422 { 1270 1423 PS_ASSERT_POLY_NON_NULL(myPoly, NAN); 1271 PS_ FLOAT_COMPARE(rangeLow, rangeHigh, NAN);1424 PS_ASSERT_FLOAT_LARGER_THAN(rangeHigh, rangeLow, NAN); 1272 1425 // We ensure that the requested f(y) value, which is getThisValue, is 1273 1426 // falls within the range of y-values of the polynomial "myPoly" in the 1274 1427 // specified x-range (rangeLow:rangeHigh). 1275 psF32 fLo = psPolynomial1DEval( 1276 myPoly, 1277 rangeLow 1278 ); 1279 psF32 fHi = psPolynomial1DEval( 1280 myPoly, 1281 rangeHigh 1282 ); 1428 psF32 fLo = psPolynomial1DEval(myPoly, rangeLow); 1429 psF32 fHi = psPolynomial1DEval(myPoly, rangeHigh); 1283 1430 if (!((fLo <= getThisValue) && (fHi >= getThisValue))) { 1284 1431 psError(PS_ERR_UNKNOWN, … … 1300 1447 oldMidpoint = midpoint; 1301 1448 1302 f = psPolynomial1DEval( 1303 myPoly, 1304 midpoint 1305 ); 1449 f = psPolynomial1DEval(myPoly, midpoint); 1306 1450 if (fabs(f - getThisValue) <= FLT_EPSILON) { 1307 1451 return (midpoint); … … 1331 1475 XXX: the vectors do not have to be the same length. Must insert the proper 1332 1476 tests to ensure that binNum is within acceptable ranges for both vectors. 1477 1478 XXX: This currently assumes that the three points are monotonically increasing 1479 or decreasing: so, it works for the cumulative histogram vectors, but not for 1480 arbitrary vectors. We should probably test that condition. 1333 1481 *****************************************************************************/ 1334 1482 psF32 fitQuadraticSearchForYThenReturnX(psVector *xVec, … … 1345 1493 PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (yVec->n - 1), NAN); 1346 1494 1347 // PS_VECTOR_DECLARE_ALLOC_STATIC(x, 3, PS_TYPE_F64);1348 // PS_VECTOR_DECLARE_ALLOC_STATIC(y, 3, PS_TYPE_F64);1349 // PS_VECTOR_DECLARE_ALLOC_STATIC(yErr, 3, PS_TYPE_F64);1350 // PS_POLY_1D_DECLARE_ALLOC_STATIC(myPoly, 2, PS_POLYNOMIAL_ORD);1351 1495 psVector *x = psVectorAlloc(3, PS_TYPE_F64); 1352 1496 psVector *y = psVectorAlloc(3, PS_TYPE_F64); … … 1364 1508 y->data.F64[1] = yVec->data.F32[binNum]; 1365 1509 y->data.F64[2] = yVec->data.F32[binNum + 1]; 1510 1511 // 1512 // Ensure that the y values are monotonic. 1513 // 1514 // XXX: This routine should probably be rewritten in a more general fashion 1515 // so that the folloiwng checks are not necessary. 1516 // 1517 if (y->data.F64[0] < y->data.F64[1]) { 1518 if (!(y->data.F64[1] <= y->data.F64[2])) { 1519 psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n"); 1520 psFree(myPoly); 1521 psFree(x); 1522 psFree(y); 1523 psFree(yErr); 1524 return(NAN); 1525 } 1526 } else { 1527 if (!(y->data.F64[1] >= y->data.F64[2])) { 1528 psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n"); 1529 psFree(myPoly); 1530 psFree(x); 1531 psFree(y); 1532 psFree(yErr); 1533 return(NAN); 1534 } 1535 } 1366 1536 1367 1537 // Ensure that yVal is within the range of the bins we are using. … … 1637 1807 // code no longer produces sensible results. 1638 1808 // XXX: Since we are no longer fitting a 1-D Gaussian, we can probably 1639 // remove some of the above code that calculated the initial estimate1809 // remove some of the above code that calculated the initial estimate 1640 1810 // for the mean and sigma. 1641 1811 … … 1645 1815 int index = i - modeBinNum + dL; 1646 1816 // XXX: Should this be the natural log? 1647 y->data.F32[index] = robustHistogramVector->data.F32[i];1648 //y->data.F32[index] = logf(robustHistogramVector->data.F32[i]);1817 // y->data.F32[index] = robustHistogramVector->data.F32[i]; 1818 y->data.F32[index] = logf(robustHistogramVector->data.F32[i]); 1649 1819 x->data.F32[index] = (psF32) index; 1650 1820 } … … 1742 1912 psFree(robustHistogramVector); 1743 1913 psFree(cumulativeRobustSums); 1914 1744 1915 return(0); 1745 1916 } 1917 1918 1919 /***************************************************************************** 1920 XXX: Is there a psLib function for this? 1921 *****************************************************************************/ 1922 psVector *PsVectorDup(psVector *in) 1923 { 1924 psVector *out = psVectorAlloc(in->n, in->type.type); 1925 1926 if (in->type.type == PS_TYPE_F32) { 1927 for (psS32 i = 0 ; i < in->n ; i++) { 1928 out->data.F32[i] = in->data.F32[i]; 1929 } 1930 } else if (in->type.type == PS_TYPE_F64) { 1931 for (psS32 i = 0 ; i < in->n ; i++) { 1932 out->data.F64[i] = in->data.F64[i]; 1933 } 1934 } else { 1935 printf("XXX: Generate an error here.\n"); 1936 return(NULL); 1937 } 1938 return(out); 1939 } 1940 1941 /****************************************************************************** 1942 XXX: This function need to be written. Actually, it simply needs to be 1943 retrieved from the CVS repository, since it was written earlier, then 1944 discarded. At present, it was deleted from the CVS repository, so we might 1945 have to retrieve it from tape. 1946 *****************************************************************************/ 1947 psVector *Fit1DGaussian(psVector *x, psVector*y) 1948 { 1949 printf("XXX: Generate an error here.\n"); 1950 printf("XXX: Error: This function was previously part of psStats.c, was removed, was purged from CVS, and now needs to be retrieved from tape.\n"); 1951 return(NULL); 1952 } 1953 1954 /****************************************************************************** 1955 1956 p_psVectorRobustStatsNew(myVector, maskVector, maskVal, stats): This is the new 1957 version of the robust stats routine. 1958 1959 XXX: MUST DO: If the errors in the input values are known, then the same 1960 approach is used, except that the histograms become probability density 1961 functions (PDFs). In this case, the input values are spread out, so that they 1962 do not simply contribute a single unit to the histogram, but rather contribute 1963 a fraction of a value, equivalent to the weight. In the interests of speed, a 1964 boxcar PDF may be used to represent each input value (as opposed to a 1965 Gaussian), where the boxcar width is equal to 2p2 ln 2 times the error and 1966 each input value contributes constant area. Then the robust median and 1967 standard deviation are estimated in the same manner as above. 1968 1969 XXX: Check for errors in psLib routines that we call. 1970 *****************************************************************************/ 1971 psS32 p_psVectorRobustStatsNew(const psVector* myVector, 1972 const psVector* errors, 1973 const psVector* maskVector, 1974 psU32 maskVal, 1975 psStats* stats) 1976 { 1977 psHistogram *robustHistogram = NULL; 1978 psHistogram *cumulativeRobustHistogram = NULL; 1979 psS32 numBins = 0; 1980 psScalar *tmpScalar = psScalarAlloc(0.0, PS_TYPE_F32); 1981 tmpScalar->type.type = PS_TYPE_F32; 1982 psS32 totalDataPoints = 0; 1983 psS32 rc = 0; 1984 psVector *tmpMaskVec = PsVectorDup((psVector *) maskVector); 1985 1986 while (1) { 1987 // 1988 // Determine the bin size of the robust histogram. This is done 1989 // by computing the total range of data values and dividing by 1000.0. 1990 // 1991 psStats* tmpStatsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); 1992 rc = p_psVectorMin(myVector, tmpMaskVec, maskVal, tmpStatsMinMax); 1993 rc|= p_psVectorMax(myVector, tmpMaskVec, maskVal, tmpStatsMinMax); 1994 if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) { 1995 psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); 1996 psFree(tmpStatsMinMax); 1997 psFree(tmpMaskVec); 1998 psFree(tmpScalar); 1999 return(1); 2000 } 2001 psF32 binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / 1000.0f; 2002 2003 // 2004 // If all data points have the same value, then we set the appropiate 2005 // members of stats and return. 2006 // 2007 if (fabs(tmpStatsMinMax->max - tmpStatsMinMax->min) <= FLT_EPSILON) { 2008 if (stats->options & PS_STAT_ROBUST_MEDIAN) { 2009 stats->robustMedian = tmpStatsMinMax->min; 2010 } 2011 if (stats->options & PS_STAT_ROBUST_QUARTILE) { 2012 stats->robustUQ = tmpStatsMinMax->min; 2013 stats->robustLQ = tmpStatsMinMax->min; 2014 } 2015 // XXX: Set these to the number of unmasked data points? 2016 stats->robustNfit = 0.0; 2017 stats->robustN50 = 0.0; 2018 psFree(tmpStatsMinMax); 2019 psFree(tmpMaskVec); 2020 psFree(tmpScalar); 2021 2022 return(0); 2023 } 2024 2025 // 2026 // ADD: Step 0. 2027 // Construct the histogram with the specified bin size. 2028 // 2029 // NOTE: we can not specify the bin size precisely since the argument 2030 // to psHistogramAlloc() is the number of bins, not the binSize. 2031 // If we get here, we know that binSize != 0.0. 2032 // 2033 numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / binSize); 2034 robustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins); 2035 cumulativeRobustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins); 2036 2037 // Populate the histogram array. 2038 psVectorHistogram(robustHistogram, myVector, errors, tmpMaskVec, maskVal); 2039 2040 // 2041 // ADD: Step 1. 2042 // Construct the cumulative histogram from the specific histogram 2043 // 2044 cumulativeRobustHistogram->nums->data.F32[0] = robustHistogram->nums->data.F32[0]; 2045 for (psS32 i = 1 ; i < robustHistogram->nums->n ; i++) { 2046 cumulativeRobustHistogram->nums->data.F32[i] = cumulativeRobustHistogram->nums->data.F32[i-1] + 2047 robustHistogram->nums->data.F32[i]; 2048 } 2049 2050 // 2051 // ADD: Step 2. 2052 // Find the bin which contains the 50% data point. 2053 // 2054 totalDataPoints = cumulativeRobustHistogram->nums->data.F32[numBins - 1]; 2055 tmpScalar->data.F32 = totalDataPoints/2.0; 2056 psS32 binMedian = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); 2057 if (binMedian != 0) { 2058 psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 50% data point.\n"); 2059 psFree(tmpStatsMinMax); 2060 psFree(robustHistogram); 2061 psFree(cumulativeRobustHistogram); 2062 psFree(tmpScalar); 2063 return(1); 2064 } 2065 2066 // 2067 // ADD: Step 3. 2068 // Interpolate to the exact 50% position: this is the robust histogram median. 2069 // XXX: Check for errors here! 2070 // 2071 stats->robustMedian = fitQuadraticSearchForYThenReturnX( 2072 *(psVector* *)&robustHistogram->bounds, 2073 *(psVector* *)&robustHistogram->nums, 2074 binMedian, 2075 totalDataPoints/2.0); 2076 2077 // 2078 // ADD: Step 4. 2079 // Find the bins which contains the 15.8655% and 84.1345% data points. 2080 // 2081 tmpScalar->data.F32 = totalDataPoints * 0.158655f; 2082 psS32 binLo = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); 2083 tmpScalar->data.F32 = totalDataPoints * 0.841345f; 2084 psS32 binHi = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); 2085 if ((binLo != 0) || (binHi != 0)) { 2086 psError(PS_ERR_UNKNOWN, false, "Failed to calculate the15.8655% and 84.1345% data point\n"); 2087 psFree(tmpStatsMinMax); 2088 psFree(robustHistogram); 2089 psFree(cumulativeRobustHistogram); 2090 psFree(tmpScalar); 2091 return(1); 2092 } 2093 2094 // 2095 // ADD: Step 4b. 2096 // Interpolate Sigma to find these two positions exactly: these are the 1sigma positions. 2097 // 2098 psF32 binLoF32 = fitQuadraticSearchForYThenReturnX( 2099 *(psVector* *)&robustHistogram->bounds, 2100 *(psVector* *)&robustHistogram->nums, 2101 binLo, 2102 totalDataPoints * 0.158655f); 2103 psF32 binHiF32 = fitQuadraticSearchForYThenReturnX( 2104 *(psVector* *)&robustHistogram->bounds, 2105 *(psVector* *)&robustHistogram->nums, 2106 binHi, 2107 totalDataPoints * 0.841345f); 2108 2109 // 2110 // ADD: Step 5. 2111 // Determine SIGMA as 1/2 of the distance between these positions. 2112 // 2113 psF32 sigma = (binHiF32 - binLoF32) / 2.0; 2114 2115 // 2116 // ADD: Step 6. 2117 // If the measured SIGMA is less than 2 times the bin size, exclude 2118 // points which are more than 25 bins from the median, 2119 // recalculate the bin size, and perform the algorithm again. 2120 // 2121 if (sigma < (2 * binSize)) { 2122 psF32 medianLo = robustHistogram->bounds->data.F32[binMedian - 25]; 2123 psF32 medianHi = robustHistogram->bounds->data.F32[binMedian + 25]; 2124 for (psS32 i = 0 ; i < myVector->n ; i++) { 2125 if ((myVector->data.F32[i] < medianLo) || 2126 (myVector->data.F32[i] > medianHi)) { 2127 tmpMaskVec->data.U8[i] = 1; 2128 } 2129 } 2130 } else { 2131 // 2132 // ADD: Step 7. 2133 // Find the bins which contains the 25% and 75% data points. 2134 // 2135 tmpScalar->data.F32 = totalDataPoints * 0.25f; 2136 psS32 binLo25 = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); 2137 tmpScalar->data.F32 = totalDataPoints * 0.75f; 2138 psS32 binHi25 = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); 2139 if ((binLo25 != 0) || (binHi25 != 0)) { 2140 psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 25% and 75% data points\n"); 2141 psFree(tmpStatsMinMax); 2142 psFree(robustHistogram); 2143 psFree(cumulativeRobustHistogram); 2144 psFree(tmpScalar); 2145 return(1); 2146 } 2147 2148 // 2149 // ADD: Step 8. 2150 // Interpolate to find these two positions exactly: these are the upper 2151 // and lower quartile positions. 2152 // XXX: Check for errors. 2153 // 2154 psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX( 2155 *(psVector* *)&robustHistogram->bounds, 2156 *(psVector* *)&robustHistogram->nums, 2157 binLo25, 2158 totalDataPoints * 0.25f); 2159 psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX( 2160 *(psVector* *)&robustHistogram->bounds, 2161 *(psVector* *)&robustHistogram->nums, 2162 binHi25, 2163 totalDataPoints * 0.75f); 2164 2165 stats->robustLQ = binLo25F32; 2166 stats->robustUQ = binHi25F32; 2167 // XXX: No idea how to calculate stats->stdev 2168 2169 // PS_BIN_MIDPOINT(robustHistogram, modeBinNum); 2170 2171 // XXX: I think sumNfit == sumN50 here. 2172 stats->robustNfit = -1; 2173 stats->robustN50 = -1; 2174 2175 // 2176 // Perform the Robust Histogram Statistics algorithm above 2177 // 2178 2179 // 2180 // Smooth the resulting histogram with a Gaussian with SIGMA_s = 1 2181 // bin. 2182 // 2183 // XXX: SIGMA_s is defined nowhere in the ADD. 2184 // 2185 psF32 SIGMA_S = 1.0; 2186 p_psVectorSmoothHistGaussian(robustHistogram, SIGMA_S); 2187 2188 // 2189 // Find the bin with the peak value in the range 2 SIGMA of the 2190 // robust histogram median. 2191 // 2192 // XXX: SIGMA is defined nowhere in the ADD. 2193 // 2194 psF32 SIGMA = 2.0; 2195 psS32 binMin = binMedian - (SIGMA * PS_GAUSS_WIDTH); 2196 if (binMin < 0) { 2197 binMin = 0; 2198 } 2199 psS32 binMax = binMedian + (2 * PS_GAUSS_WIDTH); 2200 if (binMin > (robustHistogram->nums->n - 1)) { 2201 binMin = (robustHistogram->nums->n - 1); 2202 } 2203 psS32 binNum = binNum; 2204 psF32 binMaxNums = robustHistogram->nums->data.F32[binNum]; 2205 for (psS32 i = binNum+1 ; i <= binMax ; i++) { 2206 if (robustHistogram->nums->data.F32[i] > binMaxNums) { 2207 binNum = i; 2208 binMaxNums = robustHistogram->nums->data.F32[i]; 2209 } 2210 } 2211 2212 // 2213 // Fit a Gaussian to the bins in the range 2 SIGMA of the robust 2214 // histogram median. 2215 // 2216 // XXX: SIGMA is defined nowhere in the ADD. 2217 // 2218 psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); 2219 psVector *x = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); 2220 psS32 j = 0; 2221 for (psS32 i = binNum ; i <= binMax ; i++) { 2222 y->data.F32[j] = robustHistogram->nums->data.F32[i]; 2223 x->data.F32[j] = PS_BIN_MIDPOINT(robustHistogram, i); 2224 } 2225 // 2226 // XXX: This function need to be written. Actually, it simply 2227 // needs to be retrieved from the CVS repository, since it was 2228 // written earlier, then discarded. At present, it was deleted 2229 // from the CVS repository, so we might have to retrieve it from 2230 // tape. 2231 // 2232 psVector *params = Fit1DGaussian(x, y); 2233 2234 // 2235 // The robust mean mean_r is derived directly from the fitted 2236 // Gaussian mean. 2237 // 2238 stats->robustMean = params->data.F32[0]; 2239 2240 // 2241 // The robust standard deviation, SIGMA_r is determined by 2242 // subtracting the smoothing scale in quadrature: SIGMA_r^2 = SIGMA^2 2243 // - SIGMA_s^2 2244 // 2245 // XXX: SIGMA and SIGMA_s are defined nowhere in the ADD. We must figure 2246 // out what they are. 2247 // 2248 stats->robustStdev = sqrt(PS_SQR(SIGMA) - PS_SQR(SIGMA_S)); 2249 2250 psFree(tmpStatsMinMax); 2251 psFree(robustHistogram); 2252 psFree(cumulativeRobustHistogram); 2253 psFree(tmpScalar); 2254 psFree(params); 2255 2256 return(0); 2257 } 2258 2259 psFree(tmpStatsMinMax); 2260 psFree(robustHistogram); 2261 psFree(cumulativeRobustHistogram); 2262 } 2263 return(1); 2264 } 2265 2266 2267 2268 2269 2270 2271 2272 1746 2273 1747 2274 /*****************************************************************************/
Note:
See TracChangeset
for help on using the changeset viewer.
