Index: /trunk/psLib/src/math/psStats.c
===================================================================
--- /trunk/psLib/src/math/psStats.c	(revision 7869)
+++ /trunk/psLib/src/math/psStats.c	(revision 7870)
@@ -16,6 +16,6 @@
  * use ->min and ->max (PS_STAT_USE_RANGE)
  *
- *  @version $Revision: 1.176 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-06-30 02:20:06 $
+ *  @version $Revision: 1.177 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-07-12 03:46:11 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -172,4 +172,7 @@
     long count = 0;                     // Number of points contributing to this mean
 
+    psF32 *data = myVector->data.F32;   // Dereference
+    int numData = myVector->n;          // Number of data points
+
     // If PS_STAT_USE_RANGE is requested, then we enter a slightly different loop.
     if (!errors) {
@@ -177,40 +180,38 @@
             if (maskVector) {
                 // No errors, check range, check mask
-                for (long i = 0; i < myVector->n; i++) {
+                psU8 *maskData = maskVector->data.U8;
+                for (long i = 0; i < numData; i++) {
                     // Check if the data is with the specified range
-                    if (!(maskVal & maskVector->data.U8[i]) &&
-                            (stats->min <= myVector->data.F32[i]) &&
-                            (myVector->data.F32[i] <= stats->max)) {
-                        mean += myVector->data.F32[i];
+                    if (!(maskVal & maskData[i]) && stats->min <= data[i] && data[i] <= stats->max) {
+                        mean += data[i];
                         count++;
                     }
                 }
-
             } else {
                 // No errors, check range, no mask
-                for (long i = 0; i < myVector->n; i++) {
-                    if ((stats->min <= myVector->data.F32[i]) &&
-                            (myVector->data.F32[i] <= stats->max)) {
-                        mean += myVector->data.F32[i];
+                for (long i = 0; i < numData; i++) {
+                    if (stats->min <= data[i] && data[i] <= stats->max) {
+                        mean += data[i];
+                    }
+                }
+                count = numData;
+            }
+        } else {
+            if (maskVector) {
+                // No errors, no range check, check mask
+                psU8 *maskData = maskVector->data.U8;
+                for (long i = 0; i < numData; i++) {
+                    if (!(maskVal & maskData[i])) {
+                        mean += data[i];
                         count++;
                     }
                 }
-            }
-        } else {
-            if (maskVector) {
-                // No errors, no range check, check mask
-                for (long i = 0; i < myVector->n; i++) {
-                    if (!(maskVal & maskVector->data.U8[i])) {
-                        mean += myVector->data.F32[i];
-                        count++;
-                    }
-                }
 
             } else {
                 // No errors, no range check, no mask
-                for (long i = 0; i < myVector->n; i++) {
-                    mean += myVector->data.F32[i];
-                }
-                count = myVector->n;
+                for (long i = 0; i < numData; i++) {
+                    mean += data[i];
+                }
+                count = numData;
             }
         }
@@ -224,13 +225,13 @@
     } else {
         psF32 sumWeights = 0.0;         // The sum of the weights
+        psF32 *errorsData = errors->data.F32;
         if (stats->options & PS_STAT_USE_RANGE) {
-            if (maskVector != NULL) {
-                for (long i = 0; i < myVector->n; i++) {
+            if (maskVector) {
+                psU8 *maskData = maskVector->data.U8;
+                for (long i = 0; i < numData; i++) {
                     // Check if the data is with the specified range
-                    if (!(maskVal & maskVector->data.U8[i]) &&
-                            (stats->min <= myVector->data.F32[i]) &&
-                            (myVector->data.F32[i] <= stats->max)) {
-                        float weight = 1.0 / PS_SQR(errors->data.F32[i]);
-                        mean += myVector->data.F32[i] * weight;
+                    if (!(maskVal & maskData[i]) && stats->min <= data[i] && data[i] <= stats->max) {
+                        float weight = 1.0 / PS_SQR(errorsData[i]);
+                        mean += data[i] * weight;
                         sumWeights += weight;
                         count++;
@@ -239,8 +240,7 @@
             } else {
                 for (long i = 0; i < myVector->n; i++) {
-                    if ((stats->min <= myVector->data.F32[i]) &&
-                            (myVector->data.F32[i] <= stats->max)) {
-                        float weight = 1.0 / PS_SQR(errors->data.F32[i]);
-                        mean += myVector->data.F32[i] * weight;
+                    if (stats->min <= data[i] && data[i] <= stats->max) {
+                        float weight = 1.0 / PS_SQR(errorsData[i]);
+                        mean += data[i] * weight;
                         sumWeights += weight;
                         count++;
@@ -249,9 +249,10 @@
             }
         } else {
-            if (maskVector != NULL) {
+            if (maskVector) {
+                psU8 *maskData = maskVector->data.U8;
                 for (long i = 0; i < myVector->n; i++) {
-                    if (!(maskVal & maskVector->data.U8[i])) {
-                        float weight = 1.0 / PS_SQR(errors->data.F32[i]);
-                        mean += myVector->data.F32[i] * weight;
+                    if (!(maskVal & maskData[i])) {
+                        float weight = 1.0 / PS_SQR(errorsData[i]);
+                        mean += data[i] * weight;
                         sumWeights += weight;
                         count++;
@@ -260,6 +261,6 @@
             } else {
                 for (long i = 0; i < myVector->n; i++) {
-                    float weight = 1.0 / PS_SQR(errors->data.F32[i]);
-                    mean += myVector->data.F32[i] * weight;
+                    float weight = 1.0 / PS_SQR(errorsData[i]);
+                    mean += data[i] * weight;
                     sumWeights += weight;
                 }
@@ -731,20 +732,8 @@
 
     if (errors) {
-        #if 1
-        // XXX: The ADD specifies this as the definition of the standard
-        // deviation if the errors are known.  Verify this with IfA: none of
-        // the data points in the vector are used.  Verify that the masks and
-        // data ranges are used correctly.
         stats->sampleStdev = (1.0 / sqrtf(errorDivisor));
-        #else
-
-        stats->sampleStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / (float)count)) /
-                                  (float)(count - 1));
-        #endif
-
     } else {
         stats->sampleStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / (float)count)) /
                                   (float)(count - 1));
-
     }
     psTrace(__func__, 4, "---- %s() end ----\n", __func__);
@@ -789,12 +778,14 @@
     psStats *statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
 
-    // We've already copied the mask vector, so it's safe to simply increment the reference counter.
-    psVector *tmpMask = NULL;
+    // We copy the mask vector, to preserve the original
+    psVector *tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
+    tmpMask->n = tmpMask->nalloc;
+    psVectorInit(tmpMask, 0);
     if (maskVector) {
-        tmpMask = psMemIncrRefCounter(maskVector);
-    } else {
-        tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
-        tmpMask->n = tmpMask->nalloc;
-        psVectorInit(tmpMask, 0);
+        for (long i = 0; i < myVector->n; i++) {
+            if (maskVector->data.U8[i] & maskVal) {
+                tmpMask->data.U8[i] = 0xff;
+            }
+        }
     }
 
@@ -840,5 +831,5 @@
         if (errors) {
             for (long j = 0; j < myVector->n; j++) {
-                if (! tmpMask->data.U8[j] & maskVal &&
+                if (!tmpMask->data.U8[j] &&
                         fabsf(myVector->data.F32[j] - clippedMean) > stats->clipSigma * errors->data.F32[j]) {
                     tmpMask->data.U8[j] = 0xff;
@@ -851,5 +842,5 @@
         } else {
             for (long j = 0; j < myVector->n; j++) {
-                if (! tmpMask->data.U8[j] & maskVal &&
+                if (!tmpMask->data.U8[j] &&
                         fabsf(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) {
                     tmpMask->data.U8[j] = 0xff;
@@ -1131,12 +1122,13 @@
     tmpScalar.type.type = PS_TYPE_F32;
 
-    // We need a mask for these statistics, so if the one we're given is NULL, make our own.
-    psVector *mask = NULL;              // The actual mask we will use
+    psVector *mask = psVectorAlloc(myVector->n, PS_TYPE_MASK); // The actual mask we will use
+    mask->n = myVector->n;
+    psVectorInit(mask, 0);
     if (maskVector) {
-        mask = psMemIncrRefCounter(maskVector);
-    } else {
-        mask = psVectorAlloc(myVector->n, PS_TYPE_U8);
-        mask->n = myVector->n;
-        psVectorInit(mask, 0);
+        for (long i = 0; i < myVector->n; i++) {
+            if (maskVector->data.U8[i] & maskVal) {
+                mask->data.U8[i] = 0xff;
+            }
+        }
     }
 
@@ -1440,5 +1432,5 @@
     long N50 = 0;
     for (long i = 0 ; i < myVector->n ; i++) {
-        if (!(mask->data.U8[i] & maskVal) &&
+        if (!mask->data.U8[i] &&
                 (binLo25F32 <= myVector->data.F32[i]) && (binHi25F32 >= myVector->data.F32[i])) {
             N50++;
@@ -2047,16 +2039,8 @@
     psVector *maskU8 = NULL;            // Input mask vector, U8 version
     if (mask) {
-        // We want a copy of the mask, since we may change values
-        maskU8 = psVectorCopy(NULL, mask, PS_TYPE_U8);
-    } else {
-        maskU8 = psVectorAlloc(in->n, PS_TYPE_U8);
-        maskU8->n = in->n;
-        memset(maskU8->data.U8, 0, in->n * sizeof(psU8));
-    }
-
-    // Mask bad data
-    for (long i = 0; i < inF32->n; i++) {
-        if (!isfinite(inF32->data.F32[i])) {
-            maskU8->data.U8[i] = 0xff;
+        if (mask->type.type == PS_TYPE_MASK) {
+            maskU8 = psMemIncrRefCounter((psPtr)mask);
+        } else {
+            maskU8 = psVectorCopy(NULL, mask, PS_TYPE_MASK);
         }
     }