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

Timestamp:

Mar 29, 2005, 12:34:59 PM (21 years ago)

Author:

gusciora

Message:

Fixes bugs 338. Deeals with robust stats, and min->lastDelta.

File:

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

trunk/psLib/src/math/psStats.c

@@ -9 +9 @@
  *  @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 +1489 @@
     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 +1501 @@
         }
         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 +1512 @@
 
     // 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 +1519 @@
         }
     }
-    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 +1531 @@
     // 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 +1620 @@
     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;i<robustHistogramVector->n;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;i<robustHistogramVector->n;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 +1713 @@
     }
 
-    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 +2203 @@
  *****************************************************************************/
 psStats* psVectorStats(psStats* stats,
-                       psVector* in,
-                       psVector* errors,
-                       psVector* mask,
+                       const psVector* in,
+                       const psVector* errors,
+                       const psVector* mask,
                        psU32 maskVal)
 {
@@ -2197 +2226 @@
     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;
     }