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

Timestamp:

Mar 1, 2005, 11:20:51 AM (21 years ago)

Author:

gusciora

Message:

Fixed bugs reported in 316. Not tested yet.

File:

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

trunk/psLib/src/math/psStats.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.115 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-03-01 21:03:25 $
+ *  @version $Revision: 1.116 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-03-01 21:20:51 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -1003 +1003 @@
     -2: warning
  
-XXX: Create a new psStats object for the sample mean and sample stdev.
- 
 XXX: Do we really need to calculate median and mean?
  
@@ -1015 +1013 @@
                              psStats* stats)
 {
-    psS32 i = 0;                  // Loop index variable
-    psS32 j = 0;                  // Loop index variable
     psF32 clippedMean = 0.0;    // self-explanatory
     psF32 clippedStdev = 0.0;   // self-explanatory
-    psVector* tmpMask = NULL;   // Temporary vector
-    psStats *statsTmp = NULL;
+    psVector* tmpMask = NULL;   // Temporary vector for masks during iterations.
+    psStats *statsTmp = NULL;   // Temporary psStats struct.
+    psS32 rc = 0;               // Return code.
 
     // Ensure that stats->clipIter is within the proper range.
@@ -1032 +1029 @@
                        PS_CLIPPED_SIGMA_UB, -1);
 
+    // Allocate a psStats structure for calculating the mean, median, and
+    // stdev.
     if (statsTmp == NULL) {
         statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
@@ -1045 +1044 @@
     // mask vector with those values.
     if (maskVector != NULL) {
-        for (i = 0; i < tmpMask->n; i++) {
+        for (psS32 i = 0; i < tmpMask->n; i++) {
             tmpMask->data.U8[i] = maskVector->data.U8[i];
         }
     }
+
     // 1. Compute the sample median.
     p_psVectorSampleMedian(myVector, maskVector, maskVal, statsTmp);
@@ -1069 +1069 @@
     // 3. Use the sample median as the first estimator of the mean X.
     clippedMean = statsTmp->sampleMean;
+
     // 4. Use the sample stdev as the first estimator of the mean stdev.
     clippedStdev = statsTmp->sampleStdev;
 
-    // 5. Repeat N times:
-    for (i = 0; i < stats->clipIter; i++) {
+    // 5. Repeat N (stats->clipIter) times:
+    for (psS32 iter = 0; iter < stats->clipIter; iter++) {
         // a) Exclude all values x_i for which |x_i - x| > K * stdev
         if (errors != NULL) {
-            for (j = 0; j < myVector->n; j++) {
+            for (psS32 j = 0; j < myVector->n; j++) {
                 if (fabs(myVector->data.F32[j] - clippedMean) >
                         (stats->clipSigma * errors->data.F32[j])) {
@@ -1083 +1084 @@
             }
         } else {
-            for (j = 0; j < myVector->n; j++) {
+            for (psS32 j = 0; j < myVector->n; j++) {
                 if (fabs(myVector->data.F32[j] - clippedMean) >
                         (stats->clipSigma * clippedStdev)) {
@@ -1095 +1096 @@
         p_psVectorSampleStdev(myVector, errors, tmpMask, maskVal, statsTmp);
 
-        // If the new mean and stdev are not NAN, then we use them.
-        // Otherwise, keep the old ones, and exit the loop.
-        if (! (isnan(statsTmp->sampleMean) || isnan(statsTmp->sampleStdev))) {
+        // If the new mean and stdev are NAN, we must exit the loop.
+        // Otherwise, use the new results and continue.
+        if (isnan(statsTmp->sampleMean) || isnan(statsTmp->sampleStdev)) {
+            // Exit loop.  XXX: Should we throw an error or warning here?
+            iter = stats->clipIter;
+            rc = -1;
+        } else {
             clippedMean = statsTmp->sampleMean;
             clippedStdev = statsTmp->sampleStdev;
-        } else {
-            // Exit loop
-            // Should we throw an error or warning here?
-            i = stats->clipIter;
         }
     }
@@ -1117 +1118 @@
 
     psFree(tmpMask);
-    return(0);
+    return(rc);
 }