Changeset 2340 for trunk/psLib/src/math/psStats.c

Timestamp:

Nov 12, 2004, 9:08:17 AM (22 years ago)

Author:

gusciora

Message:

...

File:

• trunk/psLib/src/math/psStats.c (modified) (34 diffs)

trunk/psLib/src/math/psStats.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.90 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-11-11 20:13:57 $
+ *  @version $Revision: 1.91 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-11-12 19:08:17 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -29 +29 @@
 #include "psVector.h"
 #include "psTrace.h"
+#include "psLogMsg.h"
 #include "psError.h"
 #include "psStats.h"
@@ -142 +143 @@
 /******************************************************************************
 p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the
-mean of the input vector.
-Inputs
-    myVector
-    maskVector
-    maskVal
-    stats
-Returns
-    NULL
- 
+mean of the input vector.  If there was a problem with the mean calculation,
+this routine sets stats->sampleMean to NAN.
  *****************************************************************************/
 void p_psVectorSampleMean(const psVector* restrict myVector,
@@ -220 +214 @@
 /******************************************************************************
 p_psVectorSampleMax(myVector, maskVector, maskVal, stats): calculates the
-max of the input vector.
-Inputs
-    myVector
-    maskVector
-    maskVal
-    stats
-Returns
-    NULL
+max of the input vector.  If there was a problem with the max calculation,
+this routine sets stats->max to NAN.
  *****************************************************************************/
 void p_psVectorMax(const psVector* restrict myVector,
@@ -289 +277 @@
 /******************************************************************************
 p_psVectorSampleMin(myVector, maskVector, maskVal, stats): calculates the
-minimum of the input vector.
-Inputs
-    myVector
-    maskVector
-    maskVal
-    stats
-Returns
-    NULL
+minimum of the input vector.  If there was a problem with the min calculation,
+this routine sets stats->min to NAN.
  *****************************************************************************/
 void p_psVectorMin(const psVector* restrict myVector,
@@ -357 +339 @@
 
 /******************************************************************************
+p_psVectorNValues(myVector, maskVector, maskVal, stats): This routine returns
+"true" if the inputPsVector as 1 or more valid elements (a valid element is an
+unmasked element within the specifined min/max range).  Otherwise, return
+"false".
+ *****************************************************************************/
+bool p_psVectorCheckNonEmpty(const psVector* restrict myVector,
+                             const psVector* restrict maskVector,
+                             psU32 maskVal,
+                             psStats* stats)
+{
+    PS_VECTOR_CHECK_NULL(myVector, -1);
+    PS_PTR_CHECK_NULL(stats, -1);
+    psS32 i = 0;                // Loop index variable
+
+    if (stats->options & PS_STAT_USE_RANGE) {
+        if (maskVector != NULL) {
+            for (i = 0; i < myVector->n; i++) {
+                if (!(maskVal & maskVector->data.U8[i]) &&
+                        (stats->min <= myVector->data.F32[i]) &&
+                        (myVector->data.F32[i] <= stats->max)) {
+                    return(true);
+                }
+            }
+        } else {
+            for (i = 0; i < myVector->n; i++) {
+                if ((stats->min <= myVector->data.F32[i]) &&
+                        (myVector->data.F32[i] <= stats->max)) {
+                    return(true);
+                }
+            }
+        }
+    } else {
+        if (maskVector != NULL) {
+            for (i = 0; i < myVector->n; i++) {
+                if (!(maskVal & maskVector->data.U8[i])) {
+                    return(true);
+                }
+            }
+        } else {
+            if (myVector->n > 0) {
+                return(true);
+            }
+        }
+    }
+    return(false);
+}
+
+
+/******************************************************************************
 p_psVectorNValues(myVector, maskVector, maskVal, stats): calculates the
 number of non-masked pixels in the vector that fall within the min/max
 range, if specified.
-Inputs
-    myVector
-    maskVector
-    maskVal
-    stats
-Returns
-    int
  *****************************************************************************/
 psS32 p_psVectorNValues(const psVector* restrict myVector,
@@ -412 +436 @@
 /******************************************************************************
 p_psVectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the
-median of the input vector.
-Inputs
-    myVector
-    maskVector
-    maskVal
-    stats
-Returns
-    NULL
+median of the input vector.  Returns true on success (including if there were
+no valid input vector elements).
  *****************************************************************************/
 bool p_psVectorSampleMedian(const psVector* restrict myVector,
@@ -429 +447 @@
     psVector* sortedVector = NULL;      // Temporary vector
     psS32 i = 0;                        // Loop index variable
-    psS32 count = 0;                    // # of points in this mean?
-    psS32 nValues = 0;                  // # of points in vector
+    psS32 count = 0;
+    psS32 nValues = 0;
 
     // Determine how many data points fit inside this min/max range
-    // and are not masked, if the maskVector is not NULL>
+    // and are not masked, if the maskVector is not NULL.
     nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
 
@@ -439 +457 @@
     // sampleMedian to?  Should we generate an error?
     if (nValues <= 0) {
-        printf("WARNING: p_psVectorSampleMedian(): no valid elements in input vector.\n");
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMedian(): no valid elements in input vector.\n");
         return(true);
-        stats->sampleMedian = 0;
+        stats->sampleMedian = NAN;
     }
 
@@ -569 +587 @@
         } else if (x.data.F32 >= lastBound) {
             jMin = robustHistogram->bounds->n - 1;
-
         }
 
@@ -575 +592 @@
         if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
             jMax = p_psVectorBinDisect(robustHistogram->bounds, &x);
+            if (jMax < 0) {
+                psError(PS_ERR_UNEXPECTED_NULL,
+                        false,
+                        PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM);
+                return(NULL);
+            }
         } else if (x.data.F32 <= firstBound) {
             jMax = 0;
         } else if (x.data.F32 >= lastBound) {
             jMax = robustHistogram->bounds->n - 1;
-
         }
 
@@ -606 +628 @@
     NULL
  *****************************************************************************/
-void p_psVectorSampleQuartiles(const psVector* restrict myVector,
+bool p_psVectorSampleQuartiles(const psVector* restrict myVector,
                                const psVector* restrict maskVector,
                                psU32 maskVal,
@@ -616 +638 @@
     psS32 count = 0;              // # of points in this mean.
     psS32 nValues = 0;            // # data points
-    float rangeMin = 0.0;       // Exclude data below this
-    float rangeMax = 0.0;       // Exclude date above this
 
     // Determine how many data points fit inside this min/max range
-    // and are not maxed, IF the maskVector is not NULL>
+    // and are not maxed, IF the maskVector is not NULL.
     nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
 
     // Allocate temporary vectors for the data.
     unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
-    sortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
 
     // Determine if we must only use data points within a min/max range.
     if (stats->options & PS_STAT_USE_RANGE) {
-        rangeMin = stats->min;
-        rangeMax = stats->max;
         // Store all non-masked data points within the min/max range
         // into the temporary vectors.
@@ -637 +654 @@
             for (i = 0; i < myVector->n; i++) {
                 if (!(maskVal & maskVector->data.U8[i]) &&
-                        (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
+                        (stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) {
                     unsortedVector->data.F32[count++] = myVector->data.F32[i];
                 }
@@ -643 +660 @@
         } else {
             for (i = 0; i < myVector->n; i++) {
-                if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
+                if ((stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) {
                     unsortedVector->data.F32[count++] = myVector->data.F32[i];
                 }
@@ -665 +682 @@
 
     // Sort the temporary vectors.
-    psVectorSort(sortedVector, unsortedVector);
+    sortedVector = psVectorSort(sortedVector, unsortedVector);
+    if (sortedVector == NULL) {
+        psError(PS_ERR_UNEXPECTED_NULL,
+                false,
+                PS_ERRORTEXT_psStats_STATS_SAMPLE_MEDIAN_SORT_PROBLEM);
+        return(false);
+    }
 
     // Calculate the quartile points exactly.
+    // XXX: We should probably do a bin-midpoint if the quartile is not an
+    // integer.
     stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)];
     stats->sampleLQ = sortedVector->data.F32[nValues / 4];
@@ -674 +699 @@
     psFree(unsortedVector);
     psFree(sortedVector);
+    return(true);
 }
 
@@ -700 +726 @@
     float sumSquares = 0.0;     // temporary variable
     float sumDiffs = 0.0;       // temporary variable
-    float rangeMin = 0.0;       // Exclude data below this
-    float rangeMax = 0.0;       // Exclude date above this
 
     // This procedure requires the mean.  If it has not been already
@@ -708 +732 @@
         p_psVectorSampleMean(myVector, maskVector, maskVal, stats);
     }
+    // If the mean is NAN, then generate a warning and set the stdev to NAN.
+    if (0 != isnan(stats->sampleMean)) {
+        stats->sampleStdev = NAN;
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): p_psVectorSampleMean() reported a NAN mean.\n");
+        return;
+    }
+
     mean = stats->sampleMean;
-
     if (stats->options & PS_STAT_USE_RANGE) {
         if (maskVector != NULL) {
             for (i = 0; i < myVector->n; i++) {
                 if (!(maskVal & maskVector->data.U8[i]) &&
-                        (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
+                        (stats->min <= myVector->data.F32[i]) &&
+                        (myVector->data.F32[i] <= stats->max)) {
                     diff = myVector->data.F32[i] - mean;
                     sumSquares += (diff * diff);
@@ -723 +754 @@
         } else {
             for (i = 0; i < myVector->n; i++) {
-                if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
+                if ((stats->min <= myVector->data.F32[i]) &&
+                        (myVector->data.F32[i] <= stats->max)) {
                     diff = myVector->data.F32[i] - mean;
                     sumSquares += (diff * diff);
@@ -752 +784 @@
         }
     }
-    if (countInt < 2) {
+    if (countInt <= 1) {
         stats->sampleStdev = NAN;
-        // XXX PS WARNING
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): no valid psVector elements.  Setting stats->sampleStdev = NAN.\n");
     } else {
         countFloat = (float)countInt;
         #ifdef DARWIN
 
-        stats->sampleStdev = (float)sqrt((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
+        stats->sampleStdev = (float)sqrt((sumSquares -
+                                          (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
         #else
 
-        stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
+        stats->sampleStdev = sqrtf((sumSquares -
+                                    (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
         #endif
 
@@ -779 +813 @@
 Returns
     0 for success.
+    -1: error
+    -2: warning
  *****************************************************************************/
 int p_psVectorClippedStats(const psVector* restrict myVector,
@@ -794 +830 @@
 
     // Ensure that stats->clipIter is within the proper range.
-    PS_INT_CHECK_RANGE(stats->clipIter,PS_CLIPPED_NUM_ITER_LB,PS_CLIPPED_NUM_ITER_UB,-1);
+    PS_INT_CHECK_RANGE(stats->clipIter,
+                       PS_CLIPPED_NUM_ITER_LB,
+                       PS_CLIPPED_NUM_ITER_UB, -1);
 
     // Ensure that stats->clipSigma is within the proper range.
-    PS_INT_CHECK_RANGE(stats->clipSigma,PS_CLIPPED_SIGMA_LB,PS_CLIPPED_SIGMA_UB,-1);
+    PS_INT_CHECK_RANGE(stats->clipSigma,
+                       PS_CLIPPED_SIGMA_LB,
+                       PS_CLIPPED_SIGMA_UB, -1);
 
     // We allocate a temporary mask vector since during the iterative
@@ -803 +843 @@
     // However, we do no want to modify the original mask vector.
     tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
-    tmpMask->n = myVector->n;
 
     // If we were called with a mask vector, then initialize the temporary
@@ -814 +853 @@
     // 1. Compute the sample median.
     p_psVectorSampleMedian(myVector, maskVector, maskVal, stats);
+    if (isnan(stats->sampleMean)) {
+        psLogMsg(__func__, PS_LOG_WARN, "Call to p_psVectorSampleMean returned NAN\n");
+        stats->clippedMean = NAN;
+        stats->clippedStdev = NAN;
+        return(-2);
+    }
 
     // 2. Compute the sample standard deviation.
     p_psVectorSampleStdev(myVector, maskVector, maskVal, stats);
+    if (isnan(stats->sampleStdev)) {
+        psLogMsg(__func__, PS_LOG_WARN, "Call to p_psVectorSampleStdev returned NAN\n");
+        stats->clippedMean = NAN;
+        stats->clippedStdev = NAN;
+        return(-2);
+    }
 
     // 3. Use the sample median as the first estimator of the mean X.
@@ -836 +887 @@
                 tmpMask->data.U8[i] = 0xff;
             }
-            // b) compute new mean and stdev
-            p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats);
-            p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats);
-
-            // c) Use the new mean for x
+        }
+
+        // b) compute new mean and stdev
+        p_psVectorSampleMean(myVector, tmpMask, maskVal, stats);
+        p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats);
+
+        // If the new mean and stdev are not NAN, then we use them.
+        // Otherwise, keep the old ones, and exit the loop.
+        if (! (isnan(stats->sampleMean) || isnan(stats->sampleStdev))) {
             clippedMean = stats->sampleMean;
-
-            // d) Use the new stdev for stdev
             clippedStdev = stats->sampleStdev;
-        }
-
-    }
+
+            // Exit loop
+            i = stats->clipIter;
+        }
+    }
+
     stats->sampleMean = oldStanMean;
     stats->sampleStdev = oldStanStdev;
@@ -865 +921 @@
 
 /*****************************************************************************
- *****************************************************************************/
-void p_psNormalizeVectorRangeF32(psVector* myData,
-                                 float outLow,
-                                 float outHigh)
-{
-    float min = PS_MAX_F32;
-    float max = -PS_MAX_F32;
-    psS32 i = 0;
-
-    for (i = 0; i < myData->n; i++) {
-        if (myData->data.F32[i] < min) {
-            min = myData->data.F32[i];
-        }
-        if (myData->data.F32[i] > max) {
-            max = myData->data.F32[i];
-        }
-    }
-
-    // Ensure that max!=min before we divide by (max-min)
-    if (FLT_EPSILON < fabs(max - min)) {
-        for (i = 0; i < myData->n; i++) {
-            myData->data.F32[i] = outLow + (myData->data.F32[i] - min) *
-                                  (outHigh - outLow) / (max - min);
-        }
-    } else {
-        // XXX: PS_WARNING
-        for (i = 0; i < myData->n; i++) {
-            myData->data.F32[i] = outLow;
-        }
-    }
-
-    for (i = 0; i < myData->n; i++) {
-        myData->data.F32[i] = outLow + (myData->data.F32[i] - min) * (outHigh - outLow) / (max - min);
-    }
-}
-void p_psNormalizeVectorRangeF64(psVector* myData,
-                                 float outLow,
-                                 float outHigh)
-{
-    float min = PS_MAX_F32;
-    float max = -PS_MAX_F32;
-    psS32 i = 0;
-
-    for (i = 0; i < myData->n; i++) {
-        if (myData->data.F32[i] < min) {
-            min = myData->data.F64[i];
-        }
-        if (myData->data.F32[i] > max) {
-            max = myData->data.F64[i];
-        }
-    }
-
-    // Ensure that max!=min before we divide by (max-min)
-    if (FLT_EPSILON < fabs(max - min)) {
-        for (i = 0; i < myData->n; i++) {
-            myData->data.F64[i] = outLow + (myData->data.F64[i] - min) * (outHigh - outLow) / (max - min);
-        }
-    } else {
-        // XXX: PS_WARNING
-        for (i = 0; i < myData->n; i++) {
-            myData->data.F64[i] = outLow;
-        }
-    }
-}
-
-/*****************************************************************************
-psNormalizeVectorRange(myData, low, high): this is a private function which
+psNormalizeVectorRange##TYPE(myData, low, high): this is a private function which
 normalizes the elements of a vector to a range between low and high.
- *****************************************************************************/
-void psNormalizeVectorRange(psVector* myData,
-                            float low,
-                            float high)
-{
-    if (myData->type.type == PS_TYPE_F32) {
-        p_psNormalizeVectorRangeF32(myData, low, high);
-    } else if (myData->type.type == PS_TYPE_F64) {
-        p_psNormalizeVectorRangeF64(myData, low, high);
-    } else {
-        char* strType;
-        PS_TYPE_NAME(strType,myData->type.type);
-        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
-                PS_ERRORTEXT_psStats_NOT_F32_F64,
-                strType);
-    }
-}
+ 
+XXX: Should we make the arguments psScalars?
+ *****************************************************************************/
+#define PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(TYPE) \
+void p_psNormalizeVectorRange##TYPE(psVector* myData, \
+                                    ps##TYPE outLow, \
+                                    ps##TYPE outHigh) \
+{ \
+    ps##TYPE min = (ps##TYPE) PS_MAX_##TYPE; \
+    ps##TYPE max = (ps##TYPE) -PS_MAX_##TYPE; \
+    psS32 i = 0; \
+    \
+    for (i = 0; i < myData->n; i++) { \
+        if (myData->data.TYPE[i] < min) { \
+            min = myData->data.TYPE[i]; \
+        } \
+        if (myData->data.TYPE[i] > max) { \
+            max = myData->data.TYPE[i]; \
+        } \
+    } \
+    \
+    /* Ensure that max!=min before we divide by (max-min) */ \
+    if (max != min) { \
+        for (i = 0; i < myData->n; i++) { \
+            myData->data.TYPE[i] = (outLow + (myData->data.TYPE[i] - min) * \
+                                    (outHigh - outLow) / (max - min)); \
+        } \
+    } else { \
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: (max==min).  Setting all elements to min.\n");
+        for (i = 0; i < myData->n; i++)
+        {
+            \
+            myData->data.TYPE[i] = outLow;
+            \
+        } \
+    } \
+} \
+
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U8)
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U16)
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U32)
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U64)
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S8)
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S16)
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S32)
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S64)
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F32)
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F64)
 
 /******************************************************************************
@@ -969 +990 @@
                        float getThisValue)
 {
+    PS_FLOAT_COMPARE(rangeLow, rangeHigh, NAN);
     psS32 numIterations = 0;
     float midpoint = 0.0;
@@ -1005 +1027 @@
  
 XXX: After you fit the polynomial, solve for X analytically.
+ 
+XXX: Check for errors in psLib routines that we call.
 *****************************************************************************/
 float fitQuadraticSearchForYThenReturnX(psVector *xVec,
@@ -1085 +1109 @@
 Returns
     0 on success.
+ 
+XXX: Check for errors in psLib routines that we call.
 *****************************************************************************/
 int p_psVectorRobustStats(const psVector* restrict myVector,
@@ -1256 +1282 @@
     psVector *y = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
 
-    psNormalizeVectorRange(robustHistogramVector, 0.0, 1.0);
+    psNormalizeVectorRangeF32(robustHistogramVector, 0.0, 1.0);
     for (i=0;i<robustHistogramVector->n;i++) {
         myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
@@ -1282 +1308 @@
     if (stats->options & PS_STAT_ROBUST_MEAN) {
         if (fabs((myParams->data.F32[0] - myMean)/myMean) > 0.1) {
-            printf("WARNING: the fitted Gaussian has more than 10% error for the mean.\n");
-            printf("Using the calculated mean (calc, fit) is (%f, %f)\n", myMean, myParams->data.F32[0]);
+            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 (calc, fit) is (%f, %f)\n",
+                     myMean, myParams->data.F32[0]);
             stats->robustMean = myMean;
         } else {
@@ -1296 +1325 @@
     if (stats->options & PS_STAT_ROBUST_STDEV) {
         if (fabs((myParams->data.F32[1] - myStdev)/myStdev) > 0.1) {
-            printf("WARNING: the fitted Gaussian has more than 10% error for the stdev.\n");
-            printf("Using the calculated stdev (calc, fit) is (%f, %f)\n", myStdev, myParams->data.F32[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 (calc, fit) is (%f, %f)\n",
+                     myMean, myParams->data.F32[0]);
             stats->robustStdev = myStdev;
         } else {
@@ -1691 +1723 @@
             psError(PS_ERR_UNKNOWN, false,
                     PS_ERRORTEXT_psStats_STATS_FAILED);
+            // XXX: Set to NAN
         }
     }
 
     if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
-        if (0 != p_psVectorClippedStats(inF32, mask, maskVal, stats)) {
+        if (-1 == p_psVectorClippedStats(inF32, mask, maskVal, stats)) {
             psError(PS_ERR_UNKNOWN, false,
-                    "Failed to calculate statistics for input psVector.");
+                    "Failed to calculate clipped statistics for input psVector.");
+            stats->clippedMean = NAN;
+            stats->clippedStdev = NAN;
+        } else if (-2 == p_psVectorClippedStats(inF32, mask, maskVal, stats)) {
+            psLogMsg(__func__, PS_LOG_WARN, "Failed to calculate clipped statistics for input psVector.");
+            stats->clippedMean = NAN;
+            stats->clippedStdev = NAN;
         }
     }