Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 3540) +++ trunk/psLib/src/math/psStats.c (revision 3547) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.122 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-29 19:41:56 $ + * @version $Revision: 1.123 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-29 22:34:59 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -1489,11 +1489,11 @@ if (fabs(binSize) <= FLT_EPSILON) { if (stats->options & PS_STAT_ROBUST_MEAN) { - stats->robustMean = stats->clippedMean; + stats->robustMean = tmpStats->clippedMean; } if (stats->options & PS_STAT_ROBUST_MEDIAN) { - stats->robustMedian = stats->clippedMean; + stats->robustMedian = tmpStats->clippedMean; } if (stats->options & PS_STAT_ROBUST_MODE) { - stats->robustMode = stats->clippedMean; + stats->robustMode = tmpStats->clippedMean; } if (stats->options & PS_STAT_ROBUST_STDEV) { @@ -1501,6 +1501,6 @@ } if (stats->options & PS_STAT_ROBUST_QUARTILE) { - stats->robustUQ = stats->clippedMean; - stats->robustLQ = stats->clippedMean; + stats->robustUQ = tmpStats->clippedMean; + stats->robustLQ = tmpStats->clippedMean; } // XXX: Set these to the number of unmasked data points? @@ -1512,6 +1512,6 @@ // Determine minimum and maximum values in the data vector. - if (isnan(stats->min)) { - if (0 != p_psVectorMin(myVector, maskVector, maskVal, stats)) { + if (isnan(tmpStats->min)) { + if (0 != p_psVectorMin(myVector, maskVector, maskVal, tmpStats)) { psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorMin(): p_psVectorMin() reported a NAN mean.\n"); @@ -1519,6 +1519,6 @@ } } - if (isnan(stats->max)) { - if (0 != p_psVectorMax(myVector, maskVector, maskVal, stats)) { + if (isnan(tmpStats->max)) { + if (0 != p_psVectorMax(myVector, maskVector, maskVal, tmpStats)) { psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorMin(): p_psVectorMax() reported a NAN mean.\n"); @@ -1531,6 +1531,6 @@ // bins, not the binSize. Also, if we get here, we know that // binSize != 0.0. - numBins = (psS32)((stats->max - stats->min) / binSize); - robustHistogram = psHistogramAlloc(stats->min, stats->max, numBins); + numBins = (psS32)((tmpStats->max - tmpStats->min) / binSize); + robustHistogram = psHistogramAlloc(tmpStats->min, tmpStats->max, numBins); // Populate the histogram array. @@ -1620,50 +1620,92 @@ myStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); - // Using the above (myMean, myStdev) as initial estimates, we fit a - // Gaussian to the robustHistogramVector. - psMinimization *min = psMinimizationAlloc(100, 0.1); - psVector *myParams = psVectorAlloc(2, PS_TYPE_F32); - psArray *myCoords = psArrayAlloc(robustHistogramVector->n); - psVector *y = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); - p_psNormalizeVectorRangeF32(robustHistogramVector, 0.0, 1.0); - for (i=0;in;i++) { - myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32); - ((psVector *) (myCoords->data[i]))->data.F32[0] = (psF32) i; - y->data.F32[i] = robustHistogramVector->data.F32[i]; - } - - myParams->data.F32[0] = myMean; - myParams->data.F32[1] = myStdev; - psImage *covar = NULL; - psMinimizeLMChi2(min, - NULL, - myParams, - NULL, - myCoords, - y, - NULL, - (psMinimizeLMChi2Func) psMinimizeLMChi2Gauss1D); - psFree(covar); - psFree(min); - psFree(myParams); - psFree(myCoords); - psFree(y); + + if ((stats->options & PS_STAT_ROBUST_MEAN) || + (stats->options & PS_STAT_ROBUST_STDEV)) { + + // Using the above (myMean, myStdev) as initial estimates, we fit a + // Gaussian to the robustHistogramVector. + psImage *covar = NULL; + psMinimization *min = psMinimizationAlloc(100, 0.1); + psVector *myParams = psVectorAlloc(2, PS_TYPE_F32); + psArray *myCoords = psArrayAlloc(robustHistogramVector->n); + psVector *y = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); + + + myParams->data.F32[0] = myMean; + myParams->data.F32[1] = myStdev; + + /* XXX: old. Remove this. + for (i=0;in;i++) { + myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32); + ((psVector *) (myCoords->data[i]))->data.F32[0] = (psF32) i; + y->data.F32[i] = robustHistogramVector->data.F32[i]; + } + */ + + for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) { + int index = i - modeBinNum + dL; + myCoords->data[index] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32); + ((psVector *) (myCoords->data[index]))->data.F32[0] = PS_BIN_MIDPOINT(robustHistogram, i); + y->data.F32[index] = robustHistogramVector->data.F32[i]; + } + + bool rc = psMinimizeLMChi2(min, + NULL, + myParams, + NULL, + myCoords, + y, + NULL, + (psMinimizeLMChi2Func) psMinimizeLMChi2Gauss1D); + if (rc == false) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: failed to minimize with psMinimizeLMChi2().\n"); + } + + // XXX: Verify this with IfA + // XXX: The check on the minimization is better than the difference from myMean. + // Do they still want this code? + + if (stats->options & PS_STAT_ROBUST_MEAN) { + if (fabs((myParams->data.F32[0] - myMean)/myMean) > 0.1) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: the fitted Gaussian has more than 10%% error for the mean.\n"); + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: Using the calculated mean instead of Gaussian-fitted mean."); + // "WARNING: Using the calculated mean instead of Gaussian-fitted mean.(calc, fit) is (%f, %f)\n", + // myMean, myParams->data.F32[0]); + stats->robustMean = myMean; + } else { + stats->robustMean = myParams->data.F32[0]; + } + } + + + if (stats->options & PS_STAT_ROBUST_STDEV) { + if (fabs((myParams->data.F32[1] - myStdev)/myStdev) > 0.1) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: the fitted Gaussian has more than 10%% error for the stdev.\n"); + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev."); + // "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev.(calc, fit) is (%f, %f)\n", + // myStdev, myParams->data.F32[1]); + stats->robustStdev = myStdev; + } else { + stats->robustStdev = myParams->data.F32[1]; + } + } + psFree(covar); + psFree(min); + psFree(myCoords); + psFree(y); + psFree(myParams); + } + + /************************************************************************** Set the appropriate members in the output stats struct. **************************************************************************/ - if (stats->options & PS_STAT_ROBUST_MEAN) { - if (fabs((myParams->data.F32[0] - myMean)/myMean) > 0.1) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: the fitted Gaussian has more than 10%% error for the mean.\n"); - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: Using the calculated mean instead of Gaussian-fitted mean."); - // "WARNING: Using the calculated mean instead of Gaussian-fitted mean.(calc, fit) is (%f, %f)\n", - // myMean, myParams->data.F32[0]); - stats->robustMean = myMean; - } else { - stats->robustMean = myParams->data.F32[0]; - } - } if (stats->options & PS_STAT_ROBUST_MODE) { @@ -1671,17 +1713,4 @@ } - if (stats->options & PS_STAT_ROBUST_STDEV) { - if (fabs((myParams->data.F32[1] - myStdev)/myStdev) > 0.1) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: the fitted Gaussian has more than 10%% error for the stdev.\n"); - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev."); - // "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev.(calc, fit) is (%f, %f)\n", - // myStdev, myParams->data.F32[1]); - stats->robustStdev = myStdev; - } else { - stats->robustStdev = myParams->data.F32[1]; - } - } // To determine the median, we fit a quadratic y=f(x) to the three bins @@ -2174,7 +2203,7 @@ *****************************************************************************/ psStats* psVectorStats(psStats* stats, - psVector* in, - psVector* errors, - psVector* mask, + const psVector* in, + const psVector* errors, + const psVector* mask, psU32 maskVal) { @@ -2197,10 +2226,10 @@ inF32 = p_psConvertToF32((psVector *) in); if (inF32 == NULL) { - inF32 = in; + inF32 = (psVector *) in; mustFreeVectorIn = 0; } errorsF32 = p_psConvertToF32((psVector *) errors); if (errorsF32 == NULL) { - errorsF32 = errors; + errorsF32 = (psVector *) errors; mustFreeVectorErrors = 0; }