Changeset 15370
- Timestamp:
- Oct 24, 2007, 9:08:24 AM (19 years ago)
- File:
-
- 1 edited
-
trunk/psLib/src/math/psStats.c (modified) (35 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/math/psStats.c
r15126 r15370 13 13 * use ->min and ->max (PS_STAT_USE_RANGE) 14 14 * 15 * @version $Revision: 1.21 8$ $Name: not supported by cvs2svn $16 * @date $Date: 2007- 09-29 23:58:55$15 * @version $Revision: 1.219 $ $Name: not supported by cvs2svn $ 16 * @date $Date: 2007-10-24 19:08:24 $ 17 17 * 18 18 * Copyright 2006 IfA, University of Hawaii … … 152 152 // static prototypes: 153 153 static psF32 minimizeLMChi2Gauss1D(psVector *deriv, const psVector *params, const psVector *coords); 154 static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal); 154 static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, psVector *yVec, psS32 binNum, 155 psF32 yVal); 155 156 156 157 /****************************************************************************** … … 309 310 310 311 // Allocate temporary vectors for the data. 311 psVector *outVector = psVectorAllocEmpty(inVector->n, PS_TYPE_F32); // Vector containing theunmasked values312 psVector *outVector = psVectorAllocEmpty(inVector->n, PS_TYPE_F32); // Vector containing unmasked values 312 313 psF32 *output = outVector->data.F32; // Dereference the vector 313 314 … … 401 402 // If the mean is NAN, then generate a warning and set the stdev to NAN. 402 403 if (isnan(stats->sampleMean)) { 403 psTrace(TRACE, 5, "WARNING: vectorSampleStdev(): sample mean is NAN. Setting stats->sampleStdev = NAN.\n"); 404 psTrace(TRACE, 5, "WARNING: vectorSampleStdev(): sample mean is NAN. " 405 "Setting stats->sampleStdev = NAN.\n"); 404 406 stats->sampleStdev = NAN; 405 407 return true; … … 443 445 // Assume that the user knows what he's doing when he masks out everything --> no error. 444 446 stats->sampleStdev = NAN; 445 psTrace(TRACE, 5, "WARNING: vectorSampleStdev(): no valid psVector elements (%ld). Setting stats->sampleStdev = NAN.\n", count); 447 psTrace(TRACE, 5, "WARNING: vectorSampleStdev(): no valid psVector elements (%ld). " 448 "Setting stats->sampleStdev = NAN.\n", count); 446 449 return true; 447 450 } 448 451 if (count == 1) { 449 452 stats->sampleStdev = 0.0; 450 psTrace(TRACE, 5, "WARNING: vectorSampleStdev(): only one valid psVector elements (%ld). Setting stats->sampleStdev = 0.0.\n", count); 453 psTrace(TRACE, 5, "WARNING: vectorSampleStdev(): only one valid psVector elements (%ld). " 454 "Setting stats->sampleStdev = 0.0.\n", count); 451 455 return true; 452 456 } … … 623 627 // a) Exclude all values x_i for which |x_i - x| > K * stdev 624 628 if (errors) { 625 // XXXX if we convert errors to variance, this should square the other terms (A*A faster than sqrt(A)) 629 // XXXX if we convert errors to variance, this should square the other terms (A*A faster than 630 // sqrt(A)) 626 631 for (long j = 0; j < myVector->n; j++) { 627 632 if (!tmpMask->data.U8[j] && … … 749 754 // Iterate to get the best bin size 750 755 for (int iterate = 1; iterate > 0; iterate++) { 751 psTrace(TRACE, 6, "-------------------- Iterating on Bin size. Iteration number %d --------------------\n", iterate); 756 psTrace(TRACE, 6, 757 "-------------------- Iterating on Bin size. Iteration number %d --------------------\n", 758 iterate); 752 759 753 760 // Get the minimum and maximum values … … 832 839 binMedian, totalDataPoints/2.0); 833 840 if (isnan(stats->robustMedian)) { 834 psError(PS_ERR_UNKNOWN, false, "Failed to fit a quadratic and calculate the 50-percent position.\n"); 841 psError(PS_ERR_UNKNOWN, false, 842 "Failed to fit a quadratic and calculate the 50-percent position.\n"); 835 843 goto escape; 836 844 } 837 845 psTrace(TRACE, 6, "Current robust median is %f\n", stats->robustMedian); 838 846 839 // ADD step 4: Find the bins which contains the 15.8655% (-1 sigma) and 84.1345% (+1 sigma) data points840 // also find the 30.8538% (-0.5 sigma) and 69.1462% (+0.5 sigma) points847 // ADD step 4: Find the bins which contains the 15.8655% (-1 sigma) and 84.1345% (+1 sigma) data 848 // points also find the 30.8538% (-0.5 sigma) and 69.1462% (+0.5 sigma) points 841 849 PS_BIN_FOR_VALUE(binLo, cumulative->nums, totalDataPoints * 0.158655f, 0); 842 850 PS_BIN_FOR_VALUE(binHi, cumulative->nums, totalDataPoints * 0.841345f, 0); … … 858 866 } 859 867 860 // ADD step 4b: Interpolate Sigma (linearly) to find these two positions exactly: these are the 1sigma positions. 868 // ADD step 4b: Interpolate Sigma (linearly) to find these two positions exactly: these are the 1sigma 869 // positions. 861 870 psTrace(TRACE, 6, "binLo is %ld. Nums at that bin and the next are (%.2f, %.2f)\n", 862 871 binLo, cumulative->nums->data.F32[binLo], cumulative->nums->data.F32[binLo+1]); … … 868 877 // (extrapolation should not be needed and will result in errors) 869 878 float binLoF32, binHiF32, binL2F32, binH2F32, binL4F32, binH4F32; 870 PS_BIN_INTERPOLATE (binLoF32, cumulative->nums, cumulative->bounds, binLo, totalDataPoints * 0.158655f); 871 PS_BIN_INTERPOLATE (binHiF32, cumulative->nums, cumulative->bounds, binHi, totalDataPoints * 0.841345f); 872 PS_BIN_INTERPOLATE (binL2F32, cumulative->nums, cumulative->bounds, binL2, totalDataPoints * 0.308538f); 873 PS_BIN_INTERPOLATE (binH2F32, cumulative->nums, cumulative->bounds, binH2, totalDataPoints * 0.691462f); 874 PS_BIN_INTERPOLATE (binL4F32, cumulative->nums, cumulative->bounds, binL4, totalDataPoints * 0.022481f); 875 PS_BIN_INTERPOLATE (binH4F32, cumulative->nums, cumulative->bounds, binH4, totalDataPoints * 0.977519f); 879 PS_BIN_INTERPOLATE (binLoF32, cumulative->nums, cumulative->bounds, binLo, 880 totalDataPoints * 0.158655f); 881 PS_BIN_INTERPOLATE (binHiF32, cumulative->nums, cumulative->bounds, binHi, 882 totalDataPoints * 0.841345f); 883 PS_BIN_INTERPOLATE (binL2F32, cumulative->nums, cumulative->bounds, binL2, 884 totalDataPoints * 0.308538f); 885 PS_BIN_INTERPOLATE (binH2F32, cumulative->nums, cumulative->bounds, binH2, 886 totalDataPoints * 0.691462f); 887 PS_BIN_INTERPOLATE (binL4F32, cumulative->nums, cumulative->bounds, binL4, 888 totalDataPoints * 0.022481f); 889 PS_BIN_INTERPOLATE (binH4F32, cumulative->nums, cumulative->bounds, binH4, 890 totalDataPoints * 0.977519f); 876 891 877 892 // report +/- 1 sigma points … … 891 906 float sigma4 = (binH4F32 - binL4F32) / 4.0; 892 907 893 // take the smallest of the three: if we have a clump with wide outliers, sigma2 and894 // sigma4 will be biased high; if we have a bi-modal distribution, sigma1 and sigma2895 // will be biased high.896 sigma = PS_MIN (sigma1, PS_MIN (sigma2, sigma4));908 // take the smallest of the three: if we have a clump with wide outliers, sigma2 and 909 // sigma4 will be biased high; if we have a bi-modal distribution, sigma1 and sigma2 910 // will be biased high. 911 sigma = PS_MIN (sigma1, PS_MIN (sigma2, sigma4)); 897 912 898 913 psTrace(TRACE, 6, "The 2x sigma is %f.\n", sigma1); … … 1328 1343 1329 1344 if (poly->coeff[2] >= 0.0) { 1330 psTrace(TRACE, 6, "Parabolic fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); 1345 psTrace(TRACE, 6, "Parabolic fit results: %f + %f x + %f x^2\n", 1346 poly->coeff[0], poly->coeff[1], poly->coeff[2]); 1331 1347 psError(PS_ERR_UNKNOWN, false, "fit did not converge\n"); 1332 1348 psFree(x); … … 1344 1360 done = true; 1345 1361 } else { 1346 psTrace(TRACE, 6, "Parabolic fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); 1362 psTrace(TRACE, 6, "Parabolic fit results: %f + %f x + %f x^2\n", 1363 poly->coeff[0], poly->coeff[1], poly->coeff[2]); 1347 1364 psTrace(TRACE, 6, "The new guess mean is %f.\n", guessMean); 1348 1365 psTrace(TRACE, 6, "The new guess stdev is %f.\n", guessStdev); … … 1490 1507 valPeak = histogram->nums->data.F32[binPeak]; 1491 1508 } 1492 psTrace (TRACE, 6, "(%f = %.0f) ", histogram->bounds->data.F32[i], histogram->nums->data.F32[i]);1493 } 1494 psTrace (TRACE, 6, "\n");1509 psTrace (TRACE, 6, "(%f = %.0f) ", histogram->bounds->data.F32[i], histogram->nums->data.F32[i]); 1510 } 1511 psTrace (TRACE, 6, "\n"); 1495 1512 1496 1513 // assume a reasonably well-defined gaussian-like population; run from peak out until val < 0.25*peak 1497 1514 1498 1515 psTrace(TRACE, 6, "The clipped numBins is %ld\n", binMax - binMin); 1499 psTrace(TRACE, 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin); 1500 psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax - 1), binMax - 1); 1501 psTrace(TRACE, 6, "The clipped peak is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binPeak), binPeak); 1516 psTrace(TRACE, 6, "The clipped min is %f (%ld)\n", 1517 PS_BIN_MIDPOINT(histogram, binMin), binMin); 1518 psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", 1519 PS_BIN_MIDPOINT(histogram, binMax - 1), binMax - 1); 1520 psTrace(TRACE, 6, "The clipped peak is %f (%ld)\n", 1521 PS_BIN_MIDPOINT(histogram, binPeak), binPeak); 1502 1522 psTrace(TRACE, 6, "The clipped peak value is %f\n", histogram->nums->data.F32[binPeak]); 1503 1523 … … 1513 1533 } 1514 1534 } 1515 psTrace(TRACE, 6, "Lower bound for lower half: %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binS), binS); 1516 psTrace(TRACE, 6, "Upper bound for lower half: %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binE), binE); 1535 psTrace(TRACE, 6, "Lower bound for lower half: %f (%ld)\n", 1536 PS_BIN_MIDPOINT(histogram, binS), binS); 1537 psTrace(TRACE, 6, "Upper bound for lower half: %f (%ld)\n", 1538 PS_BIN_MIDPOINT(histogram, binE), binE); 1517 1539 1518 1540 psVector *y = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of coordinates … … 1523 1545 continue; 1524 1546 x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i); 1525 y->data.F32[j] = log(histogram->nums->data.F32[i]); // note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2) 1547 // note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2) 1548 y->data.F32[j] = log(histogram->nums->data.F32[i]); 1526 1549 j++; 1527 1550 } … … 1544 1567 1545 1568 if (poly->coeff[2] >= 0.0) { 1546 psTrace(TRACE, 6, "Failed parabolic fit: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); 1569 psTrace(TRACE, 6, "Failed parabolic fit: %f + %f x + %f x^2\n", 1570 poly->coeff[0], poly->coeff[1], poly->coeff[2]); 1547 1571 psFree(poly); 1548 1572 psFree(histogram); 1549 1573 psFree(statsMinMax); 1550 1574 1551 // sometimes, the guessStdev is much too large. in this case, the entire real population tends to be found 1552 // in a single bin. make one attempt to recover by dropping the guessStdev down by a jump and trying again 1575 // sometimes, the guessStdev is much too large. in this case, the entire real population 1576 // tends to be found in a single bin. make one attempt to recover by dropping the guessStdev 1577 // down by a jump and trying again 1553 1578 if (iteration == 0) { 1554 1579 guessStdev = 0.25*guessStdev; … … 1568 1593 done = true; 1569 1594 } 1570 psTrace(TRACE, 6, "Parabolic Lower fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); 1595 psTrace(TRACE, 6, "Parabolic Lower fit results: %f + %f x + %f x^2\n", 1596 poly->coeff[0], poly->coeff[1], poly->coeff[2]); 1571 1597 psTrace(TRACE, 6, "The lower mean is %f.\n", guessMean); 1572 1598 psTrace(TRACE, 6, "The lower stdev is %f.\n", guessStdev); … … 1575 1601 } 1576 1602 1577 // for test, measure the same result for the upper section1603 // for test, measure the same result for the upper section 1578 1604 { 1579 1605 // fit the upper half of the distribution … … 1587 1613 } 1588 1614 } 1589 psTrace(TRACE, 6, "Lower bound for upper half: %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binS), binS); 1590 psTrace(TRACE, 6, "Upper bound for upper half: %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binE), binE); 1615 psTrace(TRACE, 6, "Lower bound for upper half: %f (%ld)\n", 1616 PS_BIN_MIDPOINT(histogram, binS), binS); 1617 psTrace(TRACE, 6, "Upper bound for upper half: %f (%ld)\n", 1618 PS_BIN_MIDPOINT(histogram, binE), binE); 1591 1619 1592 1620 psVector *y = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of coordinates … … 1597 1625 continue; 1598 1626 x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i); 1599 y->data.F32[j] = log(histogram->nums->data.F32[i]); // note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2) 1627 // note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2) 1628 y->data.F32[j] = log(histogram->nums->data.F32[i]); 1600 1629 j++; 1601 1630 } … … 1618 1647 1619 1648 // calculate upper mean & stdev from parabolic fit -- ignore this value 1620 #ifndef PS_NO_TRACE1621 1649 float upperStdev = sqrt(-0.5/poly->coeff[2]); 1622 1650 float upperMean = poly->coeff[1]*PS_SQR(upperStdev); 1623 psTrace(TRACE, 6, "Parabolic Upper fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); 1651 #ifndef PS_NO_TRACE 1652 psTrace(TRACE, 6, "Parabolic Upper fit results: %f + %f x + %f x^2\n", 1653 poly->coeff[0], poly->coeff[1], poly->coeff[2]); 1624 1654 psTrace(TRACE, 6, "The upper mean is %f.\n", upperMean); 1625 1655 psTrace(TRACE, 6, "The upper stdev is %f.\n", upperStdev); 1626 1656 #endif 1627 1657 1628 // if the resulting value is outside of the range binMin - binMax, use the upper value1629 if (done && (guessMean > PS_BIN_MIDPOINT(histogram, binMax - 1))) {1630 guessMean = upperMean;1631 guessStdev = upperStdev;1632 }1658 // if the resulting value is outside of the range binMin - binMax, use the upper value 1659 if (done && (guessMean > PS_BIN_MIDPOINT(histogram, binMax - 1))) { 1660 guessMean = upperMean; 1661 guessStdev = upperStdev; 1662 } 1633 1663 1634 1664 psFree (poly); … … 1718 1748 // Calculate the number of bins. 1719 1749 // XXX can we calculate the binMin, binMax **before** building this histogram? 1720 // if the range is too absurd, adjust numBins & binSize1750 // if the range is too absurd, adjust numBins & binSize 1721 1751 long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize)); 1722 binSize = (max - min) / (float) numBins;1752 binSize = (max - min) / (float) numBins; 1723 1753 psTrace(TRACE, 6, "The new min/max values are (%f, %f).\n", min, max); 1724 1754 psTrace(TRACE, 6, "The new bin size is %f.\n", binSize); … … 1774 1804 valPeak = histogram->nums->data.F32[binPeak]; 1775 1805 } 1776 psTrace (TRACE, 6, "(%f = %.0f) ", histogram->bounds->data.F32[i], histogram->nums->data.F32[i]);1777 } 1778 psTrace (TRACE, 6, "\n");1806 psTrace (TRACE, 6, "(%f = %.0f) ", histogram->bounds->data.F32[i], histogram->nums->data.F32[i]); 1807 } 1808 psTrace (TRACE, 6, "\n"); 1779 1809 1780 1810 // assume a reasonably well-defined gaussian-like population; run from peak out until val < 0.25*peak … … 1782 1812 psTrace(TRACE, 6, "The clipped numBins is %ld\n", binMax - binMin); 1783 1813 psTrace(TRACE, 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin); 1784 psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax - 1), binMax - 1); 1814 psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", 1815 PS_BIN_MIDPOINT(histogram, binMax - 1), binMax - 1); 1785 1816 psTrace(TRACE, 6, "The clipped peak is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binPeak), binPeak); 1786 1817 psTrace(TRACE, 6, "The clipped peak value is %f\n", histogram->nums->data.F32[binPeak]); … … 1804 1835 } 1805 1836 } 1806 psTrace(TRACE, 6, "Lower bound for lower half: %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binS), binS); 1807 psTrace(TRACE, 6, "Upper bound for lower half: %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binE), binE); 1837 psTrace(TRACE, 6, "Lower bound for lower half: %f (%ld)\n", 1838 PS_BIN_MIDPOINT(histogram, binS), binS); 1839 psTrace(TRACE, 6, "Upper bound for lower half: %f (%ld)\n", 1840 PS_BIN_MIDPOINT(histogram, binE), binE); 1808 1841 1809 1842 psVector *y = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of coordinates … … 1814 1847 continue; 1815 1848 x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i); 1816 y->data.F32[j] = log(histogram->nums->data.F32[i]); // note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2) 1849 // note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2) 1850 y->data.F32[j] = log(histogram->nums->data.F32[i]); 1817 1851 j++; 1818 1852 } … … 1835 1869 1836 1870 if (poly->coeff[2] >= 0.0) { 1837 psTrace(TRACE, 6, "Failed parabolic fit: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); 1871 psTrace(TRACE, 6, "Failed parabolic fit: %f + %f x + %f x^2\n", 1872 poly->coeff[0], poly->coeff[1], poly->coeff[2]); 1838 1873 psFree(poly); 1839 1874 psFree(histogram); 1840 1875 psFree(statsMinMax); 1841 1876 1842 // sometimes, the guessStdev is much too large. in this case, the entire real population tends to be found 1843 // in a single bin. make one attempt to recover by dropping the guessStdev down by a jump and trying again 1877 // sometimes, the guessStdev is much too large. in this case, the entire real population 1878 // tends to be found in a single bin. make one attempt to recover by dropping the guessStdev 1879 // down by a jump and trying again 1844 1880 if (iteration == 0) { 1845 1881 guessStdev = 0.25*guessStdev; … … 1859 1895 done = true; 1860 1896 } 1861 psTrace(TRACE, 6, "Parabolic Lower fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); 1897 psTrace(TRACE, 6, "Parabolic Lower fit results: %f + %f x + %f x^2\n", 1898 poly->coeff[0], poly->coeff[1], poly->coeff[2]); 1862 1899 psTrace(TRACE, 6, "The lower mean is %f.\n", guessMean); 1863 1900 psTrace(TRACE, 6, "The lower stdev is %f.\n", guessStdev); … … 1866 1903 } 1867 1904 1868 // if we converge on a solution outside the range binMin - binMax, use a more conservative range1869 float minValue = PS_BIN_MIDPOINT(histogram, binMin);1870 float maxValue = PS_BIN_MIDPOINT(histogram, binMax - 1);1871 1872 if (done && ((guessMean < minValue) || (guessMean > maxValue))) {1905 // if we converge on a solution outside the range binMin - binMax, use a more conservative range 1906 float minValue = PS_BIN_MIDPOINT(histogram, binMin); 1907 float maxValue = PS_BIN_MIDPOINT(histogram, binMax - 1); 1908 1909 if (done && ((guessMean < minValue) || (guessMean > maxValue))) { 1873 1910 psTrace(TRACE, 6, "Inconsistent result, re-trying the fit\n"); 1874 1911 … … 1890 1927 } 1891 1928 } 1892 psTrace(TRACE, 6, "Lower bound for symmetric range: %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binS), binS); 1893 psTrace(TRACE, 6, "Upper bound for symmetric range: %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binE), binE); 1929 psTrace(TRACE, 6, "Lower bound for symmetric range: %f (%ld)\n", 1930 PS_BIN_MIDPOINT(histogram, binS), binS); 1931 psTrace(TRACE, 6, "Upper bound for symmetric range: %f (%ld)\n", 1932 PS_BIN_MIDPOINT(histogram, binE), binE); 1894 1933 1895 1934 psVector *y = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of coordinates … … 1900 1939 continue; 1901 1940 x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i); 1902 y->data.F32[j] = log(histogram->nums->data.F32[i]); // note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2) 1941 // note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2) 1942 y->data.F32[j] = log(histogram->nums->data.F32[i]); 1903 1943 j++; 1904 1944 } … … 1921 1961 1922 1962 // calculate upper mean & stdev from parabolic fit -- ignore this value 1923 #ifndef PS_NO_TRACE1924 1963 guessStdev = sqrt(-0.5/poly->coeff[2]); 1925 1964 guessMean = poly->coeff[1]*PS_SQR(guessStdev); 1926 psTrace(TRACE, 6, "Parabolic Symmetric fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); 1965 #ifndef PS_NO_TRACE 1966 psTrace(TRACE, 6, "Parabolic Symmetric fit results: %f + %f x + %f x^2\n", 1967 poly->coeff[0], poly->coeff[1], poly->coeff[2]); 1927 1968 psTrace(TRACE, 6, "The symmetric mean is %f.\n", guessMean); 1928 1969 psTrace(TRACE, 6, "The symmetric stdev is %f.\n", guessStdev); 1929 1970 #endif 1930 1971 1931 // if we converge on a solution outside the range binMin - binMax, use a more conservative range1932 float minValueSym = PS_BIN_MIDPOINT(histogram, binS);1933 float maxValueSym = PS_BIN_MIDPOINT(histogram, binE - 1);1934 1935 // saturate on min or max value1936 if (guessMean < minValueSym) {1937 guessMean = minValueSym;1938 psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean);1939 }1940 1941 // saturate on min or max value1942 if (guessMean > maxValueSym) {1943 guessMean = maxValueSym;1944 psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean);1945 }1972 // if we converge on a solution outside the range binMin - binMax, use a more conservative range 1973 float minValueSym = PS_BIN_MIDPOINT(histogram, binS); 1974 float maxValueSym = PS_BIN_MIDPOINT(histogram, binE - 1); 1975 1976 // saturate on min or max value 1977 if (guessMean < minValueSym) { 1978 guessMean = minValueSym; 1979 psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean); 1980 } 1981 1982 // saturate on min or max value 1983 if (guessMean > maxValueSym) { 1984 guessMean = maxValueSym; 1985 psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean); 1986 } 1946 1987 1947 1988 psFree (poly); … … 1996 2037 // However, it is also not used anywhere, yet. 1997 2038 // 1998 psWarning("WARNING: vectorSmoothHistGaussian() on non-uniform histograms has not been tested or used.\n"); 2039 psWarning("WARNING: vectorSmoothHistGaussian() on non-uniform " 2040 "histograms has not been tested or used.\n"); 1999 2041 2000 2042 for (long i = 0; i < numBins; i++) { … … 2607 2649 if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) { 2608 2650 psError(PS_ERR_UNEXPECTED_NULL, false, 2609 _("Failed to fit a 1-dimensional polynomial to the three specified data points. Returning NAN.")); 2651 _("Failed to fit a 1-dimensional polynomial to the three specified data points. " 2652 "Returning NAN.")); 2610 2653 psFree(x); 2611 2654 psFree(y);
Note:
See TracChangeset
for help on using the changeset viewer.
