Index: /trunk/psLib/src/dataManip/psStats.c
===================================================================
--- /trunk/psLib/src/dataManip/psStats.c	(revision 4100)
+++ /trunk/psLib/src/dataManip/psStats.c	(revision 4101)
@@ -9,6 +9,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.129 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-05-25 20:26:55 $
+ *  @version $Revision: 1.130 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-06-03 22:30:20 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -1628,86 +1628,37 @@
     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;
-        /* XXX: Old, remove.
-                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;
-                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];
-                }
-        */
-
-        psMinimization *min = psMinimizationAlloc(100, 0.01);
-        psVector *myParams = psVectorAlloc(2, PS_TYPE_F32);
-        psArray *myCoords = psArrayAlloc(2 * dL + 1);
+        // We fit a 1-D polynomial to the data.
+        // XXX: Since we are no longer fitting a 1-D Gaussian, we can probably
+        // remove some of theabove code that calculated the initial estimate
+        // for the mean and sigma.
+
         psVector *y = psVectorAlloc(2 * dL + 1, PS_TYPE_F32);
-        myParams->data.F32[0] = myMean;
-        myParams->data.F32[1] = myStdev;
-
         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);
+            // XXX: Should this be the natural log?
             y->data.F32[index] = robustHistogramVector->data.F32[i];
-        }
-
-        // XXX: Must fit a 2-D polynomial, not a Gaussian.  (bug 366)
-        psBool rc = psMinimizeLMChi2(min,
-                                     NULL,
-                                     myParams,
-                                     NULL,
-                                     myCoords,
-                                     y,
-                                     NULL,
-                                     (psMinimizeLMChi2Func) psMinimizeLMChi2Gauss1D);
-        if (rc == false) {
+            //            y->data.F32[index] = logf(robustHistogramVector->data.F32[i]);
+        }
+
+        psPolynomial1D *tmpPoly = psPolynomial1DAlloc(3, PS_POLYNOMIAL_ORD);
+        // XXX: What about the NULL x argument?
+        tmpPoly = psVectorFitPolynomial1D(tmpPoly, NULL, y, NULL);
+        if (tmpPoly == NULL) {
             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?
+                     "WARNING: failed fit a 1D polynomial.\n");
+        }
+        psF32 polyFitSigma = PS_SQRT_F32(-0.5 / tmpPoly->coeff[2]);
+        psF32 polyFitMean = tmpPoly->coeff[1] * PS_SQR(polyFitSigma);
+        // psF32 polyFitNorm = exp(tmpPoly->coedd[0] + PS_SQR(polyFitMean) / (2.0 * PS_SQR(polyFitSigma)));
 
         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];
-            }
+            stats->robustMean = polyFitMean;
         }
 
         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);
+            stats->robustStdev = polyFitSigma;
+        }
         psFree(y);
-        psFree(myParams);
+        psFree(tmpPoly);
     }
 
Index: /trunk/psLib/src/math/psStats.c
===================================================================
--- /trunk/psLib/src/math/psStats.c	(revision 4100)
+++ /trunk/psLib/src/math/psStats.c	(revision 4101)
@@ -9,6 +9,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.129 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-05-25 20:26:55 $
+ *  @version $Revision: 1.130 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-06-03 22:30:20 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -1628,86 +1628,37 @@
     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;
-        /* XXX: Old, remove.
-                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;
-                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];
-                }
-        */
-
-        psMinimization *min = psMinimizationAlloc(100, 0.01);
-        psVector *myParams = psVectorAlloc(2, PS_TYPE_F32);
-        psArray *myCoords = psArrayAlloc(2 * dL + 1);
+        // We fit a 1-D polynomial to the data.
+        // XXX: Since we are no longer fitting a 1-D Gaussian, we can probably
+        // remove some of theabove code that calculated the initial estimate
+        // for the mean and sigma.
+
         psVector *y = psVectorAlloc(2 * dL + 1, PS_TYPE_F32);
-        myParams->data.F32[0] = myMean;
-        myParams->data.F32[1] = myStdev;
-
         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);
+            // XXX: Should this be the natural log?
             y->data.F32[index] = robustHistogramVector->data.F32[i];
-        }
-
-        // XXX: Must fit a 2-D polynomial, not a Gaussian.  (bug 366)
-        psBool rc = psMinimizeLMChi2(min,
-                                     NULL,
-                                     myParams,
-                                     NULL,
-                                     myCoords,
-                                     y,
-                                     NULL,
-                                     (psMinimizeLMChi2Func) psMinimizeLMChi2Gauss1D);
-        if (rc == false) {
+            //            y->data.F32[index] = logf(robustHistogramVector->data.F32[i]);
+        }
+
+        psPolynomial1D *tmpPoly = psPolynomial1DAlloc(3, PS_POLYNOMIAL_ORD);
+        // XXX: What about the NULL x argument?
+        tmpPoly = psVectorFitPolynomial1D(tmpPoly, NULL, y, NULL);
+        if (tmpPoly == NULL) {
             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?
+                     "WARNING: failed fit a 1D polynomial.\n");
+        }
+        psF32 polyFitSigma = PS_SQRT_F32(-0.5 / tmpPoly->coeff[2]);
+        psF32 polyFitMean = tmpPoly->coeff[1] * PS_SQR(polyFitSigma);
+        // psF32 polyFitNorm = exp(tmpPoly->coedd[0] + PS_SQR(polyFitMean) / (2.0 * PS_SQR(polyFitSigma)));
 
         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];
-            }
+            stats->robustMean = polyFitMean;
         }
 
         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);
+            stats->robustStdev = polyFitSigma;
+        }
         psFree(y);
-        psFree(myParams);
+        psFree(tmpPoly);
     }