IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 26, 2006, 1:49:11 PM (20 years ago)
Author:
gusciora
Message:

Made the requested changes to the psStats structure and added the fitted
statistics.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/math/psStats.c

    r6204 r6215  
    1414 *      stats->binsize
    1515 *
    16  *  @version $Revision: 1.162 $ $Name: not supported by cvs2svn $
    17  *  @date $Date: 2006-01-26 21:10:22 $
     16 * XXX: Must do
     17 * nSubsample points
     18 * use ->min and ->max (PS_STAT_USE_RANGE)
     19 * use ->binsize (PS_STAT_USE_BINSIZE)
     20 *
     21 *
     22 *
     23 *
     24 *  @version $Revision: 1.163 $ $Name: not supported by cvs2svn $
     25 *  @date $Date: 2006-01-26 23:49:11 $
    1826 *
    1927 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    97105        return true;
    98106
    99     case PS_STAT_ROBUST_MEAN:
    100         *value = stats->robustMean;
    101         psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
    102         return true;
    103 
    104107    case PS_STAT_ROBUST_MEDIAN:
    105108        *value = stats->robustMedian;
    106         psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
    107         return true;
    108 
    109     case PS_STAT_ROBUST_MODE:
    110         *value = stats->robustMode;
    111109        psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
    112110        return true;
    … …  
    449447unmasked element within the specified min/max range).  Otherwise, return
    450448"false".
     449 
     450XXX: Can you use psVectorCountPixelMask here?
    451451 *****************************************************************************/
    452452bool p_psVectorCheckNonEmpty(const psVector* myVector,
    … …  
    503503number of non-masked pixels in the vector that fall within the min/max
    504504range, if specified.
     505 
     506XXX: Can you use psVectorCountPixelMask here?
    505507 *****************************************************************************/
    506508psS32 p_psVectorNValues(const psVector* myVector,
    … …  
    766768        PS_VECTOR_PRINT_F32(smooth);
    767769    }
    768     psTrace(__func__, 4, "---- %s(psVector) end ----\n", __func__);
     770    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
    769771    return(smooth);
    770772}
    … …  
    858860
    859861/******************************************************************************
    860 p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
     862p_psVectorSampleStdevOLD(myVector, maskVector, maskVal, stats): calculates the
    861863stdev of the input vector.
    862864Inputs
    … …  
    868870    NULL
    869871 
     872XXX: remove this
    870873 *****************************************************************************/
    871874void p_psVectorSampleStdevOLD(const psVector* myVector,
    … …  
    11261129        p_psMemSetPersistent(statsTmp, true);
    11271130    } else {
    1128         // EAM : initialize structure if already allocated
     1131        // Initialize structure if already allocated
    11291132        statsTmp->sampleMean = NAN;
    11301133        statsTmp->sampleMedian = NAN;
    … …  
    11321135        statsTmp->sampleUQ = NAN;
    11331136        statsTmp->sampleLQ = NAN;
    1134         statsTmp->robustMean = NAN;
    11351137        statsTmp->robustMedian = NAN;
    1136         statsTmp->robustMode = NAN;
    11371138        statsTmp->robustStdev = NAN;
    11381139        statsTmp->robustUQ = NAN;
    11391140        statsTmp->robustLQ = NAN;
    1140         statsTmp->robustN50 = -1;            // XXX: This is never used
    1141         statsTmp->robustNfit = -1;
     1141        statsTmp->robustN50 = -1;
     1142        statsTmp->fittedMean = NAN;
     1143        statsTmp->fittedStdev = NAN;
     1144        statsTmp->fittedNfit = -1;
    11421145        statsTmp->clippedMean = NAN;
    11431146        statsTmp->clippedStdev = NAN;
    1144         statsTmp->clippedNvalues = -1;     // XXX: This is never used
     1147        statsTmp->clippedNvalues = -1;
    11451148        statsTmp->clipSigma = 3.0;
    11461149        statsTmp->clipIter = 3;
    11471150        statsTmp->min = NAN;
    11481151        statsTmp->max = NAN;
    1149         statsTmp->binsize = NAN;          // XXX: This is never used
     1152        statsTmp->binsize = NAN;
     1153        statsTmp->nSubsample = 100000;
    11501154    }
    11511155
    … …  
    16881692        }
    16891693    }
    1690     psBool iterate = true;
     1694    psS32 iterate = 1;
    16911695    psF32 sigma;
    16921696    psStats *tmpStatsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX);
    16931697
    1694     while (iterate) {
     1698    while (iterate > 0) {
    16951699        psTrace(__func__, 6, "Iterating on Bin size.\n");
    1696         //
    1697         // Determine the bin size of the robust histogram.  This is done
    1698         // by computing the total range of data values and dividing by 1000.0.
    1699         //
    1700         rc = p_psVectorMin(myVector, tmpMaskVec, 1, tmpStatsMinMax);
    1701         rc|= p_psVectorMax(myVector, tmpMaskVec, 1, tmpStatsMinMax);
    1702         if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) {
    1703             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
    1704             psFree(tmpStatsMinMax);
    1705             psFree(tmpMaskVec);
    1706             psFree(tmpScalar);
    1707             psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
    1708             return(1);
    1709         }
    1710         psTrace(__func__, 6, "Data min/man is (%.2f, %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
    1711         psF32 binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / 1000.0f;
    1712         psTrace(__func__, 6, "Robust bin size is %.2f\n", binSize);
     1700        psF32 binSize = 0.0;
     1701        if ((iterate == 1) && (stats->options & PS_STAT_USE_BINSIZE)) {
     1702            // Set initial bin size to the specified value.
     1703            binSize = stats->binsize;
     1704            psTrace(__func__, 6, "Setting initial robust bin size to %.2f\n", binSize);
     1705        } else {
     1706            // Determine the bin size of the robust histogram.  This is done
     1707            // by computing the total range of data values and dividing by 1000.0.
     1708
     1709            rc = p_psVectorMin(myVector, tmpMaskVec, 1, tmpStatsMinMax);
     1710            rc|= p_psVectorMax(myVector, tmpMaskVec, 1, tmpStatsMinMax);
     1711            if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) {
     1712                psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
     1713                psFree(tmpStatsMinMax);
     1714                psFree(tmpMaskVec);
     1715                psFree(tmpScalar);
     1716                psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
     1717                return(1);
     1718            }
     1719            psTrace(__func__, 6, "Data min/man is (%.2f, %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
     1720            binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / 1000.0f;
     1721            psTrace(__func__, 6, "Robust bin size is %.2f\n", binSize);
     1722        }
     1723        psTrace(__func__, 6, "Initial robust bin size is %.2f\n", binSize);
    17131724
    17141725        //
    … …  
    17261737            }
    17271738            // XXX: Set these to the number of unmasked data points?
    1728             stats->robustNfit = 0.0;
    1729             stats->robustN50 = 0.0;
     1739            stats->robustN50 = 0;
    17301740            psFree(tmpStatsMinMax);
    17311741            psFree(tmpMaskVec);
    … …  
    19001910        sigma = (binHiF32 - binLoF32) / 2.0;
    19011911        psTrace(__func__, 6, "The current sigma is %f.\n", sigma);
     1912        stats->robustStdev = sigma;
    19021913
    19031914        //
    … …  
    19211932            psFree(robustHistogram);
    19221933            psFree(cumulativeRobustHistogram);
    1923 
     1934            iterate++;
    19241935        } else {
    19251936            psTrace(__func__, 6, "*************: No more iteration.  sigma is %f\n", sigma);
    1926             iterate = false;
     1937            iterate = 0;
    19271938        }
    19281939    }
    … …  
    19811992    psTrace(__func__, 6, "The 25 and 75 percent data point exact positions are (%f, %f).\n", binLo25F32, binHi25F32);
    19821993
    1983     //
    1984     //
    1985     // New algorithm for calculated mean and stdev.
    1986     //
    1987     //
    19881994    psS32 N50 = 0;
    19891995    for (psS32 i = 0 ; i < myVector->n ; i++) {
     1996        // XXX: use maskVal here?
    19901997        if ((0 == tmpMaskVec->data.U8[i]) &&
    19911998                (binLo25F32 <= myVector->data.F32[i]) &&
    … …  
    19972004    psTrace(__func__, 6, "The robustN50 is %d.\n", N50);
    19982005
    1999     psF32 dN = (psF32) (0.17 * N50);
    2000     if (dN < 1.0) {
    2001         dN = 1.0;
    2002     } else if (dN > 4.0) {
    2003         dN = 4.0;
    2004     }
    2005     psF32 newBinSize = sigma / dN;
    2006 
    2007     rc = p_psVectorMin(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax);
    2008     rc|= p_psVectorMax(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax);
    2009     if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) {
    2010         psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
    2011         psFree(tmpStatsMinMax);
    2012         psFree(robustHistogram);
    2013         psFree(cumulativeRobustHistogram);
    2014         psFree(tmpScalar);
    2015         psFree(tmpMaskVec);
    2016         psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
    2017         return(1);
    2018     }
    2019 
    2020     numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / newBinSize);
    2021     psTrace(__func__, 6, "The new min/max values are (%f, %f).\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
    2022     psTrace(__func__, 6, "The new bin size is %f.\n", newBinSize);
    2023     psTrace(__func__, 6, "The numBins is %d\n", numBins);
    2024 
    2025     psHistogram *newHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins);
    2026     newHistogram = psVectorHistogram(newHistogram, myVector, errors, tmpMaskVec, maskVal|1);
    2027     if (psTraceGetLevel(__func__) >= 8) {
    2028         PS_VECTOR_PRINT_F32(newHistogram->nums);
    2029     }
    2030 
    2031     //
    2032     // Smooth the resulting histogram with a Gaussian with sigma_s = 1
    2033     // bin.
    2034     //
    2035     psF32 sigma_s = newBinSize;
    2036 
    2037     psVector *newHistogramSmoothed = p_psVectorSmoothHistGaussian(newHistogram, sigma_s);
    2038     if (psTraceGetLevel(__func__) >= 8) {
    2039         PS_VECTOR_PRINT_F32(newHistogramSmoothed);
    2040     }
    2041 
    2042     //
    2043     // Find the bin with the peak value in the range 2 sigma of the
    2044     // robust histogram median.
    2045     //
    2046 
    2047     psS32 binMin = 0;
    2048     psS32 binMax = 0;
    2049     tmpScalar->data.F32 = stats->robustMedian - (2.0 * sigma);
    2050     if (tmpScalar->data.F32 <= newHistogram->bounds->data.F32[0]) {
    2051         binMin = 0;
    2052     } else {
    2053         binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
    2054     }
    2055 
    2056     tmpScalar->data.F32 = stats->robustMedian + (2.0 + sigma);
    2057     if (tmpScalar->data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) {
    2058         binMax = newHistogram->bounds->n-1;
    2059     } else {
    2060         binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
    2061     }
    2062     if ((binMin < 0) || (binMax < 0)) {
    2063         psError(PS_ERR_UNKNOWN, false, "Failed to calculate the +- 2.0 * sigma bins\n");
    2064         psFree(tmpMaskVec);
    2065         psFree(robustHistogram);
    2066         psFree(cumulativeRobustHistogram);
    2067         psFree(tmpStatsMinMax);
    2068         psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
    2069         return(1);
    2070     }
    2071 
    2072     psS32 binNum = binMin;
    2073     psF32 binMaxNums = newHistogramSmoothed->data.F32[binNum];
    2074     for (psS32 i = binMin+1 ; i <= binMax ; i++) {
    2075         if (newHistogramSmoothed->data.F32[i] > binMaxNums) {
    2076             binNum = i;
    2077             binMaxNums = newHistogramSmoothed->data.F32[i];
    2078         }
    2079     }
    2080     psTrace(__func__, 6, "The peak bin is %d, with %f data.n", binNum, binMaxNums);
    2081 
    2082     //
    2083     // Fit a Gaussian to the bins in the range 20 sigma of the robust
    2084     // histogram median.
    2085     //
    2086     tmpScalar->data.F32 = stats->robustMedian - (20.0 * sigma);
    2087     if (tmpScalar->data.F32 < tmpStatsMinMax->min) {
    2088         binMin = 0;
    2089     } else {
    2090         binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
    2091         // XXX: check for errors here.
    2092     }
    2093     tmpScalar->data.F32 = stats->robustMedian + (20.0 * sigma);
    2094     if (tmpScalar->data.F32 > tmpStatsMinMax->max) {
    2095         binMax = newHistogramSmoothed->n - 1;
    2096     } else {
    2097         binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
    2098         // XXX: check for errors here.
    2099     }
    2100     psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32);
    2101     psVector *xTmp = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32);
    2102     psArray *x = psArrayAlloc((1 + (binMax - binMin)));
    2103     stats->robustNfit = 0;
    2104     psS32 j = 0;
    2105 
    2106     for (psS32 i = binMin ; i <= binMax ; i++) {
    2107         y->data.F32[j] = newHistogramSmoothed->data.F32[i];
    2108         x->data[j] = (psPtr *) psVectorAlloc(1, PS_TYPE_F32);
    2109         ((psVector *) x->data[j])->data.F32[0] = PS_BIN_MIDPOINT(newHistogram, i);
    2110         xTmp->data.F32[j] = PS_BIN_MIDPOINT(newHistogram, i);
    2111 
    2112         stats->robustNfit+= newHistogramSmoothed->data.F32[i];
    2113         j++;
    2114     }
    2115     if (psTraceGetLevel(__func__) >= 8) {
    2116         // XXX: Print the x array somehow.
    2117         PS_VECTOR_PRINT_F32(y);
    2118     }
    2119 
    2120     // XXX: Use the min/max routines for this
    2121     psF32 minY = FLT_MAX;
    2122     psF32 maxY = -FLT_MAX;
    2123     for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) {
    2124         if (y->data.F32[i] > maxY) {
    2125             maxY = y->data.F32[i];
    2126         }
    2127         if (y->data.F32[i] < minY) {
    2128             minY = y->data.F32[i];
    2129         }
    2130     }
    2131     //
    2132     // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0])
    2133     // XXX: Use the normalize routines for this.
    2134     //
    2135     for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) {
    2136         y->data.F32[i]= (y->data.F32[i] - minY) / (maxY - minY);
    2137     }
    2138 
    2139     //
    2140     psMinimization *min = psMinimizationAlloc(100, 0.01);
    2141     psVector *params = psVectorAlloc(2, PS_TYPE_F32);
    2142     // Initial guess for the mean ([0]) and standard dev.
    2143     params->data.F32[0] = stats->robustMedian;
    2144     params->data.F32[1] = sigma;
    2145     if (psTraceGetLevel(__func__) >= 8) {
    2146         PS_VECTOR_PRINT_F32(params);
    2147         PS_VECTOR_PRINT_F32(y);
    2148     }
    2149     psFree(xTmp);
    2150     rcBool = psMinimizeLMChi2(min, NULL, params, NULL, x, y, NULL, psMinimizeLMChi2Gauss1D);
    2151 
    2152     if (rcBool != true) {
    2153         psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
    2154         psFree(tmpStatsMinMax);
    2155         psFree(robustHistogram);
    2156         psFree(cumulativeRobustHistogram);
    2157         psFree(tmpScalar);
     2006    // ************************************************************************
     2007    if ((stats->options & PS_STAT_FITTED_MEAN) ||
     2008            (stats->options & PS_STAT_FITTED_STDEV)) {
     2009        //
     2010        // New algorithm for calculated mean and stdev.
     2011        //
     2012
     2013        // XXX: Where is this documented?
     2014        psF32 dN = (psF32) (0.17 * N50);
     2015        if (dN < 1.0) {
     2016            dN = 1.0;
     2017        } else if (dN > 4.0) {
     2018            dN = 4.0;
     2019        }
     2020        psF32 newBinSize = sigma / dN;
     2021
     2022        rc = p_psVectorMin(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax);
     2023        rc|= p_psVectorMax(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax);
     2024        if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) {
     2025            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
     2026            psFree(tmpStatsMinMax);
     2027            psFree(robustHistogram);
     2028            psFree(cumulativeRobustHistogram);
     2029            psFree(tmpScalar);
     2030            psFree(tmpMaskVec);
     2031            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
     2032            return(1);
     2033        }
     2034
     2035        numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / newBinSize);
     2036        psTrace(__func__, 6, "The new min/max values are (%f, %f).\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
     2037        psTrace(__func__, 6, "The new bin size is %f.\n", newBinSize);
     2038        psTrace(__func__, 6, "The numBins is %d\n", numBins);
     2039
     2040        psHistogram *newHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins);
     2041        newHistogram = psVectorHistogram(newHistogram, myVector, errors, tmpMaskVec, maskVal|1);
     2042        if (psTraceGetLevel(__func__) >= 8) {
     2043            PS_VECTOR_PRINT_F32(newHistogram->nums);
     2044        }
     2045
     2046        //
     2047        // FITTED STATISTICS HERE
     2048        //
     2049        // Smooth the resulting histogram with a Gaussian with sigma_s = 1
     2050        // bin.
     2051        //
     2052        psF32 sigma_s = newBinSize;
     2053
     2054        psVector *newHistogramSmoothed = p_psVectorSmoothHistGaussian(newHistogram, sigma_s);
     2055        if (psTraceGetLevel(__func__) >= 8) {
     2056            PS_VECTOR_PRINT_F32(newHistogramSmoothed);
     2057        }
     2058
     2059        //
     2060        // Find the bin with the peak value in the range 2 sigma of the
     2061        // robust histogram median.
     2062        //
     2063
     2064        psS32 binMin = 0;
     2065        psS32 binMax = 0;
     2066        tmpScalar->data.F32 = stats->robustMedian - (2.0 * sigma);
     2067        if (tmpScalar->data.F32 <= newHistogram->bounds->data.F32[0]) {
     2068            binMin = 0;
     2069        } else {
     2070            binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
     2071        }
     2072
     2073        tmpScalar->data.F32 = stats->robustMedian + (2.0 + sigma);
     2074        if (tmpScalar->data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) {
     2075            binMax = newHistogram->bounds->n-1;
     2076        } else {
     2077            binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
     2078        }
     2079        if ((binMin < 0) || (binMax < 0)) {
     2080            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the +- 2.0 * sigma bins\n");
     2081            psFree(tmpMaskVec);
     2082            psFree(robustHistogram);
     2083            psFree(cumulativeRobustHistogram);
     2084            psFree(tmpStatsMinMax);
     2085            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
     2086            return(1);
     2087        }
     2088
     2089        psS32 binNum = binMin;
     2090        psF32 binMaxNums = newHistogramSmoothed->data.F32[binNum];
     2091        for (psS32 i = binMin+1 ; i <= binMax ; i++) {
     2092            if (newHistogramSmoothed->data.F32[i] > binMaxNums) {
     2093                binNum = i;
     2094                binMaxNums = newHistogramSmoothed->data.F32[i];
     2095            }
     2096        }
     2097        psTrace(__func__, 6, "The peak bin is %d, with %f data.n", binNum, binMaxNums);
     2098
     2099        //
     2100        // Fit a Gaussian to the bins in the range 20 sigma of the robust
     2101        // histogram median.
     2102        //
     2103        tmpScalar->data.F32 = stats->robustMedian - (20.0 * sigma);
     2104        if (tmpScalar->data.F32 < tmpStatsMinMax->min) {
     2105            binMin = 0;
     2106        } else {
     2107            binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
     2108            // XXX: check for errors here.
     2109        }
     2110        tmpScalar->data.F32 = stats->robustMedian + (20.0 * sigma);
     2111        if (tmpScalar->data.F32 > tmpStatsMinMax->max) {
     2112            binMax = newHistogramSmoothed->n - 1;
     2113        } else {
     2114            binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
     2115            // XXX: check for errors here.
     2116        }
     2117        psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32);
     2118        psVector *xTmp = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32);
     2119        psArray *x = psArrayAlloc((1 + (binMax - binMin)));
     2120        psS32 j = 0;
     2121
     2122        for (psS32 i = binMin ; i <= binMax ; i++) {
     2123            y->data.F32[j] = newHistogramSmoothed->data.F32[i];
     2124            x->data[j] = (psPtr *) psVectorAlloc(1, PS_TYPE_F32);
     2125            ((psVector *) x->data[j])->data.F32[0] = PS_BIN_MIDPOINT(newHistogram, i);
     2126            xTmp->data.F32[j] = PS_BIN_MIDPOINT(newHistogram, i);
     2127            j++;
     2128        }
     2129        if (psTraceGetLevel(__func__) >= 8) {
     2130            // XXX: Print the x array somehow.
     2131            PS_VECTOR_PRINT_F32(y);
     2132        }
     2133
     2134        // XXX: Use the min/max routines for this
     2135        psF32 minY = FLT_MAX;
     2136        psF32 maxY = -FLT_MAX;
     2137        for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) {
     2138            if (y->data.F32[i] > maxY) {
     2139                maxY = y->data.F32[i];
     2140            }
     2141            if (y->data.F32[i] < minY) {
     2142                minY = y->data.F32[i];
     2143            }
     2144        }
     2145        //
     2146        // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0])
     2147        // XXX: Use the normalize routines for this.
     2148        //
     2149        for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) {
     2150            y->data.F32[i]= (y->data.F32[i] - minY) / (maxY - minY);
     2151        }
     2152
     2153        //
     2154        psMinimization *min = psMinimizationAlloc(100, 0.01);
     2155        psVector *params = psVectorAlloc(2, PS_TYPE_F32);
     2156        // Initial guess for the mean ([0]) and standard dev.
     2157        params->data.F32[0] = stats->robustMedian;
     2158        params->data.F32[1] = sigma;
     2159        if (psTraceGetLevel(__func__) >= 8) {
     2160            PS_VECTOR_PRINT_F32(params);
     2161            PS_VECTOR_PRINT_F32(y);
     2162        }
     2163        psFree(xTmp);
     2164        rcBool = psMinimizeLMChi2(min, NULL, params, NULL, x, y, NULL, psMinimizeLMChi2Gauss1D);
     2165
     2166        if (rcBool != true) {
     2167            psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
     2168            psFree(tmpStatsMinMax);
     2169            psFree(robustHistogram);
     2170            psFree(cumulativeRobustHistogram);
     2171            psFree(tmpScalar);
     2172            psFree(newHistogram);
     2173            psFree(x);
     2174            psFree(y);
     2175            psFree(min);
     2176            psFree(params);
     2177            psFree(tmpMaskVec);
     2178            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
     2179            return(1);
     2180        }
     2181        if (psTraceGetLevel(__func__) >= 8) {
     2182            PS_VECTOR_PRINT_F32(params);
     2183        }
     2184
     2185        //
     2186        // The fitted mean mean_r is derived directly from the fitted
     2187        // Gaussian mean.
     2188        //
     2189        stats->fittedMean = params->data.F32[0];
     2190        psTrace(__func__, 6, "The fitted mean is %f.\n", params->data.F32[0]);
     2191
     2192        //
     2193        // The fitted standard deviation, SIGMA_r is determined by
     2194        // subtracting the smoothing scale in quadrature:
     2195        //     SIGMA_r^2 = SIGMA^2 - sigma_s^2
     2196        //
     2197        stats->fittedStdev = sqrt(PS_SQR(params->data.F32[1]) - PS_SQR(sigma_s));
     2198        psTrace(__func__, 6, "The fitted stdev is %f.\n", stats->fittedStdev);
     2199
     2200        psFree(newHistogramSmoothed);
    21582201        psFree(newHistogram);
    21592202        psFree(x);
    … …  
    21612204        psFree(min);
    21622205        psFree(params);
    2163         psFree(tmpMaskVec);
    2164         psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
    2165         return(1);
    2166     }
    2167     if (psTraceGetLevel(__func__) >= 8) {
    2168         PS_VECTOR_PRINT_F32(params);
    2169     }
    2170 
    2171     //
    2172     // The robust mean mean_r is derived directly from the fitted
    2173     // Gaussian mean.
    2174     //
    2175     stats->robustMean = params->data.F32[0];
    2176     psTrace(__func__, 6, "The robust mean is %f.\n", params->data.F32[0]);
    2177 
    2178     //
    2179     // The robust standard deviation, SIGMA_r is determined by
    2180     // subtracting the smoothing scale in quadrature:
    2181     //     SIGMA_r^2 = SIGMA^2 - sigma_s^2
    2182     //
    2183     stats->robustStdev = sqrt(PS_SQR(params->data.F32[1]) - PS_SQR(sigma_s));
    2184     psTrace(__func__, 6, "The robust stdev is %f.\n", stats->robustStdev);
    2185 
    2186     psFree(newHistogramSmoothed);
     2206    }
     2207
    21872208    psFree(tmpStatsMinMax);
    21882209    psFree(cumulativeRobustHistogram);
    21892210    psFree(tmpScalar);
    2190     psFree(newHistogram);
    2191     psFree(x);
    2192     psFree(y);
    2193     psFree(min);
    2194     psFree(params);
    21952211    psFree(tmpMaskVec);
    21962212    psFree(robustHistogram);
    … …  
    22362252    newStruct->sampleUQ = NAN;
    22372253    newStruct->sampleLQ = NAN;
    2238     newStruct->robustMean = NAN;
    22392254    newStruct->robustMedian = NAN;
    2240     newStruct->robustMode = NAN;
    22412255    newStruct->robustStdev = NAN;
    22422256    newStruct->robustUQ = NAN;
    22432257    newStruct->robustLQ = NAN;
    22442258    newStruct->robustN50 = -1;            // XXX: This is never used
    2245     newStruct->robustNfit = -1;
     2259    newStruct->fittedMean = NAN;
     2260    newStruct->fittedStdev = NAN;
     2261    newStruct->fittedNfit = -1;
    22462262    newStruct->clippedMean = NAN;
    22472263    newStruct->clippedStdev = NAN;
    22482264    newStruct->clippedNvalues = -1;     // XXX: This is never used
    2249     // XXX: Where do these values come from?
    22502265    newStruct->clipSigma = 3.0;
    22512266    newStruct->clipIter = 3;
    … …  
    22532268    newStruct->max = NAN;
    22542269    newStruct->binsize = NAN;          // XXX: This is never used
     2270    newStruct->nSubsample = 100000;
    22552271    newStruct->options = options;
    22562272
    2257     psTrace(__func__, 3, "---- %s(psStats) end  ----\n", __func__);
     2273    psTrace(__func__, 3, "---- %s() end  ----\n", __func__);
    22582274    return (newStruct);
    22592275}
    … …  
    23182334    newHist->uniform = true;
    23192335
    2320     psTrace(__func__, 3, "---- %s(psHistogram) end  ----\n", __func__);
     2336    psTrace(__func__, 3, "---- %s() end  ----\n", __func__);
    23212337    return (newHist);
    23222338}
    … …  
    23632379    newHist->uniform = false;
    23642380
    2365     psTrace(__func__, 3, "---- %s(psHistogram) end  ----\n", __func__);
     2381    psTrace(__func__, 3, "---- %s() end  ----\n", __func__);
    23662382    return (newHist);
    23672383}
    … …  
    25942610    }
    25952611
    2596     psTrace(__func__, 3, "---- %s(psHistogram) end  ----\n", __func__);
     2612    psTrace(__func__, 3, "---- %s() end  ----\n", __func__);
    25972613    return (out);
    25982614}
    … …  
    26742690    }
    26752691
    2676     psTrace(__func__, 4,"---- %s(psVector) end  ----\n", __func__);
     2692    psTrace(__func__, 4,"---- %s() end  ----\n", __func__);
    26772693    return (tmp);
    26782694}
    … …  
    27082724        PS_ASSERT_VECTOR_TYPE(errors, in->type.type, stats);
    27092725    }
     2726    // XXX: Assert that "in" is F64, F32, U16, or S8
    27102727
    27112728    psVector* inF32 = NULL;
    … …  
    27292746    }
    27302747
     2748    if ((stats->options & PS_STAT_USE_BINSIZE) && (stats->min >= stats->max)) {
     2749        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(stats->binsize, 0.0, stats);
     2750    }
     2751
    27312752    // ************************************************************************
    27322753    if (stats->options & PS_STAT_SAMPLE_MEAN) {
    27332754        if (0 != p_psVectorSampleMean(inF32, errorsF32, mask, maskVal, stats)) {
    2734             psLogMsg(__func__, PS_LOG_WARN,
    2735                      "WARNING: psVectorStats(): p_psVectorSampleMean() returned an error.\n");
    2736         }
    2737     }
     2755            psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMean() returned an error.\n");
     2756            stats->sampleMean = NAN;
     2757        }
     2758    }
     2759
    27382760    // ************************************************************************
    27392761    if (stats->options & PS_STAT_SAMPLE_MEDIAN) {
    27402762        if (false == p_psVectorSampleMedian(inF32, mask, maskVal, stats)) {
    2741             psLogMsg(__func__, PS_LOG_WARN,
    2742                      "WARNING: psVectorStats(): p_psVectorSampleMedian() returned an error.\n");
    2743         }
    2744     }
     2763            psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMedian() returned an error.\n");
     2764            stats->sampleMedian = NAN;
     2765        }
     2766    }
     2767
    27452768    // ************************************************************************
    27462769    if (stats->options & PS_STAT_SAMPLE_STDEV) {
    27472770        if (0 != p_psVectorSampleMean(inF32, errorsF32, mask, maskVal, stats)) {
    2748             psLogMsg(__func__, PS_LOG_WARN,
    2749                      "WARNING: psVectorStats(): p_psVectorSampleMean() returned an error.\n");
    2750         }
    2751         p_psVectorSampleStdev(inF32, errorsF32, mask, maskVal, stats);
    2752     }
     2771            psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMean() returned an error.\n");
     2772            stats->sampleMean = NAN;
     2773        } else {
     2774            p_psVectorSampleStdev(inF32, errorsF32, mask, maskVal, stats);
     2775        }
     2776    }
     2777
    27532778    // ************************************************************************
    27542779    if (stats->options & PS_STAT_SAMPLE_QUARTILE) {
    27552780        if (false == p_psVectorSampleQuartiles(inF32, mask, maskVal, stats)) {
    2756             psLogMsg(__func__, PS_LOG_WARN,
    2757                      "WARNING: psVectorStats(): p_psVectorSampleQuartiles() returned an error.\n");
    2758         }
    2759     }
    2760     // Since the various robust stats quantities share much computation, they
    2761     // are grouped together in a single private function:
    2762     // p_psVectorRobustStats()
    2763     if ((stats->options & PS_STAT_ROBUST_MEAN) ||
    2764             (stats->options & PS_STAT_ROBUST_MEDIAN) ||
    2765             (stats->options & PS_STAT_ROBUST_MODE) ||
     2781            psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleQuartiles() returned an error.\n");
     2782            stats->sampleLQ = NAN;
     2783            stats->sampleUQ = NAN;
     2784        }
     2785    }
     2786
     2787    // ************************************************************************
     2788    if ((stats->options & PS_STAT_ROBUST_MEDIAN) ||
    27662789            (stats->options & PS_STAT_ROBUST_STDEV) ||
    2767             (stats->options & PS_STAT_ROBUST_QUARTILE)) {
     2790            (stats->options & PS_STAT_ROBUST_QUARTILE) ||
     2791            (stats->options & PS_STAT_FITTED_MEAN) ||
     2792            (stats->options & PS_STAT_FITTED_STDEV)) {
    27682793        if (0 != p_psVectorRobustStats(inF32, errorsF32, mask, maskVal, stats)) {
    2769             psError(PS_ERR_UNKNOWN, false,
    2770                     PS_ERRORTEXT_psStats_STATS_FAILED);
    2771             // XXX: Set to NAN
    2772             // XXX: Is this the right thing to do?
    2773             // XXX: If so, do it for other funcs?
     2794            psError(PS_ERR_UNKNOWN, false, PS_ERRORTEXT_psStats_STATS_FAILED);
    27742795            psFree(stats);
    27752796            psTrace(__func__, 3,"---- %s(NULL) end  ----\n", __func__);
    … …  
    27782799    }
    27792800
    2780     // XXX: Different conditions for return -1 and -2?
     2801    // ************************************************************************
    27812802    if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
    27822803        psS32 rc = p_psVectorClippedStats(inF32, errorsF32, mask, maskVal, stats);
    2783         if (-1 == rc) {
    2784             psError(PS_ERR_UNKNOWN, false,
    2785                     "Failed to calculate clipped statistics for input psVector.\n");
     2804        if (rc < 0) {
     2805            psError(PS_ERR_UNKNOWN, false, "Failed to calculate clipped statistics for input psVector.\n");
    27862806            stats->clippedMean = NAN;
    27872807            stats->clippedStdev = NAN;
    2788         } else if (-2 == rc) {
    2789             psLogMsg(__func__, PS_LOG_WARN, "Failed to calculate clipped statistics for input psVector.");
    2790             stats->clippedMean = NAN;
    2791             stats->clippedStdev = NAN;
    2792         }
    2793     }
     2808        }
     2809    }
     2810
    27942811    // ************************************************************************
    27952812    if (stats->options & PS_STAT_MAX) {
    27962813        if (0 != p_psVectorMax(inF32, mask, maskVal, stats)) {
    2797             psError(PS_ERR_UNKNOWN, false,
    2798                     "Failed to calculate vector maximum");
     2814            psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector maximum");
    27992815            stats->max = NAN;
    28002816        }
    28012817    }
     2818
    28022819    // ************************************************************************
    28032820    if (stats->options & PS_STAT_MIN) {
    28042821        if (0 != p_psVectorMin(inF32, mask, maskVal, stats)) {
    2805             psError(PS_ERR_UNKNOWN, false,
    2806                     "Failed to calculate vector minimum");
     2822            psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector minimum");
    28072823            stats->min = NAN;
    28082824        }
    … …  
    28152831        psFree(errorsF32);
    28162832    }
    2817     psTrace(__func__, 3,"---- %s(psStats) end  ----\n", __func__);
     2833    psTrace(__func__, 3,"---- %s() end  ----\n", __func__);
    28182834    return (stats);
    28192835}
Note: See TracChangeset for help on using the changeset viewer.